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 novel dynamics and scaling effects beyond the scope of existing theories at superfluid-insulator quantum phase transitions and highlight spin populations as a new 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.

Physical modeling of three-dimensional (3D) strongly correlated many-body systems, especially their critical dynamics across quantum phase transitions, is a fast-moving research frontier with immediate applications spanning from adiabatic quantum state preparations to developing novel materials [1][2][3][4][5][6][7][8][9]. Impeded by the intrinsic complexity of these systems and limitations of existing numerical techniques, theoretical studies have been performed mainly in homogeneous systems and lower dimensions [1][2][3][10][11][12][13][14][15][16]. Here we experimentally demonstrate that such complex dynamics can be efficiently studied in 3D spinor Bose-Hubbard model quantum simulators, consisting of lattice-confined spinor condensates. Possessing spin degrees of freedom and exhibiting magnetic order and superfluidity, these 3D quantum simulators are highly programmable [3][4][5][6][17][18][19][20]. We find novel dynamics and scaling effects beyond the scope of existing theories at superfluid-insulator quantum phase transitions, and highlight spin populations as a new observable to probe these dynamics. We also conduct numerical simulations in lower dimensions using time-dependent Gutzwiller approximations, which qualitatively describe our observations.
The 3D spinor quantum simulators realized in this work, directly implementing the important Bose-Hubbard (BH) model, are well-isolated quantum systems consisting of antiferromagnetic spinor Bose-Einstein condensates (BECs) in cubic optical lattices [17]. When lattice-trapped bosons are suddenly quenched from the superfluid (SF) phase to the Mott-insulator (MI) phase, the BH model is nearly integrable, yielding non-thermal steady states preserving memory of the initial state [9,10]. One distinctive 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 [4,5,[16][17][18][19]. 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,[17][18][19][20]. The 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-6, 17, 18].
In this paper, we experimentally explore these quantum critical dynamics by demonstrating asymmetric impulse-adiabatic regimes [12,21] as spinor gases are quenched across superfluid to spinor Mott-insulator quantum phase transitions at different speeds. 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,[22][23][24][25]. Our observations, going beyond the adiabatic-impulse-adiabatic regimes studied in many publications [9,12,[22][23][24][25][26][27], 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 [12,21]. We also find novel 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. In addition, our data present advantages of using spin populations as a new observable for probing nonequilibrium dynamics, as they are less sensitive to latticeinduced heating than observables (e.g., visibility) used in scalar bosons. While numerical simulations in 3D multicomponent systems are prohibitively difficult to perform for further verification of our results, we conduct theoretical simulations in lower dimensions and find qualitative agreements with the experiment.
We start each experimental cycle with a F = 1 antiferromagnetic sodium spinor BEC at its SF ground state, the longitudinal polar state, in which ρ 0 = 1 and m = 0. Here ρ mF is the fractional population of atoms in the spin-m F state and m = ρ −1 + ρ +1 is the magnetization. We then linearly quench the depth u L of cubic lattices, which results in an exponential increase in the ratio of the spin-independent interaction 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,[16][17][18][19]. The static and dynamic properties of the spin-1 BH Hamiltonian are computed using the Gutzwiller approximation (see Methods). Figure 1 compares density profiles of various spin components, derived from numerical simulations in systems Predicted density profiles of all atoms, the mF = 0 atoms, and the mF = ±1 atoms after a slow quench from a shallow lattice of uL = 2ER to a deep lattice at 40ER, respectively. These are derived from our numerical simulations for 2D lattice-trapped sodium spinor gases of n peak = 6 at q/h = 42 Hz using the BH 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,e,f) Similar to Panels (a,b,c) but after a fast lattice quench.
similar to the experimental system with fewer dimensions, after the lattice depth is quenched from 2E R to 40E R at different speeds. Here E R is the recoil energy [5,17]. For slow quenches, Fig. 1(a)-(c) clearly display the formation of Mott insulating shells and spin structures. Whereas for fast quenches, Fig. 1(d)-(f) show 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 Information). 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. 2(a). 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. 2(a)) 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 7E R /ms curve in Fig. 2(a), 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. 2(a)). Qualitatively similar results are obtained from our 2D numerical simulations (see Fig. 2(b)), 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 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), due to extremely high power requirements on the lattice beams.
We apply a single ρ 0 cut-off value to mark the edge of the impulse regime. For example, Fig. 2(c) 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. 2(c)). For relatively slow quenches (v ramp 13E R /ms), the impulse regime does not strongly depend on quench speed. For faster quenches, the length of the impulse regime 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 lattice-trapped 3D scalar bosons [9]. The extracted scaling exponents show no strong dependence on the ρ 0 cut-off value as long as the cutoff is sufficiently larger than the minimum observed ρ 0 values. We find a qualitative agreement between the 2D simulation (dashed line in Fig. 2(c)) and our data, although the scaling exponents are hard to match between theory and experiment due to the presence of inhomogeneous systems and difference in dimensionality. 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 non-universal scaling relations. Our numerical simulations have taken this feature into account, while omitting any finite-temperature or heating effects.
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. A typical example at an even n and an odd n is shown in Fig. 2(d). To emphasize a key advantage of using spin populations as a new observable to probe SF-MI transitions, we further analyze the 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,28]. The inset of Fig. 3(a) illustrates how we extract φ. Typical evolutions of φ for the three spin com-ponents during a quench at an intermediate speed are presented in Fig. 3(a). 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 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. 3(b).
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. 3(c). These regions are divided at roughly the same critical v ramp (≈ 13E R /ms) as the two regions found in Fig. 2(c). For fast (slow) quenches where v ramp is faster (slower) than 13E R /ms, Fig. 3(c) shows that the extracted U 0 /J linearly increases (decreases) with v ramp , instead of the power law dependence found in Fig. 2(c). This apparent discrepancy together with large disagreements between the numerical simulations and the data (see Fig. 3(c)) may confirm SF order parameters are poor observables for non-adiabatic  quenches across SF-MI transitions. One possible explanation is the dynamic heating rate effects upon SF order parameters found in rapid lattice quenches [29,30]. This, coupled with the coexistence of SF and MI regions in inhomogeneous systems, makes it difficult to compare the impulse regimes derived from SF order parameters across varying quench speeds.
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 for our system, the lack of a significant q dependence is consistent with our theoretical understanding of these transitions. From these fittings, 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. 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,[16][17][18][19], Fig. 4 implies the scaling exponents are independent of the nature of the SF-MI phase transitions.
In conclusion, we have experimentally demonstrated novel quench dynamics and scaling effects outside the scope of existing theories in 3D spinor Bose-Hubbard model quantum simulators. Our data indicate the scaling exponents are independent of the nature of SF-MI phase transitions. We have also presented the advantages of using spin populations as a new observable to probe dynamics of SF-MI 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
Each experimental sequence begins by creating a spin-1 antiferromagnetic spinor BEC, containing up to 1.5 × 10 5 sodium ( 23 Na) atoms, in free space 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 |F = 1, m F = ±1 states. This leaves us with a BEC in its superfluid 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 coils 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 the ratio of U 0 /J. Here J is the hopping energy and U 0 is the spinindependent interaction. The atoms are then abruptly released and ballistically expand, before we measure ρ 0 and the superfluid order parameter with a two-stage microwave imaging method [6,18]. 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. 2(a) and Fig. 3(a) of the main text. Each experimental data point represents the average of around 15 shots at a given u final L and t ramp , with error bars representing one standard deviation of the shots. Each complete set at a given v ramp and q contains 14 or more points at representative u final L to allow fitting with either a sigmoid function or a piecewise function given by x > x 0 . The fitting function used for each set was determined by the chi-square of the fitted functions and was then used to calculate the critical lattice depth and its associated error.
The BH Hamiltonian of lattice-trapped F =1 spinor gases can be expressed as [17], wheren i = mF a † i,mF a i,mF is the total atom number at site-i, a † mF (a mF ) are boson creation (destruction) operators with m F being ±1 and 0 for spin-1 atoms, the total spin is F = mF ,m ′ F a † mF F mF ,m ′ F a m ′ F , the standard spin-1 matrices are F mF ,m ′ F , 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 Hamil- Here z is the number of nearest neighbors, ψ mF = a † i,mF = a i,mF is the superfluid order parameter, and we absorb the V T term in the chemical potential. The side peak fraction studied in Fig. 3 is given by | â i,mF | 2 for any site i and is thus proportional to the superfluid order parameter.
We compute the static and dynamic properties of the spin-1 BH Hamiltonian using the Gutzwiller approximation [13,16], 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 = L i=1 |φ i with |φ i = ni,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.