Prethermal gauge structure and surface growth in $\mathbb{Z}_2$ lattice gauge theories

Abstract

Universal aspects of thermalization in interacting many-body systems are typically challenging to derive microscopically, yet provide a powerful framework for understanding emergent phenomena. Here, we numerically study the mean-field dynamics of a $(2+1)$D spin system with thousands of spins and show that experimentally-feasible two-body Ising interactions can stabilize a prethermal $\mathbb{Z}_2$ lattice gauge structure with dynamical matter, manifested by a gauge-invariant plateau with exponentially long lifetime. Eventually, the metastable prethermal $\mathbb{Z}_2$ gauge structure breaks down via a proliferation of Gauss' law defects, similar to bubble formation in false vacuum decay. In this regime, we discover spatio-temporal correlations described by a non-linear surface growth consistent with the $(1+1)$D Kardar-Parisi-Zhang (KPZ) universality class. We benchmark our results in small systems against semi-classical discrete time Wigner approximation (DTWA) and exact diagonalization (ED), where the breakdown of DTWA signals the emergence of an extensive number of local symmetries that strongly influence the thermalization pathway. Our model provides a testbed for quantum simulators and is directly implementable in large-scale arrays of Rydberg atoms.

Figures (6)

• Figure 1: Spin model and emergent gauge structure.(a) Spins are located on the sites (matter) and links (electric/gauge fields) of a honeycomb lattice. We define a $\mathbb{Z}_2$ Gauss' law, see Eq. \ref{['eq:Gauss-law']}, together with a two-body Hamiltonian $\hat{H}_V$ that energetically protects a metastable sector at energy $E=0$ containing the gauge-invariant configurations $G_j=+1$ (inset). (b), (c) The mean-field dynamics describe the evolution of classical spins, i.e., unit vectors on a sphere. After the system is prepared in a random gauge-invariant state, the perturbation $\hat{H}_\Omega$ rotates the spins about the $x$ axis, leading to the dynamical violation of the local Gauss' law constraints at later times $t$. (d) We evolve gauge-invariant initial states for various protections strengths $V/\Omega=4,6,10,20,40,80$ (light to dark color) averaging over $200$ random realizations in a $20\times 20$ plaquette system. The gauge-symmetry breaking errors remain small and controlled for a long time, leading to prethermal plateau with an emergent $\mathbb{Z}_2$ gauge structure. (e) The protection strength $V/\Omega$ controls the magnitude (top) and critical time (bottom) of the prethermal plateau; black lines indicate fit functions. (f) Local Gauss' laws of one exemplary trajectory for multiple times after the plateau phase, see yellow markers in panel (d). Local nucleation regions followed by a proliferation of Gauss' law defects imply a rich spatio-temporal correlation structure.
• Figure 2: Mean-field phase transition.(a) Potential $F(G)$ for the translationally-invariant model without dynamical matter, for protection strengths $V/\Omega = 2,..,5$ (bright to dark) and $(V/\Omega)_c \approx 3.33$ (orange). The dynamics is governed by the kinetic equation $(\dot{G})^2 + F(G) = 0$ with initial conditions $G(0)=+1$ and $\dot{G}(0)=0$. (b) For a protection strength below the critical value, the error grows to its maximum in the translationally-invariant setting (red curve). As the protection is increased, the error is suppressed since the root in the potential $F(G)$ moves towards $G=1$. Inset: The Gauss' law shows large gauge-symmetry breaking oscillations for $(V/\Omega)=3.3 < (V/\Omega)_c$ (red curve), and constrained dynamics for $(V/\Omega)=3.4 > (V/\Omega)_c$ (light gray curve) with smaller oscillations for increasing protection, e.g., $V/\Omega=5$ (dark gray curve).
• Figure 3: Surface growth.(a) We investigate the dynamics starting from an initial configuration with Gauss' law defects in the bottom row. We find defect and defect-free regions separated by a surface (red line, right) that we quantify by the height function $h_n(X,t)$, where $n$ labels a trajectory for a set of random protection strengths $\{\delta_j \}$; here we set $V/\Omega=10$. (b) We analyze the mean height $h$ and width $\delta h$, Eq. \ref{['eq:height_function']}, showing ballistic spreading and a non-linear roughening, respectively; statistical error bars of the numerical data are smaller than the solid lines. The gray lines show power-law fits to the data for $50 < t \Omega < 1000$, where the width $\delta h$ does not yet crossover to the asymptotic Family-Vicsek behavior. (c) Top: Width $\delta h$ for increasing linear system size $L=40,80,160,320$ (light to dark). Bottom: The dynamical scaling analysis reveals an excellent collapse of the data for a roughening exponent $\alpha=1/2$ and dynamical exponent $z=3/2$.
• Figure 4: Quantum dynamics.(a) The model without dynamical matter is naturally realized in Rydberg tweezer arrays with nearest-neighbor density-density interactions of strength $V$ between excited state atoms $\ket{r}$ and a laser drive $\Omega$ between $\ket{g}$ and $\ket{r}$. By choosing the detuning of the drive in the facilitation-like regime, $\Delta=V$, an even number of excitations sector is stabilized. The energetic separation of even and odd excitation number sectors stabilizes the $\mathbb{Z}_2$ symmetry. (b,c) We compare the dynamics in small systems using classical mean-field dynamics, DTWA, and ED for models with (b) and without (c) dynamical matter, for $J=0$. In both cases, the mean-field dynamics has excellent agreement with ED, whereas the DTWA only marginally captures the prethermal plateau due to the emergent local $\mathbb{Z}_2$ symmetry.
• Figure 5: Numerical simulations. We solve the set of differential equations \ref{['eq:DiffEq']} for a system with dynamical matter on a $L=2\times 2$ plaquette lattice and periodic boundary conditions; the number of spins is $N_{\rm spins} = 20$. We set the protection strength $V/\Omega=10$ and $J = \Omega = 1$ and we consider a single trajectory in a random potential with $\delta V=0.1V$ as well as a random initial configuration in the Gauss' law sector with $G_j=+1$ for all $j$. For the trotterization procedure, we choose the timestep to be $\delta t = 0.01$. For the ODE solver, we use the MATLAB function ode45. (a) We plot the time-averaged Gauss' law error $\varepsilon(t)$ solving the differential equations with an ODE solver and using the trotterization procedure. Only at late times ($t\Omega > 10^2$) we observe slight deviations. (b) The equations-of-motion conserve energy $E(t)=E(t=0)$ and hence we consider the energy density $e(t) = E(t)/N_{\rm spins}$ as another observable to benchmark the validity of the trotterization procedure. For an initial gauge-invariant state, the energy density is zero, $e(0)=0$. While the ODE solver exactly conserves energy, the trotterized time evolution shows small oscillations $|e(t)| < 10^{-2} \cdot V$ at later times.