Visualizing dynamics of charges and strings in (2 + 1)D lattice gauge theories

Abstract

Lattice gauge theories (LGTs)^1,2,3,4 can be used to understand a wide range of phenomena, from elementary particle scattering in high-energy physics to effective descriptions of many-body interactions in materials^5,6,7. Studying dynamical properties of emergent phases can be challenging, as it requires solving many-body problems that are generally beyond perturbative limits^8,9,10. Here we investigate the dynamics of local excitations in a ${{\mathbb{Z}}}_{2}$ LGT using a two-dimensional lattice of superconducting qubits. We first construct a simple variational circuit that prepares low-energy states that have a large overlap with the ground state; then we create charge excitations with local gates and simulate their quantum dynamics by means of a discretized time evolution. As the electric field coupling constant is increased, our measurements show signatures of transitioning from deconfined to confined dynamics. For confined excitations, the electric field induces a tension in the string connecting them. Our method allows us to experimentally image string dynamics in a (2+1)D LGT, from which we uncover two distinct regimes inside the confining phase: for weak confinement, the string fluctuates strongly in the transverse direction, whereas for strong confinement, transverse fluctuations are effectively frozen^11,12. We also demonstrate a resonance condition at which dynamical string breaking is facilitated. Our LGT implementation on a quantum processor presents a new set of techniques for investigating emergent excitations and string dynamics.

Main

Present models for fundamental forces are formulated as gauge theories. The common element of these theories is a local symmetry action and its corresponding gauge field that mediates interaction between matter particles⁵. Gauge theories are not limited to high-energy physics but can also capture emergent phenomena in condensed matter physics^6,7 and have seen applications in quantum information¹³. One of the earliest examples of the interplay between these research fields was the development of LGT, in which space is discretized to a lattice^1,2,3,4. In particular, a motivation for introducing quantum LGTs was to describe a mechanism for confinement of quarks in quantum chromodynamics². Within this framework, confined matter particles are the open ends of a string with finite tension. The discrete nature of LGTs has also been important in forming a framework for numerical calculations of equilibrium properties, for instance, using Monte Carlo or tensor-network-based methods⁹.

Understanding the non-equilibrium dynamics of string excitations in LGTs is of fundamental importance in various disciplines, ranging from transport properties of the quark–gluon plasma to spectral properties in correlated quantum materials. However, theoretical approaches to this problem face substantial obstacles: non-equilibrium dynamics is beyond perturbative treatments, numerical methods based on Monte Carlo run into sign problems and tensor-network approaches work only as long as entanglement remains sufficiently low^8,9,10. Quantum devices have been proposed as a viable alternative for the study of LGTs (refs. ^{14,15,16,17,18,19,20} for early works and reviews); their experimental implementations, on the other hand, have been limited to one spatial dimension or small scales, which limits the ability to examine string fluctuations^{21,22,23,24,25,26,27,28,29,30,31,32,33,34,35}. Because conventional LGT Hamiltonians have a constrained structure dictated by the local symmetry action, directly simulating their dynamics on quantum processors requires the ability to perform evolution generated by specific multibody local terms.

Here we realize a two-dimensional LGT on a superconducting quantum processor and use this platform to investigate and visualize the string dynamics. We consider an LGT in which the interaction between matter fields (filled circles in Fig. 1a), placed on the vertices of a square lattice, is mediated by ${{\mathbb{Z}}}_{2}$ gauge fields, located on the links that connect them³ (green diamonds in Fig. 1a). This structure is a simplification of quantum electrodynamics, in which both space and the gauge group are discretized: space becomes a lattice and the U(1) gauge group is discretized to ${{\mathbb{Z}}}_{2}$. We make use of the gauge redundancy to eliminate matter fields^36,37. In the resulting ‘matter-free’ LGT, the motion and interaction of matter fields are captured by the ${{\mathbb{Z}}}_{2}$ gauge fields with the Hamiltonian:

$${\mathcal{H}}=-{J}_{{\rm{E}}}\sum _{v}{A}_{v}-{J}_{{\rm{M}}}\sum _{p}{B}_{p}-{h}_{{\rm{E}}}\sum _{{\rm{links}}}{Z}_{l}-\lambda \sum _{{\rm{links}}}{X}_{l}.$$

(1)

Fig. 1: An LGT and its phase diagram.

a, A full two-dimensional LGT (top left) can be realized by placing charged matter (grey circles) on vertices of a square lattice and gauge fields on the links between them (green diamonds). The local gauge structure can be used to eliminate the matter field and arrive at an effective theory involving gauge fields only (right). The presence/absence of charge excitations (red/blue) or magnetic fluxes (yellow/purple) is then sensed through the links. b, Zero-temperature phase diagram of the LGT in equation (1). c, In the deconfined phase, charges move freely. In the confined phase, charges oscillate around an equilibrium configuration. We can picture an elastic string connecting them that fluctuates in both longitudinal and transverse directions, limiting their motion.

The vertex operators A_v = ∏_i∈vZ_i are products of local Z operators on link qubits emanating from a vertex v and represent the electric charge (red or blue tiles in Fig. 1a). The plaquette operators B_p = ∏_i∈pX_i are products of Pauli-X operators on link qubits encircling a plaquette and represent the presence or absence of magnetic flux (yellow or purple tiles). We consider vertex and plaquette operators of equal strength, which sets the unit of energy, J_E = J_M = 1. The h_E terms denote an electric field on each link that creates magnetic flux excitations. The coupling terms λ generate hopping of matter fields located at adjacent vertices, as mediated by a gauge field on their connecting link.

Since the foundational work of Fradkin and Shenker, it has been known that the zero-temperature phase diagram of ${\mathcal{H}}$ has two distinct phases^38,39,40,41 (Fig. 1b). One phase is the deconfined and topologically ordered phase that exists near h_E, λ ≈ 0. The quantum phase transition along the λ = 0 line can be understood by a duality mapping to the transverse-field Ising model¹, in which domain walls of the Ising model correspond to closed strings in ${\mathcal{H}}$. For small but non-zero λ, the duality breaks down: we must include contributions to dynamics from open strings, which cannot correspond directly to domain walls. Crossing this transition into the confining phase leads to a condensation of magnetic excitations and confinement of electric excitations. The deconfinement to confinement transition can be seen in the non-equilibrium dynamics of a pair of charge excitations. In the deconfined case, the excitations move freely, whereas in the confined case, the string between them acquires a tension and restricts their motion (Fig. 1c).

For h_E = λ = 0, ${\mathcal{H}}$ reduces to the celebrated toric code Hamiltonian¹³, which underlies several quantum computing error-correction codes. In that limit, all terms in ${\mathcal{H}}$ commute with each other, [A_v, B_p] = 0 ∀ v, p; hence ${\mathcal{H}}$ is exactly solvable. The efficient preparation of the toric code ground state is well studied and can be achieved with circuits that scale linearly with the shorter dimension of the lattice^42,43,44. In the limit h_E, λ ≫ 1, the ground state is a product state of the qubits with all qubits pointing in the same direction, which can be prepared with single-qubit operations. A key aspect of ${\mathcal{H}}$ at the system sizes we study is the existence of an efficient algorithm to prepare states at energy densities low enough to resolve characteristic dynamics throughout the phase diagram. We make use of a variational ansatz based on a parameterization of the gate sequence used to generate the toric code wavefunction, which we call the weight-adjustable loop ansatz (WALA)^45,46 (Fig. 2a,b and Methods). To implement this ansatz, we use a grid of qubits with fourfold connectivity (diamonds) and ancilla qubits (circles; Fig. 2a) at the centre of each plaquette of the link qubits. All qubits begin in the |0⟩ state. The state-preparation sequence starts with a single-qubit rotation R_y(θ) on each of the ancilla qubits at the centre of the plaquettes. The rest of the gate sequence does not have any adjustable parameters and is composed of controlled NOT (C-NOT) gates that generate entanglement between the qubits, starting with the centre columns of plaquettes and spreading to the edges of the lattice. The final C-NOT gate disentangles the ancilla qubits, returning them to the |0⟩ state. We use a classical computer to find the optimal angle θ that minimizes the initial state energy, as a function of h_E (Fig. 2b and Supplementary Information Section II). The WALA circuit is equivalent to a mean-field ansatz for the dual Ising model and has the advantage of being very easy to optimize classically even for larger system sizes, in contrast to the actual ground state⁴⁵ (Methods). The resulting quantum state, |ψ₀⟩, is then used as a low-energy-density initial state for simulating the dynamics. In classical simulations, the dynamics are hard to capture owing to the fast entanglement growth (Supplementary Information Section V).

