Interrelated Thermalization and Quantum Criticality in a Lattice Gauge Simulator

Since quantum gauge theories are computationally intractable in the non-perturbative regime, the idea of formulating gauge theories on discretized space-time lattices led to LGTs, and the developments of LGTs enable numerical simulation of gauge theories with various classical algorithms in the past decades Rothe:Book . Recently, a new trend is realizing LGTs in quantum simulators with ultracold atoms or trapped ions Martinez:2016rt ; Kokail:2019sv ; Schweizer:2019fa ; Mil:2020as ; Yang:2020oo ; Zhou:2021td . By taking the quantum advantage, these quantum simulation platforms offer the promise of a better understanding of LGTs than classical simulations.

One potential advantage of quantum simulation for studying LGTs lies in the non-equilibrium dynamics, such as quantum thermalization. In ultracold atom systems, the essential parameters in the gauge theory are tunable, allowing accessing different phases and the quantum criticality in between. Quantum criticality refers to a qualitative change of the ground state and the universal behavior of low-energy physics. On the contrary, thermalization dynamics usually involves highly excited states. Remarkably, it has been recently proposed that the thermalization dynamics of the $\ket{\mathbb{Z}_{2}}$ state can also signal the quantum critical regime in the U(1) LGT Yao:2022qm . However, it is challenging to study these physics due to the lack of techniques for manipulating and detecting local matter and gauge fields in the previous experiments.

In this work, we report on an experimental study of both equilibrium and thermalization dynamics in the quantum critical regime of the U(1) LGT. We integrate the techniques of single-site manipulation of matter and gauge fields and atom-number-resolved detection into an updated LGT simulator. With these technical advances, we can monitor the local Gauss’s law and then perform the post-selection, which eliminates the processes violating local gauge symmetries. We can measure the order parameter with various system sizes to perform a proper finite-size scaling, which overcomes the finite-size effects and determines the critical point accurately. We can also deterministically prepare the $\ket{\mathbb{Z}_{2}}$ state by addressing single atoms in a programmed manner. Equipped with these capabilities, we observe the nontrivial thermalization dynamics across the quantum critical regime Yao:2022qm .

Our experimental system is shown in Fig. 1(a), which realizes a U(1) LGT using bosons in optical lattices, and the protocol is the same as that reported in the previous work Yang:2020oo . We combine a short lattice and a long lattice with twice the lattice spacing to create a one-dimensional superlattice, as shown in Fig. 1(b) Zhang:2022fb . When the deep lattice sites are doubly occupied, the on-site energy is $U-2\delta$, where $U$ is the on-site interaction energy and $\delta$ is the energy offset between deep and shallow lattices. When $U\approx 2\delta$, the on-site energy of a doubly occupied site is nearly degenerate with an empty site, but is off-resonant with a singly occupied state. Thus, we only retain the empty and the doubly occupied states, and they are denoted by $\ket{\downarrow}$ and $\ket{\uparrow}$, respectively. These deep lattice sites are viewed as gauge sites, and the shallow lattice sites are viewed as matter sites. This realizes the LGT written as Yang:2020oo

where $l$ labels the matter sites, and $\hat{\psi}_{l}$ are bosonic operators for the matter fields. Spin-$1/2$ operators $\hat{S}_{l,l+1}$ are defined in the deep lattices. Here $m=\delta-U/2$ is the mass of the matter field. $\tilde{t}$ is generated by a second-order hopping process because the long-range direct hoppings are suppressed by applying a tilting potential.

This model possesses local gauge symmetries because the Hamiltonian Eq. (1) is invariant under the local gauge transformation $\hat{\psi}_{l}\rightarrow e^{i\theta_{l}}\hat{\psi}_{l}$, $\hat{S}^{+}_{l,l+1}\rightarrow e^{-i\theta_{l}}\hat{S}^{+}_{l,l+1}$ and $\hat{S}^{+}_{l-1,l}\rightarrow e^{-i\theta_{l}}\hat{S}^{+}_{l-1,l}$ Yang:2020oo . This local gauge symmetry gives rise to a set of local conserved quantities $G_{l}=S^{z}_{l-1,l}+S^{z}_{l,l+1}+n_{l}$. By introducing the electric field $E_{l-1,l}=(-1)^{l}S^{z}_{l-1,l}$ and the physical charge $Q_{l}=(-1)^{l}n_{l}$, the conservation $G_{l}=0$ can be written as

which is exactly the Gauss’s law of U(1) LGTs Yang:2020oo . A configuration of the electric field and Gauss’s law is also schematically shown in Fig. 1(b).

