Universal space-time scaling symmetry in the dynamics of bosons across a quantum phase transition

The dynamics of many-body systems spanning condensed matter, cosmology, and beyond is hypothesized to be universal when the systems cross continuous phase transitions. The universal dynamics is expected to satisfy a scaling symmetry of space and time with the crossing rate, inspired by the Kibble-Zurek mechanism. We test this symmetry based on Bose condensates in a shaken optical lattice. Shaking the lattice drives condensates across an effectively ferromagnetic quantum phase transition. After crossing the critical point, the condensates manifest delayed growth of spin fluctuations and develop anti-ferromagnetic spatial correlations resulting from sub-Poisson generation of topological defects. The characteristic times and lengths scale as power-laws of the crossing rate, yielding the temporal exponent 0.50(2) and the spatial exponent 0.26(2), consistent with theory. Furthermore, the fluctuations and correlations are invariant in scaled space-time coordinates, in support of the scaling symmetry of quantum critical dynamics.

The dynamics of many-body systems spanning condensed matter, cosmology, and beyond is hypothesized to be universal when the systems cross continuous phase transitions. The universal dynamics is expected to satisfy a scaling symmetry of space and time with the crossing rate, inspired by the Kibble-Zurek mechanism. We test this symmetry based on Bose condensates in a shaken optical lattice. Shaking the lattice drives condensates across an effectively ferromagnetic quantum phase transition. After crossing the critical point, the condensates manifest delayed growth of spin fluctuations and develop anti-ferromagnetic spatial correlations resulting from sub-Poisson generation of topological defects. The characteristic times and lengths scale as power-laws of the crossing rate, yielding the temporal exponent 0.50 (2) and the spatial exponent 0.26 (2), consistent with theory. Furthermore, the fluctuations and correlations are invariant in scaled space-time coordinates, in support of the scaling symmetry of quantum critical dynamics.
Critical phenomena near a continuous phase transition reveal fascinating connections between seemingly disparate systems that can be described via the same universal principles. Examples exist in superfluid helium [1], liquid crystals [2], biological cell membranes [3], the early universe [4], and cold atoms [5,6]. An important prediction is the power-law scaling of the topological defect density with the rate of crossing a critical point, as first discussed by T. Kibble in cosmology [4] and extended by W. Zurek in the context of condensed matter [1]. Their theory, known as the Kibble-Zurek mechanism, has been the subject of intense experimental study that has largely supported the scaling laws [7]. Recent theoretical works further propose the universality hypothesis that the collective dynamics across a critical point should be invariant in the space and time coordinates that scale with the Kibble-Zurek power-law [8][9][10].
In this report we study the critical dynamics of Bose condensates in a shaken optical lattice crossing an effectively ferromagnetic quantum phase transition. The transition occurs when we ramp the shaking amplitude across a critical value, causing the atomic population to bifurcate into two pseudo-spinor ground states [28]. We measure the growth of spin fluctuations and the spatial spin correlations for ramping rates varied over two orders of magnitude. Beyond the critical point we observe delayed development of spin domains with antiferromagnetic correlations, a distinctive feature of the non-equilibrium dynamics. The times and lengths char-acterizing the critical dynamics agree excellently with the scaling predicted by the Kibble-Zurek mechanism. We further observe that the measured fluctuations and correlations collapse onto single curves in scaled space and time coordinates, supporting the universality hypothesis.
Our experiments utilize Bose-Einstein condensates (BECs) of cesium atoms. We optically confine the condensates with trap frequencies of (ω x , ω y , ω z ) = 2π × (12, 30, 70) Hz, where the long (x ) and short (y ) axes are oriented at 45 o with respect to the x and y coordinates (Fig. 1A). The tight confinement along the vertical z-axis suppresses non-trivial dynamics in that direction, which is also the optical axis of our imaging system. We adiabatically load the condensates into a one-dimensional (1D) optical lattice [11] along the x-axis with a lattice spacing of 532 nm and a depth of 8.86 E R , where E R = h × 1.33 kHz is the recoil energy and h is Planck's constant.
To induce the ferromagnetic quantum phase transition, we modulate the phase of the lattice beam to periodically translate the lattice potential by ∆x(t) = (s/2)sin(ωt), where s is the shaking amplitude and the modulation frequency ω is tuned to mix the ground and first excited lattice bands [28,30]. The hybridized ground band dispersion can be modelled for small quasimomentum q = (q x , q y , q z ) by where m is the atomic mass, and the coefficients of its quadratic (α) and quartic (β) terms depend on the shaking amplitude (Fig. 1B). For shaking amplitudes below the critical value the coefficient α is positive and the BEC occupies the lone ground state at momentum q = 0. The quantum phase transition occurs when the quadratic term crosses zero at s = s c , where α = 0 and β > 0. Stronger shaking converts the dispersion into a  double-well with α < 0, yielding two degenerate ground states with q x = ±q * . Repulsively-interacting bosons with this double-well dispersion are effectively ferromagnetic, leaving two degenerate many-body ground states with all atoms either pseudo-spin up (q x = q * ) or down (q x = −q * ). Notably, transitioning to one of these two ground states requires the system to spontaneously break the symmetry of its Hamiltonian. Describing the dynamics across the critical point presents a major challenge due to the divergence of the correlation length and relaxation time (critical slowing).
The Kibble-Zurek mechanism provides a powerful insight into quantum critical dynamics. According to this theory, when the time remaining to reach the critical point inevitably becomes shorter than the relaxation time, the system becomes effectively frozen, see Fig. 1C. The system only unfreezes at a delay time t KZ after passing the critical point, when relaxation becomes faster than the ramp. At this time topological defects form, and the typical distance d KZ between neighboring defects is determined by the equilibrium correlation length. The Kibble-Zurek mechanism predicts that t KZ and d KZ depend on the quench rateṡ as where z and ν are the equilibrium dynamical and correlation length exponents given by the universality class of the phase transition.
For slow ramps t KZ and d KZ diverge and become separated from other scales in the system, making them the dominant scales for characterizing the collective critical dynamics [8][9][10]. This idea motivates the universality hypothesis, which can be expressed as indicating that the critical dynamics of any collective observable f obeys the scaling symmetry and can be described by a universal function F of the scaled coordinates x/d KZ and t/t KZ . The only effect of the quench rate is to modify the length and time scales. We test the scaling symmetry of time by monitoring the emergence of quasimomentum fluctuations at different quench rates. After loading the condensates into the lattice, we ramp the shaking amplitude linearly from s = 0 to values well above the critical amplitude s c = 13.1 nm [31] and interrupt the ramps at various times to perform a brief time-of-flight (TOF) before detection. After TOF the quasimomentum distribution of the sample can be extracted from the deviation δn(r) in the density difference between the +1 and −1 Bragg diffraction peaks [31]. This detection method is particularly sensitive when the quasimomentum just starts deviating from zero, indicating the emergence of fluctuations in the ferromagnetic phase where the ground states have non-zero quasimomentum.
Over a wide range of quench rates the evolution of quasimomentum fluctuation can be described in three phases ( Fig. 2A). First, below the critical point, quasimomentum fluctuation does not exceed its baseline level. Second, just after passing the critical point, critical slowing keeps the system "frozen", and fluctuation remains low. Finally, the system unfreezes and quasimomentum fluctuation quickly increases and saturates, indicating the emergence of ferromagnetic domains. We quantify this progression by investigating the fluctuation of contrast ∆C = δn 2 /n 2 (Fig. 2B) that tracks quasimomentum fluctuation in our condensates [31], where n is the to-tal density and the angle brackets denote averaging over space and over multiple images. We find empirically that the growth of fluctuations is well fit by the function where the time t is defined relative to when the system crosses the critical point at t = 0, t d characterizes the delay time when the system unfreezes, and t f is the formation time over which the fluctuation grows. The measurement of fluctuation over time provides a critical test for both the Kibble-Zurek scaling and the universality hypothesis. First, both t d and t f exhibit clear power-law scaling with the quench rateṡ varied over more than two orders of magnitude (Fig. 2C). Power-law fits yield the exponents of a d = 0.50(2) and a f = 0.50(6), respectively. The nearly equal exponents are suggestive of the universality hypothesis, which requires all times to scale identically. Indeed, the growth of contrast fluctuation ∆C follows a universal curve when time is scaled by t d (Fig. 2D), strongly supporting the universality hypothesis (Eq. 4).
We next test the spatial scaling symmetry based on the structures of pseudo-spin domains that emerge after the system unfreezes. Here, we cross the critical point with two different protocols: the first is a linear ramp starting from s = 0, while the second begins with a jump to s = s c , followed by a linear ramp. We detect domains near the time t = 1.4 t d in the spin density distribution j z (r) = n + (r) − n − (r) based on the density n +/− of atoms with spin up/down [31]. At this time the spin domains are fully-formed and clearly separated by topological defects (domain walls), shown in Fig. 3A. We characterize the domain distribution with the spin correlation function [17,28], averaged over multiple images, see Fig. 3B. Both ramping protocols lead to similar correlation functions, suggesting that the formation of topological defects does not depend on dynamics below the critical point. The spin correlations reveal rich domain structure that strongly depends on the quench rate. For slower rampsṡ < 1.3 nm/ms the structures are predominantly one-dimensional and the density of topological defects increases with the quench rate. When the quench rate exceeds 1.3 nm/ms, defects start appearing along the y-axis, and the domain structures become multidimensional. We attribute this dimensional crossover to the unfreezing time becoming too short to establish correlation in the non-lattice directions. For the remainder of this work we focus on the slower quenches and investigate the spin correlations along the lattice direction.
We examine the one-dimensional correlations using line cuts of the density-weighted correlation functions g(r) = G(r)/ n(R + r)n(R)dR [17,28]. The results exhibit prominent decaying oscillation, shown in Fig. 3C. We extract two essential length scales from the correlation functions: the average domain size d, or equivalently the distance between neighboring topological defects, and the correlation length ξ, indicating the width of the envelope function. These two scales are determined from the position and width of the peak in the Fourier transform of g(x) [31].
These length scales enable us to test the spatial scaling symmetry. The lengths d and ξ both display powerlaw scaling consistent with the Kibble-Zurek mechanism, see Fig. 3D, with fits yielding exponents b d = 0.26 (2) for the domain size and b ξ = 0.26(5) for the correlation length. Furthermore, the correlations, measured at the same scaled time, collapse to a single curve in spatial coordinates scaled by the domain size d, see Fig. 3E. This result strongly supports the spatiotemporal scaling from the universality hypothesis (Eq. 4). An empirical curve, well fits the universal correlation function, yielding σ = 1.01(1) and γ = 1.04(1), indicating that the width of the envelope is close to the typical domain size. Furthermore, the Gaussian envelope can be characterized by a thermal length of λ s = h/ √ 2πmk B T s = √ 4πσd, where k B is the Boltzmann constant and the effective spin temperature T s ∝ṡ 2b d scales with the quench rate and reaches 20 pK for our slowest ramps ofṡ = 0.04 nm/ms.
The most striking feature of the universal correlation function is the emergence of oscillatory, antiferromagnetic order in the ferromagnetic phase. In thermal equilibrium, ferromagnets are expected to have a finite correlation length but no anti-correlation. We attribute the appearance of anti-ferromagnetic order in our system to the preferential generation of domains with a certain size during the quantum critical dynamics. A statistical analysis of the topological defect distribution reveals that the domain sizes are bunched with their standard deviation σ d = 0.31(2)d well below their mean, indicating that the topological defects are created by a sub-Poisson process [31].
Finally, the combined scaling exponents of space and time allow us to extract the equilibrium critical exponents based on the Kibble-Zurek mechanism, see Fig. 3F. Solving Eqs. 2 and 3, we obtain the dynamical exponent z = 1.9(2) and correlation length exponent ν = 0.52 (5), which agree with the mean-field values z = 2 and ν = 1/2 within our experimental uncertainty. Note that the dynamical critical exponent z = 2 results from the unique quartic dispersion = βq 4 x of our system at the critical point [31].
In summary, our experiment reveals a universal, spatiotemporal scaling symmetry of the dynamics across a quantum critical point. The observed scaling laws are in excellent agreement with the prediction from the Kibble-Zurek mechanism. Furthermore, the universal correlations exhibit intriguing anti-ferromagnetic order which would not be expected in equilibrium. Direct identification of the domain walls enables us to show that the antiferromagnetic correlations are connected to sub-Poisson generation of topological defects. The scaling of the correlation functions suggests that the anti-ferromagnetic order may be a shared feature of quantum critical dynamics for phase transitions in the same universality class, meriting future experiments.
We SUPPLEMENTARY MATERIAL Experiment Setup. Our condensates form in an optical dipole trap at the crossing of three lasers with wavelength λ = 1064 nm. After evaporation the condensates are nearly pure, consisting of 30 000 to 40 000 atoms with temperatures less than 10 nK. We adiabatically load the condensates into an optical lattice by retro-reflecting the trapping beam along the x-axis. The apparatus which enables us to shake the lattice is described in Ref. (28).
Our experiments rely on a careful choice of the parameters governing the shaking optical lattice. We set the shaking frequency ω = 2π × 8.00 kHz slightly above the zero-momentum band gap h × 7.14 kHz, such that shaking raises the energy near the center of the ground band. During shaking we reduce the scattering length to a = 2.1 nm using a Feshbach resonance to lower the heating rate (32). Finally, immediately before time-offlight (TOF) we reduce the scattering length to a = 0 to prevent collisions while the atoms separate into distinct Bragg peaks.
Based on the lattice depth and shaking frequency, we calculate the critical shaking amplitude s c = 13.1 nm using Floquet theory (28), above which the system acquires a double-well dispersion. We base our calculation of a d = 0.50(2) on this theoretical critical point. To ver-ify the critical value s c = 13.1 nm, we allow the critical shaking amplitude to be a free parameter in a power law fit of t d while fixing the exponent to its theoretical value of 0.5. The fit yields s c = 13.8 (6) nm which is consistent with the calculated critical amplitude.