Fig. 2: WALA.

a, WALA gate sequence used for a two-dimensional grid of 35 qubits, consisting of 17 link qubits (diamonds) and 18 ancilla qubits (circles). The sequence begins with applying R_y(θ) to ancilla qubits of each plaquette, followed by applying C-NOT gates to qubit pairs, starting at the centre columns and moving outwards. b, Optimized θ angle used in WALA. The green curve is based on numerical calculations for a 35-qubit grid and the grey curve shows the thermodynamic limit. c, Energy error, compared with the exact ground state, of three ansatzes: (1) WALA (green); (2) toric code, θ = π/2 (blue); and (3) product state, |0⟩^⊗N (yellow), for λ = 0.25. Solid lines correspond to circuit simulations and filled circles are extracted from our experiment after readout error mitigation (Supplementary Information Section III.A). d, Experimentally measured expectation values of plaquette, vertex and Pauli-Z operators, for λ = 0.25 and h_E ∈ {0, 0.3, 0.6, 1.0}, from WALA. We post-select the measured data on the ancilla being in the expected |0⟩ state to mitigate decoherence of the device for this and all other figures of the main text (Methods).

In Fig. 2c, we show the energy error for λ = 0.25 as a function of h_E using the optimized angles in the WALA (green markers/line). The energy error is small for all values of h_E for the WALA initial state (Supplementary Information Section III.A). Preparing instead the toric code ground state (θ = π/2) gives good overlap with the true ground state only for h_E ≪ 1. Away from this limit, the energy error grows rapidly (blue markers/line). In the opposite limit, h_E ≫ 1, the polarized state is the ground state of ${\mathcal{H}}$. Considering this ansatz yields acceptable energy errors for large values of h_E (orange markers/line), but when reducing h_E, the energy error becomes large as well. The good performance of the WALA relies on the finite size of the system. In the thermodynamic limit, the WALA reduces to the toric code ground state for h_E ∈ [0, h_mf], in which h_mf = 0.25 is the mean-field transition point (grey line in Fig. 2b). To characterize the WALA pulse sequence, we measure the expectation values of A_v and B_p and also Pauli-Z on individual qubits (diamonds in Fig. 2d). These local observables show changes of the various terms of ${\mathcal{H}}$, that is, the local energy density. Owing to the form of our ansatz, all ⟨A_v⟩ (=1) and ⟨X_l⟩ (=0, not shown) terms remain constant as θ changes with increasing h_E. Variations in the energy density arise from the decrease of the magnetic parity values (B_p terms) and the emergence of the Pauli-Z polarization values (Z_l terms). The non-uniform Z_l expectation values result from the distinct connectivity of the boundary qubits.

Having designed a circuit that approximates the ground state, we next study charge confinement by measuring the dynamics of a pair of electric excitations (Fig. 3). By using ancilla qubits at each vertex and plaquette centre, we are able to implement an efficient Suzuki–Trotter expansion of the time evolution operator generated by equation (1). Each time step has eight distinct layers, consisting of single-qubit rotations and controlled Z (CZ) gates, totalling 116 CZ gates per time step for the grid of 35 qubits (Methods). We prepare a pair of electric excitations on neighbouring sites in the centre of the system by applying a single Pauli-X on top of the WALA state (Fig. 3a). This operator creates a non-equilibrium initial state that has a notable overlap with the elementary excitations and allows us to examine the deconfinement–confinement dynamics of the charges. By measuring the average separation and the spatially resolved average position of these excitations as a function of time (Fig. 3a,b), we find that they show qualitatively distinct dynamical signatures as the electric field is tuned. Although for weak electric fields the excitations spread quickly across the whole system, at strong electric fields, the two charges stay together, as indicated by the small average separation; this observation constitutes a dynamical signature of charge confinement. Although the excitation separation increases for all values of h_E, we compare with the case when we evolve under the pure toric code Hamiltonian, for which the separation is exactly stationary in theory (grey region in Fig. 3a). The increase of separation in this case indicates that decoherence of the quantum state is pushing the system towards the maximally mixed state, which has an expected separation of 7/3. By also adding depolarization mitigation (Methods), we obtain quantitative agreement with numerical results and even reveal oscillations about an average separation that is much smaller than the system size when h_E = 2.0, indicative of a confining potential. These dynamical signatures support the onset of a confining potential near the Ising critical point, which—in the thermodynamic limit—is located between h_E ≈ 0.33 and  0.34, depending on λ as obtained from numerical ground-state studies^40,41,47,48.

Fig. 3: Confinement of electric excitations.

a, A Pauli-X gate applied to the WALA initial state creates two electric excitations on adjacent vertices (red tiles). The coupling induces dynamics to the excitations and is set to λ = 0.25 for all data in this figure. After post-selecting bitstrings that correspond to two electric excitations, the separation of the excitations is monitored as a function of time for different electric fields h_E. The grey area is bounded by the separation measured when evolving under the pure toric code Hamiltonian (λ = h_E = 0). The lower panel shows the rescaled data, assuming a global depolarizing noise channel, and compares it with exact circuit simulation. b, Average density of electric excitations as measured by ⟨A_v⟩ for h_E = 0 (left) and h_E = 2.0 (right). c, Two superposition states in which the excitations interfere constructively (|ψ₊⟩) and destructively (|ψ₋⟩) at short distances, respectively. d, The separation of the excitations as a function of time for different electric fields h_E. Lines represent guides for the eye. All error bars in this work show standard deviation of the mean and are smaller than the markers when not visible. e, Spatial maps of the probability that an electric vertex A_i,j is excited, conditioned on the electric vertex being excited at time t = 3.5: P(i, j|). The different columns correspond to different , indicated by the boxed symbols. The unconditioned probability that is excited is written below the symbol. The colour scale represents P(i, j|) for all i and j. The top and bottom rows show results for h_E = 0 and 2.0, respectively. We implement dynamical decoupling and randomized compiling to mitigate control errors, as well as idle dephasing (Methods).

To further accentuate the difference between confined and deconfined dynamics, we consider two other initial configurations, |ψ₊⟩ and |ψ₋⟩ (Fig. 3c). These are positive and negative superpositions of a pair of excitations at a lattice distance of two. The intuition for choosing these initial states comes from approximations to different angular momentum eigenstates. Their dynamics can be understood from quantum interference: at short times, hopping of electric excitations that brings the pair closer together interferes constructively for |ψ₊⟩ and destructively for |ψ₋⟩. In the deconfined phase, this leads to the excitations initially moving closer together for |ψ₋⟩ and further apart for |ψ₋⟩, as observed in Fig. 3d for h_E = 0. By contrast, when h_E = 2.0, the string tension dominates and the excitations remain close to their initial separation for both initial states. The excitation separation initialized in this state has the added feature of being robust against Trotter error, allowing us to increase the Trotter step from dt = 0.3 to dt = 0.5 and reach later times (Supplementary Information Section IV.C), increasing the signal strength. The confinement signatures observed in excitation separations can be further corroborated by analysing the probability of finding an excitation at a given site, conditioned on measuring another excitation somewhere else in the lattice (Fig. 3e). The data shown are for a fixed time t = 3.5 for the |ψ₋⟩ configuration (see Methods for |ψ₊⟩). For h_E = 0, the probabilities are spread across the system, with a higher probability of the charges being found further apart than their initial separation. For h_E = 2.0, there is only a substantial probability observed at separations 1, 2 and 3, indicating that the excitations tend to stay close to their initial position or hop together in a correlated fashion and demonstrating confinement of pairs of electric charges.