When focusing on the gauge sector with all $G_{l}=0$, there are two phases in the Hamiltonian Eq. (1) as tuning $m/\tilde{t}$. This can be easily seen by looking at two limits with $m/\tilde{t}=\pm\infty$. When $m/\tilde{t}\rightarrow+\infty$, the matter field favors $n_{l}=0$ and therefore $S^{z}_{l-1,l}+S^{z}_{l,l+1}=0$, leading to an antiferromagnetic state $\ket{\mathbb{Z}_{2}}=\ket{\uparrow\downarrow\uparrow\downarrow\dots}$ or $\ket{\downarrow\uparrow\downarrow\uparrow\dots}$ at gauge sites. This leads to two-fold degenerate ground state. When $m/\tilde{t}\rightarrow-\infty$, the matter field favors $n_{l}=1$ and therefore $S^{z}_{l-1,l}=S^{z}_{l,l+1}=-1/2$. This state does not break the original lattice translational symmetry. Hence the transition is a $Z_{2}$ symmetry breaking transition of the Ising type, and detailed theoretical calculations predict the critical point at $m/\tilde{t}\approx 0.655$ Sachdev:2002mi ; Fendley:2004cd ; Rico:2014tn . Signatures of two different phases have also been observed in the previous experiment Yang:2020oo .

To accurately determine the quantum critical point, we first prepare a $\ket{\mathbb{Z}_{2}}$ state with high fidelity utilizing the ability of single-site addressing (see Methods) at $\delta=447(1)$ Hz, $U=732(2)$ Hz, and $\tilde{t}=23.5(2)$ Hz, corresponding to $m/\tilde{t}=3.44(7)$ Weitenberg:2011ss . This $\ket{\mathbb{Z}_{2}}$ state is an antiferromagnetic state $\ket{\uparrow\downarrow\uparrow\downarrow\dots}$ in the gauge sites and empty sites in the matter sites. Here we prepare $5\sim 10$ copies of identical chains along the $x$ direction, as shown in Fig. 1(a) [See also Fig. S1(f) in Methods]. Each chain contains $L$ gauge sites and $L+1$ matter sites, with a total of $L+1$ atoms in each chain ($L=5,7,9$ as shown below). Then, we adiabatically ramp $\delta$ to the targeted critical regime around $373$ Hz, with a constant ramping speed $\dot{\delta}=2.1$ Hz/$\mathrm{ms}$. Such change of $\delta$ leaves $U$ and $\tilde{t}$ nearly unaffected, and it changes the value of $m/\tilde{t}$ to the critical regime, as shown in Fig. 2. Finally, the system is suddenly frozen and the atom number at each site is read out with the single-site-resolved microscope. Note that usually the fluorescence imaging cannot distinguish two atoms from zero atom at the same site Bakr:2009aq ; Sherson:2010sa . Here, before detection, we split two atoms into two sites of a double well along the $y$ direction if there are two atoms at the same site. This allows us to resolve atom number two from zero at a single site (see Methods). To mitigate the undesired effects from processes beyond the LGT, we post-select our data based on two rules: (i) the total atom number remains the same as that of the initial state, and (ii) the Gauss’s law Eq. (2) has to be obeyed for all sites. With the adiabatic ramping and the post-selection, we ensure that the ground state of the LGT at the targeted $m/\tilde{t}$ is reached. The results of single-site-resolved measurement of atom number distribution is shown in Fig. 2(a) for a range of $m/\tilde{t}$.

To locate the critical point, we measure the order parameter $L^{-1}\sum_{l}(-1)^{l}S^{z}_{l-1,l}$. This order parameter is the staggered magnetization, and given the definition of the electric field, this order parameter is also the spatial averaged electric field strength, denoted by $\bar{E}$. We plot $|\bar{E}|$ as a function of $m/\tilde{t}$ with $L=7$ in Fig. 2(c), which exhibits a rapid change when $m/\tilde{t}$ is in the range $\sim[0,1]$. To more accurately determine the critical point, we zoom in this range and perform a finite-size scaling. We plot $|\bar{E}|L^{\beta/\nu}$ for different system size $L=5,7,9$, and for the two-dimensional Ising universality class, the order parameter critical exponent $\beta=1/8$ and the correlation length critical exponent $\nu=1$ Francesco:Book . The crossing point of these curves locates the critical point (see Methods for numerical simulation). Here we indeed observe the crossing of three curves. However, since each data point has an error bar, it is not clear to visualize the exact crossing point. Hence, we consider each point as a normalized Gaussian distribution $p_{i}(x)=\mathcal{N}(\mu_{i},\sigma_{i}^{2})$ where $x=|\bar{E}|L^{1/8}$, $\mu_{i}$ and $\sigma_{i}$ are respectively the value and the error bar of each data point. We calculate the averaged Kullback–Leibler (KL) divergence of the three Gaussian distributions at each $m/\tilde{t}$, which quantifies the degree of overlap between three Gaussian distributions. The averaged KL divergence is defined as Sgarro:1981id

We plot the $\text{D}_{\text{KL}}$ as a function of $m/\tilde{t}$ in Fig. 4(a), and fitting the data yields a peak position at $m/\tilde{t}=0.59(8)$. This value is consistent with theoretical predictions Sachdev:2002mi ; Fendley:2004cd ; Rico:2014tn . Given the energy scale of our $\tilde{t}$, the deviation from the theoretical critical value is less than $2$ Hz.

