Thermalization and criticality on an analogue–digital quantum simulator

Abstract

Understanding how interacting particles approach thermal equilibrium is a major challenge of quantum simulators^1,2. Unlocking the full potential of such systems towards this goal requires flexible initial state preparation, precise time evolution and extensive probes for final state characterization. Here we present a quantum simulator comprising 69 superconducting qubits that supports both universal quantum gates and high-fidelity analogue evolution, with performance beyond the reach of classical simulation in cross-entropy benchmarking experiments. This hybrid platform features more versatile measurement capabilities compared with analogue-only simulators, which we leverage here to reveal a coarsening-induced breakdown of Kibble–Zurek scaling predictions³ in the XY model, as well as signatures of the classical Kosterlitz–Thouless phase transition⁴. Moreover, the digital gates enable precise energy control, allowing us to study the effects of the eigenstate thermalization hypothesis^5,6,7 in targeted parts of the eigenspectrum. We also demonstrate digital preparation of pairwise-entangled dimer states, and image the transport of energy and vorticity during subsequent thermalization in analogue evolution. These results establish the efficacy of superconducting analogue–digital quantum processors for preparing states across many-body spectra and unveiling their thermalization dynamics.

Main

The advent of quantum simulators in various platforms^{8,9,10,11,12,13,14} has opened a powerful experimental avenue towards answering the theoretical question of thermalization^5,6, which seeks to reconcile the unitarity of quantum evolution with the emergence of statistical mechanics in constituent subsystems. A particularly interesting setting is that in which a quantum system is swept through a critical point^15,16,17,18, as varying the sweep rate can allow for accessing markedly different paths through phase space and correspondingly distinct coarsening behaviour. Such effects have been theoretically predicted to cause deviations^19,20,21,22 from the celebrated Kibble–Zurek (KZ) mechanism, which states that the correlation length ξ of the final state follows a universal power-law scaling with the ramp time t_r (refs. ^3,23,24,25).

Whereas tremendous technical advancements in quantum simulators have enabled the observation of a wealth of thermalization-related phenomena^{26,27,28,29,30,31,32,33,34,35}, the analogue nature of these systems has also imposed constraints on the experimental versatility. Studying thermalization dynamics necessitates state characterization beyond density–density correlations and preparation of initial states across the entire eigenspectrum, both of which are difficult without universal quantum control³⁶. Although digital quantum processors are in principle suitable for such tasks, implementing Hamiltonian evolution requires a high number of digital gates, making large-scale Hamiltonian simulation infeasible under current gate errors.

In this work, we present a hybrid analogue–digital^37,38 quantum simulator comprising 69 superconducting transmon qubits connected by tunable couplers in a two-dimensional (2D) lattice (Fig. 1a). The quantum simulator supports universal entangling gates with pairwise interaction between qubits, and high-fidelity analogue simulation of a U(1) symmetric spin Hamiltonian when all couplers are activated at once. The low analogue evolution error, which was previously difficult to achieve with transmon qubits due to correlated cross-talk effects, is enabled by a new scalable calibration scheme (Fig. 1b). Using cross-entropy benchmarking (XEB)³⁹, we demonstrate analogue performance that exceeds the simulation capacity of known classical algorithms at the full system size.

Fig. 1: Analogue–digital simulation with high-precision calibration.

a, Our platform combines analogue evolution with digital gates for extensive state preparation and characterization. b, Schematic of new scalable analogue calibration scheme. Swap (blue) and single-photon (red) spectroscopy is used to extract dressed coupling rates ($\{\widetilde{g}\}$) and qubit frequencies ($\{{\widetilde{\omega }}_{qi}\}$) of two-qubit analogue evolution (U_A), which are converted to bare qubit and coupler frequencies ({ω_qi}, {ω_cj}) through detailed device modelling. The bare frequencies allow for establishing the device Hamiltonian of the full system, which is finally projected to a spin Hamiltonian, H_s.

Leveraging these capabilities, we prepare and characterize states of a 2D XY magnet with broadly tunable energy density, allowing us to study the interplay between quantum and classical critical behaviour in the rich phase diagram of our system. Specifically, we observe finite-size signatures of the Kosterlitz–Thouless topological phase transition—including the emergence of algebraically decaying correlations with exponent near the expected universal value of $\frac{1}{4}$—and demonstrate a breakdown of the KZ mechanism. Our study takes advantage of extensive measurement capabilities to characterize, for example, entanglement entropy for subsystems up to 12 qubits, multi-qubit vortex correlators and energy fluctuations. We also leverage our hybrid analogue–digital scheme (Fig. 1a) to prepare entangled initial states, allowing us to spatially tailor the energy density and vorticity, and investigate the subsequent thermalization dynamics and energy transport.

Operating coupled transmons as a high-fidelity analogue quantum simulator requires precise knowledge of the many-body spin Hamiltonian H_s, which depends on the ‘bare’ frequencies, {ω_qi} and {ω_cj}, of qubits qi and couplers cj. However, experimental calibration is only capable of resolving ‘dressed’ frequencies that—unlike the bare frequencies—change from local (isolated) calibrations to full-scale experiments due to hybridization with neighbouring qubits and couplers. Given this difficulty, past experimental studies^30,31 either suffered from large errors or resorted to multi-parameter optimization protocols that are difficult to scale up. Sophisticated Hamiltonian learning techniques^40,41 can circumvent these issues, but still have potential vulnerabilities to Hamiltonian ramps, noise and errors in state preparation and measurement (SPAM).

In this work, we present a scalable calibration protocol that achieves low error by explicitly calibrating the bare frequencies. As illustrated in Fig. 1b, the protocol begins with two-qubit calibration measurements (single-photon and swap spectroscopy, which is robust to ramps and SPAM errors) to determine the effective coupling $\widetilde{g}$ and dressed qubit frequencies $\{{\widetilde{\omega }}_{qi}\}$ of every qubit pair. Next, we use extensive modelling of the underlying device physics to convert the dressed quantities to the bare frequencies {ω_qi}, {ω_cj}. Finally, a projection technique is applied to approximate our high-dimension device Hamiltonian, H_d({ω_qi}, {ω_cj}), into a spin Hamiltonian, H_s:

$${H}_{{\rm{s}}}=\sum _{i}{\omega }_{i}{n}_{i}+\sum _{\langle i,\,j\rangle }{g}_{ij}({X}_{i}{X}_{j}+{Y}_{i}{Y}_{j})/2+{\mathcal{O}}({g}^{2}/\eta )$$

(1)

where ω_i and ∣g_ij∣ ≈ g are tunable on-site potentials and nearest-neighbour couplings, respectively. The latter is notably smaller than the qubit anharmonicity η ≫ g. This restricts the photon occupation numbers to n_i = 0, 1 and X_i, Y_i are Pauli operators acting in this subspace. The Hamiltonian in equation (1) is in the universality class of an XY model with on-site z-fields. A natural consequence of the hybridization in our system is that H_s contains not only nearest-neighbour hopping, but also density–density interactions and next-nearest-neighbour terms, which scale as order ${\mathcal{O}}$(g²/η) and are typically five to ten times smaller than g (see further details in Methods).