In the confined regime, electric charges are located at the ends of an elastic string, which—in our two-dimensional setting—can vibrate transversely, akin to a violin string. We generate the string by applying X gates that traverse the system on top of the WALA circuit from an auxiliary qubit on the left to another one on the right (Fig. 4a). By performing a Trotterized time evolution, which excludes field terms on these extra edge qubits, the so-created electric charges will remain pinned at the edges, whereas the string itself can evolve dynamically. To investigate the vibrational dynamics of the string, we measure a two-time correlator in the Z basis:

$${{\mathcal{S}}}_{ZZ}(t)={\rm{Re}}[\langle Z(t)Z(0)\rangle ]\times \langle Z(0)\rangle $$

(2)

for each qubit. We measure ${{\mathcal{S}}}_{ZZ}(t)$ using a Hadamard test with an auxiliary circuit (Fig. 4b and Methods). This correlation function is a product of two terms. The first term is sensitive to whether the presence of the string has changed compared with its initial value at time zero, that is, it captures the stiffness of the string. The second term measures whether a string has been created on top of the WALA sequence initially, which is only possible in the confined regime. The combined correlation function, ${{\mathcal{S}}}_{ZZ}(t)$, allows us to determine the string dynamics. Note that, although for λ = 0 strings correspond to Ising domain walls, for finite fields λ ≠ 0, there exists no direct mapping to domain walls.

Fig. 4: Dynamics of the string connecting two spatially localized electric charges.

a, Schematic of the initial state preparation. Starting from the WALA initial state as vacuum, we create a pair of separated electric excitations by applying a string of X gates spanning from an extra qubit on the left (leftmost diamond) to one on the right (rightmost diamond). By not applying the local field terms of the time evolution on those two extra qubits, the excitations remain pinned, whereas the string itself can evolve dynamically. b, Circuit for measuring the unequal-time correlation function Re[⟨Z(t)Z(0)⟩]. The P gate is H(S^†)^b, in which b = 0/1 corresponds to measuring the real/imaginary part of ⟨Z(t)Z(0)⟩. The CZ gate entangling the auxiliary qubit (Q_b) to the gauge qubit of interest (Q_l) may need to be mediated by means of swaps through an ancillary qubit (Q_a) and another gauge qubit (Q_g). c, Spatial maps of ${{\mathcal{S}}}_{ZZ}(t)={\rm{Re}}[\langle Z(t)Z(0)\rangle ]\times \langle Z(0)\rangle$ for varying times and confining field h_E, at λ = 0.25 and dt = 0.3 (the same for all data in the figure). The extra qubits on either side used for state preparation are not shown. d, Re[⟨Z(t)Z(0)⟩] and ${{\mathcal{S}}}_{ZZ}(t)$ for qubits Q₁ and Q₂ in the centre-top and centre-bottom respectively, as labelled in panel c. Lines are guides for the eye. The grey regions in these plots correspond to the region limited by decoherence and is bounded by $| \langle Z(t)Z(0)\rangle {| }_{\lambda ={h}_{{\rm{E}}}=0}$.

Our measurements of ${{\mathcal{S}}}_{ZZ}(t)$ reveal three distinct regimes (Fig. 4c,d). (1) In the deconfined phase, h_E = 0.1, applying an X-string on top of the WALA sequence does not create a string excitation. The correlation function Re[⟨Z(t)Z(0)⟩] experiences coherent oscillations from the energy difference between the WALA state and the first excited state, whereas the small expectation value of ⟨Z(0)⟩ suppresses ${{\mathcal{S}}}_{ZZ}(t)$. (2) In the intermediate electric field regime, h_E = 0.6, the dynamics are already confining but the string tension is not too large. Thus, changes of the string length from higher-order dynamical processes are energetically accessible. Our measurements show a clear initial string along the path we prepared. Already after a short temporal evolution, the ${{\mathcal{S}}}_{ZZ}(t)$ correlations of the qubits both at the bottom and at the top of the grid quickly decay to zero and the string is equally likely found on either side of the system. In this regime, the string is floppy and fluctuates strongly, even though charges remain confined^11,12. (3) Deep in the confined regime, for large h_E = 1.4, we see dynamics of the initial bump at the top of the grid but very little probability that the string moves to the bottom on the timescales of the simulation. This can be understood from the large string tension in the deeply confined regime, which suppresses the flopping to the other side within experimentally accessible timescales. However, the string is still moving freely around the top qubits despite the large string tension. These length-scale-preserving moves result from the lattice discretization of our LGT, as the plaquette terms B_p of the Hamiltonian can deform the string without changing its length.

Having visualized the vibrations of the string connecting two electric charges in the confined regime, we now investigate string-breaking and pair-creation dynamics. The coupling λ can dynamically create pairs of electric excitations. When this process occurs in the ground state, the energy is increased as a result of the cost of creating an excitation pair (Fig. 5a). When, by contrast, this process occurs on a string, there is a string-energy gain that competes with the energy cost of the pair creation.

Fig. 5: String breaking.

a, Schematics for pair creation from vacuum fluctuation and string breaking. b, Difference in the charge excitation values in the presence and absence of the string ⟨A_v⟩_string − ⟨A_v⟩_vacuum, for λ ∈ {0, 0.25, 0.50} at h_E = 1.4 and t = 2.7, with dt = 0.3. c, Probability of a vertex excitation P(A_v) on three distinct vertices A₁ (gold), A₂ (green) and A_vac (black) for λ ∈ {0, 0.25, 0.50} and h_E = 1.4. The grey ‘decoherence-limited’ region is defined by the average of P(A_v) over all vertices having evolved the initial state with the X string for λ = 0, h_E = 1.4. d, Dependence of P(A_v) on h_E, acquired at t = 2 (dt = 0.2), for λ = 0 (pluses), λ = 0.25 (crosses) and λ = 0.50 (pentagons).

To investigate string breaking in our experiment, we first measure the electric excitations through ⟨A_v⟩ in the presence and absence of an initial string excitation in the confined regime h_E = 1.4. In Fig. 5b, we show the difference of ⟨A_v⟩ for these two initial states at time t = 2.7 for different strengths of string breaking set by the coupling λ. When λ = 0, the charge density in the presence of the string is comparable with the one in the absence of the string, ⟨A_v⟩_string − ⟨A_v⟩_vacuum ≃ 0. However, for finite λ values, the time-evolved string initial state has considerably more charges. We track the dynamics of ⟨A_v⟩ for electric excitations at the top and at the bottom of the system. As demonstrated, for this confinement field h_E, the string stays mainly at the top qubits within the accessible timescales. For λ = 0, both vertex operators A₁ on the top and A₂ at the bottom show the same trend as the vacuum state A_vac (Fig. 5c), which can be understood from decoherence of the device (grey region). However, as λ is increased to 0.25 and 0.50, the electric charge A₁ on the top side, at which the string has been created, shows a much higher number of excitations compared with A₂, which remains indistinguishable from the vacuum. This differential measurement is evidence for excitation creation from string breaking.

In the lowest-order string-breaking process, the coupling reduces the length of the string by one, leading to an energy gain of 2h_E and, at the same time, two electric excitations are created with a cost of 4J_E. Therefore, we predict that this energy trade-off could enhance the probability of string breaking near h_E = 2J_E. We measure the probability of electric charge creation P(A_v) on A₁ as a function of h_E (Fig. 5d). For finite couplings, we observe a maximum charge creation in the vicinity of h_E ≈ 2, demonstrating that string breaking is facilitated at the resonance condition.