Analysis of Quasimomentum Fluctuation
We use the density fluctuation in the Bragg peaks to detect quasimomentum fluctuation in the gas. Nonzero local quasimomentum q changes the local density difference between Bragg peaks. To detect this change we perform a brief TOF with duration t TOF = 5 ms, which is long enough to separate the Bragg peaks but short enough that spatial information is preserved. From our images we calculate the density difference ∆n(r) = n −1 (r)−n 1 (r) where n i (r) is the density of the i'th Bragg peak. To remove the small offset in ∆n which exists at momentum q = 0, we calculate the density deviation δn(r) = ∆n(r) − ∆n(r) , where the angle brackets denote averaging over multiple images. The shift δn is nearly proportional to the local quasimomentum regardless of the shaking amplitude (Fig. 4). Finally, we calculate the contrast fluctuation ∆C = δn 2 /n 2 which closely tracks quasimomentum fluctuation in our condensates, where n(r) is the total density. The angle brackets denote averaging over many images and over the position within each sample.
In order to remove spurious sources of fluctuation such as photon and atom shot noise, we subtract the baseline value of ∆C for each ramp rate, which is given by the average of the three measurements at the earliest times taken below the critical point. Furthermore, even though quasimomentum fluctuation should continue to grow as q * 2 with increasing shaking amplitude, where ±q * are the quasimomenta of the ground states, we find that ∆C appears to saturate to a nearly constant value for times well beyond t d . We attribute saturation to the typical displacement q * t TOF /m during TOF becoming larger than the correlation length, such that fluctuation is dominated by the motion rather than the quasimomentum. We normalize ∆C to its saturated value at each ramp rate for convenient comparison. We determine the saturated value by the average of the latest three measured values, which are taken well beyond the delay time t d .