A computationally challenging problem and useful benchmark for the quantum simulator is the thermalization dynamics of an initial Z-basis product state at half-filling, which has high temperature with respect to H_s and hosts many quasiparticles (Fig. 2a). When subject to the (photon number conserving) time evolution operator ${{\rm{e}}}^{-i{H}_{{\rm{s}}}t/\hbar }$ where ħ is the reduced Planck constant (set to 1 hereafter), interactions between quasiparticles are expected to drive the system into a chaotic state. To explore these dynamics, we perform a rapid (6 ns) ramp of the couplings g_ij/2π from 0 to 10 MHz. Quantum chaotic behaviour is then diagnosed by means of Z-basis measurements at different times, yielding a set of probability distributions p_meas(x, t) where {x} represents the set of D ‘bitstrings’ with the same number of photons as the initial state. Figure 2b shows the distribution Pr(p) of p_meas(x, t) for reduced system sizes up to N_q = 25 at t = 5.5/g. In each case, Pr(p) shows a clear exponential decay known as the Porter–Thomas distribution, signalling thermalization to a quantum chaotic state^39,42. By contrast, past studies have found substantial deviations from the Porter–Thomas distribution in other models of analogue dynamics^43,44.

Fig. 2: Fast thermalization dynamics and beyond-classical capabilities in the high-temperature regime.

a, Schematic representation of the experiment: N_q qubits are initialized in a half-filling state, evolved under a Hamiltonian H_s over time t with four instances of disorder in {ω_i}, and finally measured in the Z-basis. b, Distribution Pr(p) of bitstring probabilities p from experiment (coloured bars) at t = 96 ns ≈ 5.5/g and ideal Porter–Thomas distribution Pr(p) = De^−pD (dashed lines). Inset, convergence of the self-XEB with time. c, Time-dependent XEB fidelity for system sizes up to N_q = 35. Inset, system size dependence of ϵ (error per qubit per evolution time of 1/g) from exponential fits. d, Mixed-state entanglement proxy, ${{\mathcal{E}}}_{{\rm{P}}}$, obtained in this and previous studies, plotted against the effective system size ${N}_{\,{\rm{q}}}^{{\rm{eff}}}$ (with respect to entanglement of a fully chaotic state; Supplementary Information) of the respective platforms. Blue pentagons, Sycamore processor in the digital regime^45,48; diamonds, Zuchongzhi processor^46,47; circle, neutral atom analogue simulator⁴⁴; green pentagon, present experiment. ${N}_{\,{\rm{q}}}^{{\rm{eff}}}$ is equal to the actual N_q in the digital experiments, whereas analogue platforms are subject to U(1) conservation (this work) or constraints from Rydberg blockade⁴⁴. Inset, ${{\mathcal{E}}}_{{\rm{P}}}$ as a function of N_q computed from experimental data, including the linear fit used for extrapolation to 69 qubits.

Characterizing the thermalization dynamics through the second moment of the bitstring distribution, also called the self-XEB³⁹, $D{\sum }_{x}{p}_{{\rm{m}}{\rm{e}}{\rm{a}}{\rm{s}}}^{2}(x,t)-1$, we observe its fast convergence to the Porter–Thomas (PT) value of 1 within t_PT ≈ 60 ns (roughly 4/g) for all system sizes (Fig. 2b inset, see Supplementary Information for similar saturation rate of entanglement entropy). The observed fast scrambling dynamics are due to the simultaneously activated couplers, and—compared to an equivalent digital circuit—allow for less shift towards the decohered distribution, Pr(p) = D⁻¹ with self-XEB = 0.

To also characterize the coherent errors from imperfect calibration of H_s, we consider the linear XEB fidelity, $F(t)=D{\sum }_{x}{p}_{{\rm{m}}{\rm{e}}{\rm{a}}{\rm{s}}}(x,t){p}_{{\rm{s}}{\rm{i}}{\rm{m}}}(x,t)-1$, where p_sim are exactly simulated probabilities³⁹. The results, shown in Fig. 2c, show exponential decay after times roughly t_PT, where F accurately describes the state fidelity (Supplementary Information). Fitting the decay, we obtain an error rate of ϵ = 0.10 ± 0.02% per qubit per evolution time of 1/g (one cycle). ϵ is nearly independent of system size up to the largest exactly simulated system, N_q = 35 (inset of Fig. 2c). This indicates the scalability of our calibration protocol and allows extrapolation to the full system size of N_q = 69. Approximate matrix product state (MPS) simulations with bond dimension up to χ = 1,024 were found to be ineffective beyond exactly simulatable system sizes, due to the fast entanglement growth and 2D geometry of our system (Supplementary Information).

The combination of the observed fast dynamics and high fidelity enables quantum simulation of computationally complex states. A representative metric of this capability is the mixed-state entanglement proxy, ${{\mathcal{E}}}_{{\rm{P}}}=$ S^Rényi-1/2_ent$+{\log }_{2}F$, which lower bounds the mixed-state entanglement by accounting for the effects of infidelity on the pure-state Rényi-1/2 entropy⁴⁴. Figure 2d compares the estimated ${{\mathcal{E}}}_{{\rm{P}}}$ of our work and other recent state-of-the-art experiments^{44,45,46,47,48}, in which the proximity to the diagonal (ideal) line measures fidelity, indicating that our platform offers new possibilities for high-accuracy study of highly entangled states. In particular, we estimate that simulations with the level of our experimental fidelity requires more than 10⁶ years on the Frontier supercomputer (Supplementary Information).

Having explored the thermalization dynamics in the high-temperature regime, we next turn to the rest of the rich phase diagram in the XY model (equation (1)), which is expected to show both a quantum phase transition in the ground state and a classical Kosterlitz–Thouless phase transition at finite temperature⁴. To prepare low-energy states of an antiferromagnetic XY magnet, we apply a staggered z-field of magnitude h/(2π) = 30 MHz and initialize the qubits in the Z-basis Neel state, maximizing the energy with respect to the first term in equation (1). We then ramp down the staggered field while simultaneously turning on ferromagnetic couplings of magnitude g_m/(2π) = 20 MHz over a duration t_r (Fig. 3a). Under such a protocol⁴⁹, the system evolution is equivalent to that of an antiferromagnetic XY model with staggered field, initialized in the ground state. This ramp crosses a quantum phase transition between a paramagnetic phase with unbroken U(1) symmetry and the XY-ordered phase at a time t_c ≈ 0.45t_r when h_c/g_c ≈ 1.8(6) (Supplementary Information). The transition, analogous to the 2D Mott insulator–superfluid transition⁵⁰, is in the universality class of a three-dimensional (3D) XY model, with the correlation length and dynamical critical exponents ν ≈ 0.67 and z = 1, respectively. Following the ramp, we rapidly return back to the idle frequencies within 3 ns and perform measurements of correlation functions.

Fig. 3: Critical coarsening in the XY model.

a, Left, experimental schematic of qubit frequencies (blue) and coupling (yellow). Right, phase diagram. Dynamics become diabatic (dashed black) with increased temperature (T) when inverse remaining time (red) exceeds gap (green; Δ ∝ ∣g − g_c∣^νz). QCP and CCP denote the quantum and classical critical phases, respectively. b, The final energy density approaches the ground state value (ε_gs, grey) and Kosterlitz–Thouless transition value (ε_KT, black) as t_r is increased. Red circles and squares indicate single-qubit (SQ) and Bell basis measurements, respectively. Blue, MPS simulation. Purple shading indicates where classical critical behaviour is expected. c, Average correlation, $\bar{G}({\bf{r}})$ (found from averaging (⟨X_iX_j + Y_iY_j⟩ − ⟨X_i⟩⟨X_j⟩ − ⟨Y_i⟩⟨Y_j⟩)/2 over all pairs i, j separated by r) measured at various t_r. d, Decay of radially averaged correlations. Green and purple curves show examples of exponential and power-law fits, respectively, performed up to a distance of 6 to avoid finite-size effects at longer distances. Error bars represent one standard deviation estimated from bootstrapping (N_reps = 5 × 10⁴ repetitions). e, Ratio between r.m.s. errors from power-law and exponential fits (${{\epsilon }}_{{\rm{pow}}}/{{\epsilon }}_{\exp }$) decreases for g_mt_r > 15. f, Power-law exponent, γ, approaches expected value at Kosterlitz–Thouless transition (1/4; black line). g, Vortex density proxy, n_V, decreases to minimum of 2 × 10⁻². h, Correlation length increases with t_r. Both simulation results (blue) and experimental data (red) show substantially more superlinear growth than KZ predictions (dashed black). Diamonds, correlation lengths extracted at expected freezing point (i). i, Correlation length during the ramp, shown with and without rescaling (main and inset, respectively) and two-sided logarithmic axes for ∣τ∣ > τ_KZ. ξ_KZ is found from fitting $\xi (\tau ={\tau }_{{\rm{KZ}}})={\xi }_{0}{({\tau }_{{\rm{KZ}}}/{t}_{r})}^{-\beta }$ with β = 0.9(1) (difference from β = ν = 0.67 expected to be due to finite-size effects). Dashed coloured lines show the theoretically expected scaling, f(x) = ∣x∣^−ν with ν = 0.67 for x < −1 (purple) and a heuristic f(x) = x^η with η = 1 for x > 1 (teal). The increase in ξ beyond the freezing point (diamond) causes deviation from KZ predictions.