In this work, we imaged the dynamics of deconfined and confined excitations in a (2+1)D ${{\mathbb{Z}}}_{2}$ LGT and measured the vibrations of the string connecting them. Our work demonstrates the potential for quantum processors to study the dynamics of emergent excitations in correlated quantum matter, which are prohibitively hard to predict theoretically owing to their non-perturbative nature. The dynamical observables measured in this work are related to those measured in more conventional scattering experiments but examine complementary regimes; our approach provides insight into the microscopic nature of non-equilibrium dynamics of quasiparticles, as compared with the measurement of their spectral properties. Our space-resolved and time-resolved measurements provide a new visualization approach for characterizing the dynamics of interacting emergent excitations.

Methods

Experimental techniques and device characterization

Gate implementation

All experiments in this work can be carried out on a grid of 45 qubits with square connectivity and were implemented on a 72-qubit Google Sycamore processor, as used in ref. ⁴⁹ (Extended Data Fig. 1). Dominant errors come from CZ entangling gates⁵⁰ and final readout⁵¹. Qubit, coupler and readout parameters are optimized using the Snake Optimizer^52,53. A smaller contribution to the total error comes from the single-qubit microwave gates. Gates are calibrated using Google’s Optimus calibration tools^54,55.

To mitigate the effects of coherent noise, we implement randomized compiling^56,57,58. For all observables, we average over 30 compiling instances of randomly chosen single-qubit Pauli gates sandwiching each CZ gate, such that the resulting ideal unitary is unchanged. We record approximately 400 shots per instance after post-selection. All sequential single-qubit gates are then combined into a single phased-X Z gate, such that the structure of the circuit is always alternating single-qubit and two-qubit gate layers. The use of randomized compiling to convert coherent errors to incoherent ones also supports our choice to apply simple depolarizing noise mitigation to compare experimental results with numerical simulations,

During the projective state preparation and Hadamard test experiments presented in this work, there are ancilla qubits that must sit unchanged while the rest of the system undergoes up to ten Trotter cycles, which constitute 80 single-qubit layers and 80 two-qubit layers. These long evolution times place a strict constraint on the ancilla qubit to remain highly coherent. Therefore, it is important to mitigate ancilla dephasing to ensure the highest fidelity experiments. To this end, we implement dynamical decoupling whenever a qubit would otherwise be idle. Our approach is to use a simple echo sequence of X gates during each single-qubit gate layer^42,59.

In implementing circuits formulated in terms of non-native C-NOT and SWAP gates, each C-NOT gate is converted to a CZ gate, sandwiched by Hadamard gates on the target qubit. The extra Hadamard gates then get combined with other sequential single-qubit gates into a single phased-X Z gate. SWAP gates, used to initialize central qubits for the Hadamard test in Fig. 4, are broken into three C-NOT gates.

Suzuki–Trotter circuit

In implementing the WALA state and time evolution under ${\mathcal{H}}$, the fourfold connectivity of Google’s Sycamore quantum processor allows for a configuration without ancilla qubits⁴² or with an ancilla qubit at the centre of each vertex and plaquette⁴⁹. Although the first configuration allows for a denser packing of vertex sites, it complicates the Trotter evolution because all plaquette and vertex terms cannot be executed in parallel, increasing errors from T₁ and T₂ decay. Instead, we use the configuration with ancilla qubits. The grid of qubits and the corresponding vertices and plaquettes are shown in Extended Data Fig. 2. The diamonds represent physical qubits on the gauge sites and the circles indicate ancilla qubits. The blue/purple tiles correspond to vertices/plaquettes, respectively.

To carry out the time evolution, we developed a circuit to implement the first-order Suzuki–Trotter expansion (Trotter errors discussed in Supplementary Information Section IV.C). The first operation of the Trotter circuit consists of a local field unitary operator:

$${U}_{{\rm{Fields}}}=\exp \left[-{\rm{i}}\left(-\lambda \sum _{l}{X}_{l}-{h}_{{\rm{E}}}\sum _{l}{Z}_{l}\right){\rm{d}}t\right]$$

(3)

which can be implemented by a single phased-X Z gate on each physical gauge qubit. The second operator of the Trotter circuit involves the vertex and plaquette terms:

$${U}_{{\rm{Plaquettes}}}=\exp \left[-{\rm{i}}\left(-{J}_{{\rm{E}}}\sum _{v}{A}_{v}-{J}_{{\rm{M}}}\sum _{p}{B}_{p}\right){\rm{d}}t\right]$$

(4)

which can be implemented in parallel for all vertices and plaquettes in eight entangling layers. With four layers of entangling gates, the commuting A_v/B_p operators are transformed into single-qubit operators on the ancilla qubits, which are then rotated by an angle  −2J_Edt/−2J_Mdt about the Z axis to invoke the time evolution. The transformation of the vertex and plaquette operators is then reversed with another four layers of entangling gates, which returns the state to the physical qubits and disentangles the ancilla qubits. This algorithm gives a gate count of 16L_xL_y − 12(L_x + L_y) + 8 per Trotter cycle for a rectangular grid with L_x × L_y vertex sites (116 entangling gates per Trotter cycle for our experimental set-up in Extended Data Fig. 2). Executing state preparation and nine Trotter steps, our experiment uses  about 1,075 total entangling gates. Because we are focusing on local observables, we measure appreciable signals in these circuits, whose depth is slightly larger than state-of-the-art random circuit sampling experiments, which focus on cross-entropy benchmarking⁶⁰.

Post-selection

Decoherence is unavoidable on noisy intermediate-scale quantum (NISQ) processors. A common technique to deal with decoherence, which does not depend on any detailed knowledge of the device error model, is post-selection. Any observable that is known to be conserved by the quantum circuit is a good candidate for post-selection criteria, as long as they can be measured concurrently with the final observable of physical interest. In our case, we extract expectation values by measuring the physical qubits. However, we also measure the ancilla qubits concurrently, which would all remain in the |0⟩ state if no errors occurred (Extended Data Fig. 3). This is because circuits for state preparation and Trotterization entangle the ancilla qubits with the physical qubits to carry out the quantum operation, but always disentangle the ancilla and ideally return it to its original state. However, errors during a Trotter step can quickly propagate across the chip. Our measurements show increasing numbers of these errors as the number of Trotter steps increases (Extended Data Fig. 3a,b). Therefore, to mitigate these errors, we post-select all of our measurements presented in the main text and supplement such that all ancilla qubits are in the |0⟩ state.

Even after this first round of post-selection, residual experimental errors, Trotter errors and deviation from perfect overlap between the A_v operators and dressed physical particle operators contribute to experimental measurements of a different number of excitations than initialized. Although the number of vertex parity violations is not an exactly conserved quantity, it is expected to be an approximate one. Therefore, in measurements that scrutinize the properties of a set number of electric excitations, we post-select shots that have the same number of vertex parity flips as the initialized state to further mitigate experimental errors, while also eliminating some of the spurious effects of the Trotter error (Fig. 3a,d,e, Extended Data Figs. 4 and 8 and Supplementary Figs. 3, 4 and 9). An added benefit from post-selecting on the two-charge sector is the ability to unambiguously assign the distance between two charges (Fig. 3a,d).

Although these post-selection techniques increase the accuracy of our results, they come at an exponential cost (Extended Data Fig. 3c). We find that, to obtain a constant number of post-selected shots after applying both ancilla qubit and excitation number post-selection on our standard grid of 17 physical qubits and 18 ancilla qubits, the number of shots taken on the hardware scales as N_shots ≈ 2.5ⁿ, for n Trotter steps. This procedure reduces the rate that observables of a stationary state trend towards the maximally mixed value (Extended Data Fig. 3d).

Global depolarizing channel mitigation

