Out-of-equilibrium quantum magnetism and thermalization in a spin-3 many-body dipolar lattice system

Understanding quantum thermalization through entanglement build up in isolated quantum systems addresses fundamental questions on how unitary dynamics connects to statistical physics. Spin systems made of long-range interacting atoms offer an ideal experimental platform to investigate this question. Here, we study the spin dynamics and approach towards local thermal equilibrium of a macroscopic ensemble of S = 3 chromium atoms pinned in a three dimensional optical lattice and prepared in a pure coherent spin state, under the effect of magnetic dipole–dipole interactions. Our isolated system thermalizes under its own dynamics, reaching a steady state consistent with a thermal ensemble with a temperature dictated from the system’s energy. The build up of quantum correlations during the dynamics is supported by comparison with an improved numerical quantum phase-space method. Our observations are consistent with a scenario of quantum thermalization linked to the growth of entanglement entropy.

U ltra-cold atomic systems featuring long-range interactions are becoming ideal platforms for probing strongly correlated out-of-equilibrium quantum behavior and, in particular, the phenomenon of quantum magnetism, where magnetic moments with quantized energy levels (spins) interact with one another [1][2][3] . Their appeal stems from the fact that they feature internal levels that can be initialized in pure states and coherently evolved with controllable long-range interactions even under frozen conditions. Recently great advances have been accomplished, but, so far have been mostly limited to small systems (hundreds or fewer particles) [4][5][6][7][8][9][10][11][12] , or to dilute disordered molecular ensembles 13,14 .
Magnetic quantum dipoles featuring sizable magnetic moments offer unique opportunities since magnetic interactions can directly happen in an enlarged set of low-lying hyperfine Zeeman levels and are not forbidden by parity and time-reversal symmetry as is the case with electric dipoles 15 . They offer untapped opportunities as a quantum resource since S > 1/2 spin models have more complexity and cost exponentially more resources to classically simulate 16,17 . In fact, the exploration of the complex non-equilibrium dynamics of dipolar-coupled S > 1/ 2 spin models remains a fascinating territory which only starts to be explored 18,19 .
Here we compare our experimental observations of the seven Zeeman populations of an initial spin coherent state made of S = 3 spin particles, with different models. We find that our data compare well with exact short time calculations and semiclassical simulations based on a discrete Monte Carlo sampling in phase space [20][21][22] . We show that the steady state reached at long times is captured by a statistical ensemble with nonzero thermodynamic entropy, by deriving a simple analytical formula, which compares well to both data and semiclassical simulations. We finally study the growth of entanglement by computing the Renyi entropy associated with a local spin. Our studies confirm a scenario of quantum thermalization as a result of the entanglement accumulated during the dynamics.

