Quantum critical dynamics in a spinor Hubbard model quantum simulator

Three-dimensional (3D) strongly correlated many-body systems, especially their dynamics across quantum phase transitions, are prohibitively difficult to be numerically simulated. We experimentally demonstrate that such complex many-body dynamics can be efficiently studied in a 3D spinor Bose–Hubbard model quantum simulator, consisting of antiferromagnetic spinor Bose–Einstein condensates confined in cubic optical lattices. We find dynamics and scaling effects beyond the scope of existing theories at superfluid–insulator quantum phase transitions, and highlight spin populations as a good observable to probe the quantum critical dynamics. Our data indicate that the scaling exponents are independent of the nature of the quantum phase transitions. We also conduct numerical simulations in lower dimensions using time-dependent Gutzwiller approximations, which qualitatively describe our observations. Three-dimensional (3D) strongly correlated many-body systems and their dynamics across quantum phase transitions pose a challenge when it comes to numerical simulations. The authors experimentally demonstrated that such many-body dynamics can be efficiently studied in a 3D spinor Bose–Hubbard model quantum simulator, and observed dynamics and scaling effects beyond the scope of existing theories at superfluid–insulator quantum phase transitions.

P hysical modeling of three-dimensional (3D) strongly correlated many-body systems, especially their critical dynamics across quantum phase transitions, is a fastmoving research frontier with immediate applications spanning from adiabatic quantum state preparations to developing novel materials [1][2][3][4][5][6][7][8][9] . Dynamics of phase transitions have been a rich but challenging research theme. One paradigm of dynamical studies is the Kibble-Zurek mechanism (KZM), which characterizes dynamics into two regimes: the impulse regime, where a system is frozen due to diminished energy gaps and slowed equilibrations in the vicinity of phase transition points, and the adiabatic regime, where a system with considerable energy gaps progresses adiabatically [9][10][11][12][13] . Impeded by the intrinsic complexity of the many-body systems and limitations of existing numerical techniques, theoretical studies of these systems have been performed mainly in homogeneous systems and lower dimensions [1][2][3][14][15][16][17][18][19][20] . It has been suggested, however, that such complicated many-body dynamics could be efficiently studied in 3D spinor Bose-Hubbard (BH) model quantum simulators, consisting of antiferromagnetic spinor Bose-Einstein condensates (BECs) in cubic optical lattices 3,4,6,21 . Possessing spin degrees of freedom and exhibiting magnetic order and superfluidity, these 3D spinor quantum simulators are well-isolated quantum systems that are highly programmable with a remarkable degree of control over many parameters, such as temperature, spin, density, and dimensionality [3][4][5][6][21][22][23][24] . When lattice-trapped bosons are suddenly quenched from the superfluid (SF) phase to the Mottinsulator (MI) phase, the BH model is nearly integrable, yielding nonthermal steady states preserving memory of the initial state 9,14 . Another distinctive and useful feature of these 3D quantum simulators is their ability to tune the nature of SF-MI quantum phase transitions, i.e., spinor gases can cross first-order (second-order) SF-MI transitions when the quadratic Zeeman energy q is smaller (larger) than the spin-dependent interaction energy U 2 in even Mott lobes (see Supplementary Note 2) 4,5,20-23 . U 2 also enables spin-mixing oscillations among multiple spin components and generates ferromagnetic and antiferromagnetic order along with various magnetic states including spin singlets and nematic states 4,[21][22][23][24][25] . These 3D quantum simulators have thus been suggested as an ideal platform for studying quantum coherence, long-range order, magnetism, quantum phase transitions, and symmetry breaking [4][5][6]21,22 .
In this paper, we experimentally demonstrate the strength of 3D spinor BH model quantum simulators and explore the quantum critical dynamics by revealing asymmetric impulse-adiabatic regimes 16,26 as spinor gases are quenched across SF to spinor MI quantum phase transitions at different speeds. Our observations, going beyond the adiabatic-impulse-adiabatic regimes studied in many publications [9][10][11][12][13]16,27,28 , confirm the impulse regime is enlarged to cover the entire SF phase and part of the MI phase due to negligible excitations within the SF phase 16,26 . We find quench dynamics and scaling effects at SF-MI quantum phase transitions, which are outside the framework of existing theories including KZM. The observed power law scaling exponents are independent of the nature of SF-MI transitions. We also compare the dynamics observed using both spin populations and SF order parameters as observables and discuss the benefits and limitations of each. While numerical simulations in 3D multicomponent systems are prohibitively difficult to perform for further verification of our results, we conduct theoretical simulations in 2D and find qualitative agreements with the experiment.