Taking bitstrings directly from the experimental measurements and computing the various observables in this work results in deviations from expectations. For example, results on a 2 × 3 vertex system are shown in Extended Data Fig. 4a. Although there seems to be clear separation between the h_E = 0 data points (red) and the h_E = 2.25 data points (blue), the excitations seem to be drifting apart quickly for all values of h_E.

To estimate the effect of device noise, we perform numerical simulations with a local two-qubit depolarizing noise channel following every CZ gate in our circuit. We estimate the depolarizing probability to be 0.7%, corresponding to the mean error rate attained from cross-entropy benchmarking. The results are shown by the dashed lines in Extended Data Fig. 4a. We see that the local depolarizing error model captures the trend of the data very well for all values of h_E.

To simplify the interpretation and analysis, we consider the possibility that a global depolarizing channel may capture the behaviour of our data phenomenologically^61,62,63,64. To determine the correct error probability to mitigate, we consider Trotterized Hamiltonian evolution under parameters that leave the observable unchanged. For example, when measuring the vertex parity values, time evolution with λ = 0 results in no dynamics in the ideal case, because all of the vertex operators commute with the Hamiltonian. Then, by monitoring the measured charge separation as a function of Trotter cycle, the expectation value drifts from the expectation for the initial state to that of the depolarized state (Extended Data Fig. 4b). The effective depolarizing probability, p_eff, can be extracted for each cycle:

$${p}_{{\rm{eff}}}=\frac{{\langle \widehat{{\mathcal{O}}}\rangle }_{{\rm{measured}}}-{{\mathcal{O}}}_{{\rm{initial}}}}{{{\mathcal{O}}}_{{\rm{depolarized}}}-{{\mathcal{O}}}_{{\rm{initial}}}}$$

(5)

with ${\langle \widehat{{\mathcal{O}}}\rangle }_{{\rm{measured}}}$ being the measured expectation value, ${{\mathcal{O}}}_{{\rm{initial}}}$ being the initial value of the observable, set by the state preparation, and ${{\mathcal{O}}}_{{\rm{depolarized}}}$ being the expectation value in the completely depolarized state, that is, the maximally mixed state. Such a p_eff for the data shown in Extended Data Fig. 4a is shown in Extended Data Fig. 4c. The solid lines in Extended Data Fig. 4a correspond to the application of the global depolarizing channel with the effective depolarizing probability in Extended Data Fig. 4c. The agreement between experiment and the global depolarizing model is reasonable.

To compare experimental measurements of an observable $\widehat{{\mathcal{O}}}$ with noiseless simulations, we use this global depolarizing model to rescale the experimental data:

$${\langle \widehat{{\mathcal{O}}}\rangle }_{{\rm{rescaled}}}=\frac{{\langle \widehat{{\mathcal{O}}}\rangle }_{{\rm{measured}}}-{p}_{{\rm{eff}}}{{\mathcal{O}}}_{{\rm{depolarized}}}}{1-{p}_{{\rm{eff}}}}$$

(6)

Although our experimental data innately show robust signatures of confinement and string dynamics without any rescaling, to compare with numerical simulations, we show rescaled data (always alongside its unscaled counterpart) in Fig. 3a, Extended Data Figs. 9 and 10c and Supplementary Figs. 2, 3, 4c, 5d, 6, 7 and 8.

Two-time Pauli string correlator Hadamard test

In this work, we present two distinct two-time Pauli string correlators of the form ⟨ψ|Z_l(t)Z_l(t₀)|ψ⟩ (Fig. 4) and $\langle \psi | ({X}_{{Q}_{1}}{X}_{{Q}_{2}}\ldots {X}_{{Q}_{j}})(t){X}_{{Q}_{1}}({t}_{0})| \psi \rangle$ (Supplementary Information Section III.C.2). Theoretical works have outlined schemes for measuring such quantities with quantum circuits^65,66 and, recently, experiments have used Hadamard tests with two controlled operations to measure two-time correlators⁶⁷. Because a generic controlled-P operation, in which P is an arbitrary Pauli string, is not necessarily native to our quantum hardware, we measure correlation functions with a version of the Hadamard test with only a single controlled operation, C-A at time t₀ (ref. ⁶⁸). In fact, because both of these correlators are of the form ⟨ψ|B(t)A(t₀)|ψ⟩ with simple operators A = Z_l and $A={X}_{{Q}_{1}}$, respectively, the implementation is straightforward and only requires C-NOT and CZ gates as controlled operations.

To measure the two-time Pauli string correlator, ${\mathcal{C}}(A({t}_{0}),B(t))\,=$ $\langle \psi | B(t)A({t}_{0})| \psi \rangle$, we first prepare the ancilla in the state:

$$| \eta ({\vartheta },\phi )\rangle =\cos \frac{{\vartheta }}{2}| 0\rangle +\sin \frac{{\vartheta }}{2}{{\rm{e}}}^{{\rm{i}}\phi }| 1\rangle $$

(7)

using an arbitrary single-qubit gate. The controlled operator we want to apply to the system to measure ${\mathcal{C}}$ is C-A. Thus, having the grid of qubits starting in state |ψ⟩, we apply C-A to the system, controlled by the ancilla, and have the following resulting state for the entire system (grid plus ancilla):

$$\begin{array}{l}| \Psi \rangle \,=\,\cos \frac{{\vartheta }}{2}| \psi \rangle \otimes | 0\rangle +\sin \frac{{\vartheta }}{2}{{\rm{e}}}^{{\rm{i}}\phi }A| \psi \rangle \otimes | 1\rangle \\ \,\,=\,\frac{1}{\sqrt{2}}\left(\cos \frac{{\vartheta }}{2}{\mathbb{1}}+\sin \frac{{\vartheta }}{2}{{\rm{e}}}^{{\rm{i}}\phi }A\right)| \psi \rangle \otimes | +\rangle \\ \,\,+\,\frac{1}{\sqrt{2}}\left(\cos \frac{{\vartheta }}{2}{\mathbb{1}}-\sin \frac{{\vartheta }}{2}{{\rm{e}}}^{{\rm{i}}\phi }A\right)| \psi \rangle \otimes | -\rangle \end{array}$$

(8)

Then, it can be shown that measuring the expectation value of the Pauli string B × X_a, in which a is the ancilla qubit, results in:

$$\begin{array}{l}\langle \Psi | B(t){X}_{a}(t)| \Psi \rangle =\sin ({\vartheta })\cos (\phi ){\rm{Re}}[\langle \psi | B(t)A({t}_{0})| \psi \rangle ]\\ \,\,\,\,\,\,\,\,\,-\sin ({\vartheta })\sin (\phi ){\rm{Im}}[\langle \psi | B(t)A({t}_{0})| \psi \rangle ]\end{array}$$

(9)

By choosing the initial state of the ancilla to be $| \eta (\frac{\pi }{2},0)\rangle$, we get ⟨ψ|B(t)X_a(t)|ψ⟩ = Re[⟨ψ|B(t)A(t₀)|ψ⟩]. By choosing instead the ancilla in the state $| \eta (\frac{\pi }{2},\frac{-\pi }{2})\rangle$, we get ⟨ψ|B(t)X_a(t)|ψ⟩ = Im[⟨ψ|B(t)A(t₀)|ψ⟩]. These initial states of the ancilla correspond to the usual version of the Hadamard test in which a Hadamard gate applied to the ancilla measures the real part of the controlled unitary and a Hadamard gate times the S-adjoint gate measures the imaginary part.

Symmetry

For the measurement of ⟨Z_l(t)Z_l(0)⟩ on qubit l, we take advantage of the symmetry of the initial state. For this correlator, separate circuits must be run for each qubit. Thus, it is a much more demanding experiment than the measurements of ⟨Z_l⟩, ⟨A_v⟩ or ⟨B_p⟩. To expedite the acquisition of the data shown in Fig. 4c, we only perform the measurement on 10 out of 17 physical qubits and use the vertical mirror symmetry plane to assign the values of the other seven qubits. This is reasonable because the initial state and dynamics respect this mirror symmetry and all qubits are used in the time evolution circuits for each qubit measured. Symmetrization is not used in any other figure.