Results
Realization of an XXZ Heisenberg spin model. In our system the spin degree of freedom is encoded in the Zeeman levels of the purely electronic S = 3 ground state of 52 Cr atoms. The experiment starts with the production of a spin-3 Bose-Einstein condensate (BEC) of~4 × 10 4 atoms in the m S = −3 state, following the procedure described in ref. 23 . We then adiabatically load the BEC into a three dimensional (3D) optical lattice made by laser beams at 532 nm 18 . The lattice structure is rectangular in the horizontal plane, and uses a standard retro-reflecting scheme on the vertical axis (see Methods). After loading the atoms into a deep optical lattice, the sample forms a Mott insulator consisting of a core with doubly occupied sites n ¼ 2 ð Þ, surrounded by a 3D shell of singly occupied sites n ¼ 1 ð Þ, see Fig. 1. The experimental procedure to induce spin dynamics is shown in Fig. 1. We initialize the system in a well-characterized state consisting of a macroscopic array of long-lived singly occupied sites close to unit filling, by performing first a filtering protocol. It relies on dipolar relaxation 18 to empty all doubly occupied sites within the n ¼ 2 Mott core after the application of a π rf pulse that promotes the atoms to the most energetic spin state m S = 3.
The filtering protocol takes about 7 ms (see Methods). To trigger the spin dynamics we then apply a second rf pulse. This rotates the coherent spin state, such that it forms an angle θ with respect to the external magnetic field which sets the quantization axis (see Fig. 1). This prepares a tilted spin coherent state. The spin dynamics is studied by monitoring the time evolution of the population of the different Zeeman states, using absorption imaging after a Stern-Gerlach separation procedure 23 .
A unit-filled array of frozen magnetic dipoles in a lattice interact via dipolar interactions. In the presence of an external magnetic B field strong enough to generate Zeeman splittings larger than nearest-neighbor dipolar interactions, only those processes that conserve the total magnetization are energetically allowed and the dynamics is described by the following secular Hamiltonian 18 (with B along the z axis): where the sum runs over all pairs of atoms (i,j  iŜ À j þ h:cÞ terms. c A first π pulse is used to promote all atoms to the most excited spin state. Dipolar relaxation empties doubly occupied sites. Once this is achieved, a second pulse collectively rotates all spins by an angle θ from the external magnetic field B which sets the quantization direction. We then study spin dynamics due to intersite dipole-dipole interactions by registering the relative populations of the different Zeeman states of vacuum, g ' 2 is the Landé factor, and μ B the Bohr magneton. r (i,j) is the distance between atoms, and ϕ (i,j) the angle between their inter-atomic axis and the external magnetic field. The Hamiltonian is given in terms of spin-3 angular momentum operators, b S i ¼ fŜ x i ;Ŝ y i ;Ŝ z i g, associated to atom i. An important feature is that the dynamical redistribution of populations can happen for large spins (S > 1/2), even though both the total particle number N and the collective magnetization M ¼ hŜ z i are conserved quantities (withŜ x;y;z ¼ P N j¼1Ŝ x;y;z i ). The magnetic dipolar interaction energy between S = 3 spins is 36 times larger than the one for S = 1/2 alkali atoms, allowing us to probe such population dynamics at milliseconds time scales, as seen in Fig. 2.
Perturbation theory. We will first introduce the expected basic dynamical features according to time-dependent perturbation theory. We will focus on the main differences when assuming a classical behavior, or when taking into account quantum correlations. The simplest possible picture for the population dynamics relies on a mean-field treatment (i.e., neglecting quantum correlations), where each atom undergoes Larmor precession around an effective dipolar field created by all the other spins, 2 fhŜ x j i; hŜ y j i; À2hŜ z j ig. Time-dependent perturbation theory yields the following equation for p m S , the relative population of Zeeman level m S = −3, …, 3: We give in the Methods section exact formulas for α m S ðθÞ. For instance, α m S ¼fÀ3;À2;À1;0;1;2;3g ðπ=2Þ ¼ 135=512 f1; 2; À1; À4; À1; 2; 1g.
For a homogeneous gas, B dih vanishes, and the population dynamics with it. This behavior remains valid at all times given that, by preparation, all spins point along the same direction initially, they precess around the same classical dipolar field, and thus evolve identically. Therefore, the local magnetization hS z i i ¼ M=N remains constant for each spin, cancelling population dynamics altogether. On the other hand, in a trapped gas, the inhomogeneous dipolar field introduces a differential precession rate between spins, which results in population dynamics. Note that B dih is determined by border effects and also that K d ðtÞ itself is time dependent (∝t 2 ). Therefore, classically, population redistribution is a slow t 4 process. We emphasize that B dih is proportional to cos θ and thus vanishes when θ = π/2 where no mean-field dynamics takes place at all.
Quantum fluctuations can drastically modify this behavior and induce much faster population dynamics even for a homogeneous gas 24 . Second-order time-dependent perturbation theory on the exact Hamiltonian 25,26 in Eq. (1) yields: In contrast to the mean-field case, the dynamics grows as t 2 and is driven by V eff ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi P N i;j≠i V 2 ij =N q . We emphasize the relatively fast decay of V 2 eff with interparticle distance r (as r −6 ), which makes the short time evolution mainly determined by the nearest-neighbor interactions. As V eff is independent of θ the tipping angle θ provides a way to study an out-of-equilibrium magnetism increasingly determined by quantum correlations as θ → π/2.
In the experiment, external systematics such as quadratic Zeeman fields, B Q , generated by tensorial light shifts induced by the lattice lasers-with eigenenergies B Q m 2 S -or inhomogeneities associated with magnetic field gradients, Δ ij = B i − B j , need to be accounted for. Their role in the short time dynamics can be understood using perturbation theory (see Methods  K d ðtÞ À 4=3Q 2 and V 2 eff ! V 2 eff À 4=3Q 2 in the classical and quantum cases, respectively, with Q 2 1 Thus dipolar inhomogeneities and magnetic field gradients are in direct competition. In the quantum case magnetic field gradients also enter as t 4 but in this case they play a subdominant role since the leading dipolar dynamics is significantly faster (∝t 2 ). Numerical methods to model the spin dynamics. Although perturbation theory allows emphasizing some of the main qualitative differences in the classical and in the quantum regime, to accurately describe the population dynamics we need to go beyond perturbation theory. To accomplish that we parameterize each spin i by a generalized Bloch vector,λ ½i . In contrast to spin-1/2 systems, this vector is a 48-dimensional object that determines all independent elements of the 7 × 7 (= (2S + 1) 2 ) individual spin-3 density matrices,ρ i ðλ ½i Þ 27 . Inserting the product-state ansatz of the system density matrix,ρ ¼ Q N i¼1ρi ðλ ½i Þ, into the von-Neumann equation, dρ=dt ¼ ðÀi= hÞ½Ĥ;ρ, yields N × 48 independent non-linear mean-field equations, in which each generalized Bloch vector evolves in the field of the others. The mean-field "classical" results are obtained by numerically integrating these equations of motion (see Methods).
To capture the build up of quantum correlations we developed a generalization of a semiclassical method (generalized discrete truncated Wigner approximation, GDTWA) based on a discrete Monte Carlo sampling in phase space originally derived in the framework of the so-called truncated Wigner approximation (TWA) 28 . It describes the initial state in terms of a probability distribution. Initial spin coherent states are ideal since they can be fully described by a positive discrete probability distribution. For spin-1/2 systems, randomly sampling this initial "Wigner function", leads to the discrete truncated Wigner approximation (DTWA) [20][21][22] , an approximation that has been remarkably successful and can capture complex quantum aspects of spin dynamics. In contrast to the spin-1/2 case, here (GDTWA) the discrete probabilities are not provided by the eigenvalues of the three Pauli matrices, but instead by the eigenvalues of the corresponding 48 generalized SU(7) generators (see Methods). This semiclassical approach is benchmarked by comparison with exact diagonalization predictions (see Supplementary Note 1 and Supplementary Fig. 1).

Comparisons between experiment and numerical simulations.
We now describe how our data compare with simulations for different values of θ. In Fig. 2, we show our data and the comparisons to both the classical and the GDTWA models. The theoretical models take into account the 3D lattice structure and the measured magnetic field gradients along all three directions. We also include the weak quadratic Zeeman field present in the experiment. Since we could not measure it directly we allow it to be a fitting parameter (see Supplementary Note 2 and Supplementary Fig. 2). For each of the four tilting angles used for the measurements we plot the evolution of fractional populations in different Zeeman states. We only plot the most relevant Zeeman states (the most populated for most of the angles, and only the negative Zeeman states for the symmetric case π/2-see Supplementary Fig. 3 for extensive data). Experimentally, we find that the amplitude of spin dynamics (i.e., the amplitude of the variations of the populations in the different spin states) is stronger when the angle increases.
At small angles the experimental data are qualitatively reproduced by both classical and GDTWA simulations. As can be seen in Fig. 2, both simulations then yield similar results, but nevertheless show systematic differences. This shows that even at the smallest angles that we have probed, beyond mean-field effects are in principle already at play. However, given the signalto-noise ratio in the experimental data, it is difficult to quantify the contribution of beyond mean-field effects to the spin dynamics at weak rotations. When increasing the angle, it becomes increasingly clear that only the beyond-mean-field simulation accounts for the observed dynamics, both at short times (t < 20 ms, see Fig. 2) and at long times (t > 40 ms, see Fig. 3a).
We have performed a systematic study in our numerical simulations, by varying the size of the system which we use. We find that a good agreement between experiment and beyond mean-field theory is only reached, provided the number of interacting spins in the simulation is larger than about 60 (see Supplementary Note 4 and Supplementary Fig. 4). Taken together, these data show that spin dynamics after the initial quench is inherently many-body, and beyond the grasp of meanfield models. As can also be seen in Fig. 2, our experimental data at short times are also in excellent agreement with the exact dynamics calculated within the framework of second-order time-dependent perturbation theory (see Eq. (3)). Also in good agreement with this equation, we find the dynamics at short times to be roughly independent of the magnetic field gradient applied to the sample (up to values >30 MHz m −1 ; see Supplementary Note 5 and Supplementary Fig. 5). In contrast, we point out that the experimental data at short times systematically show faster dynamics than predicted in the classical picture, whose initial t 4 dependence (see Eq. (2)) fails to reproduce the experimental observations. Finally, we checked that the effect of imperfections in the lattice preparation on GDTWA predictions is small (see Supplementary Note 6 and Supplementary Fig. 6).
Models for thermalization at long times. For an isolated system, entanglement build up after a quench into a nonequilibrium situation is tied to the scenario of quantum thermalization. To support the relevance of quantum correlations during dynamics, we thus analyze the long-time behavior of the populations. For all tipping angles θ, we observe that the experimental system approaches a steady state, which is in agreement with predictions of closed system quantum thermalization, given, e.g., by the eigenstate thermalization hypothesis (ETH) 29,30 . In particular, we find that the long-time average populations are very well described by the effective thermal distributionρ cT ðβ; μÞ ¼ e ÀβĤ T ÀμŜ z tr½e ÀβĤ T ÀμŜ z where the chemical potential μ and inverse temperature β = 1/k B T are set by the energy and magnetization of the initial pure state: which are conserved throughout the evolution. HereĤ T 1 H þ P i B Q ðŜ z i Þ 2 is the total Hamiltonian. As shown in Fig. 2, the steady-state populations approach the ones (indicated by the arrows for all tipping angles) dictated by the thermal ensemble when simply setting β = 0 (in which case the maximum-entropy state only depends on magnetization). For angles close to π/2, however, where quantum effects are most significant, we find a deviation compared with this simplistic prediction. We therefore proceed to study this interesting regime.  Fig. 3 we show dynamics up to longer times, and confirm that a steady state is indeed reached after 40 ms for the π/2 case, a feature that the GDTWA simulation reproduces, while classical simulations predict an oscillatory behavior for a much longer duration. This qualitative difference between the classical and the quantum behavior is associated with the different origin of thermalization in both pictures: while quantum-mechanically thermalization is tied to the growth of entanglement, classically, reaching a steady state in a system of frozen particles is a consequence of the single-particle dephasing induced by field inhomogeneities. The inhomogeneous fields arise either from external fields, or from effective fields generated on one particle from the mean-field interactions with the surrounding particles. This behavior differs from the typical thermalization scenario in mobile particles where collisions can classically change both motional and internal degrees of freedom while redistributing energies and momenta 29,31-33 . Our observations clearly rule out a simple mean-field (classical) behavior. Most interestingly, we find a very good agreement between the experimental data points taken at long times, i.e. after the system has reached its steady state, if instead of setting β = 0 we account for the corrections generated by the quadratic Zeeman field B Q , and the finite but small energy of the initial state in Eq. (4). Using a simple perturbative approach (see Methods) we obtain that the β (0) = μ (0) = 0 solutions should be replaced by leading to: with V 1=N P i>j V ij . We find V=h % À0:57 Hz and V eff =h % 6:13 Hz for our lattice geometry. As V<0, negative temperatures are expected for low enough B Q (as allowed for a system whose maximum energy is bounded). Figure 3b shows a very good agreement between the equilibrium data and the analytical model (see Supplementary Note 7 and Supplementary Fig. 8 for extensive equilibrium data). For this comparison, there is no free parameter, since we use the value of B Q for which the dynamical evolution of the spins is best reproduced by GDTWA simulations. This good agreement confirms the scenario that the coherent Hamiltonian evolution of the many-body system drives it toward a strongly entangled pure state for which the observables display thermal-like behavior. The agreement between the analytical model and the GDTWA at long times shown in Fig. 3b also indicates that the GDTWA not only captures the short term dynamics (as previously known from the theoretical point of view), but also the approach to equilibrium.
To compare the analytical formula to the data we have ignored magnetic field gradients in Eq. (6) (See Methods). In principle, magnetic field gradients should lead to an equilibrium state where a spatial texture of magnetization develops. However, such a texture requires long-range interactions between remote parts of the cloud, which only occurs for an extremely long timescale for dipolar interactions. We have verified (see Supplementary Note 7 and Supplementary Fig. 7) that indeed magnetic field gradients can be neglected to evaluate the quasi-steady state populations reached at 100 ms. This shows that a local equilibrium is first reached, well before the full many-body system may reach true equilibrium with maximum entropy, where all populations in the different Zeeman states would be equally populated.
Study of quantum correlations. To quantify the importance of quantum correlations in the spin dynamics as a function of the tilting angle θ, we analyze from the theoretical point of view the properties of the reduced density matrix for each spin,ρ i ðλ ½i Þ. In our simulations those density matrices are readily available from the generalized Bloch vectors. To minimize finite size and boundary effects we focus on the density matrix of the central spin of our simulated blockρ 0 . Even when, as in our simulations, the quantum state of the full systemρ is pure, the reduced singlespin density matrices can assume a mixed character due to the buildup of entanglement between the spins. This mixed character is quantified by a reduced purity, trðρ 2 0 Þ < 1 and thus an increased entropy, which we compute in terms of the second-order Renyi entropy, S ð2Þ 0 ¼ Àlog 2 ½trðρ 2 0 Þ. If the state of the full system is pure, trðρ 2 Þ ¼ 1, the Renyi entropy is a measure of entanglement 34 : it is zero for product states, and reaches the maximum value of S ð2Þmax 0 ¼ log 2 ½7 (the value for a fully mixed state of a spin-3 particle) for many-body states where the quantum information encoded in an individual spin is completely scrambled due to entanglement with other ones.
As can be seen in Fig. 4, the quantum evolution leads to a growth of S ð2Þ 0 already for the smallest investigated angle θ = 0.2π. The dynamical growth of entanglement increases significantly for larger tilt angles. At θ = π/2, we find that S . Although we cannot perform a full state-tomography from the experimental data, we can compare our experimental data to a "diagonal entropy" computed in terms of the diagonal part of the averaged single-particle density matrixρ S ¼ ð1=NÞ P N i¼1ρi ðλ ½i Þ. Note that for an homogeneous systemρ S ¼ρ 0 . In this case we can define this entropy as S D 2 ¼ Àlog 2 ftr½diagðρ S Þ 2 g, which can be readily accessed from the population data, assuming homogeneity: S D 2 ¼ Àlog 2 f P p 2 m S g. This diagonal entropy is not an entanglement witness, but it provides an upper bound of the entanglement entropy, S D 2 ! S ð2Þ 0 . In a translationally invariant system it increases as quantum correlations build up with time and approaches the full entropy as the single-spin density matrices decohere due to entanglement. However, in our finite system, boundary effects can obscure this behavior. For example, at small angles, diagonal entropy shows a slight reduction as a function of time (see Fig. 4 for θ ≤ 0.3π). On the other hand, we do observe it increases with time as the system thermalizes for large values of θ, in which case the non-trivial growth of the experimental diagonal entropy is in excellent agreement with our theoretical estimates, provided quantum fluctuations are taken into account. Moreover it also approaches S ð2Þ 0 for θ = π/2.

Discussion
In summary, our study demonstrates the dominant role of quantum correlations in the out-of-equilibrium dynamics of an initially uncorrelated spin coherent state, when the angle it makes with the external magnetic field is close to π/2. We have shown that our long-range interacting many-particle isolated spin system internally thermalizes through entanglement buildup, and develops an effective thermal-like behavior through a mechanism which is purely quantum and conservative. The comparison between experiment and theory shows that the GDTWA simulations can be trusted for studying the dynamics in a complex quantum many-body system, provided a sufficient number of atoms is included in the simulation. Thus, our experiment provides a test-bed for a theoretical method based on the GDTWA, for systems of large spins, and in a many-body regime where simulations based on exact diagonalization techniques are intractable with current computational resources. In turn, our study can be used as a benchmark of a quantum simulator of the spin-3 XXZ Heisenberg model and opens a path toward the study of open problems in quantum many-body physics. For example, by operating the experiment at smaller lattice depths, where tunneling is allowed, we will have the exciting opportunity to study itinerant magnetism, whose description is typically unaccessible to theory, but which is believed to be at the heart of the physics behind high-temperature superconductivity 35 .

Methods
Description of the 3D lattice. The 3D lattice is made with five laser beams at 532 nm. On the horizontal plane, three beams with the same frequency define a rectangular pattern, with respective directions u H 1 ¼ cosðαÞu x À sinðαÞu y , . Two other beams, contra-propagating, with a frequency offset by 30 MHz compared with the beams in the horizontal plane, with directions u V 1 ¼ ÀcosðβÞu z þ sinðβÞu x ¼ Àu V 2 , β = 7/180π, form an independent light pattern. Calibration of the lattice is performed by standard matter wave diffraction pattern analysis after pulsing lattice beams onto the BEC, with the three pairs of beams (H 1 , H 2 ), (H 1 , H 3 ), and (V 1 , V 2 ). The laser powers are chosen so that these three couples of beams induce almost equal lattice depths, larger than 25 recoil energy. For these lattice depths, the tunneling time is typically 100 ms, and tunneling events can safely be neglected during dynamics.
Preparation of a lattice with only singly occupied sites. To prepare a lattice of atoms at unit filling, we first slowly load the BEC into a 3D optical lattice, to reach a Mott-insulating state. For our experimental parameters, there exists a core with only doubly occupied sites, surrounded by a 3D shell of atoms at unit filling. We empty the doubly occupied sites by performing a rf pulse to promote all atoms from the lowest energy Zeeman state m s = −3 into the state m s = 3, which triggers dipolar relaxation. We perform our experiment in presence of an external magnetic field which is large enough that dipolar relaxation can be considered as a shortrange process 36 . Thus, only atoms in doubly occupied sites undergo dipolar relaxation, and each dipolar relaxation event empties one doubly occupied lattice site. We estimate the probability of secondary collisions during this filtering procedure to be below 0.05. After 7 ms, all doubly occupied sites are empty, with about 10,000 remaining atoms.
The spin dynamics experiment is then performed using the atoms remaining in the shell with unit occupancy. Because the sample during dynamics consists of a 3D shell of atoms with unit occupancy within the lattice, border effects might not be fully negligible during dynamics. Indeed, our estimates is that about 20 percent of the atoms within the shell of singly occupied sites are close to the boundary. It is likely that spin dynamics is slower for these atoms lying close to the frontier of the shell.
Note that the experiment could not be performed at arbitrarily high magnetic field intensities. As a consequence, some of the atoms which underwent dipolar relaxation remain trapped in very highly excited states of the combined latticedipole trap potentials. This translates into losses affecting the sample with unit filling. After 40 ms, from 20 to up to 40 percent of the atoms are typically missing, depending on the magnetic field strength. This phenomenon does not seem to impact the agreement of our spin dynamics data with GDTWA theory as long as losses are below 30 percent.
Atom number calibration. The number of atoms in different spin states is estimated using standard absorption imaging, after spin separation using an applied magnetic field gradient during the free fall of atoms, following a Stern-Gerlach procedure. The cross section for absorption of resonant light strongly depends on the m s states, through Clebsch-Gordan coefficients. Therefore, we calibrate the relative sensitivity of the imaging system for the different spin states by comparing the measured populations just after the rf pulse to the theoretically expected values. This calibration depends on the external magnetic field direction during spin dynamics, as eddy currents do not allow to rapidly set its direction during imaging.
For the specific case of θ = π/2, we employ a slightly different method to calibrate the different sensitivities. Indeed, the number of atoms in m s = +3 is then very small just after the rf pulse and the detectivity of this Zeeman state is the lowest, due to unfavorable Clebsch-Gordan coefficients. For this specific dataset, we thus enforce that the m s = −3 and m s = 3 average atom number after spin dynamics is identical. This choice is motivated by the fact that the Hamiltonian preserves magnetization (as experimentally verified for all other datasets), and by the initially symmetric theoretical populations in the different Zeeman states. For example, for the π/2 data in Fig. 2  Short-time analysis of population dynamics. Using time-dependent perturbation theory we analyze the contribution of the different terms in the Hamiltonian at short times.
For our system the initial population is given by Generalized Bloch vectors and the GDTWA. A generic density matrix for a discrete system with D states on site i takes the formρ i ¼ P D α¼1;β¼1 c α;β jαihβj. For a spin-3 atom D = 7, and to the states α ¼ 1; 2; 3; ; 6; 7 j i we may associate the spin states m S ¼ 3; 2 ; À2; À3 j i . Since ðρ i Þ y ¼ρ i and trðρ i Þ ¼ 1 a total of D 2 − 1 real numbers are needed to describe an arbitrary state. Those numbers can be expressed as expectation values of D 2 − 1 orthogonal observables:Λ The mean-field equations can be written as (D 2 − 1) × N coupled non-linear equations for the expectation values of λ ½i μ ¼ hΛ ½i μ i. The λ ½i μ can be interpreted as components of a D 2 − 1 dimensional Bloch vector via the expansion Here the tensor f μ,κ,η is defined via ½Λ ½i μ ;Λ ½i κ ¼ if μ;κ;ηΛ ½i η , whose elements are the structure constants of the SU(D) group. The full mean-field equations for the generalized Bloch vector at site i are κ is the expansion of the single-site Hamiltonians containing all local terms (field gradients, quadratic Zeeman fields, etc.) into GGMs. It is straightforward to construct the equations for arbitrary Hamiltonians containing single-and two-site terms numerically, as well as to evolve the generalized Bloch vectors in time.
In the numerical mean-field simulations, the quantum state is represented by N time-dependent generalized Bloch vectors, λ ½i μ ðtÞ. We evolve the vectors for the initial state In contrast, in the GDTWA approach we describe the initial state not by a generalized Bloch vector, but instead by a probability "Wigner" distribution, p ½i μ;a μ , for certain discrete configurations of Bloch vector elements, λ ½i μ;a μ . Initially, the probabilities and configurations are chosen in such a way that on average λ ½i μ ðt ¼ 0Þ ¼ P a μ p ½i μ;a μ λ ½i μ;a μ λ ½i μ;a μ . In practice, the initial multi-spin configurations are selected via a random sampling of p In particular, as discrete set of initial configurations, fλ ½i μ;a μ g, we use a set which is inspired by a "projective measurement, of the GGMs": for each λ ½i μ ðt ¼ 0Þ, we choose a set of initial configurations given by the eigenvalues of each GGM. Quantum thermalization. It is generally believed that the unitary quantum evolution of a complex quantum system leads to an apparent maximum-entropy state that can be described by thermodynamical ensembles that properly account for the conserved quantities. In our systems those are the energy and magnetization. We thus postulate that the steady-state properties of local observables, such as the relative population of Zeeman levels, can be described in our system by the thermal distributionρ cT ðβ; μÞ ¼ e ÀβĤ T ÀμŜ z tr½e ÀβĤ T ÀμŜ z where μ and β = 1/k B T are the chemical potential and inverse temperature set by the energy and magnetization, respectively, accordingly to Eq. (4). While the determination of β and μ can be a challenging task for a complex many-body system, the anisotropic character of the dipolar interactions facilitates an analytic high-temperature expansion around β = 0 since V is small (see main text). Under this assumption, the chemical potential to leading order, is set by hŜ z i ¼ tr½ρ cT ð0; μ ð0Þ ÞŜ z ¼  Fig. 2. The case θ = π/2 is particularly simple since hŜ z i ¼ 0 and thus μ ð0Þ ¼ 0 and p ð0Þ m S ðt SS Þ ¼ 1=7. This solution, however, shows deviations with the observed long-time dynamics indicating that finite β corrections are relevant. To first order in β (1) the chemical potential can be written as μ ð1Þ ¼ μ ð0Þ þ β ð1Þ δν and the solutions of Eq. (4) are described by the relations: A. Solutions of those equations, taking into account the intial magnetization and total energy in the intial product state, are particularly simple for the θ = π/2 case, where μ (0) = 0, ê S z ¼ 0, Δ ĝ H TŜ z ¼ 0 and Δ ĝ S zŜz ¼ NI 2 , f H T ¼ NB Q I 2 and Δ ĝ H TĤT ¼ NðB 2 Q ðI 4 À I 2 2 Þ þ 3=4V 2 eff I 2 2 Þ with I r ¼ ð P 3 m¼À3 m r Þ=7 (thus I 2 = 4 and I 4 = 28). Those yield the expressions for β (1) and μ (1) quoted in Eq. (5). In the presence of linear gradientsĤ T !Ĥ T þ P N i¼1 B iŜ z i , under the assumption that PN i¼1 B i ¼ 0, the inverse temperature equation in Eq. (6) for the case of θ = π/2 should be replaced by β ð1Þ ¼

Data availability
The experimental data supporting the findings of this study are available within the paper and its Supplementary Material. Additional numerical data and computer codes used in this study are available from the corresponding author upon request.