Reconstruction and Analysis of Pseudo-spin Domains
Reconstructing spin density. Reconstruction enables us to study the in-situ spin density distribution, including the domain structure. After linearly ramping the shaking amplitude at ramp rateṡ until t = 1.4 t d when domains have fully formed, we rapidly increase the shaking amplitude by as much as a factor of 15 over only 0.5 ms (Fig. 5A). The sudden increase in shaking amplifies the signal by exciting the two spin states to predominantly occupy different Bragg peaks. Depending on the shaking amplitude at t = 1.4 t d , we adjust the timing and the shaking amplitude immediately before release in order to maximize the distinguishability of the two spin states. We use images of condensates with uniform spin (Fig. 5B) to calibrate the projection of each spin state onto Bragg peaks (Fig. 5C). In comparison to the procedure without enhanced shaking used in Ref. (28), the amplification stage improves the fraction of atoms which distinguish the spin states from 23% to about 71%, corresponding to an increase in signal by more than a factor of three. After the enhanced shaking period, we perform a 3 ms TOF and measure the density in each Bragg peak (Fig. 5D), from which we can reconstruct the spin density distribution (Fig. 5E) using an algorithm similar to that described in Ref. (28).
Minimizing bias of the total polarization. The effectively ferromagnetic quantum phase transition can be biased by a nonzero initial velocity of the condensate relative to the lattice (28). In order to focus our study on the dynamics across an unbiased quantum phase transition, we test the total spin polarization P = j z ( R)d R/ n( R)d R of each reconstructed image, which is expected to be close to zero for unbiased samples. Indeed, under most conditions (0.16 ≤ṡ < 1.0 nm/ms) we find that more than 90% of images have total polarization |P | < 0.3. The correlation analysis excludes the remaining biased images with |P | > 0.3. For very slow ramps (ṡ < 0.16 nm/ms) starting from s = 0, we find that many samples are biased, likely due to increased susceptibility to a small, uncontrolled velocity between the condensate and the lattice. We have excluded data from these conditions to avoid poor statistics.
Removing systematic effects due to finite imaging resolution. We study the one-dimensional domain structures along the lattice direction (x) by taking cuts g meas (x) along the long-axis of the measured, normalized correlation functions g(r). Since the domain walls are predominantly oriented along the non-lattice direction (y), long axis cuts maximize the range of the correlation functions that we can evaluate but still reflect the structure along the lattice direction.
To obtain the physical spin correlations we must remove the systematic effects of our finite imaging resolution. Since the correlation functions depend on the spin density at both ends of the displacement vector, the measured correlations g meas (x) are the physical correlations g(x) convolved with the point spread function P (x) twice (33). We calculate the Fourier transform of the deconvolved correlation functiong(k) =g meas (k)/P 2 (k) from the Fourier transforms of the measured correlations g meas (k) and of the point spread functionP (k). Inverting the Fourier transform produces the correlation functions g(x) shown in Fig. 3C. Furthermore, from the peak position k p ing(k) we extract the typical domain size d = π/k p , and from the full width at half maximum ∆k of the peak we extract the correlation length ξ = π/∆k.

Sub-Poisson generation of domain walls
To better understand the process which generates domain walls we calculate the domain size distribution from our images. We identify domain walls by integrating the spin density along the y-direction, filtering noise at the single pixel scale (0.6 µm) which is below our resolution limit, and locating where the spin density changes sign. We calculate the domain sizes from the distances between neighboring walls. Since the correlation functions in scaled space are invariant with quench rate we focus on a single ratė s = 0.08 nm/ms, for which the domains are relatively large (d = 6.0 µm) and easy to resolve. The resulting domain size distribution (Fig. 6A) is tightly bunched around its mean. This bunching would not be expected for a Poisson process, which should exhibit an exponential distribution due to the constant probability of forming a domain wall at any location. Similarly, Poisson generation of domain walls would lead to exponentially decaying correlations that are qualitatively distinct from the oscillatory correlations observed in our experiments (Fig. 6B).

Equilibrium critical exponents ν and z
The partition function Z of our system near the critical point where α → 0 and β > 0 can be written in the path integral form as