Variational quantum circuit

In this section, we discuss the variational circuit used in the main text. In recent years, there have been other theoretical proposals based on variational methods for realizing the ground state of LGTs with shallow quantum circuits^69,70,71. Here we show that our protocol is equivalent to a mean-field ansatz for the dual Ising model, for which the expectation values of the local terms of the Hamiltonian can be evaluated analytically as a function of the rotation angle θ (ref. ⁴⁵). Thus, the energy optimization of the variational circuit reduces to finding the minimum of a simple polynomial of trigonometric functions, which can be solved efficiently for arbitrary system sizes. The variational state, proposed in refs. ^45,46, is of the form:

$$| \psi \rangle =\prod _{p}\left(\cos \left(\frac{\theta }{2}\right){\mathbb{1}}+\sin \left(\frac{\theta }{2}\right){B}_{p}\right)| 0\rangle ,$$

(10)

in which the rotation angle θ is the only variational parameter. The idea behind the ansatz is to create a weight-adjustable loop gas: starting from the initial state |0⟩^⊗N, applying B_p flips all spins around this plaquette and creates a closed loop of spins in the |1⟩ state. Applying the operator $\cos \left(\frac{\theta }{2}\right){\mathbb{1}}+\sin \left(\frac{\theta }{2}\right){B}_{p}$ on a plaquette creates a weighted superposition of a closed loop and no loop around that plaquette. In the special case in which both possibilities have the same weight, that is, $\theta =\frac{\pi }{2}$, the circuit prepares the toric code ground state, which is an equal-weight superposition of all closed-loop configurations. Decreasing the angle to tune away from the toric code gives less weight to the configurations with flipped plaquettes—in particular, configurations with large loops are now exponentially suppressed, scaling with the area of the loop (that is, the number of enclosed plaquettes). Such a suppression of large loops is what we might expect when adding an onsite field term such as −h_E∑_lZ_l with h_E > 0 to the toric code Hamiltonian, because now long loops of flipped spins have an energy cost. However, such a suppression should scale with the perimeter of the loop, not with its area. This overpenalization of long loops leads to the fact that the ansatz cannot support topological order in the thermodynamic limit when we tune the angle away from the toric code fixed point at $\theta =\frac{\pi }{2}$. We will see this explicitly in the next section, by mapping the variational circuit to a mean-field ansatz for the dual transverse-field Ising model, in which tuning the angle away from $\theta =\frac{\pi }{2}$ means explicitly breaking the symmetry of the Ising model and thus tuning from the symmetric into the symmetry-broken phase.

Variational circuit as a mean-field ansatz of an Ising model

In the main text, we considered a circuit using ancilla qubits to prepare the variational ansatz because the use of ancillas then reduced the depth of the circuit needed for the time evolution. Here we are only interested in the action of the circuit on the physical system, so we can ignore the ancillas. An alternative (but equivalent) circuit for constructing the state using only the physical qubits is given in Extended Data Fig. 5a,b (refs. ^42,46). Extended Data Fig. 5a shows the repeating circuit element that is applied to every plaquette. It consists of two parts: first, on the top qubit of each plaquette, a y-rotation gate R_y(θ) = exp(−iθY/2) is applied. Then, three C-NOT gates are applied, in which the top qubit acts as the control qubit and the remaining qubits are the target qubits. Extended Data Fig. 5b shows one possible choice of the order of applying the repeating circuit element: the element is applied to each plaquette sequentially, starting from the bottom right, traversing each row right to left, before moving up to the next row. Such an ordering ensures that the rotation gate always acts on a qubit in the |0⟩ state and, thus, the circuit effectively implements the operator $\cos \left(\frac{\theta }{2}\right){\mathbb{1}}+\sin \left(\frac{\theta }{2}\right){B}_{p}$ on each plaquette, which yields the state in equation (10). Within a row, the circuit elements commute and we could choose different orderings. In fact, the presented ordering is not optimal but will simplify the calculations in the following. An optimal ordering is given in ref. ⁴².

To better understand the state prepared by this circuit, we can transform the Hamiltonian in the main text with only the C-NOT gates of the circuit and consider the resulting transformed Hamiltonian. The expectation value of a single Pauli-X operator evaluated in the variational state is always zero, which means that the state is completely insensitive to adding an X-field and it is enough to consider the toric code Hamiltonian with Z-field only. This is because the variational state is a superposition of closed loops only and applying a single Pauli-X operator creates open-ended loops in each state of the superposition. A schematic of the Hamiltonian is shown in Extended Data Fig. 5c. There we show a 4 × 4 lattice of vertex operators. The terms in the Hamiltonian corresponding to the vertex operators are shown in orange, the plaquette terms are shown in blue and the onsite Z-field is coloured in green. Conjugating this Hamiltonian by only the C-NOT gates of the circuit, as we will show below, leads to the Hamiltonian in Extended Data Fig. 5d. There (in the bulk) the vertex terms have shrunk to two-site Ising interactions, the plaquette terms have shrunk to onsite Pauli-X operators and the onsite Z-fields have extended to two-site or three-site Pauli-Z terms. On all sites for which there are no Pauli-X operators, the transformed Hamiltonian commutes with a single-site Pauli-Z operator, so the eigenstates of the Hamiltonian are labelled by product states of |0⟩ or |1⟩ on those sites. The lowest-energy states are given by |0⟩ on all of those sites, as the field term in the Hamiltonian in the main text has a negative prefactor. Note that this is precisely the state of those qubits in the variational circuit before applying the C-NOT gates. Focusing on this subspace, in which all qubits except the top qubit on each plaquette are in the |0⟩ state, we get the Hamiltonian in Extended Data Fig. 5e, which now only lives on the sites on which the y-rotation gates act in the variational circuit. We can see that, on those sites, the Hamiltonian is a two-dimensional transverse-field Ising model. The remaining variational circuit, that is, the variational circuit without the C-NOT gates, simply describes a product state with spins rotated in the x–z plane. Such a state is a mean-field ansatz for the Ising model. In particular, if the rotation angle is $\theta =\frac{\pi }{2}$, all spins are aligned in the x direction, which corresponds to the ground state of the transverse-field Ising model in the limit at which the strength of the Ising interaction is taken to zero. In that case, the ground state is symmetric under the global spin-flip symmetry of the Ising model, that is, under the simultaneous application of a Pauli-X gate on every site. However, if we tune the angle away from $\theta =\frac{\pi }{2}$, the spins are no longer oriented along the x direction and the state is no longer symmetric under spin flips—there is a mean-field phase transition from the symmetric phase into the symmetry-broken phase. (Technically, the single-site Pauli-Z operators at the boundary explicitly break the global spin-flip symmetry of the Ising model in this case; however, in the thermodynamic limit, their contribution to the ground-state energy becomes negligible compared with the bulk and a transition occurs. The only effect of these boundary terms will then be to always favour the same direction of the symmetry-breaking, namely, towards the all |0⟩ state). In the original model of the toric code in a field, this phase transition corresponds to a (mean-field) transition from the topological toric code ground state to the trivial paramagnet. More generally, there is a known duality transformation between the toric code with a Z-field and the two-dimensional transverse-field Ising model on the dual lattice, which relates the topological phase of the toric code to the symmetric phase of the Ising model, and the trivial phase of the toric code with a large Z-field to the symmetry-broken phase of the Ising model^1,4.