Figure 3b shows the ramp time dependence of the average energy density, $\varepsilon ={n}_{{\rm{B}}}^{-1}{\sum }_{\langle i,j\rangle }\langle {X}_{i}{X}_{j}+{Y}_{i}{Y}_{j}\rangle /2$ averaged over n_B = 110 bonds (N_q = 65) and corrected for readout errors (Methods). As t_r increases and the dynamics become more adiabatic, we observe a decrease in energy density towards the theoretically predicted ground state value of ε_gs = −0.56, as well as the predicted Kosterlitz–Thouless (KT) transition energy density, ε_KT = −0.53 ± 0.01 (grey and black lines, respectively). As demonstrated below, the final states are thermalized to a strong extent, so ε can be used to evaluate the final effective temperature. To correct for photon decay errors, we apply digital entangling gates at the end of the circuit to convert each pair of qubits to the Bell basis (Methods). This allows for postselecting with respect to photon number conservation (red squares), which yields an improved value of ε = −0.53 ± 0.01, roughly equal to the Kosterlitz–Thouless transition point. The remaining discrepancy from ε_gs is attributed to dephasing effects, which are not corrected by this technique.

As the energy itself does not reveal the effects of thermalization, we next turn to correlations at longer distances and consider the average correlation, $\bar{G}({\bf{r}})$, between pairs of qubits separated by r, shown in Fig. 3c. We observe antiferromagnetic ordering, with the range and magnitude of correlations increasing notably with ramp time, as expected for states with decreasing energy. We next compute the radial average, $\bar{G}(| {\bf{r}}| )$, and fit the resulting decay profiles with exponential fits to extract the correlation length, ξ, as well as with power-law fits to evaluate the type of distance-scaling (Fig. 3d). At short ramp times, the correlations are found to decay exponentially, as theoretically expected for states above the Kosterlitz–Thouless transition, in which freely proliferating vortices preclude long-range order. At longer ramp times, on the other hand, the decay behaviour is better described by power-law fits, as shown in Fig. 3e; specifically, we observe a marked decrease in the ratio between the root-mean-square (r.m.s.) errors of power-law and exponential fits to well below 1 near g_mt_r = 25, where the energy is also close to its minimum value. This behaviour is consistent with that expected in the classical critical regime, where free vortices become entropically unfavourable and are replaced by bound vortex–antivortex pairs, leading to algebraically decaying correlations. (We note that finite-size scaling analysis of the Kosterlitz–Thouless transition is challenging, owing to characteristic rapid growth of the correlation length, and is not attempted here.) In the region with good power-law agreement, we extract a power-law exponent of γ = −0.29 (Fig. 3f), close to the theoretically expected universal value of $-\frac{1}{4}$ at the Kosterlitz–Thouless transition⁵¹. To further substantiate our interpretation, we also measure four-qubit correlators to construct the Swendsen proxy for the vortex density⁵², given by ${n}_{{\rm{V}}}=\frac{1}{4{N}_{{\rm{P}}}}{\sum }_{{\rm{i=1}}}^{{N}_{{\rm{P}}}}(1-{X}_{i1}{X}_{i3}-{Y}_{i1}{Y}_{i3})(1-{X}_{i2}{X}_{i4}-{Y}_{i2}{Y}_{i4})$ for plaquettes i = 1, ‥, N_P with vertices {i1, i2, i3, i4}. Indeed, we find a rapid decrease in n_V as t_r is increased (Fig. 3g), to a minimum value of 2 × 10⁻² in the low-energy regime.

Having studied the classical critical behaviour, we next explore the scaling of the correlation length with the duration t_r over which we sweep through the quantum critical point (Fig. 3h). The correlation length rises to a maximum of ξ ≈ 10 at g_mt_r = 25, which is equal to the longest dimension of our system. At long ramp times, we observe a slight decrease in ξ, attributed to qubit decoherence, as well as periodic oscillations. The latter are also observed in MPS simulations and expected to be caused by finite-size effects. Focusing on shorter ramp times for which these additional effects are absent and the correlations show a more clear exponential decay, we observe strong deviation from the power-law scaling with exponent ν/(1 + νz) = 0.4 predicted by KZ theory. Specifically, ξ grows substantially more superlinearly, and clear discrepancies from power-law scaling are observed in both experiment and simulation. We attribute the observed breakdown of KZ scaling to coarsening beyond the expected freezing point^19,21.

To demonstrate this more explicitly, we measure the correlation length along the Hamiltonian ramp (Fig. 3i). The KZ prediction assumes that the dynamics freeze at t_KZ, when the inverse gap exceeds ∣t − t_c∣. By contrast, we find that ξ continues to increase, suggesting that the system is instead able to further thermalize, thus causing a different correlation length at the end of the ramp. To illuminate this point, we plot the experimentally measured correlation lengths at t_KZ in Fig. 3h and find better agreement with the KZ prediction. Notably, by rescaling to ξ/ξ_KZ and (t − t_c)/∣t_KZ − t_c∣ ≡ τ/τ_KZ, we find that the curves collapse to a common f(τ/τ_KZ), consistent with predictions of universal coarsening dynamics^19,21,22. The collapse extends well beyond the quantum critical regime, −τ_KZ < τ < τ_KZ, indicative of dynamical universality driven by coarsening. We observe behaviour similar to the theoretically predicted f(x) ∝ ∣x∣^−ν for x < −1 (small deviations probably related to effects of finite size and short ξ), and f(1)/f(−1) = 2.3 ± 0.1. We heuristically find approximately $f(x)\propto x$ for x > 1, in which the interplay of gapped and gapless modes is expected to cause different behaviour from quantum Ising models.

Thus far, we have tuned the final energy through the ramp rate of the Hamiltonian. To further study thermalization, as well as the scaling relations near the Kosterlitz–Thouless transition, we prepare a variable number of excitations, n₀ (pairs of spin flips in randomized locations) in the initial state⁵³. Whereas we find that the final average energy density depends linearly on n₀ (Fig. 4a), the behaviour of the correlation length is more intricate (Fig. 4b) and is best understood by plotting ξ versus energy density for all n₀ and g_mt_r > 5 (Fig. 4c). Notably, the points show a collapse (also for n_V in inset), suggesting that the final states are well thermalized, such that the energy density determines ξ and n_V, as expected from the eigenstate thermalization hypothesis (ETH)^5,6. Barring ξ near the system size, we find that $\log \,\xi$ is nearly linear in ∣ε − ε_KT∣^−0.5, as expected near the Kosterlitz–Thouless transition. This is incompatible with naive KZ scaling, and thus further corroborates its breakdown.

