Quantum computing universal thermalization dynamics in a (2 + 1)D Lattice Gauge Theory

Abstract

Simulating non-equilibrium phenomena in strongly-interacting quantum many-body systems, including thermalization, is a promising application of near-term and future quantum computation. By performing experiments on a digital quantum computer consisting of fully-connected optically-controlled trapped ions, we study the role of entanglement in the thermalization dynamics of a Z₂ lattice gauge theory in 2+1 spacetime dimensions. Using randomized-measurement protocols, we efficiently learn a classical approximation of non-equilibrium states that yields the gap-ratio distribution and the spectral form factor of the entanglement Hamiltonian. These observables exhibit universal early-time signals for quantum chaos, a prerequisite for thermalization. Our work, therefore, establishes quantum computers as robust tools for studying universal features of thermalization in complex many-body systems, including in gauge theories.

Probing quantum many-body systems while undergoing thermalisation is challenging, especially when looking for signatures of ergodicity and quantum chaos. Here, the authors study a lattice gauge theory in 2+1 dimensions using a trapped-ion-based universal digital quantum computer, unveiling the role of entanglement in the thermalization dynamics.

Introduction

Thermalization of isolated quantum many-body systems in e.g., ultra-cold atomic gases^{1, 2, 3, 4, 5, 6, 7–8}, trapped ions^{9, 10, 11–12}, condensed matter physics^13,14, cosmology¹⁵, and nuclear and high-energy physics^16,17, remains a vibrant frontier. Most quantum systems thermalize according to the Eigenstate Thermalization Hypothesis (ETH)^18,19, which posits that, under broadly applicable conditions, the long-time average of certain observables agrees with microcanonical predictions^20,21. Yet, probing large quantum many-body systems undergoing thermalization is inherently challenging due to the non-equilibrium nature of the process, which precludes using first-principle classical computational techniques. Recent advancements in quantum information theory and experiment^{22, 23, 24, 25–26} have brought this topic within immediate experimental reach, potentially enabling, via quantum simulation^{27, 28–29}, verification of important thermalization paradigms, and discovery of novel principles grounded in quantum information science.

While it is commonly posited that quantum chaos and ergodicity are prerequisites for thermalization^30,31, their demonstration in the context of quantum many-body systems remains somewhat elusive. Indicators of chaos and ergodicity involve measures associated with the eigenvalues of a given Hamiltonian^{32, 33–34}, or the properties of its eigenstates³⁵. A connection has been drawn recently between the entanglement of quantum states and quantum chaos, via the so-called entanglement Hamiltonians (EHs)^36,37. This connection can, in principle, be leveraged to experimentally probe thermalization via quantum simulation. Most studies to date, nonetheless, have stayed in the realm of theoretical exploration^{38, 39, 40, 41, 42–43}.

Analog quantum simulation allows one to monitor a quantum system continuously in time. However, current analog quantum simulators have limited programmability and are restricted to probing specific physical models and a limited set of observables^{9, 10, 11–12}. In contrast, universal digital quantum computers allow, in principle, the probing of the dynamics of any physical model. The universal control in digital quantum computers enables the use of tomographic techniques to extract a wide range of observables, including entanglement. Digitized time evolution, via Trotterization^44,45 or other schemes^46,47, limits the near-term simulations to times that are short compared to those of thermalization. Nonetheless, there exist universal phenomena that are indicative of quantum chaos and ergodicity at earlier times. Probing these phenomena, therefore, is an attractive near-term opportunity for digital quantum computers.

Gauge theories and their lattice formulations are among prime physical models whose simulations will benefit from quantum-computing technology^{29,48, 49, 50–51}. Gauge theories are key in high-energy and nuclear physics^52,53, condensed and synthetic quantum matter^{54, 55, 56–57}, local fermion-to-qubit mappings^{58, 59–60}, and quantum-error correction^{61, 62, 63, 64–65}. Studying thermalization dynamics of gauge theories, e.g., in the early universe and in high-energy particle collisions, remains challenging using first-principles simulation methods¹⁷. As a first step in experimentally probing thermalization dynamics of gauge theories, we study a Z₂ lattice gauge theory (LGT) in 2 + 1 spacetime dimensions^66,67 using a digital trapped-ion quantum computer^{67, 68, 69–70}. We use a chain of fifteen ¹⁷¹Yb⁺ ions to realize a general-purpose fully-connected digital quantum computer with twelve qubits, schematically shown in Fig. 1a, and use this computer to natively and accurately encode the system’s initial state, evolve this state in time, and measure observables.

[See PDF for image]

Fig. 1

Overview.