To show that the toric code Hamiltonian with a Z-field indeed transforms into an Ising model under the action of the C-NOT gates of the circuit, we first consider the action of the C-NOT gates on a single plaquette, before we then transform the full Hamiltonian plaquette by plaquette. Extended Data Fig. 6a shows all relevant operator transformations on a single plaquette. The transformation of the vertex operators on a plaquette is shown in the left column of the table. Generally, the vertex operator also includes Pauli-Z operators on neighbouring plaquettes; however, they are unaffected by the C-NOT gates applied to the selected plaquette. The right column of the table in the figure shows the transformation of the onsite Z-field. The last line at the bottom of the table shows the transformation of the plaquette operator under the C-NOT gates. Equipped with these transformation rules, we can transform the full Hamiltonian. Note that the first set of C-NOT gates that acts on the Hamiltonian is the last one that is applied to the circuit. So, to transform the Hamiltonian plaquette by plaquette, we need to proceed in the opposite order of how we constructed the circuit in Extended Data Fig. 5b. The process of the transformation is graphically depicted in Extended Data Fig. 6b. At each step, the plaquette highlighted in yellow is transformed next and each plaquette is transformed according to the rules in Extended Data Fig. 6a. The first row in Extended Data Fig. 6b shows the transformation of each plaquette in the top row of the lattice individually. In the second row, because the C-NOT gates applied on plaquettes in the same row commute, we transform a whole row in one step. The final result is the last diagram in the bottom-right corner, which is the same as the Hamiltonian in Extended Data Fig. 5d.

Classically optimizing the variational circuit

The mapping of the variational circuit to the two-dimensional Ising model also helps us to optimize the variational parameter θ in the ansatz, as the expectation values of all terms in the Hamiltonian can be evaluated analytically as a function of θ. Minimizing the energy to find the optimal θ then reduces to finding the minimum of a simple polynomial of trigonometric functions of θ.

As computed in the previous section, the vertex operators map to Ising-type interactions that only have support on sites at which the variational ansatz before the C-NOT layers is in the state |0⟩—see also Extended Data Fig. 6b. Thus, the expectation values of all vertex operators in the variational state is:

$$\langle {A}_{s}\rangle =1.$$

(11)

The plaquette operators map to single-site Pauli-X operators, with support on those sites at which the variational circuit has the y-rotation gates before the C-NOT layers. Thus, for the expectation values of the plaquette operators, we have:

$$\langle {B}_{p}\rangle =\langle 0| {{\rm{e}}}^{{\rm{i}}\theta Y/2}X{{\rm{e}}}^{-{\rm{i}}\theta Y/2}| 0\rangle =\sin (\theta ).$$

(12)

There are two different cases for the Z-field in the original Hamiltonian; in the bulk, it maps to an Ising interaction between two qubits that are acted on with a y-rotation gate and at the boundary, it maps to a single-site Pauli-Z operator on a qubit with a y rotation. Thus, we have:

$$\begin{array}{l}\langle {Z}_{{\rm{bulk}}}\rangle =(\langle 0| {{\rm{e}}}^{{\rm{i}}\theta Y/2}\otimes \langle 0| {{\rm{e}}}^{{\rm{i}}\theta Y/2})(Z\otimes Z)({{\rm{e}}}^{-{\rm{i}}\theta Y/2}| 0\rangle \otimes {{\rm{e}}}^{-{\rm{i}}\theta Y/2}| 0\rangle )\\ \,\,\,\,=\,{\cos }^{2}(\theta )\end{array}$$

(13)

and:

$$\langle {Z}_{{\rm{boundary}}}\rangle =\langle 0| {{\rm{e}}}^{{\rm{i}}\theta Y/2}Z{{\rm{e}}}^{-{\rm{i}}\theta Y/2}| 0\rangle =\cos (\theta ).$$

(14)

As discussed in the previous section, the expectation value of single-site Pauli-X operators in the variational circuit is zero because the state is a superposition of closed loops only.

The energy of the system is then given by the sum of all terms, which for a lattice of L_x × L_y vertex operators is:

$$\begin{array}{l}E\,=\,-{J}_{{\rm{E}}}{L}_{x}{L}_{y}\\ \,\,-{J}_{{\rm{M}}}({L}_{x}-1)({L}_{y}-1)\sin (\theta )\\ \,\,-{h}_{{\rm{E}}}(({L}_{x}-2)({L}_{y}-1)+({L}_{x}-1)({L}_{y}-2)){\cos }^{2}(\theta )\\ \,\,-{h}_{{\rm{E}}}2({L}_{x}-1+{L}_{y}-1)\cos (\theta ).\end{array}$$

(15)

This equation can be minimized efficiently numerically for arbitrary system sizes. For this, we use the default settings of the function scipy.optimize.minimize_scalar from the Python library SciPy⁷², which implements Brent’s algorithm⁷³.

In the limit L_x, L_y → ∞, we can minimize the energy density ${\mathcal{E}}=E/({L}_{x}{L}_{y})$ analytically, as was similarly done in ref. ⁴⁵. We only keep terms in the energy proportional to L_xL_y and obtain:

$${\mathcal{E}}=-{J}_{{\rm{E}}}-{J}_{{\rm{M}}}\sin (\theta )-2{h}_{{\rm{E}}}{\cos }^{2}(\theta ).$$

(16)

Taking the derivative with respect to θ, we have:

$$\frac{{\rm{d}}{\mathcal{E}}}{{\rm{d}}\theta }=-{J}_{{\rm{M}}}\cos (\theta )+4{h}_{{\rm{E}}}\cos (\theta )\sin (\theta )=0.$$

(17)

This equation has two solutions for $0\le \theta \le \frac{\pi }{2}$. The first is given by cos(θ) = 0 or $\theta =\frac{\pi }{2}$, which corresponds to the toric code ground state, and the second is given by $\sin (\theta )=\frac{{J}_{{\rm{M}}}}{4{h}_{{\rm{E}}}}$ or $\theta =\arcsin \,\left(\frac{{J}_{{\rm{M}}}}{4{h}_{{\rm{E}}}}\right)$. Now, we only need to check in which regime which of the two solutions is energetically favourable. For $\theta =\frac{\pi }{2}$, we find:

$${\mathcal{E}}=-{J}_{{\rm{E}}}-{J}_{{\rm{M}}}$$

(18)

and for $\theta =\arcsin \left(\frac{{J}_{{\rm{M}}}}{4{h}_{{\rm{E}}}}\right)$, we find:

$$\begin{array}{l}{\mathcal{E}}\,=\,-{J}_{{\rm{E}}}-{J}_{{\rm{M}}}\frac{{J}_{{\rm{M}}}}{4{h}_{{\rm{E}}}}-2{h}_{{\rm{E}}}\left(1-{\left(\frac{{J}_{{\rm{M}}}}{4{h}_{{\rm{E}}}}\right)}^{2}\right)\\ \,=\,-{J}_{{\rm{E}}}-\frac{{J}_{{\rm{M}}}^{2}}{8{h}_{{\rm{E}}}}-2{h}_{{\rm{E}}}.\end{array}$$

(19)

Comparing the two energies, we find that they are equal when:

$$\begin{array}{c}-{J}_{{\rm{E}}}-{J}_{{\rm{M}}}=-{J}_{{\rm{E}}}-\frac{{J}_{{\rm{M}}}^{2}}{8{h}_{{\rm{E}}}}-2{h}_{{\rm{E}}}\\ \Rightarrow 0=\frac{{J}_{{\rm{M}}}^{2}}{16{h}_{{\rm{E}}}^{2}}-\frac{{J}_{{\rm{M}}}}{2{h}_{{\rm{E}}}}+1={\left(\frac{{J}_{{\rm{M}}}}{4{h}_{{\rm{E}}}}-1\right)}^{2}\\ \Rightarrow {h}_{{\rm{E}}}=\frac{{J}_{{\rm{M}}}}{4}.\end{array}$$

(20)