Fig. 4: Tunable thermalized states through initial excitations.

a,b, Energy density (a) and correlation length (b) versus number of initial excitations, n₀, averaged over three randomizations. Error bars smaller than markers (N_reps = 6 × 10³ per data point). c, Energy dependence of ξ, demonstrating collapse when data from sweeps of t_r and n₀ are plotted together. $\log \,\xi$ is near-linear in ∣ε − ε_KT∣^−0.5 as theoretically expected. Inset, vortex density versus energy density, showing similar collapse. Error estimated from bootstrapping in c and e. d, t_r-dependence of energy density and its fluctuations (width), for various n₀. e, Energy fluctuations versus energy density, showing an absence of collapse becase σ_ε does not thermalize. f, Second Rényi entropy versus subsystem size for various n₀ at t_r = 200 ns ≈ 23/g. Increasing n₀ causes transition from area- to volume-law behaviour, also seen from the extracted contributions (inset, error bars represent one standard deviation across five excitation randomizations and four subsystems).

Although thermalization causes states created with different n₀ and t_r to have the same observables (for example, n_V and ξ) if their final energies are equal, the states themselves are not necessarily the same. This can be seen by studying observables such as the energy fluctuations, ${\sigma }_{\varepsilon }={({n}_{{\rm{B}}}{g}_{{\rm{m}}})}^{-1}\sqrt{\langle {H}_{XY}^{2}\rangle -{\langle {H}_{XY}\rangle }^{2}}$ with H_XY = ∑_⟨i,j⟩(X_iX_j + Y_iY_j)/2, which trivially commute with H_XY and are thus not thermalized under ETH. We next reconstruct σ_ε from two- and four-qubit correlators (Methods) and find that it decreases from 0.07 to 0.02 as we approach the ground state for n₀ = 0, whereas its dependence on n₀ is weaker (Fig. 4d). The low value of σ_ε compared to the tunable energy range indicates our ability to probe specific parts of the spectrum. Notably, when the full dataset across t_r and n₀ is plotted against energy density, the points do not collapse (Fig. 4e). This shows the difference in states accessed by the two tuning techniques, which was previously concealed by the thermalization of n_V and ξ.

To further characterize the degree of thermalization, we leverage the fast data acquisition rate of our platform to measure the entanglement entropy for subsystem sizes up to 12 qubits, using randomized measurements⁵⁴. At n₀ = 0, we find area-law behaviour (Fig. 4f), which, up to a subleading logarithmic contribution, is consistent with predictions for low-energy states in the XY model⁵⁵. However, tuning to higher final energies by means of larger n₀, we find a continuous crossover to volume-law behaviour (area- and volume-law components in inset), as is expected from ETH for thermalized states at finite energy density^2,36.

We have so far observed signatures of thermalization in the final state of the dynamics, but the thermalization dynamics themselves are still left unexplored. Although we have shown that t_r and n₀ are effective for realizing and studying states with a desired energy and energy fluctuations, they are limited when it comes to studying spatiotemporal dynamics; to study a state with substantial correlations (⟨XX⟩ > 0.1), a ramp time of more than roughly 1/g is required, at which point the system is typically already near equilibrium. Moreover, although these knobs allow for tuning energy density and antiferromagnetic correlations, quantities such as vorticity are out of reach.

Next, we therefore expand the capabilities of our platform by combining the analogue evolution with entangled state preparation by means of high-fidelity (digital) two-qubit gates (Fig. 5a,b). Following the preparation of the dimer state, ${(| 01\rangle -| 10\rangle )}^{\otimes {N}_{{\rm{q}}}/2}$, we rapidly turn on H_s with g/(2π) = 10 MHz and observe very fast thermalization of the energy density on a timescale of just around 1.5/g (Fig. 5c). As the system thermalizes, the range of correlations increases rapidly (Fig. 5d), converging to a correlation length of roughly 1.0 (Fig. 5e). As is expected from ETH, this is in good agreement with ξ ≈ 1.1 observed for the same energy density (−0.23g) when tuning t_r and n₀.

Fig. 5: Transport and thermalization dynamics with entangled initial states.

a, Dimer states are prepared using digital gates, and their thermalization and transport dynamics are realized with analogue evolution, before finally measuring energy, spin current and vorticity. b, We prepare dimer states with spatially tunable phase, ϕ. Energy gradients between ϕ = 0 (ε > 0) and ϕ = π (ε < 0) drive energy current, whereas ϕ = π/2 gives non-zero spin current and vorticity. c,d, Time evolution of energy density (c) and correlations (d) after dimer preparation demonstrate rapid thermalization. e, Correlations become increasingly long-ranged as the system thermalizes. Dashed line, exponential fit. f,g, Energy density (f) and energy gradients (g) after dimer preparation with ϕ = 0 and π in the left and right halves of the system, respectively, showing energy transport on much longer timescales. Colour and length scales of arrows in g and i are logarithmic. h, Time dependence of average energy density along various vertical cuts (coloured circles) and energy imbalance across x = 5 (black circles), showing very good agreement with diffusion model (dashed lines). Error bars smaller than markers. i,j, Spin current (i) and vorticity (j) for ϕ = π/2, showing rapid thermalization. k, The r.m.s. vorticity shows initial slow dynamics followed by near-exponential decay with rate Γ = 49 MHz = 0.85g (fit shown by dashed line). Error bars in e and k are estimated by bootstrapping (N_reps = 10⁴).

Next we leverage the tunability of the phases of the initial dimer states to enable study of transport (Fig. 5b). Specifically, we prepare the dimers in one half of the device in the higher-energy dimer state, $| 01\rangle +| 10\rangle$ (Fig. 5f). Now the dynamics are found to be substantially slower, with clear spatial non-uniformity remaining even after 23 cycles. We also plot the energy density gradient in Fig. 5g, which quickly establishes a uniform field in the +x direction. Figure 5h shows the time dependence of the average energy density at various vertical cuts (coloured circles), as well as the total energy transfer across x = 5 (black circles), which both show excellent agreement with a diffusion model (dashed lines). The energy transport is indeed expected to be diffusive due to the relatively high energy of the dimer state. The data allow for extracting a diffusion constant of D = 29.6 MHz = 0.52g.

The use of initial entangled states in our hybrid analogue–digital platform enables not only tailoring the initial energy landscape, but also other observables such as vorticity and spin current. We achieve this by further tuning the initial dimer phases to π/2 (Fig. 5b). This gives rise to local spin currents, ⟨X_iY_i+1 − Y_iX_i+1⟩/2 ≠ 0, and a sea of vortices and anti-vortices, quantified by the vorticity, ${V}_{i}=\frac{1}{4}({X}_{i1}{Y}_{i2}-{Y}_{i2}{X}_{i3}+{X}_{i3}{Y}_{i4}-{Y}_{i4}{X}_{i1})$ for each plaquette i with vertices {i1, i2, i3, i4}. The temporal evolution of the spin current and vorticity is presented in Fig. 5i,j, respectively, showing thermalization on a fast timescale similar to that in Fig. 5c. Specifically, after an initial super-exponential decay, the r.m.s. vorticity decays near-exponentially with a rate of Γ = 49 MHz = 0.85g (Fig. 5k).