Results and discussion
Each experimental cycle is started with a spin-1 (F = 1) sodium spinor BEC at its SF ground state, the longitudinal polar state where ρ 0 = 1 and m = 0, at a given q. Here ρ m F is the fractional population of atoms in the F ¼ 1; m F j istate, m F is the magnetic quantum number, and m = ρ −1 + ρ +1 is the magnetization. The investigated q values are of the same order of magnitude as the spin-dependent interaction U 2 . Spinor gases are characterized as ferromagnetic (antiferromagnetic) when U 2 is negative (positive) 4,23 . Sodium spinor BECs with U 2 ≃ 0.035U 0 > 0 studied in this paper are a good example of antiferromagnetic spinor gases (See Supplementary Note 1) 6 . Here, U 0 is the spin-independent interaction. We load spinor BECs into a cubic lattice by linearly quenching the lattice depth u L , which results in an exponential increase in the ratio of U 0 to the hopping energy J. Similar to scalar bosons, spinor gases can cross SF-MI transitions when U 0 /J becomes larger than critical values 5,6,20,21,[21][22][23] . We describe the lattice-confined spinor gases with the spin-1 BH Hamiltonian and compute their static and dynamic properties using the Gutzwiller approximation (see "Methods"). Note that the linear Zeeman energy, which remains constant during collisional spin interconversions, can be neglected in the Hamiltonian 4,21,23 .
Quantum critical dynamics probed via spin populations. Figure 1 compares density profiles of various spin components, derived from numerical simulations in systems similar to the experimental system but in two dimensions (2D), after the lattice depth u L is quenched from 2E R to 40E R at different speeds v ramp . Here, the unit of v ramp is E R per millisecond (E R ⋅ ms −1 ) and E R is the recoil energy 5,21 . For slow quenches, Fig. 1a-c, g clearly display the formation of Mott insulating shells and spin structures. Whereas for fast quenches, Fig. 1d-f, h shows persistence of SF type density profiles. Lattice gases studied in this paper are inhomogeneous systems confined by an overall harmonic trap. In these systems, the BH model predicts that SF and MI regions coexist and many Mott lobes are arranged in a wedding-cake structure (see Fig. 1 and Supplementary Note 2). Such inhomogeneous systems have been suggested as good candidates for adiabatic quantum state preparations, because inhomogeneous phase transitions suppress excitations during quenches 9 .
Our experimental data agree with these predictions, as shown by spin dynamic curves in Fig. 2a. During fast quenches, the observed spin-0 population ρ 0 implies that the initial SF state is frozen as the lattice is quenched across predicted SF-MI transition points (shaded gray area in Fig. 2a) and then starts to evolve inside isolated lattice sites in deep lattices. Whereas in sufficiently slow quenches, our data indicate that atoms stay in the instantaneous ground state during the dynamics and enter the MI phase at a critical U 0 /J. For intermediate quench speeds, for example the v ramp = 7E R ⋅ ms −1 curve in Fig. 2a, partial revivals in the spin-0 population may be observed after ρ 0 reaches its minimum value. Such revivals are not found at slow and fast quench speeds (see Fig. 2a). Qualitatively similar results are obtained from our 2D numerical simulations (see Fig. 2b), although it is harder to realize adiabatic phase transitions in lower dimensional systems where more low-energy states are accessible for excitations 9 . Theoretical calculations also predict that the frequency of the periodic revivals is determined by U 2 and q for a fixed atom number n, while the amplitude of the revivals depends on the quench speed v ramp . The physical limitations of our system, however, prevent us from safely observing full periods of the revivals at intermediate v ramp , as well as the minimum ρ 0 in fast quenches (v ramp ≳ 28E R ⋅ ms −1 ), due to extremely high power requirements on the lattice beams.
We define the length of the impulse regime as the length of the frozen region, where atoms cannot remain in instantaneous ground states and thus stay frozen in their initial states during quenches, and apply a single ρ 0 cut-off value to mark the edge of the impulse regime. For example, Fig. 2c is extracted from the observed dynamics by setting the end of the impulse regime at the lattice depth where ρ 0 = 0.9. When U 0 /J at this lattice depth is plotted against v ramp , two distinct dynamical regions emerge (see Fig. 2c). For relatively slow quenches (v ramp ≲ 13E R ⋅ ms −1 ), the length of the impulse regime does not strongly depend on quench speed. For faster quenches, the length obeys a power law relationship with v ramp as expected from a quantum version of the KZM, although the extracted scaling exponent 1.55(8) is much larger than the 0.67 predicted by a simple KZM for latticetrapped 3D scalar bosons 9 . We find that the 2D simulations (dashed line in Fig. 2c) and our experimental data are numerically very similar when v ramp > 18E R ⋅ ms −1 and somewhat similar when v ramp ≤ 6E R ⋅ ms −1 , while displaying large discrepancies in the region near the critical ramp speed (≈13E R ⋅ ms −1 ). These discrepancies may be due to the presence of inhomogeneous systems and difference in dimensionality. The scaling exponents extracted from the experiments and the 2D simulations show no strong dependence on the ρ 0 cut-off value as long as the cutoff is sufficiently larger than the minimum observed ρ 0 values. The quenches studied in this paper are exponential in U 0 /J, whereas the KZM usually applies to quenches linear in system parameters, so our experiments may have revealed nonuniversal scaling relations. Our numerical simulations have taken this feature into account, while omitting any finite-temperature or heating effects.
Quantum critical dynamics probed via SF order parameters. Both ρ 0 and SF order parameters are theoretically good observables for SF-MI phase transitions, due to the similarities in the topology and boundaries of the phase diagrams derived from SF order parameters and ρ 0 when n > 1 (see Fig. 2d and Supplementary Note 2). Solid lines in Fig. 2d indicate that ρ 0 in even Mott lobes changes more drastically than in the odd Mott lobes near the SF-MI transition points while within the even Mott lobes the change in ρ 0 gets smaller at higher n. We further analyze these dynamics by examining the SF order parameter, which is proportional to the side peak fraction ϕ measured in experiments. ϕ is a measure of coherence of the wavefunction and is proportional to the number of atoms at zero momentum in the spin-dependent momentum distributions (see "Methods") 5,29 . The inset of Fig. 3a illustrates how we extract ϕ. Typical evolutions of ϕ for the three spin components during a quench at an intermediate speed are presented in Fig. 3a. For the spin-0 atoms, ϕ increases steadily until it obtains its maximum around 12E R and then decreases as the lattice gets deeper and the interference pattern disappears. In Here, E R is the recoil energy, m F is the magnetic quantum number, and the unit of v ramp is E R per millisecond (E R ⋅ ms −1 ). These density profiles are derived from our numerical simulations for 2D lattice-trapped sodium spinor gases of n peak = 6 at q/h = 42 Hz using the Bose-Hubbard model and the Gutzwiller approximation (see "Methods"). Here, n peak is the peak occupation number per lattice site, q is the quadratic Zeeman energy, and h is the Planck constant. d-f Similar to panels (a-c) but after a fast lattice quench with v ramp = 50E R ⋅ ms −1 . The horizontal axes indicate individual lattice sites along the spatial xaxis and y-axis, respectively. The color scale and vertical axis indicate n, the occupation number per lattice site. g, h 2D cross sections of panels (a-c) and panels (d-f) at the 16th lattice site along the y-axis, respectively.
contrast, the observed ϕ for the spin ± 1 components always stays at zero throughout the evolution, which indicates m F = ±1 atoms are only formed in the MI phase. These observations agree with our theoretical simulations shown in Fig. 3b. By defining a cut-off ϕ at which the fitting curve for ϕ of the spin-0 atoms equals 0.2, two different regions again emerge, as shown in Fig. 3c. These regions are divided at roughly the same critical v ramp (≈13E R ⋅ ms −1 ) as the two regions found in Fig. 2c. For fast (slow) quenches where v ramp is faster (slower) than 13E R ⋅ ms −1 , Fig. 3c shows that the extracted U 0 /J linearly increases (decreases) with v ramp , instead of the power law dependence found in Fig. 2c. This apparent discrepancy together with large disagreements between the numerical simulations and the data (see Fig. 3c) may be due to various reasons. First, spin-0 population ρ 0 is a local   Here, q is the quadratic Zeeman energy, h is the Planck constant, E R is the recoil energy, m F is the magnetic quantum number, and the unit of v ramp is E R per millisecond (E R ⋅ ms −1 ). The red (blue) line is a Gaussian (linear) fit. The top horizontal axis lists the corresponding U 0 /J, where U 0 is the spin-independent interaction energy and J is the hopping energy. The inset shows a time-of-flight (TOF) absorption image of the spin-0 atoms and the corresponding density profile of the top, bottom, center peak populations (red line) and of the background population (dashed line) created by integrating along a narrow vertical slice (including only the top, center, and bottom peaks) of the TOF image and by fitting to four Gaussian curves. From these fittings we can extract ϕ by taking the area of the two side peaks and dividing by the total area. The color scale in the TOF image shows the optical density (OD). b Similar to panel a but derived from the Bose-Hubbard model for 2D lattice-trapped sodium spinor gases (see "Methods"). c Markers represent U 0 /J, where the measured ϕ of the spin-0 atoms equals 0.2, extracted from curves similar to the one shown in panel a but taken at various v ramp . The solid line is a two-line fit and the dashed line represents the 2D numerical simulation results. Error bars represent one standard deviation. observable for a system confined in a harmonic trap, while the SF order parameter is a global observable that takes into account twobody correlations of all sites by doing a Fourier transform. For a trapped system, it is thus conceivable that the SF order parameter is more affected by noise, temperature and other nonideal situations. Second, because our two-stage microwave imaging method (see "Methods") is capable of recording both condensed and thermal atoms presented in the optical dipole trap, the measurement of ρ 0 should be much less sensitive to heatings. Nonadiabatic quenches across SF-MI transitions could inherently excite/heat atoms leading to an unavoidable reduction in the visibility of the interference pattern 30,31 , which may set insurmountable technical limitations for using the SF order parameter as an observable to accurately study the nonadiabatic dynamics. On the other hand, the utilization of spin populations to investigate SF-MI transitions is limited to small q, because these transitions are no longer marked by a significant change in ρ 0 when q/h > 120 Hz 5 . Other potential factors that may contribute to the disagreements shown in Fig. 3c include nonequilibrium dynamic and thermalization effects, harmonic trapping effect, and the difference in dimensionality. These coupled with the coexistence of SF and MI regions in inhomogeneous systems make it difficult to compare the impulse regimes derived from SF order parameters across varying quench speeds. Our data indicate that the choice of the cut-off ϕ does not significantly alter the overall observed behavior of the length of the impulse regime, although slopes of the linear fits (solid lines in Fig. 3c) increase as the cutoff value decreases. Another obvious difference from the theoretical predictions for both observables is the observed critical quench speed shown in Fig. 2c and Fig. 3c, which our numerical simulations have not been able to reproduce. The underlying physics of the critical quench speed requires further study, although it could also be due to a dimensional effect, harmonic trapping effect, or experimental peculiarity related to heating effects.
Quantum critical dynamics in different magnetic fields. To test how different local magnetic fields might affect these results, we repeat these experiments at multiple quadratic Zeeman energy q within 0 < q/h < 100 Hz. We find no significant q dependence in the scaling exponent or the critical v ramp , as shown in Fig. 4. Because q ≪ U 0 , the lack of a significant q dependence is consistent with our theoretical understanding of these transitions. From these experiments, we find the scaling exponent of our system to be 1.6(1) and the critical quench speed at which the impulse regime begins to significantly vary with quench speed to be 13(2)E R ⋅ ms −1 . Because SF-MI phase transitions are first order in even Mott lobes at q < U 2 ≈ h × 50 Hz while second order in odd Mott lobes at any q 4,5,20-23 , Fig. 4 implies the scaling exponents are independent of the nature of the SF-MI phase transitions.