We then consider the quench dynamics. The same $\ket{\mathbb{Z}_{2}}$ state is prepared as the initial state. Now, instead of adiabatic ramping, we suddenly change the value of $m/\tilde{t}$ to the targeted value nearby the critical point, with a typical time scale of $1~{}\mathrm{ms}$. Here we fix $\tilde{t}=34.1(3)$ Hz and $U=676(2)$ Hz. The value of $\tilde{t}$ is larger than what we used for the adiabatic ramping process. This is because we find that even a weak spatial inhomogeneity can cause noticeable variations of $m$ from site to site, which can affect the time dynamics. A larger $\tilde{t}$ helps to reduce this variation in terms of $m/\tilde{t}$ and can significantly suppress the effects of the inhomogeneity. We record the spatial atom number distribution at different times during the real-time dynamics using a single-site-resolved microscope. We also apply the same post-selection to ensure that the dynamics is governed by the LGT. The measurements of time dynamics are shown in Fig. 3 for different values of $m/\tilde{t}$.

We fit the experimental data in Fig. 3 with a damped sinusoidal function $Ae^{-t/\tau}\sin(\omega t)+E_{\infty}$, where $A$, $\tau$, $\omega$, and $E_{\infty}$ are all fitting parameters. We obtain $E_{\infty}$ as the long-time steady value, as illustrated by the dashed lines in Fig. 3. Meanwhile, we can also theoretically extract the thermalization values $E_{\text{th}}$, provided that the system obeys the eigenstate thermalization hypothesis and the initial $\ket{\mathbb{Z}_{2}}$ state thermalizes Deutsch:1991qs ; Srednicki:1994ca ; Rigol:2008ta ; DAlessio:2016fq . $E_{\text{th}}$ is obtained by calculating $\text{Tr}\left[\rho(T)L^{-1}\sum_{l}(-1)^{l}S^{z}_{l-1,l}\right]$. Here $\rho(T)$ is an equilibrium density matrix of the LGT system, with the temperature $T$ determined by matching energy $\text{Tr}[\rho(T)\hat{H}]=\bra{\mathbb{Z}_{2}}\hat{H}\ket{\mathbb{Z}_{2}}$ (see Methods). Previous work has predicted that for the PXP model, $E_{\infty}$ matches $E_{\text{th}}$ and the initial $\ket{\mathbb{Z}_{2}}$ state thermalizes only in the quantum critical regime Yao:2022qm . These two values depart from each other away from the critical point $(m/\tilde{t})_{\text{c}}$. When $m/\tilde{t}>(m/\tilde{t})_{\text{c}}$, the ground state is two-fold degenerate and the $\ket{\mathbb{Z}_{2}}$ state has large overlap with the ground state, and the ground state is always not thermalized. When $m/\tilde{t}<(m/\tilde{t})_{\text{c}}$, especially around $m/\tilde{t}=0$, it is known that the PXP model hosts a set of many-body scar states as the system’s eigenstates Bernien:2017pm ; Turner:2018we ; Turner:2018qs . The $\ket{\mathbb{Z}_{2}}$ state also has large overlap with the scar states, preventing it from thermalization Bernien:2017pm ; Turner:2018we ; Turner:2018qs . The PXP model is equivalent to this LGT model under the local gauge constraints of $G_{l}=0$ for all $l$, and therefore, the discussion also applies to this LGT model Surace:2020lg ; ChengTA . In Fig. 4(b), we compare $E_{\infty}$ with $E_{\text{th}}$, and our measurements agree with the prediction that $E_{\infty}\approx E_{\text{th}}$ only in the quantum critical regime Yao:2022qm .

In summary, we have performed a single-site-resolved quantitative experimental study of the quantum criticality in the U(1) LGT realized with bosons in optical lattices. Our study combines both the equilibrium property and the thermalization dynamics. Our system size is up to $\sim 19$ lattice sites, and our quantitative results agree well with the numerical results of exact diagonalization. This agreement benchmarks the validity of our ultracold atom quantum simulator quantitatively, and demonstrates this simulator as a powerful platform to study non-equilibrium dynamics of the gauge theory. In the near future, when our system size is enlarged several times, it will be beyond the capability of exact diagonalization. The experimental control and detection capability developed in this work can be used to study other interesting dynamical phenomena in this system, such as string breaking Banerjee:2012aq ; Hebenstreit:2013rt , dynamical transition between quantum phases Huang:2019dq ; Zache:2019dt , the false vacuum decay, and the confinement-deconfinement transition Surace:2020lg ; ChengTA ; Hauke . The current scheme of implementing the LGT can also be extended to higher dimensions Ott:2021sc .

Acknowledgement We thank Meng-Da Li, Qian Xie, Bo Xiao, Hui Sun, Wan Lin, An Luo for their help on building the experimental setup. This work was supported by NNSFC grant 12125409, Innovation Program for Quantum Science and Technology 2021ZD0302000.

METHODS AND SUPPLEMENTARY MATERIALS