Our results demonstrate a high-fidelity quantum simulator with the capability of emulating beyond-classical chaotic dynamics, a wide range of characterization probes and versatile analogue–digital control. Leveraging these features has enabled new insights about the rich interplay of quantum and classical critical behaviour in the 2D XY model, including the Kosterlitz–Thouless transition, thermalization dynamics and a breakdown of the KZ scaling relations. The effects of the co-existing gapped longitudinal modes and gapless (finite-size limited) transverse modes, specifically on the coarsening dynamics, is of particular interest for future theoretical study. Looking ahead, the new platform presented here is expected to offer an invaluable playground for studies of classically intractable many-body quantum physics, including, for example, dynamical response functions and magnetic frustration. Finally, we note that during the preparation of this paper, we became aware of a related work studying coarsening near an Ising quantum phase transition with Rydberg atoms⁵⁶.

Methods

Device details

The experiments are performed on a superconducting quantum processor with frequency-tunable transmon qubits and couplers, with a similar design to that in ref. ⁴⁵. Extended Data Fig. 1a,b show the measured Ramsey dephasing (${T}_{2}^{* }$) and photon relaxation (T₁) times at the interaction frequency of 5.93 GHz used in our experiments, with median values of 2.0 and 18.8 μs, respectively. Characterizing our digital gate performance, we find a median Pauli error of 4.5 × 10⁻³ for combined $\sqrt{\text{iSWAP}\,}$ and single-qubit gates (Extended Data Fig. 1c), and 1.0 × 10⁻³ for single-qubit gates alone (Extended Data Fig. 1d). Finally, Extended Data Fig. 1e shows our readout errors, with a median of 1.4 × 10⁻².

Analogue calibration

In this section, we describe our new, scalable analogue calibration framework that enables roughly 0.1% cycle error per qubit. To achieve a scalable scheme, we perform pairwise calibration measurements—specifically single-photon and swap spectroscopy—which allows for accurately setting the effective coupling $\widetilde{g}$ and dressed qubit frequencies ${\widetilde{\omega }}_{qi}$ in each qubit pair. A key challenge in analogue calibration that contrasts with its digital counterpart is that these dressed quantities in the pairwise scenario change drastically when all couplers are turned on in the fully coupled global case. Therefore, we perform extensive modelling of the device physics to accurately convert them to the bare qubit and coupler frequencies, {ω_qi},{ω_cj}, which, crucially, do not change from the local calibration measurements to the full-scale experiments.

Model device Hamiltonian

We model both the qubits and couplers in our tunable coupler architecture as Kerr oscillators, with four or five levels in each transmon, depending on the number of photons involved in the Hamiltonian term of interest. Specifically, in calculations concerning one- and two-photon terms, we include four and five levels, respectively. This is done to account for effects that do not obey the rotating-wave approximation, which couple $| 1\rangle$ to $| 3\rangle$ and $| 2\rangle$ to $| 4\rangle$. To ensure high accuracy, we account for not only coupling terms between neighbouring qubits and couplers, but also diagonal pathways, including between couplers:

$$\begin{array}{c}\genfrac{}{}{0ex}{}{{\rm{S}}{\rm{i}}{\rm{n}}{\rm{g}}{\rm{l}}{\rm{e}}\,{\rm{q}}{\rm{u}}{\rm{b}}{\rm{i}}{\rm{t}}}{{H}_{{\rm{d}}}\,=\,\overbrace{\sum _{qi}{\omega }_{qi}{\hat{n}}_{qi}-{\eta }_{qi}{\hat{n}}_{qi}({\hat{n}}_{qi}-1)/2}}\\ \,\,\genfrac{}{}{0ex}{}{{\rm{S}}{\rm{i}}{\rm{n}}{\rm{g}}{\rm{l}}{\rm{e}}\,{\rm{c}}{\rm{o}}{\rm{u}}{\rm{p}}{\rm{l}}{\rm{e}}{\rm{r}}}{+\overbrace{\sum _{cj}{\omega }_{cj}{\hat{n}}_{cj}-{\eta }_{cj}{\hat{n}}_{cj}({\hat{n}}_{cj}-1)/2}}\\ \,\,\genfrac{}{}{0ex}{}{{\rm{Q}}{\rm{u}}{\rm{b}}{\rm{i}}{\rm{t}}\,-\,{\rm{q}}{\rm{u}}{\rm{b}}{\rm{i}}{\rm{t}}\,{\rm{c}}{\rm{o}}{\rm{u}}{\rm{p}}{\rm{l}}{\rm{i}}{\rm{n}}{\rm{g}}}{+\overbrace{\sum _{qi,qj}\frac{1}{2}{\mathop{k}\limits^{ \sim }}_{qi,qj}\sqrt{{\omega }_{qi}{\omega }_{qj}}{\hat{Q}}_{qi}{\hat{Q}}_{qj}}}\\ \,\,\genfrac{}{}{0ex}{}{{\rm{Q}}{\rm{u}}{\rm{b}}{\rm{i}}{\rm{t}}\,-\,{\rm{c}}{\rm{o}}{\rm{u}}{\rm{p}}{\rm{l}}{\rm{e}}{\rm{r}}\,{\rm{c}}{\rm{o}}{\rm{u}}{\rm{p}}{\rm{l}}{\rm{i}}{\rm{n}}{\rm{g}}}{+\overbrace{\sum _{qi,cj}\frac{1}{2}{\mathop{k}\limits^{ \sim }}_{qi,cj}\sqrt{{\omega }_{qi}{\omega }_{cj}}{\hat{Q}}_{qi}{\hat{Q}}_{cj}}}\\ \,\,\genfrac{}{}{0ex}{}{{\rm{C}}{\rm{o}}{\rm{u}}{\rm{p}}{\rm{l}}{\rm{e}}{\rm{r}}\,-\,{\rm{c}}{\rm{o}}{\rm{u}}{\rm{p}}{\rm{l}}{\rm{e}}{\rm{r}}\,{\rm{c}}{\rm{o}}{\rm{u}}{\rm{p}}{\rm{l}}{\rm{i}}{\rm{n}}{\rm{g}}}{+\overbrace{\sum _{ci,cj}\frac{1}{2}{\mathop{k}\limits^{ \sim }}_{ci,cj}\sqrt{{\omega }_{ci}{\omega }_{cj}}{\hat{Q}}_{ci}{\hat{Q}}_{cj}}},\end{array}$$

(2)

where $\widehat{Q}={a}^{\dagger }+a$ and the $\widetilde{k}$ are the effective coupling efficiencies between transmons, including both direct and indirect capacitive contributions (note that the indirect contributions should not be confused with contributions due to virtual exchange interactions, which are included indirectly when we project out the couplers later on). The coupling efficiencies for the various terms can be summarized as follows:

For k_qq, we include three types of qubit–qubit coupling, distinguished by the relative positioning of the qubits. Notably, the geometry of the transmons breaks the 90° rotational symmetry; specifically, the couplings differ along the northwest–southeast (NW–SE) and northeast-southwest (NE–SW) directions. To discuss the three types of coupling, we consider the four qubits on a plaquette shown in Extended Data Fig. 2 and consider examples of pairs of transmons (the formulas for the remaining pairs are given by reflection symmetry about the NW–SE and NE–SW axes, for example, ${\widetilde{k}}_{{q}_{1},{c}_{23}}={k}_{{q}_{1},{c}_{23}}+2{k}_{{q}_{1},{q}_{2}}{k}_{{q}_{2},{c}_{23}}$ infers that ${\widetilde{k}}_{{q}_{1},{c}_{34}}={k}_{{q}_{1},{c}_{34}}+2{k}_{{q}_{1},{q}_{4}}{k}_{{q}_{4},{c}_{34}}$):

• (1) Nearest-neighbours qubits, q₁ and q₂ separated by a coupler c₁₂: ${\widetilde{k}}_{{q}_{1},{q}_{2}}={k}_{{q}_{1},{q}_{2}}+{k}_{{q}_{1},c}{k}_{{q}_{2},c}$.