Conclusion
In conclusion, we have experimentally demonstrated quench dynamics and scaling effects outside the scope of existing theories in 3D spinor BH model quantum simulators. To explore these dynamics, we have utilized both spin populations and SF order parameters as observables and discuss the technical limitations and benefits of each. Our data indicate the scaling exponents are independent of the nature of SF-MI phase transitions. While it is prohibitively difficult to compute dynamics of multicomponent bosons in 3D lattices, we have performed numerical simulations in lower dimensions and found qualitative agreements with the experiment. This work may indicate the 3D quantum simulators could be more powerful than their classical counterparts, especially in studying many-body quantum dynamics and the interplay of magnetism and superfluidity, although their effectiveness and accuracy is worthy of further study.

Methods
Experimental procedures. Each experimental sequence begins by creating a spin-1 antiferromagnetic spinor BEC containing up to 1.5 × 10 5 sodium ( 23 Na) atoms in a crossed optical dipole trap of frequencies ω x,y,z /2π = 132, 132, 197 Hz at q/h = 42 Hz. Here q is the quadratic Zeeman energy and h is the Planck constant. We then apply resonant microwave and imaging pulses to eliminate the atoms in states. This leaves us with a BEC in its SF ground state, the longitudinal polar state, where ρ 0 = 1 and m = 0. After the pure state is prepared, we quench q to a desired value using magnetic fields whose strength is tuned within 50 and 600 mG, and then load the BEC into a cubic lattice. The cubic lattice is constructed by three orthogonal standing waves producing a lattice with spacing of 532 nm. The lattice beams originate from a single mode laser at 1064 nm and are frequency shifted by at least 20 MHz from each other. The lattice depth u L is calibrated using Kapitza-Dirac diffraction patterns. Once loaded we linearly ramp u L up to the desired value u final L in a given ramp time (t ramp ), which corresponds to a roughly exponential ramp in U 0 /J. Here, J is the hopping energy and U 0 is the spin-independent interaction. The atoms are then abruptly released and ballistically expand, before we measure ρ 0 and the SF order parameter with a two-stage microwave imaging method, in which the m F = 0 atoms are pumped up to the F = 2 state by resonant microwave pulses for imaging in the first stage and then the m F = ± 1 atoms are imaged in the second stage 6,22 . The estimated systematic errors of ρ 0 and the side peak fraction are 0.015 and 0.02, respectively. By holding the lattice ramp speed v ramp ¼ u final L =t ramp constant while varying the u final L and t ramp , we produce graphs similar to those found in Fig. 2a and Fig. 3a. Alternatively by varying the holding time after non-adiabatic lattice ramps, spin dynamics can be observed that reveal information about the Fock state distributions and key parameters determining the spinor physics, as demonstrated in our previous work 6 Fig. 2c as a function of q. Solid lines are linear fits. Here h is the Planck constant, the unit of v ramp is E R per millisecond (E R ⋅ ms −1 ), and E R is the recoil energy. Symbols correspond to values extracted from fittings over at least eight datasets taken at each q and error bars represent one standard deviation. and y ¼ y 0 þ A cosðf ðx À x 0 ÞÞ if x > x 0 . For the piecewise function, y 0 , A, and x 0 are fitting parameters, while f is held at an appropriate value in order to produce more reliable fittings. The fitting function used for each set was determined by the chisquare of the fitted functions and was then used to calculate the critical lattice depth and its associated error.
Theoretical models. The BH Hamiltonian of lattice-trapped F = 1 spinor gases can be expressed as 21 wheren i ¼ ∑ m F a y i;m F a i;m F is the total atom number at site-i, a y m F (a m F ) are boson creation (destruction) operators with m F being ±1 and 0 for spin-1 atoms, the total spin isF ¼ ∑ m F ;m 0 F a y m FF m F ;m 0 F a m 0 a m 0 F , the standard spin-1 matrices areF m F ;m 0 , the spin-dependent interaction is U 2 , the chemical potential is μ, the total number of lattice sites is L, and V T is the harmonic trap strength. In a mean-field approximation, the above Hamiltonian can be expressed as a sum of independent single site Hamiltonians H = ∑ i H i with Here z is the number of nearest neighbors, ψ m F ¼ ha y i;m F i ¼ ha i;m F i is the SF order parameter, and we absorb the V T term in the chemical potential. The side peak fraction studied in Fig. 3 is given by jhâ i;m F ij 2 for any site i and is thus proportional to the SF order parameter.
We compute the static and dynamic properties of the spin-1 BH Hamiltonian using the Gutzwiller approximation 17,20 , which assumes the many-body wave function in the full lattice can be written as a product of local states at L individual sites, Ψ GW j i¼ Q L i¼1 ϕ i with ϕ i ¼ ∑ n i;m F A i ðn i;1 ; n i;0 ; n i;À1 Þ n i;1 ; n i;0 ; n i;À1 . A i is the variational parameter at site i. Dynamics are derived from numerically solving the time-dependent Schrödinger equation for specific initial states. To perform the Gutzwiller calculations, we write the matrix elements of the Hamiltonian H i in the occupation number basis n i;À1 ; n i;0 ; n i;1 , and truncate the onsite Hilbert space by allowing a maximum number of particles per site, n max = 8, for which dynamics can be computed efficiently and the truncation effects are negligible.

Data availability
All relevant data of this study are available from the corresponding author upon request.