a Schematic of the trapped-ion experiment: 15 optically-controlled ions (yellow circles) in a linear trap (electrodes shown as brown rectangular pads) realize a universal digital quantum computer. Single-qubit and all-to-all two-qubit gates are implemented by an array of individually-focused laser beams (blue and purple, vertical) and a global laser beam (blue, horizontal). Two-qubit gates are performed using pairs of individual beams, such as the ones highlighted in purple. b Schematic of the randomized-measurement protocol for entanglement-Hamiltonian tomography. This protocol extracts a classical approximation of a reduced density matrix associated with a subsystem of the quantum state $\left.∣\right)\psi \left(t\right)\left(\right.⟩$, from which the presence or absence of quantum chaos is inferred. The protocol consists of measuring observables in a single-qubit randomized basis, then classically learning the entanglement Hamiltonian H_E({β_i}), which parameterizes the reduced quantum state with parameters {β_i}, so as to optimally reproduce all measurements. The statistical behavior of the eigenvalue spectrum of the entanglement Hamiltonian is then analyzed: e.g., eigenvalue repulsion indicates quantum chaos, as detailed in the main text.

Our analysis relies on EH tomography^{71, 72, 73, 74, 75–76} in combination with randomized-measurement protocols^{77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87–88} to learn representations approximating non-equilibrium states. The EH is defined as $H_{\mathrm{E}}:=-\log \left({\rho }_A\right)$, where ρ_A denotes the reduced density matrix of subsystem A formed by bipartitioning a (pure) quantum state. The utility of EH as a theoretical and experimental tool stems from the observation that, in many cases, it consists of approximately local operators^89,90. This resemblance to conventional Hamiltonians, i.e., energy operators governing system dynamics, facilitates simple theoretical analysis: If an EH is (k-)local, it can be described by parameters whose number scales polynomially with the system size, unlike the matrix ρ_A, which requires an exponential number of parameters.

Randomized-measurement protocols⁸⁸ can be used to learn the EH: after repeatedly preparing and then measuring the quantum state in a randomly chosen basis, one can fit a parameterized EH to the obtained measurement outcomes. Because the spectrum of H_E, known as the entanglement spectrum, scales logarithmically with the Schmidt eigenvalues of ρ_A, it is difficult to learn an EH and reproduce the entire entanglement spectrum with arbitrary accuracy. Randomized-measurement protocols may require resolving exponentially small probabilities even for relatively small systems, which is out of reach of current quantum devices and perhaps fault-tolerant devices, too. The present experiments can, therefore, only aim at approximating H_E and accurately reproducing its low-energy spectrum. Applications of this approach include the verification of topological phases^41,85.

Nonetheless, what if the precise quantitative structure of H_E were not crucial, but only its statistical properties were? This viewpoint is that of random matrix theory^91,92, where metrics like level distribution and spectral form factor differentiate between integrable and chaotic dynamics, while being indifferent to quantitative details. For instance, a model’s symmetries influence spectral correlations of the EH and must, therefore, be taken into account^93,94. This statistical perspective underpins our study. Guided by physical insights regarding the expected operator content of a non-equilibrium EH, we ask whether one can learn, from experimental data, a classical representation of a state to answer a simpler question: does a quantum state exhibit universal signatures of quantum chaos evident in the statistical properties of its EH? Crucially, can we discern this scenario from one where the state lacks chaotic behavior?

To answer these questions, we focus on two observables indicative of quantum ergodic and chaotic behavior: the entanglement(-Hamiltonian) gap-ratio distribution (EGRD)^{38, 39–40} and the entanglement spectral form factor (ESFF)^42,95. The EGRD is predicted to exhibit level repulsion for chaotic states, in contrast to uncorrelated levels in non-chaotic scenarios. Similarly, the ESFF is predicted to display a ramp-plateau structure in chaotic states. Strictly speaking, both quantities indicate quantum ergodicity, which usually implies quantum chaos. To experimentally test these predictions, we initialize our system in a product state and perform a quantum quench by digitally evolving the system for a variable length of time. Our protocol to access these quantities is depicted schematically in Fig. 1b. At short evolution times, we observe an EGRD that is characteristic of uncorrelated states. With increasing evolution time, we observe the onset of level repulsion and a ramp-plateau structure in the ESFF, indicative of quantum chaos. Detailed analysis, including comparison with emulated data, suggests that the observed behavior primarily arises from the thermalization dynamics of the isolated quantum system under simulation, with experimental inaccuracies and shortcomings of our ansatz playing a small, though non-negligible, role. Our work, therefore, establishes present-time quantum computers as robust tools for studying the universal feature of thermalization dynamics in complex many-body systems.

Results

We focus on a Z₂ LGT in 2 + 1 D with the Hamiltonian

1

$H={\sum }_{\square }{W}_{\square }+g{\sum }_{\ell }{\sigma }_{\ell }^z.$

Here, $W_{\square }:={\prod }_{\ell \in \square }{\sigma }_{\ell }^x$ is a magnetic-field operator, where □ denotes the elementary (square) plaquettes of the lattice, see the Methods section. ${\sigma }_{\ell }^x$ and ${\sigma }_{\ell }^z$ are Pauli operators representing gauge link and electric fields on edge ℓ, respectively. These act on (spin-$\frac{1}{2}$) hardcore bosons residing on the edges of a two-dimensional spatial square lattice. The first term in Eq. (1) corresponds to the magnetic energy, and the second to the electric energy, with coupling g controlling the relative strength of the two non-commuting contributions. In this formulation, the z basis, therefore, represents the electric basis, while the x basis corresponds to the magnetic basis. The Hamiltonian is expressed in dimensionless units, and ℏ is set to one throughout.

We consider the joint + 1 eigenspace of the Gauss-law operators $G_j:={\prod }_{\ell \in {+}_j}{\sigma }_{\ell }^z$ as the physical Hilbert space, where ℓ ∈ + _j denotes the links adjacent to a lattice site j. In the subsequent discussion, we consider a quasi-1D chain composed of L_x plaquettes along the x-direction with periodic boundary conditions. Fixed boundary conditions are applied along the short side of the chain (with L_y = 1 plaquettes) as illustrated in Fig. 2 (a). In this configuration, all Gauss laws contain three links at each site. Consider a short-ribbon operator $V_y:={\sigma }_{\alpha }^z{\sigma }_{\beta }^z$, where α and β are the two opposing links along the inner and the outer circumference shown in Fig. 2. The position of α and β is arbitrary, and V_y is unique: different choices of α and β are related to each other upon application of the Gauss-law operator under which states are invariant This operator commutes with the Hamiltonian and defines superselection sectors of the model; we work in the V_y = 1 sector (i.e., short-ribbon operators possess eigenvalue one).

[See PDF for image]

Fig. 2

Model overview.

a The spatial lattice considered in this work consisting of L_x  = 10 plaquettes along the x-direction with periodic boundary conditions and L_y = 1 plaquette along the y-direction with fixed boundary conditions. Degrees of freedom (DOFs) are (spin-$\frac{1}{2}$) hardcore bosons denoted by orange circles residing on the edges of a two-dimensional square array. The electric-field (link) operators on each edge are denoted by pink (blue) ellipses. Electric operators at site j build the Gauss-law operator G_j, and the four link operators covering edges of a plaquette build the magnetic operator W_□. The short-ribbon operator V_y, comprised of electric-field operators at the adjacent links α and β, is also shown. b The lattice in the dual formulation. Gauge-dependent operators (e.g., ${\sigma }_{{\ell }_{ĀA}}^z$) and the corresponding Gauss laws (e.g., $G_{ĀA}^{\mathrm{dual}}$ denoted by the pink ellipsoid) on the boundary between the system A and its complement $Ā$ (indicated by green thick dashed lines) remain unchanged from the original formulation. However, bulk operators within each subsystem are replaced with gauge-independent Pauli operators ${\mu }_i^{x,z}$ on each plaquette. The 12 quantum DOFs in this system are mapped to the trapped ion quantum computer, with DOFs shown as pink circles mapped to qubits 1 to 4 in the system and qubits 6 to 11 in the complement. Boundary DOFs, shown as yellow circles, are mapped to qubits 0 and 5.

For the thermalization study, we will bipartition the system to subsystem A and its complement $Ā$, resulting in two boundaries, which we call ${\ell }_{AĀ}$ and ${\ell }_{ĀA}$. To directly encode the degrees of freedom of the model described onto the qubits, one requires one qubit per link degree of freedom, or N_q = 3L_x qubits in total. To adapt the simulation to the available hardware, we employ a dual formulation of the model, described in the Methods section, that requires only N_q = L_x + 2 qubits. The dual formulation maps all bulk operators onto gauge-invariant ones (μ_i), except at the two boundaries where the original operators are retained, as shown in Fig. 2b. This mapping guarantees that the entanglement structure—the primary observable in our study—remains identical between the dual and the original lattice gauge theory. Similarly, both the original and the dual theory exhibit four symmetry sectors, which are introduced in the Methods section.

Digital quantum computation of real-time evolution

We first demonstrate the capability of our digital computer to compute the dynamics of the Z₂ LGT model, including initial-state preparation, time evolution, and subsequent measurements.

We initialize the system in a randomly chosen electric eigenstate that respects Gauss law ($\left.∣\right){\mathrm{\Psi }}_0\left(\right.⟩=\left.∣\right)\downarrow \downarrow \downarrow \uparrow \downarrow \downarrow \uparrow \uparrow \uparrow \downarrow \uparrow \uparrow \left(\right.⟩$, see the Methods section), and simulate its time evolution. Explicitly, we approximately obtain the state $\left.∣\right)\psi \left(t\right)\left(\right.⟩=U\left(t\right)\left.∣\right){\mathrm{\Psi }}_0\left(\right.⟩$ with $U\left(t\right)=e^{-it{H}^{\mathrm{dual}}}$, via a Trotterization protocol. Each Trotter step separately implements the non-commuting terms of H^dual, i.e., those diagonal in electric (z) or magnetic (x) bases. We implement four Trotter steps and set the coupling to g = 0.85. Our choice of coupling is such that magnetic and electric terms are of similar magnitude, making the model sufficiently non-integrable. Since the initial state is not an eigenstate of the full Hamiltonian, the system undergoes non-trivial time evolution when evolved under H^dual. For each evolution time t, we apply single-qubit gates to measure all qubits in either the Bloch x or the z bases and repeat each experiment N_shots = 500 times. Quantum circuits implementing initial-state preparation, Trotterized time evolution, and measurements are provided in the Methods section. Each Trotter step of evolution uses 12 variable-angle Mølmer-Søerensen (MS) gates.

The measurement results are used to compute single- and two-qubit observables in the Bloch x and z bases and are plotted in Fig. 3 for select observables. The statistical errors are determined by a standard bootstrap resampling analysis with N_boots = 1000 bootstrap samples. We have verified that our results are independent of N_boot. The expectation values of the physical operators (dark blue and red) follow the zero-parameter theory prediction (solid lines). The observed deviations are well-explained by ~10% systematic under-rotations of the XX gates, likely caused by calibration drift (see Supplementary Fig. S2). Importantly, the expectation values of gauge invariance—violating operators (i.e., those that do not commute with the Gauss laws at the two boundaries), shown in light blue, are consistent with zero within measurement errors. In Section 1 of the Supplementary Information, which includes refs. ^41,83,96,97, gauge-invariance violation is shown to be about 5% throughout the time evolution. The approach of observables to zero at late times is a coincidence associated with the time frame chosen—it neither signals thermalization nor decoherence. Additionally, we have verified robustness at larger gate angles by collecting data at even later times but using the same number of Trotter steps, see the Supplemental Figure 3.

[See PDF for image]

Fig. 3

Real-time observables in Z₂ LGT in (2 + 1)D.

Measured expectation values of single-qubit (a) and two-qubit (b) observables of qubits i = 0, …, 5 of the 10-plaquette theory with g = 0.85 following four steps of Trotterized time evolution from an initial state $\left.∣\right){\mathrm{\Psi }}_0\left(\right.⟩=\left.∣\right)\downarrow \downarrow \downarrow \uparrow \downarrow \downarrow \uparrow \uparrow \uparrow \downarrow \uparrow \uparrow \left(\right.⟩$. Each point represents data from N_shots = 500 experiments. Error bars indicate statistical uncertainties obtained using N_boot = 1000 bootstrap samples and do not include systematic errors. Red (blue) points correspond to measurements in the Bloch $z$-basis ($x$-basis); light-blue points indicate the expectation values of non-gauge—invariant observables, which are expected to be zero. Solid lines correspond to exact classically computed results.

Thermalization dynamics of Z₂ lattice gauge theory with trapped ions

To study thermalization, we randomly initialize the system in Gauss-law–respecting electric eigenstates. We then define a subsystem A consisting of L_A ≤ L_x/2 plaquettes. The reduced density matrix ρ_A(t) of A after evolution time t is related to the density matrix ρ(t) of the whole system as

2

${\rho }_A\left(t\right):={\mathrm{Tr}}_{Ā}\left[\rho \left(t\right)\right],\rho \left(t\right)=U\left(t\right){\rho }_{0}U{\left(t\right)}^{†},$where ${\rho }_0:=\mid {\mathrm{\Psi }}_0\left(\right.⟩\left.⟨\right){\mathrm{\Psi }}_0\mid$. Our main interest lies in the statistical properties of the Schmidt decomposition of the reduced density matrix of the subsystem A, given by ${\rho }_A\left(t\right)={\sum }_{\lambda }{p}_{\lambda }\left(t\right)\mid \lambda \left(t\right)\left(\right.⟩\left.⟨\right)\lambda \left(t\right)\mid$, and its associated EH^36,37

3

$H_{\mathrm{E}}\left(t\right):=-\log \left({\rho }_A\left(t\right)\right).$

The EH’s eigenspectrum, called the entanglement spectrum, is defined by $\left\{{\xi }_{\lambda }\left(t\right):=-\log \left(p_{\lambda }\left(t\right)\right)\right\}$.

The statistical properties of the entanglement spectrum directly inform the thermalization dynamics of the system. One statistical probe is the distribution, P(r), of the gap ratios, r_λ, of the entanglement spectrum³⁸,

4

$r_{\lambda }:=\frac{\min \left({\delta }_{\lambda },{\delta }_{\lambda -1}\right)}{\max \left({\delta }_{\lambda },{\delta }_{\lambda -1}\right)},$where δ_λ := ξ_λ − ξ_λ−1 ≥ 0 are the gaps between the eigenvalues ξ_λ of H_E and r ≡ {r_λ}. According to random matrix theory of Hermitian matrices, this quantity distinguishes chaotic from nonchaotic behavior depending on whether the distribution is centered away from zero (level repulsion) or centered around zero (uncorrelated levels)⁹¹.

In classical computation of a tractably small system, the state and its entanglement are readily available. In a quantum-simulation experiment, however, these must be inferred from measurements. To tackle this challenge, we constrain the operator content of the EH to obtain an approximate H_E, and utilize an EH-tomography scheme based on randomized measurements^{71, 72, 73, 74, 75–76}. Concretely, we perform a single-layer, single-qubit randomized-measurement⁹⁸ of the time-evolved state. At the end of the time evolution, this protocol applies one of $N_{\mathcal{U}}$ the different gates ${\mathcal{U}}_j:={\otimes }_{i=0,\cdots ,N_q-1}{u}_{i,j}$ for $j=1,\cdots ,N_{\mathcal{U}}$ consisting of independent single-qubit gates u_i,j sampled from a unitary 2-design⁹⁹. For each ${\mathcal{U}}_j$, N_shots experiments are performed, measuring all qubits in the Bloch z basis. The relative frequencies of the different bitstrings are then compared with the prediction based on an approximate EH inspired by the Bisognano-Wichmann theorem^89,90. Details are provided in the Methods section.

To study the evolution of the entanglement spectrum in the Z₂ LGT, we repeat the time-evolution experiments shown in Fig. 3 starting from 6 randomly chosen initial electric eigenstates that satisfy the Gauss laws. We numerically checked that the results remain valid for other randomly chosen states in this basis. For each evolution time t, symmetry sector, and initial state, we subsequently perform randomized-measurement tomography with g = 0.85, $N_{\mathcal{U}}=24$ bases, and N_shots = 750 bitstring measurements in each basis, and use these measurements to reconstruct the EH as described in the previous section. This procedure yields a set of gap ratios for each initial state, symmetry sector, and evolution time.

In Fig. 4(a), we plot in orange points the average gap ratio $⟨r⟩:={\sum }_{r}r\overline{P}\left(r\right)$ of the reconstructed EH as a function of the scaled evolution time gt. The plotted gap ratios are averaged over both the symmetry sectors and the 6 randomly chosen initial states. The black lines correspond to the predicted exact distributions following the Trotterized time evolution. The cyan lines correspond to the EH obtained from an optimal BW-inspired ansatz, where we numerically minimize the relative entropy (Kullback-Leibler divergence) between the exact state and the ansatz, corresponding to the limit of infinitely many measurements. A buildup of the level repulsion is discernible as the observed average gap ratio ⟨r⟩ evolves from  ≈ 0.4 predicted for a non-repulsive Poisson distribution (blue dashed line) towards  ≈ 0.6 characteristic of repulsive level statistics of a Gaussian Unitary Ensemble (GUE)³⁸. Three time regimes Early (E), Intermediate (I), and Late (L) are identified that correspond, respectively, to the evolution of the predicted average gap ratio toward the Poisson-distribution value, toward the GUE value, and to saturation at the GUE value. These regimes were determined before taking data, based on the numerically computed results (black solid lines in Fig. 4), hence constituting a prediction.

[See PDF for image]

Fig. 4

Statistics of the gap ratios of the spectrum of the entanglement Hamiltonian.

a Time evolution of the average gap ratio averaged over 6 randomly-drawn initial states and all symmetry sectors. The horizontal lines represent the averages for non-repulsive (Poisson, blue dotted) and repulsive (GUE, red dashed) distributions. Error bars indicate the standard deviation of the mean over the initial states and the symmetry superselection sectors of the reduced density matrix (see the Methods section). b Distribution of the entanglement-spectrum gap ratios, combined across 6 randomly-drawn initial states, all symmetry sectors, and all times in each of the regimes Early (E), Intermediate (I), and Late (L). A total of 840, 504, and 336 gaps are quantum-computed based on 10⁶ total shots, and their average is shown in orange. Simulated Bisognano-Wichmann results in the limit of infinite measurements are shown in cyan, and the exact predicted distributions in black. Blue-dotted and red-dashed curves represent Poisson and GUE distributions.

In Fig. 4b, we plot, for the three ranges of evolution times, the corresponding normalized distribution of the gap ratios, $\overline{P}\left(r\right)$, combined over 6 random initial states and 4 symmetry sectors. We observe a clear transition from early-time non-repulsion (Poisson distribution, blue dotted line) in the early regime (E) to level repulsion (GUE, red) in the late regime (L). At intermediate times, a distribution is observed between the initial absence of level repulsion and the subsequent emergence of level repulsion. Error bars and bands denote the variance resulting from averaging over symmetry sectors and initial states. The predicted exact distribution (black line) exhibits a sharper peak around zero at the earliest times due to the reduced density matrix not being of full rank. Conversely, the BW-inspired ansatz (cyan line) generally parametrizes a full-rank matrix unless couplings are very finely tuned, resulting in an overestimation of level repulsion in the initial stages. The time scale required to achieve a distribution consistent with GUE predictions is approximately the same for all randomly chosen initial states. This timescale is expected to scale with subsystem size as was first observed in ref. ³⁹ in a (1 + 1)D spin model, see also Supplemental Figure 5.

The main effects of Trotterization versus exact time evolution are two-fold. First, the evolution is governed by an effective instead of the exact Hamiltonian; however, we draw our BW parameterization ansatz from the exact Hamiltonian. Second, at later times and with large Trotter-step sizes (not shown here), a recurrence of the initial state is observed, and the gap-ratio statistic becomes ‘un-chaotic’ again, see Supplemental Figure 8.

A complementary view of the statistics of the entanglement spectrum that captures global correlations in the level distribution is afforded by the entanglement spectral form factor (ESFF)^40,91,

5

$\mathcal{F}\left(\theta ;t\right):=⟨\frac{1}{{\mathcal{R}}_A^2\left(t\right)}{\sum }_{\lambda ,{\lambda }^{\prime }}{e}^{i\theta \left[{\xi }_{\lambda }\left(t\right)-{\xi }_{{\lambda }^{\prime }}\left(t\right)\right]}⟩,$where ⟨ ⋅ ⟩ denotes the average over initial states $\left.∣\right){\mathrm{\Psi }}_0\left(\right.⟩$. Here, ${\mathcal{R}}_A:={lim}_{\alpha \to 0}\exp \left\{\frac{1}{1-\alpha }\log \left({\sum }_{\lambda }{p}_{\lambda }^{\alpha }\right)\right\}$⁴⁰ is the effective rank of H_E, whose value depends on the state and lies in the range $\left[1,2^{L_A}\right]$. Only the ${\mathcal{R}}_A$ lowest levels of the EH are included in the analysis. Small values of θ probe global correlations among eigenvalues, which are not universal, while large values probe correlations between eigenvalues that are progressively closer and where GUE predictions apply. In Fig. 5, we plot the ESSF reconstructed from the data used in Fig. 4. The three panels of the plot correspond to the three time regimes discussed in relation to the EGRD evolution. The displayed curves are averaged over initial states, symmetry sectors, and each of the three time ranges. Starting from an initial flat behavior of the ESSR as a function of θ in panel (E), our data in panels (I) and (L) clearly show the buildup of a ramp-plateau structure, indicating ergodic behavior and implying quantum chaos⁴⁰. In panel (L), we show a fit (purple dotted line) to the observed ramp, indicating a power-law behavior θ^κ where κ = 0.6 ± 0.2, with the fit error determined by changing the fit regime. The exponent depends on the definition of the ESSF in the presence of symmetries. Our definition leads to a value that is smaller than the one from RMT. It is consistent with the results obtained from the numerical analysis of a significantly larger system in Supplementary Information section. 2.

[See PDF for image]

Fig. 5

Entanglement spectral form factor.

The average entanglement spectral form factor across various initial states, symmetry sectors, and three distinct regimes Early (E), Intermediate (I), and Late (L) identified in Fig. 4. The figures display exact results (black curves), infinite-measurement outcomes (cyan), and quantum-computed experimental data (orange). In panel (L), a purple dotted line indicates a power-law fit of the ramp observed in the quantum-computed data. Our normalization ensures $⟨\mathcal{F}\left(0\right)⟩=1$, with the plateau occurring at $⟨\mathcal{F}\left(\infty \right)⟩=1/d_s$, where d_s denotes the dimension of a symmetry block, see the Methods section for details. Shaded areas indicate the standard deviation over initial states, symmetry sectors, and times. Data included in the fit are shown as circles.

Our results demonstrate that, using observables that reveal universal statistical or global features of the entanglement spectrum, the onset of chaotic behavior can be robustly traced in the entanglement dynamics of arbitrary initial states. Importantly, time digitization and experimental infidelities in current systems do not significantly impact qualitative features of thermalization dynamics.

Discussion

Using a digital trapped-ion quantum computer, we observe early-stage thermalization dynamics in a Z₂ lattice gauge theory in 2 + 1 dimensions by generating nonequilibrium quantum states and measuring their entanglement structure. As performing quantitative state tomography is extremely challenging and resource intensive, we focused on a simpler question: can a randomized-measurement–based entanglement-Hamiltonian tomography recover universal properties of the EH that indicate quantum chaos?

The term ‘universal’ refers to the fact that, on timescales of order one in units of inverse coupling, the entanglement spectrum of thermalizing systems develops correlations between neighboring levels that are described by a universal (GUE) statistical ensemble, independent of the initial conditions. For integrable systems, the gap ratio may not exhibit universal behavior. While a Poisson distribution emerges in certain cases⁴³, the distribution generally remains sub-Poissonian, may oscillate, or can be ill-defined if the rank of the reduced density matrix remains small throughout the evolution.

We experimentally determine the entanglement spectrum of time-dependent quantum states and use two properties of this spectrum, namely the entanglement gap-ratio distribution [Eq. (4)] and the entanglement spectral form factor [Eq. (5)], as indicators of the emergence of quantum chaos. Our data indicate that both quantities behave as expected: the initially non-repulsive level distribution transitions to a repulsive one. Likewise, a ramp-plateau feature of the ESFF is initially absent but builds up with time. The emergence of quantum chaos in the EH is a prerequisite for thermalization, which takes a parametrically longer time to set in. This separation is discussed in ref. ³⁹: while the timescale for the EH to become chaotic is governed by the Lieb-Robinson velocity describing local-operator spread, entanglement saturation—and thus thermalization—is much slower^100,101.

The key technique enabling our analysis is a Bisognano-Wichmann–inspired ansatz for the EH, which allows the parameterization of a nonequilibrium state using a polynomial number of parameters. While this classical parameterization optimally reproduces the observed (randomized) measurement results, in Supplementary Information 2, we show that the higher-lying part of the EH spectrum is not quantitatively recovered. Indeed, the measurement cost for precise state tomography still scales exponentially with the subsystem size. Fortunately, to detect the presence of quantum chaos, one needs to only distinguish repulsion from non-repulsion in the EGRD or identify a ramp-plateau structure in the ESFF. Our work indicates that selecting a sufficiently small subsystem of a much larger, potentially classically nonsimulable system can sufficiently constrain these observables.

Other signatures of thermalization can also be examined, such as the saturation of local observables or the entanglement entropy to their thermal values. However, beyond the longer simulation times required, several caveats must be considered when using such quantities to assess thermalization. For instance, finite-volume effects affect local observables at late times, and certain local observables do not accurately distinguish thermalization from integrable dynamics; see, e.g., ref. ¹⁰².

Our randomized-measurement procedure remains classically simulable because the lattice sizes we considered are manageable with exact diagonalization. Our procedure is also heavily tailored at minimizing the computational load on the quantum computer, avoiding deep quantum circuits and the need for error mitigation, at the cost of classical post-processing. To extend our approach to larger systems, several steps can be implemented in future work:

• Our study focuses on gauge theories. As demonstrated in the example studied, gauge theories possess a highly intricate Hilbert-space structure, shaped by an extensive set of local constraints, i.e., Gauss’s laws, that govern their dynamics. These constraints also give rise to a non-trivial entanglement structure. For example, resolving the symmetry properties of the reduced density matrix was crucial to identifying the presence of quantum chaos. Without knowledge of this structure, simply computing the gap-ratio distribution of all eigenvalues of the EH would have shown a Poisson, or more singular distribution, falsely implying the absence of quantum chaos because eigenvalues from different symmetry blocks are uncorrelated. It would be interesting to extend this study to continuous groups, as well as non-Abelian gauge theories, which can pave the way toward studying thermalization physics of the Standard Model.

• The employed single-qubit one-layer randomization strategy is symmetry ignorant, i.e., it randomizes regardless of the known symmetry structure of the subsystem density matrix ρ_A(t). Although our scheme is tomographically complete and enables state reconstruction along with its symmetries, it is inefficient. Symmetry-conscious protocols, for the LGT in this work and for other models, have been developed^85,103,104, providing tomographically complete circuits with polynomial depth, and could be adopted in future work.

• A major methodological uncertainty is the BW-inspired parameterization of the EH, which is entirely heuristic. While there is research exploring the operator content of the EH of ground, excited, and thermal states¹⁰⁵, further investigations are needed into the applicability and limitations of EH-based schemes for far-from-equilibrium states. Generally, the BW(-inspired) ansätze describe the low-energy regime of the EH well, and incorporating progressively more nonlocal terms improves convergence into the bulk⁷³, which the EGRD, and to a lesser extent the ESFF, predominantly rely on. Notably, our analysis reveals that despite the quantitative discrepancy in the bulk, the statistical distribution of the EH appears to be accurately reproduced. This suggests that it may be unnecessary to perform precise state tomography when one’s interest is solely in statistical properties.

• Trotterization corresponds to time evolution with an effective, rather than the desired, Hamiltonian. We observe that employing too few Trotter steps leads to poor convergence of our optimization procedure for the BW-inspired ansatz. In this case, more nonlocal operators must be included in the ansatz; constraining these becomes challenging, especially when relying on a limited number of measurements. Since using a large number of Trotter steps increases experimental errors, the optimal Trotterization and its interplay with the BW ansatz should be further investigated.

• Detailed comparison between the emulator and the experimental data (see Supplementary Information) highlights the influence of device errors, primarily for the time evolution and tomography steps of our algorithm. Developing a comprehensive error model to predict how errors impact the BW-EHT analysis for fully digital computations is a complex task that we have not investigated in detail. However, this question has been explored in ref. ⁷³ for an analog-digital scheme; a comprehensive error model for the specific computing platform used in our work can be found in ref. ¹⁰⁶. The dominant errors in our study are Z-flip errors and under-rotations due to mechanical motion of the ions¹⁰⁷. Our tomography protocol is especially susceptible to small single-qubit rotation errors during randomization⁹⁸. Recently demonstrated sympathetic cooling during circuit execution¹⁰⁷ would allow the needed measurement circuits to be executed with higher fidelity. The same technique could extend our study to larger systems, later times, and more Trotter steps, e.g., allowing us to directly test the eigenstate thermalization hypothesis in a digital quantum-computing setup. Once system sizes exceed classical simulability, verification becomes an important problem; for an overview, see, for instance, ref. ¹⁰⁸.

In summary, our results demonstrate that entanglement structure is a measurable quantity in present-day LGT quantum-simulation experiments, and illustrate the potential value of our approach to probe thermalization dynamics and its robust universal features in strongly coupled isolated quantum many-body systems. A compelling future direction is to extend the investigation to later times, to probe aspects such as pre-thermalization^{10,12,109, 110–111} or fluctuation-dissipation relations^{112, 113, 114–115}, once experimental capabilities permit. This would also allow for probing the applicability of the Eigenstate Thermalization Hypothesis, e.g., in systems with non-Abelian symmetries¹¹⁶ and other gauge theories^117,118. Quantum many-body scars have been identified in a slightly different lattice geometry of this model and may serve as a means of verifying quantum simulation experiments beyond the reach of classical computation¹¹⁹. Furthermore, the experimental and theoretical tools of this study can be applied in a number of other applications, including obtaining thermodynamic quantities such as work and heat exchanged during nonequilibrium processes¹²⁰, and detecting phases of matter, including topological phases⁸⁵, in quantum-simulation experiments. Ultimately, strategies such as those presented in this work can reveal the role of entanglement in QCD thermalization and may shed light on the short thermalization time scales observed in heavy-ion collisions¹⁷. Such studies, nonetheless, will require large-scale fault-tolerant quantum computers. In the meantime, studies in simpler gauge theories, such as in confining lower-dimensional models, can potentially provide the first clues.

Methods

In this section, we provide further details on the model studied in this work, our experimental trapped-ion setup, the entanglement-Hamiltonian tomography we employed, and the quantum-circuit representation of our full quantum-simulation protocol.

Generalized Ising duality

To adapt the simulation to the available hardware, we consider a dual formulation of the model that requires only N_q = L_x + 2 qubits. In the dual formulation, all operators, except those containing operators having support at the boundaries between the subsystems, labeled ${\ell }_{ĀA}$ and ${\ell }_{AĀ}$, are entirely represented using gauge-invariant (Ising dual) spin-$\frac{1}{2}$ variables which act solely within the physical Gauss-law subspace; they are denoted by ${\mu }_i^{x,z}$ Pauli matrices [see purple circles in Fig. 2b]. Concretely, ${\mu }_i^x\to W_i={\prod }_{\ell \in {\square }_i}{\sigma }_{\ell }^x$, while electric variables are represented by ${\mu }_i^z{\mu }_{i-1}^z\equiv {\sigma }_{\ell }^z$ where ℓ is the link between plaquettes i and i − 1. The corresponding “bulk” Hamiltonian for the subsystem A then reads

6

$H_A^{\mathrm{dual}}={\sum }_{i\in A_{\mathrm{bulk}}}{\mu }_i^x+g\left({\sum }_{⟨i,j⟩\in A}{\mu }_i^z{\mu }_j^z+\kappa {\sum }_{i\in A}{\mu }_i^z\right),$where κ := 1 + V_y. Here, i ∈ A_bulk in the sum of ${\mu }_i^x$ indicates that the original plaquette □_i is entirely in A (i.e., does not touch the boundaries). ⟨i, j⟩ ∈ A denotes nearest-neighbor bulk-spin pairs i, j in A. Finally, the sum of ${\mu }_i^z$ runs over all bulk-spin indices i. The terms in the Hamiltonian acting on the complement are defined identically. Importantly, at the boundaries between the subsystems, ${\ell }_{ĀA}$ and ${\ell }_{AĀ}$ [see the orange circles in Fig. 2b], the gauge-variant variables of the LGT, as described by Eq. (1), are retained. Denoting the Ising-spin index at the A side of one of the boundaries to be k, the Hamiltonian terms coupling the subsystems at this boundary are

7

$H_{ĀA}^{\mathrm{dual}}={\mu }_k^x{\sigma }_{{\ell }_{ĀA}}^x+{\mu }_{k-1}^x{\sigma }_{{\ell }_{ĀA}}^x+g{\sigma }_{{\ell }_{ĀA}}^z,$with a similar definition for the other boundary. Note that on the plaquette k containing one boundary link, ${\mu }_k^x{\sigma }_{{\ell }_{ĀA}}^x\equiv W_k$, where ${\mu }_k^x\equiv {\prod }_{\ell \in {\tilde{\square }}_k}{\sigma }_{\ell }^x$, with ${\tilde{\square }}_i$ referring to all but non-boundary links of plaquette k, and similarly for the plaquette k − 1. The Gauss laws at the boundaries are not eliminated by the duality. The two Gauss laws in the dual model that are independent (one at each boundary) are

8

$G_{ĀA}^{\mathrm{dual}}={\mu }_k^z{\sigma }_{{\ell }_{ĀA}}^z{\mu }_{k-1}^z,$and similarly for the boundary at ${\ell }_{AĀ}$. The dual Hamiltonian of the model is the sum of the terms above:

9

$H^{\mathrm{dual}}=H_A^{\mathrm{dual}}+H_{Ā}^{\mathrm{dual}}+H_{ĀA}^{\mathrm{dual}}+H_{AĀ}^{\mathrm{dual}}.$

This formulation ensures that all gauge-invariant variables of the LGT and its dual have the same expectation values. Importantly, by maintaining the Gauss laws at the boundaries, the entanglement properties of the dual formulation are identical to those of the LGT⁴¹. This is a subtle but crucial distinction from the standard Ising duality¹²¹, which does not preserve the entanglement structure.

The entanglement structure depends on the symmetries of the reduced density matrix, which stem from the remaining Gauss laws. Specifically, [S_j, ρ_A(t)] = 0, where in LGT variables, $S_j:={\prod }_{{\ell }_j\in A}{\sigma }_{{\ell }_j}^x$, with ℓ_j being the two links originating from a boundary site with lattice coordinates j = (j_x, j_y) within A. This is a direct consequence of the state being a Gauss-law eigenstate⁴¹. In addition, [V_x,A, ρ_A(t)] = 0 where the long-ribbon operator within A is defined as $V_{x,A}:={\prod }_{\ell \in \mathcal{P}}{\sigma }_{\ell }^x$ where ${\sigma }_{\ell }^x$ acts on the original degrees of freedom along the path $\mathcal{P}$ connecting both boundaries within A, see Fig. 6. Not all symmetry operators are independent. For the quasi-1D chain of plaquettes shown in Fig. 2, only one of the two S_j is independent at each boundary. We place one of the boundaries at sites (0, 0) and (0, 1) and the other at sites (L_A, 0) and (L_A, 1). We then choose the two independent S_j operators to be S_(0, 0) and $S_{L_A,0}$, henceforth referred to as $S_{AĀ}$ and $S_{ĀA}$, respectively. Note that, when acting on a Gauss-law–respecting state $V_{x,A}=S_{AĀ}{S}_{ĀA}$. Consequently, there are four independent symmetry blocks ρ_A(t) ≡ ⨁_s=1,⋯ ,4ρ_A,s(t) that will play a role in our analysis. In the dual formulation, these symmetry operators are simply $S_{AĀ}^{\mathrm{dual}}\equiv {\mu }_k^z{\sigma }_{{\ell }_{ĀA}}^z$ and $S_{AĀ}^{\mathrm{dual}}\equiv {\mu }_{k\prime }^z{\sigma }_{{\ell }_{ĀA}}^z$ where k and $k^{\prime }$ are the first and last plaquettes within A. We denote the eigenvalue {1, 1}, {1, − 1}, { − 1, 1}, and { − 1, − 1} sectors of $\left\{S_{ĀA}^{\mathrm{dual}},S_{AĀ}^{\mathrm{dual}}\right\}$ operators with labels s = 1, 2, 3, 4, respectively.

[See PDF for image]

Fig. 6

Symmetry operators.

Symmetry operators of the reduced density matrix of a state evolving under Gauss law—respecting dynamics. The thin dashed green line represents the path $\mathcal{P}$ along which the Pauli operators in V_x,A act on.

Experimental setup

Our experimental system consists of a linear chain of fifteen ¹⁷¹Yb⁺ ions trapped along the x-axis of a microfabricated ion trap ⁶⁶ and spaced by approximately 3.7 μm. The $\mid \downarrow \left(\right.⟩$ and $\mid \uparrow \left(\right.⟩$ spin states are respectively mapped to the first-order magnetic-field–insensitive hyperfine $\mid F=1,m_F=0\left(\right.⟩$ and $\mid F=0,m_F=0\left(\right.⟩$ states of the ²S_1/2 ground electronic manifold of the ions. Here, we use the quantum-information convention for the spin states, i.e., $\left.∣\right)\uparrow \left(\right.⟩\equiv \left.∣\right)0\left(\right.⟩$ and $\left.∣\right)\downarrow \left(\right.⟩\equiv \left.∣\right)1\left(\right.⟩$. The spin states are initialized in the $\mid \uparrow \left(\right.⟩$ state by optical pumping and are measured by coupling to the excited ²P_1/2, $\left.∣\right)F=0,m_F=0\left(\right.⟩$ state using a 369-nm laser, wherein the presence (absence) of emitted photons differentiates between the bright state $\mid \downarrow \left(\right.⟩$ and the dark state $\mid \uparrow \left(\right.⟩$ with an average 0.3% infidelity, limited by off-resonant photon scattering.

The ions are individually addressed by an equispaced array of 1 μm-diameter, 355 nm, laser beams oriented perpendicular to the trap surface (z axis), together with a global 300 μm × 30 μm beam oriented along the y-axis, which is parallel to the trap surface (see Fig. 1a and ref. ⁵). These beams drive stimulated Raman transition via the ²P_1/2 and ²P_3/2 excited electronic states, where a photon is absorbed from the global beam and emitted into an individual beam to flip the qubit from $\mid \downarrow \left(\right.⟩$ to $\mid \uparrow \left(\right.⟩$ and vice versa. All spin-manipulation light is produced by a single pulsed laser, modified to control its repetition rate to null 4-photon Stark shifts. The phase and amplitude of the 355-nm beams is controlled by single-channel and 32-channel acousto-optic modulators provided by the L3 Harris Corporation.

To minimize addressing errors, crosstalk, and stray coupling to axial motion of the ions¹⁰⁷, single-qubit operations are realized using compound SK1 pulses¹²² with Gaussian-shaped sub-segments, with typical infidelities of 0.2%. Entangling operations between any two qubits are realized via variable-angle pairwise Mølmer-Sørensen (MS) gates¹²³ employing laser-induced state-dependent forces on the  ≈ 3 MHz motional modes of the ion chain oriented along the y + z axis. Robust decoupling from the ion motion at the end of each entangling gate is accomplished using amplitude-modulated, detuning-robust pulse waveforms¹²⁴. Typical non-unitary errors of fully-entangling MS gates are 1% and consist predominantly of gate-angle errors and Z-flips of the individual qubits¹²⁵.

Entanglement Hamiltonian tomography

Our randomized-measurement scheme applies one of $N_{\mathcal{U}}$ different gates ${\mathcal{U}}_j:={\otimes }_{i=0,\cdots ,N_q-1}{u}_{i,j}$ for $j=1,\cdots ,N_{\mathcal{U}}$ consisting of independent single-qubit gates u_i,j sampled from a unitary Haar-design⁹⁹. Concretely, randomization is via single-qubit random circuits, $\mathcal{U}{:=⨂}_i{u}_i$ where u_i is the following single-qubit unitaryand for each qubit, the angles ${\gamma }_i^1$, ${\gamma }_i^2$, ${\gamma }_i^3$ are drawn according to a circular unitary ensemble (an overall phase is ignored) following refs. ^76,98. A drawback of this approach is that it does not maintain the symmetry structure of ρ_A, as it randomizes over the whole Hilbert space instead of each symmetry block of ρ_A. Such a symmetry-ignorant randomization mixes different symmetry sectors, resulting in an outcome that is inaccessible to any physical time-evolved quantum state. We choose it, nonetheless, to avoid the larger circuit depth associated with a symmetry-conscious randomization such as that proposed in ref. ⁸⁵. We note that a symmetry-conscious scheme would significantly reduce the measurement cost, as the sampling cost would only scale as the size of the symmetry block instead of the Hilbert-space size. It would thus simplify the classical optimization significantly. However, because the single-qubit scheme is tomographically complete, one can still reconstruct, approximately, the symmetry structure from the data.

The classical post-analysis of the randomized measurement results consists of comparing the relative frequencies of the different bitstrings b to the prediction based on an approximate EH inspired by the Bisognano-Wichmann (BW) theorem^89,90. Explicitly, we assume that H_E(t) is a linear combination of k-local terms ${\mathcal{O}}_i$, i.e.,

10

$H_{\mathrm{E}}\left(t;\left\{{\beta }_i\right\}\right)={\sum }_i{\beta }_i\left(t\right){\mathcal{O}}_i.$

To find a suitable set of operators {O_i}, we proceed as follows. Starting from the operator content of the physical Hamiltonian in Eq. (1), containing at most 4-local terms (2-local in the dual formulation), new operators are generated iteratively by forming non-trivial commutators from the existing set. This process can be halted after two iterations (‘commutators of commutators’), resulting in a maximum of 7-local operators within the LGT framework (or 3-local in the dual formulation).

While the BW theorem hints at an optimal choice of operators ${\mathcal{O}}_i$ for ground states, our ansatz is heuristic, as a systematic ansatz for non-equilibrium states is not known. Our criteria for the operators are that (i) they are local, involving operators with support on at most two neighboring plaquettes, (ii) they are compatible with the symmetries of the reduced state ρ_A (which will be discussed below), and (iii) they are independent, meaning they cannot be transformed into each other through Gauss laws when acting on a physical state (although they are certainly not independent in the algebraic sense).

The operators we use are pictorially represented in Fig. 7, plus all operators that are generated from those depicted upon translation. Here, for the sake of generality, we represent them in terms of Z₂ variables, albeit in practice, they are represented in the dual formulation in our algorithm. Figure 8 shows a special operator; this operator (and related examples) relates operators within the subsystem to operators in the complement via Gauss laws. Further, it commutes with ρ_A, and is hence connected to the symmetry blocks of ρ_A. All terms in Fig. 7 are obtained by considering the physical-Hamiltonian operators and by recursively commuting the physical-Hamiltonian operators, stopping at two recursions (i.e., commutators of commutators). Other selection strategies are also feasible.

[See PDF for image]

Fig. 7

Operator set for the entanglement-Hamiltonian ansatz.

Representation of the operator set in the entanglement-Hamiltonian ansatz used in our study, encompassing operators positioned inside the subsystem. All operators are Hermitian. Electric fields along the x-direction are interdependent (since V_y = 1), thus not treated as independent operators. For the specific subsystem with N_A = 4, our ansatz incorporates a total of 73 distinct operators. While this is close to the subsystem Hilbert-space dimension, for a larger system, the number of operators scales linearly with the degrees of freedom within the subsystem.

[See PDF for image]

Fig. 8

The operator related to the symmetry sectors of the ansatz state.

Because of the (remaining) Gauss laws, some operators outside of the subsystem are identical to operators within its complement. The operator shown, and a similar operator placed on the other boundary, commutes with all other operators in the BW-inspired ansatz shown in Fig. 7. Therefore, their (common) eigenspace are the symmetry sectors of the reduced density matrix.

Given the ansatz and following ref. ⁷⁴, we then minimize the functional

11

${⟨{\sum }_b{\left[P_{\mathcal{U}}\left(b\right)-\mathrm{Tr}\left[{\mathcal{U}}^{†}\mid b\left(\right.⟩\left.⟨\right)b\mid \mathcal{U}{\rho }_A\left(t;\left\{{\beta }_i\right\}\right)\right]\right]}^2⟩}_{\mathcal{U}}$where $P_{\mathcal{U}}\left(b\right)$ is the probability to measure a bitstring b in the basis determined by $\mathcal{U}$ and

12

${\rho }_A\left(t;\left\{{\beta }_i\right\}\right)=\frac{e^{-H_{\mathrm{E}}\left(t;\left\{{\beta }_i\right\}\right)}}{\mathrm{Tr}\left[e^{-H_{\mathrm{E}}\left(t;\left\{{\beta }_i\right\}\right)}\right]}.$is the normalized reduced density matrix parameterized by H_E. Here, ${⟨\cdot ⟩}_{\mathcal{U}}$ is the average over random circuits. Note that optimizing Eq. (11) effectively weighs more favorably the largest Schmidt eigenvalues of ρ_A because they, on average, contribute the most to any random observable. Because of this, the optimization more accurately reproduces the low-energy part of H_E. The optimization is performed using MATLAB’s¹²⁶ non-linear FMINCON optimization package with the ‘sqp’ algorithm¹²⁷. Convergence and uniqueness of the obtained minimum have been cross checked for several data sets using MATLAB’s GLOBALSEARCH routine with default parameters¹²⁸. All EH couplings are confined to the range β_i ∈ [ − 50, 50]. We check explicitly that the routine does not come close to the boundary of the parameter regime. The EH that is obtained in this way is projected into symmetry sectors, which we then separately analyze. For this projection, the symmetry sector can be read off directly from the row and column numbers of the EH when they are interpreted as binary. The corresponding subblock H_E,s is then selected for further analysis, where s labels one symmetry sector. For the EGRD analysis, we analyze the gap-ratio distribution separately for every sector, then combine the distributions and re-normalize the total distribution. For the ESFF analysis, we average over the sectors. Finally, EGRD and ESFF, shown in the main text and appendices, involve a combination/average over initial states. The error band and the error bars in the plots are the standard deviation for the symmetry-sector and initial-state averages. A systematic study of the performance of the EHT approach can be found in Supplementary Information sec. 3.

Circuit representation for state preparation, evolution, and measurement

In this Section, we provide circuits employed to experimentally realize state preparation, time evolution, and randomized measurement in the Z₂ LGT in (2 + 1)D within the dual Ising formulation.

Our basic circuit design is depicted in Fig. 9. We label qubits with indices q = 0, …, L − 1, with L = 12. Qubit q = 0 corresponds to the spin at ${\ell }_{ĀA}$ positioned at the interface between subsystem A and its complement, identifiable by an orange circle in Fig. 2 (b). Qubits q = 1, …, 4 represent the succeeding dual variables arranged counterclockwise at the center of plaquettes denoted by purple circles in Fig. 2 (b). Following this sequence, qubit q = 5 resides at the opposite boundary marked as ${\ell }_{AĀ}$, and so forth.

[See PDF for image]

Fig. 9

Overview of the circuits employed in the simulations of this work.

a Initial-state preparation, which involves setting a Gauss-law--respecting electric eigenstate in the Bloch z basis. The shown circuit corresponds to one such initial state. b Trotterized time evolution, consisting of magnetic interactions (in the x basis) comprised of single- and two-qubit rotations in the dual formulation. The electric part of the Hamiltonian evolution is diagonal consisting of Z and ZZ rotations. All MS gates are implemented according to Eq. (15) depending on the employed time-evolution step δt. c Randomized-measurement circuits. ${\gamma }_i^{1,2,3}$ at qubit i = 0, 1, ⋯  , 5 are drawn from a circular unitary ensemble. The qubits 0 and 5 marking the two boundaries of the periodic lattice are drawn in orange. The qubit index of the rotation gates are dropped to reduce clutter and can be deduced from the qubit(s) they act on.

Working in the electric eigenbasis and starting from an all-up spin state, we first perform single-qubit σ^x bit-flip operations to select a randomly chosen product state in the σ^z basis that is consistent with Gauss laws, i.e., $G_j\left.∣\right)\psi \left(0\right)\left(\right.⟩=\left.∣\right)\psi \left(0\right)\left(\right.⟩$ for all j. Recall that in the dual formulation, there are only two Gauss-law operators each located at one of the subsystem’s boundaries: $G_{ĀA}^{\mathrm{dual}}={\mu }_{11}^z{\sigma }_0^z{\mu }_1^z$ and $G_{AĀ}^{\mathrm{dual}}={\mu }_4^z{\sigma }_5^z{\mu }_6^z$. Importantly, any spin configuration in either bulk, A and $Ā$, is physical, as in the dual formulation, only gauge-invariant degrees of freedom are kept in the bulk.

For time evolution, we utilize a first-order Trotter scheme with the time-evolution operator

13

$U\left(t\right)=e^{-it{H}^{\mathrm{dual}}}={\prod }_{\delta t}U\left(\delta t\right).$

We take δt to be a variable time step while keeping the Trotter depth fixed. In the simulations conducted in this work, t/δt = 4, resulting in states and observables closely resembling the exact time-evolved states at early times. However, for late times, Trotter effects become more pronounced.

U(δt) in Eq. (13) can be written as U(δt) = ∏_{a∈X,Z,XX,ZZ}U_a with $U_a:=e^{-i\delta t{H}_a^{\mathrm{dual}}}$; $H_a^{\mathrm{dual}}$ are the respective 1- and 2-local operators of the dual Hamiltonian in Eq. (9), sorted in X and Z operations, as well as single- and two-qubit entangling gates. While XX, X, and Z are native operations, ZZ entangling operations are realized via basis transformation and usage of the native MS gate.

14

$R_{ij}^{zz}\left(\alpha \right):=R_i^y\left(-\frac{\pi }{2}\right)R_j^y\left(-\frac{\pi }{2}\right)R_{ij}^{xx}\left(\alpha \right)R_i^y\left(\frac{\pi }{2}\right)R_j^y\left(\frac{\pi }{2}\right).$

Here, $R_i^{x/y/z}\left(\alpha \right):=e^{-i\frac{\alpha }{2}{\sigma }_i^{x/y/z}}$ and $R_{ij}^{xx}\left(\alpha \right):=e^{-i\alpha {\sigma }_i^x{\sigma }_j^x}$. The experimental errors of our MS gates increase with the absolute value of the gate angle. To minimize the gate error, all MS gates are optimized as follows. First, all angles are mapped to the regime α ∈ [ − π, π]. Within this range, we make the following substitution

15

$R_{ij}^{xx}\left(\alpha \right)\to \left\{\begin{array}{cc}{R}_{ij}^{xx}\left(\alpha \right)\hfill & \mathrm{if}\mid \alpha \mid \le \frac{\pi }{4}\hfill \\ R_{ij}^{xx}\left(\alpha +\frac{\pi }{2}\right)R_i^x\left(\pi \right)R_j^x\left(\pi \right)\hfill & \mathrm{if}\mid \alpha \mid \in \left(\frac{\pi }{4},\frac{3\pi }{4}\right]\hfill \\ \hfill & \&\alpha <0\hfill \\ R_{ij}^{xx}\left(\alpha -\frac{\pi }{2}\right)R_i^x\left(\pi \right)R_j^x\left(\pi \right)\hfill & \mathrm{if}\mid \alpha \mid \in \left(\frac{\pi }{4},\frac{3\pi }{4}\right]\hfill \\ \hfill & \&\alpha >0\hfill \\ R_{ij}^{xx}\left(\alpha +\pi \right)\hfill & \mathrm{if}\alpha <-\frac{3\pi }{4}\hfill \\ R_{ij}^{xx}\left(\alpha -\pi \right)\hfill & \mathrm{if}\alpha >\frac{3\pi }{4}\hfill \end{array}\right)$so that all MS operations are restricted to the range $\mid \alpha \mid \le \frac{\pi }{4}$.

Acknowledgements

N.M. thanks H. Froland, A. Polkovnikov, M. Savage, M. Srednicki, X. Yao, T. Zache, P. Zoller, and the participants of the InQubator for Quantum Simulation (IQuS) workshop “Thermalization, from Cold Atoms to Hot Quantum Chromodynamics ”(https://iqus.uw.edu/events/iqus-workshop-thermalization/) at the University of Washington in September 2023 for many valuable discussions leading to this work. M.C. thanks L. Feng for valuable discussions and help with the experimental setup leading to this work. N.M. acknowledges funding by the Department of Energy (DOE), Office of Science, Office of Nuclear Physics, IQuS (https://iqus.uw.edu), via the program on Quantum Horizons: QIS Research and Innovation for Nuclear Science under Award DE-SC0020970. Z.D., M.C., and T.W. were supported by the National Science Foundation’s Quantum Leap Challenge Institute for Robust Quantum Simulation under Award OMA-2120757. Z.D. further acknowledges support by the DOE, Office of Science, Early Career Award DE-SC0020271. This work is further supported by a collaboration between the US DOE and other Agencies. This material is based upon work supported by the DOE, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator.

Author contributions

M.C., Z.D., O.K., and N.M designed the research; N.M. performed the numerical computations, designed the quantum circuits, and conducted the data analysis; M.C. and T.W. performed experiments and collected the data; M.C., Z.D., and N.M. contributed to the manuscript, with input from all authors.

Peer review

Peer review information

Nature Communications thanks the anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Data availability

The data displayed in all figures are publicly accessible at this repository. The raw data generated in this study have been deposited in Duke University’s Box cloud storage database. The raw data are available under restricted access subject to intellectual property laws, non-disclosure agreements, and trade-secret protections, in accordance with all applicable laws, regulations, agreement terms and conditions, and the policies and directives of the authors’ sponsors and affiliated institutions. Access can be obtained by contacting Marko Cetina, Zohreh Davoudi, or Niklas Mueller.

Code availability

The code and circuit designs presented in this manuscript are publicly accessible at this repository https://gitlab.com/Niklas-Mueller1988/thermalization-z2-lgt.

Competing interests

M.C. is a co-inventor on patents that are licensed from the University of Maryland to IonQ, Inc. The remaining authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-60177-7.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Eisert, J; Friesdorf, M; Gogolin, C. Quantum many-body systems out of equilibrium. Nat. Phys.; 2015; 11, 124.1:CAS:528:DC%2BC2MXhvFCgurc%3D [DOI: https://dx.doi.org/10.1038/nphys3215]

2. Schreiber, M et al. Observation of many-body localization of interacting fermions in a quasirandom optical lattice. Science; 2015; 349, 842.2015Sci..349.842S34097231:CAS:528:DC%2BC2MXhtlKltrfL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26229112][DOI: https://dx.doi.org/10.1126/science.aaa7432]

3. Schachenmayer, J; Pollet, L; Troyer, M; Daley, AJ. Thermalization of strongly interacting bosons after spontaneous emissions in optical lattices. EPJ Quantum Technol.; 2015; 2, [DOI: https://dx.doi.org/10.1140/epjqt15] 1.

4. Kaufman, AM et al. Quantum thermalization through entanglement in an isolated many-body system. Science; 2016; 353, 794.2016Sci..353.794K1:CAS:528:DC%2BC28XhtlCju73M [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27540168][DOI: https://dx.doi.org/10.1126/science.aaf6725]

5. Eigen, C et al. Universal prethermal dynamics of bose gases quenched to unitarity. Nature; 2018; 563, 221.2018Natur.563.221E1:CAS:528:DC%2BC1cXitFejsbfM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30405228][DOI: https://dx.doi.org/10.1038/s41586-018-0674-1]

6. Ueda, M. Quantum equilibration, thermalization and prethermalization in ultracold atoms. Nat. Rev. Phys.; 2020; 2, 669.1:CAS:528:DC%2BB3cXitVSqu7bK [DOI: https://dx.doi.org/10.1038/s42254-020-0237-x]

7. Le, Y; Zhang, Y; Gopalakrishnan, S; Rigol, M; Weiss, DS. Observation of hydrodynamization and local prethermalization in 1d bose gases. Nature; 2023; 618, 494.2023Natur.618.494L1:CAS:528:DC%2BB3sXhtVejsbzE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37198493][DOI: https://dx.doi.org/10.1038/s41586-023-05979-9]

8. Zhou, Z-Y et al. Thermalization dynamics of a gauge theory on a quantum simulator. Science; 2022; 377, 311.2022Sci..377.311Z44742591:CAS:528:DC%2BB38XitVWrtLfI [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35857589][DOI: https://dx.doi.org/10.1126/science.abl6277]

9. Smith, J et al. Many-body localization in a quantum simulator with programmable random disorder. Nat. Phys.; 2016; 12, 907.1:CAS:528:DC%2BC28Xpslymu7o%3D [DOI: https://dx.doi.org/10.1038/nphys3783]

10. Neyenhuis, B. et al. Observation of prethermalization in long-range interacting spin chains. Sci. Adv.3, https://doi.org/10.1126/sciadv.1700672 (2017).

11. Morong, W et al. Observation of stark many-body localization without disorder. Nature; 2021; 599, 393.2021Natur.599.393M1:CAS:528:DC%2BB3MXisFWksL7I [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34789908][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9747247][DOI: https://dx.doi.org/10.1038/s41586-021-03988-0]

12. Kyprianidis, A et al. Observation of a prethermal discrete time crystal. Science; 2021; 372, 1192.2021Sci..372.1192K43775601:CAS:528:DC%2BB3MXhtlamu77P [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34112691][DOI: https://dx.doi.org/10.1126/science.abg8102]

13. Nandkishore, R; Huse, DA. Many-body localization and thermalization in quantum statistical mechanics. Annu. Rev. Condens. Matter Phys.; 2015; 6, 15.2015ARCMP..6..15N [DOI: https://dx.doi.org/10.1146/annurev-conmatphys-031214-014726]

14. Borgonovi, F; Izrailev, FM; Santos, LF; Zelevinsky, VG. Quantum chaos and thermalization in isolated systems of interacting particles. Phys. Rep.; 2016; 626, 1.2016PhR..626..1B34789341:CAS:528:DC%2BC28XjslSjsL0%3D [DOI: https://dx.doi.org/10.1016/j.physrep.2016.02.005]

15. Micha, R; Tkachev, II. Turbulent thermalization. Phys. Rev. D; 2004; 70, 043538.2004PhRvD.70d3538M [DOI: https://dx.doi.org/10.1103/PhysRevD.70.043538]

16. Baier, R; Mueller, AH; Schiff, D; Son, DT. "bottom-up” thermalization in heavy ion collisions. Phys. Lett. B; 2001; 502, 51.2001PhLB.502..51B1:CAS:528:DC%2BD3MXhvFChtLs%3D [DOI: https://dx.doi.org/10.1016/S0370-2693(01)00191-5]

17. Berges, J; Heller, MP; Mazeliauskas, A; Venugopalan, R. QCD thermalization: Ab initio approaches and interdisciplinary connections. Rev. Mod. Phys.; 2021; 93, 035003.2021RvMP..93c5003B43371721:CAS:528:DC%2BB3MXit1yltbrK [DOI: https://dx.doi.org/10.1103/RevModPhys.93.035003]

18. Deutsch, JM. Quantum statistical mechanics in a closed system. Phys. Rev. A; 1991; 43, 2046.1991PhRvA.43.2046D1:CAS:528:DyaK3MXhsFSnsrs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9905246][DOI: https://dx.doi.org/10.1103/PhysRevA.43.2046]

19. Srednicki, M. Chaos and quantum thermalization. Phys. Rev. E; 1994; 50, 888.1994PhRvE.50.888S1:STN:280:DC%2BC2sfisFGgsA%3D%3D [DOI: https://dx.doi.org/10.1103/PhysRevE.50.888]

20. Rigol, M; Dunjko, V; Olshanii, M. Thermalization and its mechanism for generic isolated quantum systems. Nature; 2008; 452, 854.2008Natur.452.854R1:CAS:528:DC%2BD1cXkvVSis7g%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18421349][DOI: https://dx.doi.org/10.1038/nature06838]

21. D’Alessio, L; Kafri, Y; Polkovnikov, A; Rigol, M. From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics. Adv. Phys.; 2016; 65, 239.2016AdPhy.65.239D [DOI: https://dx.doi.org/10.1080/00018732.2016.1198134]

22. Grumbling, E. & Horowitz, M. (Editors) National Academies of Sciences, Engineering, and Medicine. Quantum Computing: Progress and Prospects. https://doi.org/10.17226/25196 (The National Academies Press, Washington, DC, 2019).

23. Arute, F et al. Quantum supremacy using a programmable superconducting processor. Nature; 2019; 574, 505.2019Natur.574.505A1:CAS:528:DC%2BC1MXitVagsb3I [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31645734][DOI: https://dx.doi.org/10.1038/s41586-019-1666-5]

24. Deutsch, IH. Harnessing the power of the second quantum revolution. PRX Quantum; 2020; 1, 020101. [DOI: https://dx.doi.org/10.1103/PRXQuantum.1.020101]

25. Bluvstein, D et al. Controlling quantum many-body dynamics in driven rydberg atom arrays. Science; 2021; 371, 1355.2021Sci..371.1355B42692971:CAS:528:DC%2BB3MXnsVyqtrs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33632894][DOI: https://dx.doi.org/10.1126/science.abg2530]

26. Evered, S. J. et al. High-fidelity parallel entangling gates on a neutral-atom quantum computer. Nature622, 268–272 (2023).

27. Cirac, JI; Zoller, P. Goals and opportunities in quantum simulation. Nat. Phys.; 2012; 8, 264.1:CAS:528:DC%2BC38XkvVGqs74%3D [DOI: https://dx.doi.org/10.1038/nphys2275]

28. Daley, AJ et al. Practical quantum advantage in quantum simulation. Nature; 2022; 607, 667.2022Natur.607.667D1:CAS:528:DC%2BB38XhvFKks7zE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35896643][DOI: https://dx.doi.org/10.1038/s41586-022-04940-6]

29. Bauer, CW et al. Quantum simulation for high-energy physics. PRX Quantum; 2023; 4, 027001.2023PRXQ..4b7001B [DOI: https://dx.doi.org/10.1103/PRXQuantum.4.027001]

30. Haake, F. Quantum Signatures of Chaos (Springer, 1991).

31. Goldstein, S; Lebowitz, JL; Mastrodonato, C; Tumulka, R; Zanghi, N. Normal typicality and von neumann’s quantum ergodic theorem. Proc. R. Soc. A: Math. Phys. Eng. Sci.; 2010; 466, 3203.2010RSPSA.466.3203G2726973 [DOI: https://dx.doi.org/10.1098/rspa.2009.0635]

32. Berry, MV; Tabor, M. Level clustering in the regular spectrum. Proc. R. Soc. Lond. A Math. Phys. Sci.; 1977; 356, 375.1977RSPSA.356.375B [DOI: https://dx.doi.org/10.1098/rspa.1977.0140]

33. Bohigas, O; Giannoni, M-J; Schmit, C. Characterization of chaotic quantum spectra and universality of level fluctuation laws. Phys. Rev. Lett.; 1984; 52, 1.1984PhRvL.52..1B730191 [DOI: https://dx.doi.org/10.1103/PhysRevLett.52.1]

34. Joshi, LK et al. Probing many-body quantum chaos with quantum simulators. Phys. Rev. X; 2022; 12, 011018.1:CAS:528:DC%2BB38XpsFKltb8%3D

35. Pandey, M; Claeys, PW; Campbell, DK; Polkovnikov, A; Sels, D. Adiabatic eigenstate deformations as a sensitive probe for quantum chaos. Phys. Rev. X; 2020; 10, 041017.1:CAS:528:DC%2BB3MXjtFSksLo%3D

36. Li, H; Haldane, FDM. Entanglement spectrum as a generalization of entanglement entropy: Identification of topological order in non-abelian fractional quantum hall effect states. Phys. Rev. Lett.; 2008; 101, 010504.2008PhRvL.101a0504L [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18764098][DOI: https://dx.doi.org/10.1103/PhysRevLett.101.010504]

37. Dalmonte, M; Eisler, V; Falconi, M; Vermersch, B. Entanglement hamiltonians: from field theory to lattice models and experiments. Ann. Phys.; 2022; 534, 2200064. [DOI: https://dx.doi.org/10.1002/andp.202200064]

38. Oganesyan, V; Huse, DA. Localization of interacting fermions at high temperature. Phys. Rev. B; 2007; 75, 155111.2007PhRvB.75o5111O [DOI: https://dx.doi.org/10.1103/PhysRevB.75.155111]

39. Rakovszky, T; Gopalakrishnan, S; Parameswaran, S; Pollmann, F. Signatures of information scrambling in the dynamics of the entanglement spectrum. Phys. Rev. B; 2019; 100, 125115.2019PhRvB.100l5115R1:CAS:528:DC%2BC1MXit1Shs7vJ [DOI: https://dx.doi.org/10.1103/PhysRevB.100.125115]

40. Chang, P-Y; Chen, X; Gopalakrishnan, S; Pixley, J. Evolution of entanglement spectra under generic quantum dynamics. Phys. Rev. Lett.; 2019; 123, 190602.2019PhRvL.123s0602C40327141:CAS:528:DC%2BB3cXhvFyktr8%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31765207][DOI: https://dx.doi.org/10.1103/PhysRevLett.123.190602]

41. Mueller, N; Zache, TV; Ott, R. Thermalization of gauge theories from their entanglement spectrum. Phys. Rev. Lett.; 2022; 129, 011601.2022PhRvL.129a1601M44563681:CAS:528:DC%2BB38XhvFKisrvP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35841570][DOI: https://dx.doi.org/10.1103/PhysRevLett.129.011601]

42. Ma, C-T; Wu, C-H. Quantum entanglement and spectral form factor. Int. J. Theor. Phys.; 2022; 61, 272.4522777 [DOI: https://dx.doi.org/10.1007/s10773-022-05251-2]

43. Froland, H., Zache, T. V., Ott, R. & Mueller, N. Entanglement structure of non-Gaussian states and how to measure it. Preprint at https://arxiv.org/abs/2407.12083 (2024).

44. Suzuki, M. General theory of fractal path integrals with applications to many-body theories and statistical physics. J. Math. Phys.; 1991; 32, 400.1991JMP..32.400S1088360 [DOI: https://dx.doi.org/10.1063/1.529425]

45. Childs, AM; Su, Y; Tran, MC; Wiebe, N; Zhu, S. Theory of Trotter error with commutator scaling. Phys. Rev. X; 2021; 11, 011020.1:CAS:528:DC%2BB3MXovVCntLk%3D

46. Lin, L. Lecture notes on quantum algorithms for scientific computation. Preprint at https://doi.org/10.48550/arXiv.2201.08309 (2022).

47. Childs, A. M. Lecture Notes on Quantum Algorithms. (Lecture notes at University of Maryland, 2017).

48. Bañuls, MC et al. Simulating lattice gauge theories within quantum technologies. Eur. Phys. J. D; 2020; 74, 2020EPJD..74.165B [DOI: https://dx.doi.org/10.1140/epjd/e2020-100571-8] 165.

49. Klco, N; Roggero, A; Savage, MJ. Standard model physics and the digital quantum revolution: thoughts about the interface. Rept. Prog. Phys.; 2022; 85, 064301.2022RPPh..85f4301K44255631:CAS:528:DC%2BB2cXosleqsb0%3D [DOI: https://dx.doi.org/10.1088/1361-6633/ac58a4]

50. Bauer, CW; Davoudi, Z; Klco, N; Savage, MJ. Quantum simulation of fundamental particles and forces. Nature Rev. Phys.; 2023; 5, 420.2023NatRP..5.420B [DOI: https://dx.doi.org/10.1038/s42254-023-00599-8]

51. Di Meglio, A. et al. Quantum computing for high-energy physics: State of the art and challenges. summary of the qc4hep working group. Preprint at https://doi.org/10.48550/arXiv.2307.03236 (2023).

52. Aitchison, I. J. & Hey, A. J. Gauge Theories in Particle Physics: A Practical Introduction (Taylor & Francis, 2012).

53. Quigg, C. Gauge Theories of Strong, Weak, and Electromagnetic Interactions (CRC Press, 2021).

54. Fradkin, E., Field Theories of Condensed Matter Physics (Cambridge University Press, 2013).

55. Kleinert, H. Gauge Fields in Condensed Matter (World Scientific, 1989).

56. Wen, X-G. Topological orders in rigid states. Int. J. Mod. Phys. B; 1990; 4, 239.1990IJMPB..4.239W1043302 [DOI: https://dx.doi.org/10.1142/S0217979290000139]

57. Levin, MA; Wen, X-G. String-net condensation: A physical mechanism for topological phases. Phys. Rev. B; 2005; 71, 045110.2005PhRvB.71d5110L [DOI: https://dx.doi.org/10.1103/PhysRevB.71.045110]

58. Chen, Y-A; Kapustin, A; Radičević, D. Exact bosonization in two spatial dimensions and a new class of lattice gauge theories. Ann. Phys.; 2018; 393, 234.2018AnPhy.393.234C38054201:CAS:528:DC%2BC1cXot1ygu7w%3D [DOI: https://dx.doi.org/10.1016/j.aop.2018.03.024]

59. Chen, Y-A. Exact bosonization in arbitrary dimensions. Phys. Rev. Res.; 2020; 2, 033527.1:CAS:528:DC%2BB3cXitV2gsr7E [DOI: https://dx.doi.org/10.1103/PhysRevResearch.2.033527]

60. Chen, Y-A; Xu, Y. Equivalence between fermion-to-qubit mappings in two spatial dimensions. PRX Quantum; 2023; 4, 010326.2023PRXQ..4a0326C [DOI: https://dx.doi.org/10.1103/PRXQuantum.4.010326]

61. Kitaev, AY. Fault-tolerant quantum computation by anyons. Ann. Phys.; 2003; 303, 2.2003AnPhy.303..2K19510391:CAS:528:DC%2BD3sXisl2hsg%3D%3D [DOI: https://dx.doi.org/10.1016/S0003-4916(02)00018-0]

62. Kitaev, A. Anyons in an exactly solved model and beyond. Ann. Phys.; 2006; 321, 2.2006AnPhy.321..2K22006911:CAS:528:DC%2BD28XjvVeisA%3D%3D [DOI: https://dx.doi.org/10.1016/j.aop.2005.10.005]

63. Das Sarma, S; Freedman, M; Nayak, C. Topological quantum computation. Phys. Today; 2006; 59, 32. [DOI: https://dx.doi.org/10.1063/1.2337825]

64. Nayak, C; Simon, SH; Stern, A; Freedman, M; Das Sarma, S. Non-abelian anyons and topological quantum computation. Rev. Mod. Phys.; 2008; 80, 1083.2008RvMP..80.1083N24437221:CAS:528:DC%2BD1cXhtFCku73K [DOI: https://dx.doi.org/10.1103/RevModPhys.80.1083]

65. Lahtinen, V; Pachos, JK. A short introduction to topological quantum computation. SciPost Phys.; 2017; 3, 021.2017ScPP..3..21L [DOI: https://dx.doi.org/10.21468/SciPostPhys.3.3.021]

66. Maunz, P. L. W. High Optical Access Trap 2.0. https://doi.org/10.2172/1237003 (2016).

67. Monroe, C et al. Programmable quantum simulations of spin systems with trapped ions. Rev. Mod. Phys.; 2021; 93, 025001.2021RvMP..93b5001M42909451:CAS:528:DC%2BB3MXhs1Witr7L [DOI: https://dx.doi.org/10.1103/RevModPhys.93.025001]

68. Pogorelov, I et al. Compact ion-trap quantum computing demonstrator. PRX Quantum; 2021; 2, 020343.2021PRXQ..2b0343P [DOI: https://dx.doi.org/10.1103/PRXQuantum.2.020343]

69. Yao, R et al. Experimental realization of a multiqubit quantum memory in a 218-ion chain. Phys. Rev. A; 2022; 106, 062617.2022PhRvA.106f2617Y1:CAS:528:DC%2BB3sXhtlGjtrw%3D [DOI: https://dx.doi.org/10.1103/PhysRevA.106.062617]

70. Bohnet, JG et al. Quantum spin dynamics and entanglement generation with hundreds of trapped ions. Science; 2016; 352, 1297.2016Sci..352.1297B35607621:CAS:528:DC%2BC28Xpt1Gntbk%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27284189][DOI: https://dx.doi.org/10.1126/science.aad9958]

71. Pichler, H; Zhu, G; Seif, A; Zoller, P; Hafezi, M. Measurement protocol for the entanglement spectrum of cold atoms. Phys. Rev. X; 2016; 6, 041033.

72. Dalmonte, M; Vermersch, B; Zoller, P. Quantum simulation and spectroscopy of entanglement hamiltonians. Nat. Phys.; 2018; 14, 827.1:CAS:528:DC%2BC1cXhtFGksbbL [DOI: https://dx.doi.org/10.1038/s41567-018-0151-7]

73. Kokail, C; van Bijnen, R; Elben, A; Vermersch, B; Zoller, P. Entanglement hamiltonian tomography in quantum simulation. Nat. Phys.; 2021; 17, 936.1:CAS:528:DC%2BB3MXhsVSmsb3I [DOI: https://dx.doi.org/10.1038/s41567-021-01260-w]

74. Kokail, C et al. Quantum variational learning of the entanglement hamiltonian. Phys. Rev. Lett.; 2021; 127, 170501.2021PhRvL.127q0501K1:CAS:528:DC%2BB3MXis1SmsbbP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34739272][DOI: https://dx.doi.org/10.1103/PhysRevLett.127.170501]

75. Zache, TV; Kokail, C; Sundar, B; Zoller, P. Entanglement spectroscopy and probing the Li-Haldane conjecture in topological quantum matter. Quantum; 2022; 6, 702. [DOI: https://dx.doi.org/10.22331/q-2022-04-27-702]

76. Mueller, N et al. Quantum computation of dynamical quantum phase transitions and entanglement tomography in a lattice gauge theory. PRX Quantum; 2023; 4, 030323.2023PRXQ..4c0323M [DOI: https://dx.doi.org/10.1103/PRXQuantum.4.030323]

77. Paini, M. & Kalev, A. An approximate description of quantum states. Preprint at https://doi.org/10.48550/arXiv.1910.10543 (2019).

78. Huang, H-Y; Kueng, R; Preskill, J. Predicting many properties of a quantum system from very few measurements. Nat. Phys.; 2020; 16, 1050.1:CAS:528:DC%2BB3cXht1WrtLjK [DOI: https://dx.doi.org/10.1038/s41567-020-0932-7]

79. Huang, H-Y; Kueng, R; Preskill, J. Efficient estimation of pauli observables by derandomization. Phys. Rev. Lett.; 2021; 127, 030503.2021PhRvL.127c0503H42920441:CAS:528:DC%2BB3MXhslWgu7jF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34328776][DOI: https://dx.doi.org/10.1103/PhysRevLett.127.030503]

80. Hu, H.-Y., Choi, S. & You, Y.-Z. Classical shadow tomography with locally scrambled quantum dynamics. Phys. Rev. Res.5, 023027 (2023).

81. Kunjummen, J; Tran, MC; Carney, D; Taylor, JM. Shadow process tomography of quantum channels. Phys. Rev. A; 2023; 107, 042403.2023PhRvA.107d2403K45830301:CAS:528:DC%2BB3sXhtVeku7rP [DOI: https://dx.doi.org/10.1103/PhysRevA.107.042403]

82. Levy, R; Luo, D; Clark, BK. Classical shadows for quantum process tomography on near-term quantum computers. Phys. Rev. Res.; 2024; 6, 013029.1:CAS:528:DC%2BB2cXotVChu7s%3D [DOI: https://dx.doi.org/10.1103/PhysRevResearch.6.013029]

83. Huang, H-YR et al. Demonstrating quantum advantage in learning from experiments. Science; 2021; 376, 1182.2022Sci..376.1182H [DOI: https://dx.doi.org/10.1126/science.abn7293]

84. Huang, H.-Y. Learning quantum states from their classical shadows. Nat. Rev. Phys.4, 81 (2022).

85. Bringewatt, J; Kunjummen, J; Mueller, N. Randomized measurement protocols for lattice gauge theories. Quantum; 2024; 8, 1300. [DOI: https://dx.doi.org/10.22331/q-2024-03-27-1300]

86. Elben, A; Vermersch, B; Roos, CF; Zoller, P. Statistical correlations between locally randomized measurements: A toolbox for probing entanglement in many-body quantum states. Phys. Rev. A; 2019; 99, 052323.2019PhRvA.99e2323E39578711:CAS:528:DC%2BC1MXhtFKlsLjF [DOI: https://dx.doi.org/10.1103/PhysRevA.99.052323]

87. Elben, A et al. Mixed-state entanglement from local randomized measurements. Phys. Rev. Lett.; 2020; 125, 200501.2020PhRvL.125t0501E1:CAS:528:DC%2BB3cXisVWrs7bP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33258654][DOI: https://dx.doi.org/10.1103/PhysRevLett.125.200501]

88. Elben, A et al. The randomized measurement toolbox. Nat. Rev. Phys.; 2023; 5, 9. [DOI: https://dx.doi.org/10.1038/s42254-022-00535-2]

89. Bisognano, JJ; Wichmann, EH. On the duality condition for a hermitian scalar field. J. Math. Phys.; 1975; 16, 985.1975JMP..16.985B438943 [DOI: https://dx.doi.org/10.1063/1.522605]

90. Bisognano, JJ; Wichmann, EH. On the duality condition for quantum fields. J. Math. Phys; 1976; 17, 303.1976JMP..17.303B438944 [DOI: https://dx.doi.org/10.1063/1.522898]

91. Guhr, T; Müller-Groeling, A; Weidenmüller, HA. Random-matrix theories in quantum physics: common concepts. Phys. Rep.; 1998; 299, 189.1998PhR..299.189G16284671:CAS:528:DyaK1cXivVKntb8%3D [DOI: https://dx.doi.org/10.1016/S0370-1573(97)00088-4]

92. Mehta, M. L. Random Matrices (Elsevier, 2004).

93. Rosenzweig, N; Porter, CE. repulsion of energy levels in complex atomic spectra. Phys. Rev.; 1960; 120, 1698.1960PhRv.120.1698R1:CAS:528:DyaF3MXjtFeksg%3D%3D [DOI: https://dx.doi.org/10.1103/PhysRev.120.1698]

94. Giraud, O; Macé, N; Vernier, É; Alet, F. Probing symmetries of quantum many-body systems through gap ratio statistics. Phys. Rev. X; 2022; 12, 011006.1:CAS:528:DC%2BB38XpsFKlsLg%3D

95. Dağ, CB; Mistakidis, SI; Chan, A; Sadeghpour, HR. Many-body quantum chaos in stroboscopically-driven cold atoms. Commun. Phys.; 2023; 6, 136. [DOI: https://dx.doi.org/10.1038/s42005-023-01258-1]

96. Ghosh, S; Soni, RM; Trivedi, SP. On the entanglement entropy for gauge theories. J. High Energy Phys.; 2015; 2015, 3351247 [DOI: https://dx.doi.org/10.1007/JHEP09(2015)069] 1.

97. Van Acoleyen, K et al. Entanglement of distillation for lattice gauge theories. Phys. Rev. Lett.; 2016; 117, 131602.2016PhRvL.117m1602V3636531 [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27715127][DOI: https://dx.doi.org/10.1103/PhysRevLett.117.131602]

98. Brydges, T et al. Probing rényi entanglement entropy via randomized measurements. Science; 2019; 364, 260.2019Sci..364.260B1:CAS:528:DC%2BC1MXns1CgtLs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31000658][DOI: https://dx.doi.org/10.1126/science.aau4963]

99. Dankert, C; Cleve, R; Emerson, J; Livine, E. Exact and approximate unitary 2-designs and their application to fidelity estimation. Phys. Rev. A; 2009; 80, 012304.2009PhRvA.80a2304D [DOI: https://dx.doi.org/10.1103/PhysRevA.80.012304]

100. Mezei, M; Stanford, D. On entanglement spreading in chaotic systems. J. High Energy Phys.; 2017; 2017, 3662866 [DOI: https://dx.doi.org/10.1007/JHEP05(2017)065] 1.

101. Jonay, C., Huse, D. A. & Nahum, A. Coarse-grained dynamics of operator and state entanglement. Preprint at https://doi.org/10.48550/arXiv.1803.00089 (2018).

102. Rigol, M; Dunjko, V; Yurovsky, V; Olshanii, M. Relaxation in a completely integrable many-body quantum system: an ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons. Phys. Rev. Lett.; 2007; 98, 050405.2007PhRvL.98e0405R [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17358832][DOI: https://dx.doi.org/10.1103/PhysRevLett.98.050405]

103. Hearth, S. N., Flynn, M. O., Chandran, A. & Laumann, C. R. Unitary k-designs from random number-conserving quantum circuits. Phys. Rev. X15, 021022 (2025).

104. Hearth, S. N., Flynn, M. O., Chandran, A. & Laumann, C. R. Efficient local classical shadow tomography with number conservation. Phys. Rev. Lett.133, 060802 (2024).

105. Joshi, M. K. et al. Exploring large-scale entanglement in quantum simulation. Nature624, 539–544 (2023).

106. Debroy, DM; Li, M; Huang, S; Brown, KR. Logical performance of 9 qubit compass codes in ion traps with crosstalk errors. Quantum Sci. Technol.; 2020; 5, 034002.2020QS&T..5c4002D [DOI: https://dx.doi.org/10.1088/2058-9565/ab7e80]

107. Cetina, M. et al. Control of transverse motion for quantum gates on individually addressed atomic qubits. PRX Quantum3, https://doi.org/10.1103/PRXQuantum.3.010334 (2022).

108. Carrasco, J; Elben, A; Kokail, C; Kraus, B; Zoller, P. Theoretical and experimental perspectives of quantum verification. PRX Quantum; 2021; 2, 010102. [DOI: https://dx.doi.org/10.1103/PRXQuantum.2.010102]

109. Berges, J; Borsányi, S; Wetterich, C. Prethermalization. Phys. Rev. Lett.; 2004; 93, 142002.2004PhRvL.93n2002B1:STN:280:DC%2BD2crlsVSjtg%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15524783][DOI: https://dx.doi.org/10.1103/PhysRevLett.93.142002]

110. Gring, M et al. Relaxation and prethermalization in an isolated quantum system. Science; 2012; 337, 1318.2012Sci..337.1318G1:CAS:528:DC%2BC38XhtlaitL7F [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22956685][DOI: https://dx.doi.org/10.1126/science.1224953]

111. Mori, T; Ikeda, TN; Kaminishi, E; Ueda, M. Thermalization and prethermalization in isolated quantum systems: a theoretical overview. J. Phys. B At. Mol. Opt. Phys.; 2018; 51, 112001.2018JPhB..51k2001M [DOI: https://dx.doi.org/10.1088/1361-6455/aabcdf]

112. Callen, HB; Welton, TA. Irreversibility and generalized noise. Phys. Rev.; 1951; 83, 34.1951PhRv..83..34C44778 [DOI: https://dx.doi.org/10.1103/PhysRev.83.34]

113. Kubo, R. The fluctuation-dissipation theorem. Rep. Prog. Phys.; 1966; 29, 255.1966RPPh..29.255K1:CAS:528:DyaF2sXlvVyqtA%3D%3D [DOI: https://dx.doi.org/10.1088/0034-4885/29/1/306]

114. Orioli, AP; Berges, J. Breaking the fluctuation-dissipation relation by universal transport processes. Phys. Rev. Lett.; 2019; 122, 150401. [DOI: https://dx.doi.org/10.1103/PhysRevLett.122.150401]

115. Schuckert, A; Knap, M. Probing eigenstate thermalization in quantum simulators via fluctuation-dissipation relations. Phys. Rev. Res.; 2020; 2, 043315.1:CAS:528:DC%2BB3MXhtVegu7%2FL [DOI: https://dx.doi.org/10.1103/PhysRevResearch.2.043315]

116. Murthy, C; Babakhani, A; Iniguez, F; Srednicki, M; Halpern, NY. Non-abelian eigenstate thermalization hypothesis. Phys. Rev. Lett.; 2023; 130, 140402.2023PhRvL.130n0402M45826481:CAS:528:DC%2BB3sXoslehtLo%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37084457][DOI: https://dx.doi.org/10.1103/PhysRevLett.130.140402]

117. Yao, X. SU(2) gauge theory in 2+1 dimensions on a plaquette chain obeys the eigenstate thermalization hypothesis. Phys. Rev. D; 2023; 108, L031504.2023PhRvD.108c1504Y46462571:CAS:528:DC%2BB3sXhvFyjurbL [DOI: https://dx.doi.org/10.1103/PhysRevD.108.L031504]

118. Ebner, L; Schäfer, A; Seidl, C; Müller, B; Yao, X. Eigenstate thermalization in (2 + 1)-dimensional su (2) lattice gauge theory. Phys. Rev. D; 2024; 109, 014504.2024PhRvD.109a4504E47059161:CAS:528:DC%2BB2cXjslajsLw%3D [DOI: https://dx.doi.org/10.1103/PhysRevD.109.014504]

119. Hartse, J., Fidkowski, L. & Mueller, N. Stabilizer scars. Preprint at https://doi.org/10.48550/arXiv.2411.12797 (2024).

120. Davoudi, Z. et al. Quantum thermodynamics of nonequilibrium processes in lattice gauge theories. Phys. Rev. Lett.133, 250402 (2024).

121. Wegner, F. J. Duality in Generalized Ising Models. (Oxford University Press, 2017).

122. Brown, KR; Harrow, AW; Chuang, IL. Arbitrarily accurate composite pulse sequences. Phys. Rev. A; 2004; 70, 052318.2004PhRvA.70e2318B [DOI: https://dx.doi.org/10.1103/PhysRevA.70.052318]

123. Sørensen, A; Mølmer, K. Quantum computation with ions in thermal motion. Phys. Rev. Lett.; 1999; 82, 1971.1999PhRvL.82.1971S [DOI: https://dx.doi.org/10.1103/PhysRevLett.82.1971]

124. Leung, PH et al. Robust 2-qubit gates in a linear ion crystal using a frequency-modulated driving force. Phys. Rev. Lett.; 2018; 120, 020501.2018PhRvL.120b0501L1:CAS:528:DC%2BC1MXltVyjs7k%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29376710][DOI: https://dx.doi.org/10.1103/PhysRevLett.120.020501]

125. Huang, S., Brown, K. R. & Cetina, M. Comparing Shor and Steane error correction using the Bacon-Shor code. Sci. Adv. 10, eadp2008 (2024).

126. MATLAB, version 23.2.0.2485118 (R2023a) (The MathWorks Inc., Natick, Massachusetts, 2023).

127. MathWorks, Fmincon documentation. https://www.mathworks.com/help/optim/ug/fmincon.html (2024).

128. MathWorks Global Search documentation. https://www.mathworks.com/help/gads/globalsearch.html (2024).