• (2) Diagonally separated qubits in the NW–SE direction, q₁ and q₃: ${\widetilde{k}}_{{q}_{1},{q}_{3}}={k}_{{q}_{1},{q}_{3}}+2({k}_{{q}_{1},{q}_{2}}{k}_{{q}_{2},{q}_{3}}+{k}_{{q}_{1},{q}_{4}}{k}_{{q}_{4},{q}_{3}})$.

• (3) Diagonally separated qubits in the NE–SW direction, q₂ and q₄: ${\widetilde{k}}_{{q}_{2},{q}_{4}}={k}_{{q}_{2},{q}_{4}}$.

For k_qc, we also include three types of qubit–coupler coupling:

• (1) Nearest-neighbours: ${\widetilde{k}}_{{q}_{1},{c}_{1}}={k}_{{q}_{1},{c}_{1}}$.

• (2) Diagonally separated qubit and coupler in the NW–SE direction, q₁ and c₂₃: ${\widetilde{k}}_{{q}_{1},{c}_{23}}={k}_{{q}_{1},{c}_{23}}+2{k}_{{q}_{1},{q}_{2}}{k}_{{q}_{2},{c}_{23}}$.

• (3) Diagonally separated qubit and coupler in the NE–SW direction, q₄ and c₁₂: ${\widetilde{k}}_{{q}_{4},{c}_{12}}={k}_{{q}_{4},{c}_{12}}$.

For k_cc, we consider two types of coupler–coupler coupling:

• (1) Diagonally separated couplers in the NW–SE direction c₁₂ and c₂₃: ${\widetilde{k}}_{{c}_{12},{c}_{23}}={k}_{{c}_{12},{c}_{23}}+2{k}_{{c}_{12},{q}_{2}}{k}_{{q}_{2},{c}_{23}}$.

• (2) Diagonally separated qubit and coupler in the NE–SW direction, c₁₂ and c₁₄: ${\widetilde{k}}_{{c}_{12},{c}_{14}}={k}_{{c}_{12},{c}_{14}}$.

Calibration experiments

To calibrate the bare qubit and coupler frequencies for a given set of applied biases, we perform various types of calibration measurements (Extended Data Fig. 3a):

Ramsey spectroscopy

In this measurement, we perform standard Ramsey spectroscopy for a range of applied qubit bias values, while keeping the couplers turned off and the neighbouring qubits detuned, to prevent swapping.

Swap spectroscopy

This measurement is performed on a pairwise level, in which neighbouring couplers (except the one connecting the pair) are turned off. The two qubits are prepared in the $| 10\rangle$-state and we measure the swap rate as a function of detuning between the two qubits (Extended Data Fig. 3b). The minimum swap rate tells us the effective coupling between the two qubits, $\widetilde{g}$, and the detuning at which this occurs equals the difference between the dressed frequencies of the qubits, ${\widetilde{\omega }}_{{q}_{1}}-{\widetilde{\omega }}_{{q}_{2}}$ (Extended Data Fig. 3c). Using an iterative scheme, we calibrate the coupler bias required to achieve the target effective coupling.

Single-photon spectroscopy

Whereas the swap spectroscopy provides us with the difference of the dressed frequencies, we also need to find their sum to determine the individual values, ${\widetilde{\omega }}_{{q}_{1}}$ and ${\widetilde{\omega }}_{{q}_{2}}$. We achieve this by preparing the qubits in $(| 1\rangle +| 0\rangle )| 0\rangle /\sqrt{2}$ and measuring ⟨X + iY⟩ as a function of evolution time (Extended Data Fig. 3d). The Fourier transform of the signal then reveals the eigenfrequencies of the two-qubit system, the average of which is equal to $({\widetilde{\omega }}_{{q}_{1}}+{\widetilde{\omega }}_{{q}_{2}})/2$ (Extended Data Fig. 3e).

Next, using separately calibrated coupling efficiencies, we model all the calibration experiments above with the device Hamiltonian described earlier, to find the bare qubit and coupler frequencies that give the dressed quantities observed in the calibration experiment. We model not only the two qubits and the coupler involved in pairwise experiments (single qubit involved in Ramsey), but also the neighbouring ‘padding’ qubits and couplers to account for their effects. Therefore, we start by determining the bare idle frequencies, {ω^idle}, because these must be known to represent the padding in the interaction configuration.

Projection onto computational subspace

Considering the fact that our model device Hamiltonian involves both qubits and couplers with up to five levels in each, it is computationally intractable to use it for time evolution even at small photon numbers. Moreover, in this form, it is very difficult to map its behaviour onto physically relevant systems. We therefore perform a projection technique to convert the device Hamiltonian into a spin Hamiltonian, H_s, that acts on the computational subspace. To find spin Hamiltonian terms involving n photons in a system of N_q qubits, we write ${H}^{(n)}={\sum }_{i,j}| i\rangle \langle i| {H}_{{\rm{d}}}| \,j\rangle \langle j|$, where $\{| i\rangle \}$ are our ${N}_{n}=\left(\begin{array}{c}{N}_{{\rm{q}}}\\ n\end{array}\right)$ new dressed n-photon basis states.

Let us now motivate our choice of dressed basis states, by considering a few different options. One option could have been to simply use the bare qubit states, $\{{| i\rangle }_{{\rm{bare}}}\}$; however, this would cause the spin Hamiltonian to have different eigen-energies from the low-energy spectrum of H_d. A second option would be to instead use the N_n lowest-energy n-photon eigenstates of ${H}_{{\rm{d}}},\{{| i\rangle }_{{\rm{eigen}}}\}$. In this case, the spin Hamiltonian is guaranteed to have the same N_n lowest n-photon eigen-energies as H_d. However, these basis states are highly delocalized and poorly represent our qubits. Hence, to get the best of both worlds, we turn to a third option, in which we project the bare qubit states onto the low-energy eigenspace spanned by $\{{| i\rangle }_{{\rm{eigen}}}\}$. These projections are not orthonormal, so we perform singular value decomposition and set the singular values to one to arrive at our new dressed basis states. It can be shown that this is the most localized set of states that still preserve the low-energy eigenvalues⁵⁷. These new basis states are slightly delocalized on the nearest couplers and qubits, and also have a weak overlap with states that have n + 2 and n − 2 photons due to terms beyond the rotating-wave approximation. We note that our typical coupler ramp times of more than 5 ns are sufficient to ensure adiabatic conversion between the bare qubit states (in which we perform state preparation and measurement) and the dressed basis states that are relevant under analogue evolution.

The spin Hamiltonian H⁽ⁿ⁾ found from the technique above in principle includes all terms involving ≤n photons, including very long-range interactions; however, they drop off rapidly with the photon–photon separation d (typically as (g/η)^d ~ 0.1^d). Moreover, we also find that the terms decay with the number of involved photons in a similar way. Hence, to achieve the low error demonstrated in our manuscript, it is sufficient to include only terms involving up to two photons, and where all the involved qubits are a maximum Manhattan distance of two sites apart, resulting in:

$$\begin{array}{l}H=\sum _{i}{\omega }_{i}{n}_{i}+\sum _{\langle i,\,j\rangle }{g}_{ij}({X}_{i}{X}_{j}+{Y}_{i}{Y}_{j})/2+\sum _{\langle i,\,j\rangle }{g}_{ij}^{nn}{n}_{i}{n}_{j}\\ \,+\,\sum _{\langle i,\,j,\,k\rangle }({g}_{ijk}^{XnX}{n}_{j}+{g}_{ijk}^{XIX})({X}_{i}{X}_{k}+{Y}_{i}{Y}_{k})/2\\ \,+\,\sum _{\langle i,\,j,\,k\rangle }({g}_{ijk}^{nXX}{n}_{i}({X}_{j}{X}_{k}+{Y}_{j}{Y}_{k})/2,\end{array}$$

(3)

where ${g}_{ij}^{nn},{g}_{ijk}^{XnX}$ and ${g}_{ijk}^{XIX}$ scale as g²/η, while ${g}_{ijk}^{nXX}$ scales as g³/η² and qubits i, j, k are connected (Extended Data Fig. 4). We note that there is an offset to these scaling behaviours, which arises due to the diagonal capacitive coupling. This is particularly evident for terms involving qubits along the NW–SE diagonal, because the diagonal coupling is strongest there.

Our technique requires finding the N_n lowest-energy n-photon eigenstates of H_d, which has a high computational cost for large N_q. Fortunately, for a given Hamiltonian term involving a certain set of qubits, the effect of other transmons decays quickly with distance, and we only need to include the nearest neighbouring qubits and couplers to achieve accuracies on the tens of kHz scale. To find the spin Hamiltonian terms, we therefore scan through various subsystems and perform the procedure outlined above for each of them.

Phase calibration for hybrid analogue–digital experiments

In experiments in which we prepare an entangled initial state, the frequency trajectories of the qubits lead to phase accumulation that must be characterized and corrected through phase gates, both before and after the analogue evolution (Extended Data Fig. 5a). Specifically, in the frame that rotates at the interaction frequency, the qubits in each dimer pair precess relative to each other before they reach the interaction frequency. Hence, a phase rotation ϕ_0,i must be applied to every qubit before turning on the analogue Hamiltonian to ensure that the dimer pairs have the desired phase difference when the coupling is turned on. Second, in the idle frame (in which we perform the final measurements) the qubits are precessing relative to each other while on resonance. Hence, a final phase correction ϕ_1,i + ω_it (where t is the analogue evolution time) must also be applied to every qubit before measurements. These corrections are very sensitive to timing and dispersive shifts: before the analogue evolution, a timing delay in dimer generation of only 150 ps corresponds to a 0.1-rad change in ϕ₀ for an idle frequency difference of 100 MHz. Furthermore, during the idle evolution, a 0.1% (80 kHz) change in dispersive shift leads to a 0.1-rad change in the final phase after 200 ns of analogue evolution. Hence, standard calibration techniques, such as single-qubit Ramsey spectroscopy, in which the configuration is sufficiently different from that in the actual experiment, are not accurate enough. We therefore use a set of three calibration techniques for ϕ_0,i, ϕ_1,i and ω_i that are designed to represent the configuration used in the actual experiment as well as possible:

To calibrate ϕ_0,i, we make use of the fact that the dimer state is only an eigenstate of the coupling Hamiltonian when the phase difference of the qubits is 0 or π. Hence, we sweep the phase difference and measure the population oscillations between the qubits with time. The correct phase compensation is the one that minimizes the amplitude of the population oscillations. We note two important points about this calibration step: first, as the measurements are in the Z-basis, they do not depend on the calibration of ϕ_1,i and ω_i. Second, because the phase calibrated in this step is accumulated before the couplers are turned on, it is not affected by dispersive shifts. It is therefore not a problem that neighbouring couplers are turned off during this particular step.

As mentioned previously, the calibration of ω_i is very sensitive to dispersive shifts and must therefore be performed in the exact same configuration as the actual experiment. We achieve this by performing the KZ experiment (ramp from Neel state in staggered field) with a slow ramp and leaving the analogue Hamiltonian on for a variable time (Extended Data Fig. 5d). The resultant state shows long-range XX + YY correlations, and the effect of the phase accumulation in the idle frame is to cause oscillations in the correlator between each pair i and j with a frequency ω_i − ω_j (Extended Data Fig. 5e). Hence, by measuring the frequency of oscillations of all the correlators, the full set of {ω_i} can be determined. The key advantage of this calibration measurement is that all the couplers are turned on, so that the dispersive shifts are the same as in the actual dimer experiment. However, the initial part of the KZ circuit—including the initial staggered field and the slow ramp of the couplers—is different, so the time-independent part of the phase correction, ϕ_1,i, must be calibrated separately.

Finally, to determine ϕ_1,i, we take advantage of energy conservation. Specifically, we perform the dimer experiment with single dimers while sweeping their final phase difference (Extended Data Fig. 5f). Only the correct phase compensation leads to ⟨X₁X₂⟩ = 1 and conserved energy, as can be see in Extended Data Fig. 5g. Whereas the dispersive shifts from neighbouring couplers affect the time-dependent part of the final phase ω_it and thus had to be included in the previous step, they do not have this effect on ϕ_1,i and can therefore be excluded here.

Finally, we note that for experiments not involving entangled initial states (Figs. 3 and 4), only the step for calibration of {ω_i} outlined above is required.

Readout correction and postselection schemes

Bell measurements

When measuring ⟨XX + YY⟩ correlators using standard single-qubit measurements, we cannot simultaneously get information about the number of photons measured on the pair of qubits, preventing us from postselecting our data on photon conservation. To get around this for nearest-neighbour pairs, we change our measurement basis by applying an entangling gate given by the unitary,

$$\left[\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & 1/\sqrt{2} & -1/\sqrt{2} & 0\\ 0 & 1/\sqrt{2} & 1/\sqrt{2} & 0\\ 0 & 0 & 0 & 1\end{array}\right]$$

to each pair. From these measurements, we can deduce both the nearest-neighbour correlators and the number of photons present. We use this technique to process the data labelled ‘Bell’ in Fig. 3b. We find good alignment between direct measurements of the correlators and the inferred correlators from the Bell measurements.

Bell measurements with readout corrections

Typically, one can correct for readout errors by inverting the error channel. In the case in which readout errors are uncorrelated, we can simply characterize the matrix β for each qubit

$$\beta =\left[\begin{array}{cc}{p}_{(0|0)} & {p}_{(0|1)}\\ {p}_{(1|0)} & {p}_{(1|1)}\end{array}\right]$$

where p_(i∣j) is the probability of measuring a state $| i\rangle$ given that $| \,j\rangle$ was prepared⁵⁸. In the case in which readout errors are correlated for pairs, we can similarly characterize a matrix γ for each pair

$$\gamma =\left[\begin{array}{cccc}{p}_{(00|00)} & {p}_{(00|01)} & {p}_{(00|10)} & {p}_{(00|11)}\\ {p}_{(01|00)} & {p}_{(01|01)} & {p}_{(01|10)} & {p}_{(01|11)}\\ {p}_{(10|00)} & {p}_{(10|01)} & {p}_{(10|10)} & {p}_{(10|11)}\\ {p}_{(11|00)} & {p}_{(11|01)} & {p}_{(11|01)} & {p}_{(11|11)}\end{array}\right]$$

where p_(ij∣ab) is the probability of measuring a state $| ij\rangle$ given that $| ab\rangle$ was prepared. One can compensate for the effects of readout errors on an observable by inverting these matrices and applying them to the measured distribution of bitstrings of the subsystem involved in the observable.

In a case in which we want to both correct for readout errors and postselect our data, we cannot apply the readout correction on the postselected distributions as this would overcorrect for p_(0∣1) type errors. We also cannot simply correct the distributions of subsystem bitstrings before the postselection process because we need access to the global bitstrings to postselect on photon number conservation. Instead, we use a Markov-like process in which we consider each individual bitstring, and flip pairs of spins according to the probabilities inferred from the γ matrices. We then postselect the individual bitstrings on the criteria of photon conservation and, finally, compute the quantity of interest.

To confirm the validity of this method, we classically simulate a low-temperature state of the XY model for 64 qubits (using the ground state of two disconnected sets of 32 qubits), introduce noise to the system and use the above protocol to correct for the T₁ and readout errors. In simulating the readout errors, we include a readout bias equal to that observed in experiment, namely p_(0∣1)/p_(1∣0) = 3.7. We compute the energies of the system after various correction schemes and compare to the noiseless value. The results from these simulations are shown in Extended Data Fig. 6a,b, where we evaluate the performance for a wide range of readout error and probability of photon decay, respectively. The combined technique described above is found to provide the most accurate estimate of the actual energy across a very wide parameter range, extending beyond the range relevant to our experiment (in the experiment, we have readout errors in the range 1–4% and a probability of photon decay of 3–6% for ramp times of 200–500 ns). For very high T₁ errors, we find that the error in the combined technique eventually becomes slightly higher than that of pure postselection. In the special case of very low T₁ errors, we observe an interesting effect that leads to a slight underestimate of the energy, which can be understood as follows. Whereas the stochastic compensation of readout errors perfectly re-establishes the correct distributions of subsystem bitstrings (by construction of the probabilities with which we change the two-qubit bitstrings), each individual global bitstring has a non-zero probability of having the wrong total number of photons, even in the case of zero T1 error. The lowest-energy two-qubit state, $| 10\rangle -| 01\rangle$ (converted to $| 10\rangle$ by Bell conversion) has a slightly higher chance of being postselected than other two-qubit states. The result of this is a slight underestimate of the energy, which we emphasize is very small (roughly 1%) and not relevant in the parameter range of our experiment.

Comparison of ⟨X X⟩ and ⟨Y Y⟩

The final states produced after the ramp procedures in Figs. 3 and 4 are expected to be U(1)-symmetric, and thus have equally strong XX and YY correlations. We here check this by comparing ⟨XX⟩ and ⟨YY⟩ averaged over all nearest-neighbour qubit pairs across a range of ramp times (Extended Data Fig. 7), and indeed find that the two are equal.

Diffusion model

In Fig. 5h, we fit the observed energy transport with a diffusion model, which we describe in further detail here. We define the energy density at site (i, j), e_i,j(t), as the average of the energy (⟨XX + YY⟩/2) on the bonds that include site (i, j) and model the transport using a simple discretized version of the diffusion equation:

$$\frac{{\rm{d}}{e}_{i,j}}{{\rm{d}}t}=D({e}_{i+1,j}+{e}_{i-1,j}+{e}_{i,j-1}+{e}_{i,j+1}-4{e}_{i,j})$$

(4)

where the diffusion constant, D, is the only fit parameter.

Measurements of energy density fluctuations

We use measurements of two- and four-qubit correlators to reconstruct the energy density fluctuations, ${\sigma }_{\varepsilon }={({n}_{{\rm{B}}}{g}_{{\rm{m}}})}^{-1}\sqrt{\langle {H}_{XY}^{2}\rangle -{\langle {H}_{XY}\rangle }^{2}}$, with:

$$\begin{array}{l}\frac{{H}_{XY}^{2}}{{g}_{{\rm{m}}}^{2}}={\left(\sum _{\langle i,j\rangle }({X}_{i}{X}_{j}+{Y}_{i}{Y}_{j})/2\right)}^{2}=\sum _{\langle i,j\rangle }(1-{Z}_{i}{Z}_{j})/2\\ \,\,+\,\sum _{\langle i,j\rangle }\sum _{\langle m,n\rangle }({X}_{i}{X}_{j}{X}_{m}{X}_{n}+{Y}_{i}{Y}_{j}{Y}_{m}{Y}_{n}+{X}_{i}{X}_{j}{Y}_{m}{Y}_{n}+\\ \,\,+\,{Y}_{i}{Y}_{j}{X}_{m}{X}_{n})/4+\sum _{\langle i,j\rangle ,\langle j,k\rangle }({X}_{i}{X}_{k}+{Y}_{i}{Y}_{k})/2,\end{array}$$

(5)

where ⟨i, j⟩, ⟨j, k⟩ and ⟨m, n⟩ are nearest-neighbour pairs and i, j, k, m, n are distinct (note that j is included in the last sum to count the number of length-2 paths from i to k). Almost all of these terms can be reconstructed from just three different sets of measurements, namely {X_i}, {Y_i} and {Z_i}, except the four-qubit correlators involving both X and Y. To determine these, we measure eight periodic patterns of X, Y shown in Extended Data Fig. 8a, and leverage the substantial degree of isotropy to find the remaining correlators not included in these patterns (further justification below). As shown in Extended Data Fig. 8b, the four-qubit correlators that involve both X and Y show a clear trend with the distance between the centres of mass of the two involved nearest-neighbour pairs (i, j) and (m, n), and we therefore interpolate the data obtained from these eight sets of measurements to find the remaining terms. Determining σ_ε with good relative accuracy is challenging, owing to the very small relative difference between ${\langle {H}_{XY}\rangle }^{2}$ and $\langle {H}_{XY}^{2}\rangle$. Nevertheless, we find that our technique works well, and we obtain relatively good agreement with MPS simulations (Extended Data Fig. 8c).

To further justify the use of this interpolation technique, we show the dependence of ⟨X_iX_jY_mY_n⟩ on the relative position of the centres of mass of the two involved nearest-neighbour pairs (i, j) and (m, n), showing near-isotropic distributions (Extended Data Fig. 9a). We observe a weak angular dependence with a period of π (Extended Data Fig. 9b), which becomes most pronounced when the correlation length is maximized (for example, g_mt_r = 12.3). The amplitude is only roughly ±0.01 (or roughly 5% of the signal itself) and is expected to be due to the system shape. As we are only interested in the sum of all the correlators, this small degree of isotropy has very little effect on the interpolation scheme described above. In Extended Data Fig. 9c, we compare the result of radial interpolation of ⟨XXYY⟩ at distance 5 (dashed black curve) to the actual correlators (coloured circles in main) and their average (red dashed curve in inset), and find that the difference is very small. In particular, we quantify the relative difference between the radial interpolation and the averaged actual correlators in Extended Data Fig. 9d, and find that it is on the order of a few percent, and even smaller at the long times that are most essential to our conclusions. These deviations are comparable to the statistical noise (as shown by the error bars) and do not contribute a dominant effect to the total energy fluctuations.

Correlation fitting

We here provide further details about the fitting procedures used in the main text for analysing correlations. As shown in Extended Data Fig. 10 and also in some of the curves in Fig. 3d, we observe distortions in the correlation decay at longer distances both in experiment (a) and simulation (b), which are expected to be due to the finite size of our system. Specifically, we find that the correlations drop rapidly for some ramp times and start increasing at others. If fitting up to the longest distances, these effects have a strong impact on the analysis, as can be seen from the sharp upturn in the fit error as we exceed a fit range of roughly six sites in Extended Data Fig. 10c. Informed by these findings, and the fact that the maximum distance at which such effects are still minimal is six sites, we use this as the fit-range cut-off. Note that we also observe a noise floor in the correlations around 10⁻², and we therefore do not fit data points smaller than this value.

We investigate the dependence on fit range further by plotting the r.m.s. fit errors for all ramp times and a wide range of fit-range cut-offs in Extended Data Fig. 10d,e (power-law and exponential fits, respectively). From these plots, it is again evident that the fits with distance cut-offs longer than six sites have particularly high errors (it is of course natural to see some increase in error with increasing fit range, but we are here referring to the distinct increase seen especially well in the inset of Extended Data Fig. 10d). Plotting the error ratio in Extended Data Fig. 10f, we find that all fits up to a fit range of seven sites show the same drop below one around g_mt_r = 10, and the discrepancy from KZ scaling is observed for all fit ranges (Extended Data Fig. 10g).