Thus, for h_E ≤ J_M/4, the variational state with the lowest energy is given by $\theta =\frac{\pi }{2}$, that is, the toric code ground state, whereas for h_E > J_M/4, the optimized parameter is given by $\theta =\arcsin \left(\frac{{J}_{{\rm{M}}}}{4{h}_{{\rm{E}}}}\right)$. This is the functional form of the grey line shown in Fig. 2b.

Further experiments

Absolute initial state energy

After carrying out an uncorrelated readout error mitigation procedure, the corrected values for each term in the Hamiltonian are extracted and plotted in Fig. 2d. Averaging equivalent vertices, plaquettes and links allows for a detailed understanding of how each term in ${\mathcal{H}}$ depends on h_E (Extended Data Fig. 7).

Heatmaps for the |ψ ₊⟩ state

The ability to prepare the superposition states presented in Fig. 3c–e are a distinct advantage of digital quantum simulation. The procedure to create this superposition excitation only adds five C-NOT gates and uses a single ancilla qubit, Q_b (Extended Data Fig. 8a). By first preparing a superposition between the physical system and the ancilla, the measurement of Q_b in the Z basis projects the state onto either the |ψ₊⟩ state or the |ψ₋⟩ state. This measurement of Q_b can be carried out concurrently with the measurements of all other physical and ancilla qubits.

The measurements of ⟨A_v⟩ across the grid show similar results for the |ψ₊⟩ and |ψ₋⟩ states (Extended Data Fig. 8b,c). At time t = 0, electric excitations are equally measured on the four sites touched by the initial excitation, with ⟨A_v⟩ ≈ 0 on each of these sites, consistent with the superposition. As time evolves, the excitations spread more for the deconfined case, h_E = 0, compared with the confining h_E = 2.0 measurements. The conditional location of charges show qualitatively different dynamics in the deconfined phase between the two initial states, which can be attributed to the different quantum interference, as discussed in the main text (Extended Data Fig. 8d). The confined conditional probabilities are similar between the two initial states.

Comparing the excitation separation measured on the device to numerical simulation shows qualitative but not quantitative agreement (Extended Data Fig. 8e). This is probably because of the cancellation of local errors, leading to worse rescaling using the global depolarizing model, as discussed for the distance 2 state in Supplementary Information Section III.B.

Further ⟨Z(t)Z(0)⟩ data

Although ${{\mathcal{S}}}_{ZZ}(t)$ and Re[⟨Z(t)Z(0)⟩] show the distinct behaviours of the string dynamics with the onset of confinement (Fig. 4), we may be interested in the behaviour of other closely related quantities to describe the string dynamics. Thus, in this section, we show experimental results and numerical simulations of ⟨Z(0)⟩, Im[⟨Z(t)Z(0)⟩] and ⟨Z(t)⟩.

Although ⟨Z(0)⟩ can be read from Fig. 4d as the value ${{\mathcal{S}}}_{ZZ}(0)$ (because Re[⟨Z(0)Z(0)⟩] = 1), we explicitly plot those values in Extended Data Fig. 9a. In practice, these values are extracted from the same experiment that yields Re[⟨Z(t)Z(0)⟩] by measuring ⟨Z_a⟩ for a being the ancilla qubit (standard Hadamard test).

To complement the Re[⟨Z(t)Z(0)⟩] shown in Fig. 4, we present the Im[⟨Z(t)Z(0)⟩] (Extended Data Fig. 9b). Oscillations are very apparent for h_E ∈ {0, 0.1, 0.2} in the deconfined phase. These oscillations can be understood in the toric code limit (λ, h_E = 0), in which Z couples the ground state to an excited state with  +4J higher energy. This leads to ⟨Z(t)Z(0)⟩ = e^−4iJt for a bulk qubit. However, for qubits on the boundary, such as Q₁ and Q₂, one of the magnetic excitations is outside the system and does not contribute to the energy difference, therefore ⟨Z(t)Z(0)⟩_edge = e^−2iJt. The imaginary part of ⟨Z(t)Z(0)⟩ is suppressed as h_E is increased to h_E = 1.4. This is also expected, because as h_E → ∞, the ground state becomes an eigenstate of Z. Furthermore, no strong distinction between Q₁ and Q₂ is observed across the measured range of h_E. After rescaling assuming a global depolarizing model, good agreement is observed between experiment and numerical simulations.

Motivated by the fact that the Im[⟨Z(t)Z(0)⟩] is suppressed for large h_E, we also measured ⟨Z(t)⟩ to investigate whether this observable also reveals string dynamics, because the quantum state should approach an eigenstate of Z in the large h_E limit. In this picture, ⟨Z_l(t)⟩ = −1/+1 corresponds to the presence/absence of the Wilson string on qubit l at time t. Our results show behaviour qualitatively similar to ${{\mathcal{S}}}_{ZZ}(t)$ (Extended Data Fig. 9c). We see that, for Q₁, the value of ⟨Z(t)⟩ quickly moves away from the theoretically stationary evolution (grey region) for all values of h_E, which is consistent with the string always moving away from its initial configuration regardless of the degree of confinement. However, for Q₂, in the most confining case, when h_E = 1.4, the dynamics of ⟨Z(t)⟩ are indistinguishable from the stationary evolution, which corroborates our interpretation that the string is not able to move to the bottom qubits for large h_E. After rescaling assuming a global depolarizing model, good agreement is observed between experiment and numerical simulations.

Temporal mapping of vacuum fluctuations and string breaking

Having shown that, in the strongly confining phase, when h_E = 1.4, a string with an initial bump on the top is not able to move to qubits on the bottom of the grid on the experimentally accessible timescales owing to the small matrix elements, we presented signatures of string breaking by comparing the probabilities of finding a vertex excitation at sites at the top and at the bottom of the grid (Figs. 4 and 5). In Fig. 5, we present data after an evolution time of t = 2.7. Possible contributions to this late-time charge density are, for example: (1) residual energy owing to the imperfect approximation of the WALA initial state to the true ground state; (2) the disagreement between A_v and the dressed particle operators for the full Hamiltonian; and (3) device decoherence.

Starting from the WALA initial state and evolving under a Hamiltonian with λ = 0, we expect no electric excitations to appear, because A_v commutes with the Trotterized Hamiltonian. However, in the experiment, we see excitations developing, with the highest density on bulk sites (Extended Data Fig. 10a). This is natural for an experiment on a NISQ processor, as the noise will push the system towards the maximally mixed state, in which ⟨A_v⟩ = 0 for all vertices. The qubits in the bulk take part in the most entangling gates and thus the effects of decoherence can be expected to be strongest for bulk sites. In this case of λ = 0, we observe equivalent results regardless of using an initial state with or without the string excitation, up to experimental errors (Extended Data Fig. 10b). When λ is increased to 0.25 and 0.50, the trend of increasing electric excitations takes on a faster rate for the WALA state. This indicates that the noiseless evolution begins to create pairs of electric excitations owing to reasons (1) and (2) above. In the main text, we call these ‘vacuum fluctuations’, because they are pair-creation events that spawn from evolution of the approximate ground state. When we consider non-zero λ for the string initial state, we see an increased excitation density compared with the WALA state. Indeed, the average probability of finding a vertex excitation on any site is markedly higher for the string initial state, compared with the WALA state, when λ = 0.50 (Extended Data Fig. 10c). After rescaling, assuming a global depolarizing model, we observe almost perfect agreement between our experimental data and noiseless numerical circuit simulations, indicating that this effect is not a spurious result of errors on the quantum processor.

Examining the heatmaps in Extended Data Fig. 10, we see that the extra intensity build-up is concentrated on vertices that the initial string passes through. The asymmetry between excitations on the top and on the bottom is the topic of Fig. 5c and represents strong evidence of string breaking. At the same time, the low probability of electric-charge creation in the WALA state by time t = 2.7 (about 2% for λ = 0.25 and  about 5% for λ = 0.50 after depolarization mitigation) confirms that this ansatz is a good low-energy initial state with only a small amount of intrinsic dynamics.