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

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 \(\mathcalH\), 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:

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:

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:

with \(\langle \widehat\mathcalO\rangle _\rmmeasured\) being the measured expectation value, \(\mathcalO_\rminitial\) being the initial value of the observable, set by the state preparation, and \(\mathcalO_\rmdepolarized\) 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\mathcalO\) with noiseless simulations, we use this global depolarizing model to rescale the experimental data:

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_1X_Q_2\ldots X_Q_j)(t)X_Q_1(t_)| \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, \(\mathcalC(A(t_),B(t))\,=\) \(\langle \psi | B(t)A(t_)| \psi \rangle \), we first prepare the ancilla in the state:

using an arbitrary single-qubit gate. The controlled operator we want to apply to the system to measure \(\mathcalC\) 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):

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:

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:

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)\mathbb1+\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)\mathbb1+\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:

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:

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:

and:

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:

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 \(\mathcalE=E/(L_xL_y)\) analytically, as was similarly done in ref. ⁴⁵. We only keep terms in the energy proportional to L_xL_y and obtain:

Taking the derivative with respect to θ, we have:

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 )=\fracJ_\rmM4h_\rmE\) or \(\theta =\arcsin \,\left(\fracJ_\rmM4h_\rmE\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:

and for \(\theta =\arcsin \left(\fracJ_\rmM4h_\rmE\right)\), we find:

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

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(\fracJ_\rmM4h_\rmE\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 \(\mathcalH\) 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 \(\mathcalS_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 \(\mathcalS_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 \(\mathcalS_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.