Abstract
Understanding how interacting particles approach thermal equilibrium is a major challenge of quantum simulators1,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 predictions3 in the XY model, as well as signatures of the classical Kosterlitz–Thouless phase transition4. Moreover, the digital gates enable precise energy control, allowing us to study the effects of the eigenstate thermalization hypothesis5,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.
Similar content being viewed by others
Main
The advent of quantum simulators in various platforms8,9,10,11,12,13,14 has opened a powerful experimental avenue towards answering the theoretical question of thermalization5,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 point15,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 deviations19,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 tr (refs. 3,23,24,25).
Whereas tremendous technical advancements in quantum simulators have enabled the observation of a wealth of thermalization-related phenomena26,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 control36. 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–digital37,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)39, we demonstrate analogue performance that exceeds the simulation capacity of known classical algorithms at the full system size.
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 (UA), 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, Hs.
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 Hs, 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 studies30,31 either suffered from large errors or resorted to multi-parameter optimization protocols that are difficult to scale up. Sophisticated Hamiltonian learning techniques40,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, Hd({ωqi}, {ωcj}), into a spin Hamiltonian, Hs:
where ωi and ∣gij∣ ≈ 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 ni = 0, 1 and Xi, Yi 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 Hs contains not only nearest-neighbour hopping, but also density–density interactions and next-nearest-neighbour terms, which scale as order \({\mathcal{O}}\)(g2/η) 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 Hs 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 gij/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 pmeas(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 pmeas(x, t) for reduced system sizes up to Nq = 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 state39,42. By contrast, past studies have found substantial deviations from the Porter–Thomas distribution in other models of analogue dynamics43,44.
a, Schematic representation of the experiment: Nq qubits are initialized in a half-filling state, evolved under a Hamiltonian Hs 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 Nq = 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 regime45,48; diamonds, Zuchongzhi processor46,47; circle, neutral atom analogue simulator44; green pentagon, present experiment. \({N}_{\,{\rm{q}}}^{{\rm{eff}}}\) is equal to the actual Nq in the digital experiments, whereas analogue platforms are subject to U(1) conservation (this work) or constraints from Rydberg blockade44. Inset, \({{\mathcal{E}}}_{{\rm{P}}}\) as a function of Nq 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-XEB39, \(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 tPT ≈ 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−1 with self-XEB = 0.
To also characterize the coherent errors from imperfect calibration of Hs, 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 psim are exactly simulated probabilities39. The results, shown in Fig. 2c, show exponential decay after times roughly tPT, 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, Nq = 35 (inset of Fig. 2c). This indicates the scalability of our calibration protocol and allows extrapolation to the full system size of Nq = 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}}}=\) SRényi-1/2ent\(+{\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 entropy44. Figure 2d compares the estimated \({{\mathcal{E}}}_{{\rm{P}}}\) of our work and other recent state-of-the-art experiments44,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 106 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 temperature4. 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 gm/(2π) = 20 MHz over a duration tr (Fig. 3a). Under such a protocol49, 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 tc ≈ 0.45tr when hc/gc ≈ 1.8(6) (Supplementary Information). The transition, analogous to the 2D Mott insulator–superfluid transition50, 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.
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 − gc∣ν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 tr 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 (⟨XiXj + YiYj⟩ − ⟨Xi⟩⟨Xj⟩ − ⟨Yi⟩⟨Yj⟩)/2 over all pairs i, j separated by r) measured at various tr. 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 (Nreps = 5 × 104 repetitions). e, Ratio between r.m.s. errors from power-law and exponential fits (\({{\epsilon }}_{{\rm{pow}}}/{{\epsilon }}_{\exp }\)) decreases for gmtr > 15. f, Power-law exponent, γ, approaches expected value at Kosterlitz–Thouless transition (1/4; black line). g, Vortex density proxy, nV, decreases to minimum of 2 × 10−2. h, Correlation length increases with tr. 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 nB = 110 bonds (Nq = 65) and corrected for readout errors (Methods). As tr 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 gmtr = 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 transition51. To further substantiate our interpretation, we also measure four-qubit correlators to construct the Swendsen proxy for the vortex density52, 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, ‥, NP with vertices {i1, i2, i3, i4}. Indeed, we find a rapid decrease in nV as tr is increased (Fig. 3g), to a minimum value of 2 × 10−2 in the low-energy regime.
Having studied the classical critical behaviour, we next explore the scaling of the correlation length with the duration tr over which we sweep through the quantum critical point (Fig. 3h). The correlation length rises to a maximum of ξ ≈ 10 at gmtr = 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 point19,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 tKZ, when the inverse gap exceeds ∣t − tc∣. 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 tKZ in Fig. 3h and find better agreement with the KZ prediction. Notably, by rescaling to ξ/ξKZ and (t − tc)/∣tKZ − tc∣ ≡ τ/τKZ, we find that the curves collapse to a common f(τ/τKZ), consistent with predictions of universal coarsening dynamics19,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, n0 (pairs of spin flips in randomized locations) in the initial state53. Whereas we find that the final average energy density depends linearly on n0 (Fig. 4a), the behaviour of the correlation length is more intricate (Fig. 4b) and is best understood by plotting ξ versus energy density for all n0 and gmtr > 5 (Fig. 4c). Notably, the points show a collapse (also for nV in inset), suggesting that the final states are well thermalized, such that the energy density determines ξ and nV, 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.
a,b, Energy density (a) and correlation length (b) versus number of initial excitations, n0, averaged over three randomizations. Error bars smaller than markers (Nreps = 6 × 103 per data point). c, Energy dependence of ξ, demonstrating collapse when data from sweeps of tr and n0 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, tr-dependence of energy density and its fluctuations (width), for various n0. 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 n0 at tr = 200 ns ≈ 23/g. Increasing n0 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 n0 and tr to have the same observables (for example, nV 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 HXY = ∑⟨i,j⟩(XiXj + YiYj)/2, which trivially commute with HXY 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 n0 = 0, whereas its dependence on n0 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 tr and n0 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 nV 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 measurements54. At n0 = 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 model55. However, tuning to higher final energies by means of larger n0, 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 density2,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 tr and n0 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 Hs 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 tr and n0.
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 (Nreps = 104).
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, ⟨XiYi+1 − YiXi+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 atoms56.
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. 45. Extended Data Fig. 1a,b show the measured Ramsey dephasing (\({T}_{2}^{* }\)) and photon relaxation (T1) 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−3 for combined \(\sqrt{\text{iSWAP}\,}\) and single-qubit gates (Extended Data Fig. 1c), and 1.0 × 10−3 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−2.
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:
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 kqq, 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, q1 and q2 separated by a coupler c12: \({\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, q1 and q3: \({\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, q2 and q4: \({\widetilde{k}}_{{q}_{2},{q}_{4}}={k}_{{q}_{2},{q}_{4}}\).
For kqc, 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, q1 and c23: \({\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, q4 and c12: \({\widetilde{k}}_{{q}_{4},{c}_{12}}={k}_{{q}_{4},{c}_{12}}\).
For kcc, we consider two types of coupler–coupler coupling:
-
(1)
Diagonally separated couplers in the NW–SE direction c12 and c23: \({\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, c12 and c14: \({\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, Hs, that acts on the computational subspace. To find spin Hamiltonian terms involving n photons in a system of Nq 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 Hd. A second option would be to instead use the Nn 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 Nn lowest n-photon eigen-energies as Hd. 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 eigenvalues57. 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(n) 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.1d). 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:
where \({g}_{ij}^{nn},{g}_{ijk}^{XnX}\) and \({g}_{ijk}^{XIX}\) scale as g2/η, while \({g}_{ijk}^{nXX}\) scales as g3/η2 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 Nn lowest-energy n-photon eigenstates of Hd, which has a high computational cost for large Nq. 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 ϕ0 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 ⟨X1X2⟩ = 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,
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
where p(i∣j) is the probability of measuring a state \(| i\rangle \) given that \(| \,j\rangle \) was prepared58. In the case in which readout errors are correlated for pairs, we can similarly characterize a matrix γ for each pair
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 T1 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 T1 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 T1 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), ei,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:
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:
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 {Xi}, {Yi} and {Zi}, 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 ⟨XiXjYmYn⟩ 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, gmtr = 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−2, 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 gmtr = 10, and the discrepancy from KZ scaling is observed for all fit ranges (Extended Data Fig. 10g).
Data availability
The data that support the findings in this study are available at Zenodo (https://doi.org/10.5281/zenodo.14060446)59.
References
Altman, E. et al. Quantum simulators: architectures and opportunities. PRX Quantum 2, 017003 (2021).
Abanin, D. A., Altman, E., Bloch, I. & Serbyn, M. Colloquium: many-body localization, thermalization, and entanglement. Rev. Mod. Phys. 91, 021001 (2019).
del Campo, A. & Zurek, W. H. Universality of phase transition dynamics: topological defects from symmetry breaking. Int. J. Mod. Phys. A 29, 1430018 (2014).
Kosterlitz, J. M. & Thouless, D. J. Ordering, metastability and phase transitions in two-dimensional systems. J. Phys. C. Solid State Phys. 6, 1181 (1973).
Deutsch, J. M. Quantum statistical mechanics in a closed system. Phys. Rev. A 43, 2046–2049 (1991).
Srednicki, M. Chaos and quantum thermalization. Phys. Rev. E 50, 888–901 (1994).
D’Alessio, L., Kafri, Y., Polkovnikov, A. & Rigol, M. From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics. Adv. Phys. 65, 239–362 (2016).
Daley, A. J. et al. Practical quantum advantage in quantum simulation. Nature 607, 667–676 (2022).
Bloch, I., Dalibard, J. & Zwerger, W. Many-body physics with ultracold gases. Rev. Mod. Phys. 80, 885–964 (2008).
Blatt, R. & Roos, C. F. Quantum simulations with trapped ions. Nat. Phys. 8, 277–284 (2012).
Houck, A. A., Türeci, H. E. & Koch, J. On-chip quantum simulation with superconducting circuits. Nat. Phys. 8, 292–299 (2012).
Carusotto, I. et al. Photonic materials in circuit quantum electrodynamics. Nat. Phys. 16, 268–279 (2020).
Browaeys, A. & Lahaye, T. Many-body physics with individually controlled Rydberg atoms. Nat. Phys. 16, 132–142 (2020).
King, A. D. et al. Computational supremacy in quantum simulation. Preprint at https://arxiv.org/abs/2403.00910 (2024).
Dziarmaga, J. Dynamics of a quantum phase transition: exact solution of the quantum Ising model. Phys. Rev. Lett. 95, 245701 (2005).
Zurek, W. H., Dorner, U. & Zoller, P. Dynamics of a quantum phase transition. Phys. Rev. Lett. 95, 105701 (2005).
Polkovnikov, A. Universal adiabatic dynamics in the vicinity of a quantum critical point. Phys. Rev. B 72, 161201 (2005).
Ali, A. et al. Quantum quench dynamics of geometrically frustrated Ising models. Preprint at https://arxiv.org/abs/2403.00091 (2024).
Samajdar, R. & Huse, D. A. Quantum and classical coarsening and their interplay with the Kibble–Zurek mechanism. Preprint at https://arxiv.org/abs/2401.15144 (2024).
Roychowdhury, K., Moessner, R. & Das, A. Dynamics and correlations at a quantum phase transition beyond Kibble–Zurek. Phys. Rev. B 104, 014406 (2021).
Biroli, G., Cugliandolo, L. F. & Sicilia, A. Kibble–Zurek mechanism and infinitely slow annealing through critical points. Phys. Rev. E 81, 050101 (2010).
Chandran, A., Erez, A., Gubser, S. S. & Sondhi, S. L. Kibble–Zurek problem: universality and the scaling limit. Phys. Rev. B 86, 064304 (2012).
Navon, N., Gaunt, A. L., Smith, R. P. & Hadzibabic, Z. Critical dynamics of spontaneous symmetry breaking in a homogeneous bose gas. Science 347, 167–170 (2015).
Keesling, A. et al. Quantum Kibble–Zurek mechanism and critical dynamics on a programmable Rydberg simulator. Nature 568, 207–211 (2019).
Ebadi, S. et al. Quantum phases of matter on a 256-atom programmable quantum simulator. Nature 595, 227–232 (2021).
Langen, T. et al. Experimental observation of a generalized Gibbs ensemble. Science 348, 207–211 (2015).
Prüfer, M. et al. Observation of universal dynamics in a spinor Bose gas far from equilibrium. Nature 563, 217–220 (2018).
Schreiber, M. et al. Observation of many-body localization of interacting fermions in a quasirandom optical lattice. Science 349, 842–845 (2015).
Kaufman, A. M. et al. Quantum thermalization through entanglement in an isolated many-body system. Science 353, 794–800 (2016).
Neill, C. et al. Ergodic dynamics and thermalization in an isolated quantum system. Nat. Phys. 12, 1037–1041 (2016).
Roushan, P. et al. Spectroscopic signatures of localization with interacting photons in superconducting qubits. Science 358, 1175–1179 (2017).
Braumüller, J. et al. Probing quantum information propagation with out-of-time-ordered correlators. Nat. Phys. 18, 172–178 (2022).
Zhou, Z.-Y. et al. Thermalization dynamics of a gauge theory on a quantum simulator. Science 377, 311–314 (2022).
Zhang, X., Kim, E., Mark, D. K., Choi, S. & Painter, O. A superconducting quantum simulator based on a photonic-bandgap metamaterial. Science 379, 278–283 (2023).
Scholl, P. et al. Quantum simulation of 2D antiferromagnets with hundreds of Rydberg atoms. Nature 595, 233–238 (2021).
Karamlou, A. H. et al. Probing entanglement in a 2D hard-core Bose–Hubbard lattice. Nature 629, 561–566 (2024).
Bluvstein, D. et al. A quantum processor based on coherent transport of entangled atom arrays. Nature 604, 451–456 (2022).
Lamata, L., Parra-Rodriguez, A., Sanz, M. & Solano, E. Digital-analog quantum simulations with superconducting circuits. Adv. Phys. X 3, 1457981 (2018).
Boixo, S. et al. Characterizing quantum supremacy in near-term devices. Nat. Phys. 14, 595–600 (2018).
Bairey, E., Arad, I. & Lindner, N. H. Learning a local hamiltonian from local measurements. Phys. Rev. Lett. 122, 020504 (2019).
Evans, T. J., Harper, R. & Flammia, S. T. Scalable Bayesian Hamiltonian learning. Preprint at https://arxiv.org/abs/1912.07636 (2019)
Shaw, A. L. et al. Universal fluctuations and noise learning from Hilbert-space ergodicity. Preprint at https://arxiv.org/abs/2403.11971 (2024).
Mark, D. K., Choi, J., Shaw, A. L., Endres, M. & Choi, S. Benchmarking quantum simulators using ergodic quantum dynamics. Phys. Rev. Lett. 131, 110601 (2023).
Shaw, A. L. et al. Benchmarking highly entangled states on a 60-atom analogue quantum simulator. Nature 628, 71–77 (2024).
Arute, F. et al. Quantum supremacy using a programmable superconducting processor. Nature 574, 505–510 (2019).
Wu, Y. et al. Strong quantum computational advantage using a superconducting quantum processor. Phys. Rev. Lett. 127, 180501 (2021).
Zhu, Q. et al. Quantum computational advantage via 60-qubit 24-cycle random circuit sampling. Sci. Bull. 67, 240–245 (2022).
Morvan, A. et al. Phase transitions in random circuit sampling. Nature 634, 328–333 (2024).
Chen, C. et al. Continuous symmetry breaking in a two-dimensional Rydberg array. Nature 616, 691–695 (2023).
Fisher, M. P. A., Weichman, P. B., Grinstein, G. & Fisher, D. S. Boson localization and the superfluid-insulator transition. Phys. Rev. B 40, 546–570 (1989).
Nelson, D. R. & Kosterlitz, J. M. Universal jump in the superfluid density of two-dimensional superfluids. Phys. Rev. Lett. 39, 1201–1205 (1977).
Swendsen, R. H. First-and second-order phase transitions in the d = 2xy model. Phys. Rev. Lett. 49, 1302 (1982).
Schuckert, A. et al. Observation of a finite-energy phase transition in a one-dimensional quantum simulator. Preprint at https://arxiv.org/abs/2310.19869 (2023).
Brydges, T. et al. Probing rényi entanglement entropy via randomized measurements. Science 364, 260–263 (2019).
Metlitski, M. A. & Grover, T. Entanglement entropy of systems with spontaneously broken continuous symmetry. Preprint at https://arxiv.org/abs/1112.5166 (2015).
Manovitz, T. et al. Quantum coarsening and collective dynamics on a programmable quantum simulator. Preprint at https://arxiv.org/abs/2407.03249 (2024).
Bravyi, S., DiVincenzo, D. P. & Loss, D. Schrieffer–Wolff transformation for quantum many-body systems. Ann. Phys. 326, 2793–2826 (2011).
Smith, A. W. R., Khosla, K. E., Self, C. N. & Kim, M. S. Qubit readout error mitigation with bit-flip averaging. Sci. Adv. 7, eabi8009 (2021).
Andersen, T. I. Thermalization and criticality on an analog-digital quantum simulator v.1. Zenodo https://doi.org/10.5281/zenodo.14060446 (2024).
Acknowledgements
We acknowledge useful discussions with R. Samajdar, D. A. Huse and S. Choi. A. Schuckert acknowledges support from the US Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator. J.M. acknowledges funding through SNSF Swiss Postdoctoral Fellowship, grant no. 210478. A.E. acknowledges funding by the German National Academy of Sciences Leopoldina under the grant number LPDS 2021-02 and by the Walter Burke Institute for Theoretical Physics at Caltech. Work in Grenoble is funded by the French National Research Agency through the JCJC project QRand (grant no. ANR-20-CE47-0005), Laboratoire d’excellence LANEF (grant no. ANR-10-LABX-51-01), from the Grenoble Nanoscience Foundation.
Author information
Authors and Affiliations
Contributions
T.I.A., D.A.A. and X.M. conceived the project and designed the experiments. T.I.A., A.H.K. and J. Berndtsson performed the experiments and data analysis (A.H.K., entanglement entropy measurements and J. Berndtsson, measurements of ramp time dependence). T.I.A., J.A.G., X.M. and D.A.A. developed the calibration procedures, with assistance from N.A., Y.Z., E.F., B.K., A D.P., A.R.K., I.D., A. Petukhov and L.B.I. J.M., A. Szasz, D.R. and D.A.A. performed MPS simulations and theoretical work with A.M.L. N.A. performed XEB analysis, classical complexity estimates and exact state vector simulations with assistance from T. Westerhout, V.K. and A. Szasz. A. Elben, A.R., V.V. and B. Vermersch performed theoretical work on randomized measurements. T.I.A., N.A., A.M.L., D.A.A. and X.M. wrote the manuscript. D.A.A. and X.M. led and coordinated the project. T.I.A., N.A., A.H.K., J. Berndtsson, J.M., A. Szasz, J.A.G., A. Schuckert, T. Westerhout, Y.Z., E.F., D.R., B.K., A.D.P., A.R.K., I.D., V.K., A. Petukhov, L.B.I., A. Elben, A.R., V.V., B. Vermersch, R.A., L.A.B., K. Anderson, M. Ansmann, F.A., K. Arya, A.A., J.A., B. Ballard, J.C.B., A. Bengtsson, A. Bilmes, G.B., A. Bourassa, J. Bovaird, L.B., M.B., D. A. Browne, B. Buchea, B.B.B., D. A. Buell, T.B., B. Burkett, N.B., A. Cabrera, J. Campero, H.-S.C., Z.C., B.C., J. Claes, A.Y.C., J. Cogan, R.C., P.C., W.C., A.L.C., S. Das, D.M.D., L.D.L., A.D.T.B., S. Demura, P.D., A.D., C. Earle, A. Eickbusch, A.M.E., M.E., C. Erickson, L.F., R.F., V.S.F., L.F.B., A.G.F., B.F., S.G., R. Gasca, W.G., C.G., D. Gilboa, M.G., R. Gosula, A.G.D., D. Graumann, A.G., S. Habegger, M.C.H., M.H., M.P.H., S.D.H., S. Heslin, P. Heu, G.H., M.R.H., H.-Y. Huang, T.H., A.H., W.J.H., S.V.I., E.J., Z.J., C. Jones, S.J., C. Joshi, P.J., D.K., H.K., K.K., T. Khaire, T. Khattar, M. Khezri, M. Kieferová, S.K., A.K., P.K. A.N.K., F.K., J.M.K., D. Landhuis, B.W.L., P.L., K.-M.L., L.L.G., J. Ledford, J. Lee, K. W. Lee, Y.D.L., B.J.L., W.Y.L., A.T.L., W.L., W.P.L., A. Locharla, D. Lundahl, A. Lunt, S. Madhuk, A. Maloney, S. Mandrà, L.S.M., O.M., S. Martin, C.M., J.R.M., M.M., S. Meeks, K.C.M., A. Mieszala, S. Molina, S. Montazeri, A. Morvan, R.M., C.N., A. Nersisyan, M. Newman, A. Nguyen, M. Nguyen, C.-H.N., M.Y.N., W.D.O., K.O., A. Pizzuto, R.P., O.P., L.P.P., C.Q., M.J.R., D.M.R., G.R., C.R., E.R., N.C.R., N. Saei, K.S., K.J.S., H.F.S., C.S., M.J.S., A. Shorter, N. Shutty, V. Shvarts, V. Sivak, J. Skruzny, S. Small, W.C.S., S. Springer, G.S., J. Suchard, M.S., A. Sztein, D.T., A.T., M.M.T., A.V., S.V., B. Villalonga, C.V.H., S.W., S.X.W., T. White, K.W., B.W.K.W., C.X., Z.J.Y., P.Y., B.Y., J.Y., N.Y., G.Y., A.Z., N. Zhu, N. Zobrist, H.N., R.B., S.B., J.H., E.L., A. Megrant, J.K., Y.C., V. Smelyanskiy, G.V., P.R., A.M.L., D.A.A. and X.M. contributed to revising the paper and the Supplementary Information.
Corresponding authors
Ethics declarations
Competing interests
A provisional patent application has been submitted for the analogue calibration scheme, titled ‘High-accuracy calibration of an analog quantum simulator’ (application number 63/639,509). The listed inventors are T.I.A., J.A.G., D.A.A. and X.M. The other authors declare no competing interests.
Peer review
Peer review information
Nature thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Extended data figures and tables
Extended Data Fig. 1 Device characterization.
a,b, Ramsey dephasing (\({T}_{2}^{* }\); a) and photon relaxation (T1; b) times across the qubit grid. c,d, Histogram of Pauli error for combined \(\sqrt{{\rm{iSWAP}}}\) and single qubit gates (c) and only single qubit gates (d). Red dashed lines indicate the median values. (CDF: cumulative distribution function). e, Histogram of readout errors.
Extended Data Fig. 2 Schematic of underlying coupling pathways in the device.
In addition to capacitive coupling between neighboring qubits (orange) and couplers (blue), there are also diagonal next-nearest-neighbor couplings. Asymmetry in the underlying structure of the qubits causes a difference in the couplings along the NW-SE and NE-SW diagonals.
Extended Data Fig. 3 Analogue calibration procedure.
a, Overview of calibration steps. We perform three main steps, which together allow for determining the bare frequencies of the qubits and couplers in the idle configuration in which \(\widetilde{g}=0\) (top row), as well as in the interaction configuration (bottom two rows). For each step, we model a subsystem (third column) to convert the measured dressed frequencies (fourth column) to bare frequencies (fifth column). b, Circuit schematic of swap spectroscopy. c, Top: Measured population difference, ⟨Z1 − Z2⟩, as a function of qubit detuning and time. Bottom: Extracted swap rate from Fourier transform vs qubit detuning. The position of the minimum allows for determining \(\widetilde{g}\) and the difference of the dressed qubit frequencies, \({\widetilde{\omega }}_{q1}-{\widetilde{\omega }}_{q2}\). d, Circuit schematic of single-photon spectroscopy. e, Fourier transform of the measured ⟨X⟩ + i⟨Y⟩. The average of the peak positions is equal to the average of the dressed qubit frequencies \(({\widetilde{\omega }}_{q1}+{\widetilde{\omega }}_{q2})/2\).
Extended Data Fig. 4 Higher order terms in the analogue spin Hamiltonian.
Average coupling coefficient vs nearest-neighbor hopping g for a, nini+1, b, (Xini+1Xi+2 + Yini+1Yi+2)/2, c, (XiXi+2 + YiYi+2)/2, and d, ni(Xi+1Xi+2 + Yi+1Yi+2)/2, where qubits i, i + 1, and i + 2 are placed along a connected line. Aside from an offset due to diagonal capacitive coupling, the first three terms scale as g2/η, while the fourth scales as g3η2, where η is the anharmonicity. At g = 2π × 10MHz, all higher-order terms are smaller than 1 × 2π MHz. In the three latter terms, there is asymmetry between the three possible configurations displayed in the insets (see text for details). Note that ni(Xi+1Xi+2 + Yi+1Yi+2)/2 does not differ on average from (XiXi+1 + YiYi+1)ni+2/2 in d.
Extended Data Fig. 5 Phase calibration for hybrid analogue-digital experiments.
a, Schematic of phase accumulation and correction throughout hybrid analogue-digital circuit. While we typically prepare initial dimer states, we here consider an initial state \(| ++\rangle \) for the purpose of simplified explanation. Blue and yellow lines show qubit frequency trajectories and coupling profile, respectively, while brown (beige) boxes show the relative alignment of the two spins in the idle (resonance) frame. We apply corrective phases {ϕ0,i} before the analogue circuit to ensure the correct dimer phase in the resonance frame when the analogue Hamiltonian is turned on. Additional phases {ωit + ϕ1,i} are applied after the analogue evolution in order to measure the same phase in the idle frame as was in the resonance frame. b, {ϕ0,i} are calibrated by preparing triplet states, sweeping the phase difference within each qubit pair, and measuring the population difference after a variable time t. c, Population difference after time t for an applied phase difference ϕ. Since only the dimer phases 0 and π are eigenstates of the analogue Hamiltonian, the correct ϕ0,i is determined by minimizing the population oscillations. d, {ωi} are calibrated by performing adiabatic ground state preparation with an initial staggered field and a slow (25/gm) ramp, and measuring the ⟨XX⟩ correlations a time t after the ramp. e, Top: ⟨XX⟩ after time t when applying no corrective phase after the analogue evolution. Since the low-energy final state is known to have long-range correlations, the observed oscillations can be fit to extract the time-dependent part of the corrective phase after the analogue pulse. Bottom: ⟨XX⟩ after time t when applying the corrective phase found from fitting the oscillations. The near-constant value indicates a successful correction. f, {ϕ1,i} are calibrated by preparing an initial dimer state, performing the same circuit as in the experiment with corrective pre-analogue phases {ϕ0,i} and partial post-analogue phases {ωit}, applying a variable phase ϕ to one qubit in each pair, and measuring the ⟨XX⟩ correlations a time t after the ramp. g, Top: ⟨XX⟩ after time t. Since the state is known to be the triplet state, the correct ϕ1 is found from maximizing ⟨XX⟩ correlations. Bottom: As a complementary technique, one can prepare the singlet state instead and find the ϕ that minimizes variations in ⟨XX⟩ correlations.
Extended Data Fig. 6 Correction for readout error and photon decay.
Performance of various readout and photon decay correction techniques as a function of a, readout error b, and photon decay probability. The performance is measured as the relative error between the estimated energy (Eest) and the actual ground state energy (Egs). We find that the combined technique (red) achieves the lowest relative error for a very wide range of parameters.
Extended Data Fig. 7 Comparison of XX- and YY-correlations.
Ramp time dependence of ⟨XX⟩ and ⟨YY⟩ averaged over all nearest-neighbor pairs. The two are found to be very similar, consistent with U(1)-symmetry.
Extended Data Fig. 8 Energy density fluctuations.
a, In addition to {Xi}, {Yi} and {Zi}, we measure 8 periodic patterns of XX and YY to find σε. b, ⟨XXYY⟩ has a relatively simple dependence on Euclidean distance (data from measurements shown in a), which can be interpolated to find the remaining terms. c, Energy density fluctuations, σε, displaying good agreement between experiment (red) and simulation (blue); however, at long ramp times, decoherence causes higher fluctuations in the experimental case.
Extended Data Fig. 9 Isotropy of 4-qubit correlators.
a, Dependence of ⟨XiXjYmYn⟩ on relative position between centers of mass of sites (i, j) and (m, n), showing a substantial degree of isotropy. b, Angular dependence of ⟨XiXjYmYn⟩, displaying weak π-periodic oscillations that are most pronounced in the regime with longest correlations. c, Comparison of radial interpolation (dashed black) with the actual correlators at distance 5 (colored circles in main) and their average (dashed red in inset). d, Relative difference between radial interpolation and average of actual correlators, as a function of ramp time. The error is on the order of a few percent, and is comparable to the statistical error (error bars estimated from bootstrapping with 48,000 shots at each ramp time).
Extended Data Fig. 10 Correlation fitting.
a,b, Dependence of ⟨XiXj + YiYj⟩/2 on Euclidean distance between i and j at various ramp times from experimental data (a) and simulation results (b). In both cases, we observe distortions at longer (≳ 6 sites) distances, attributed to the finite size of our system. c, Root-mean-square fit error for exponential (green) and power-law (violet) fits, as a function of the cut-off distance applied in the fits. Going beyond a distance of 6, we observe a steep increase in the fit error, arising from the effects seen in a and b. d,e, Ramp time dependence of rms error in power-law fits (d) and exponential fits (e) for various fit range cut-offs. Inset of d: Enlarged version of area indicated by red dashed square, showing abrupt increase in error for fit-range cut-offs longer than 7 sites. f,g, Ramp time dependence of rms error ratio (f) and correlation length (g). The drop in the fit error ratio below 1 is observed for all fit range cut-offs shorter or equal to 7 sites (f), and the discrepancy from KZ-scaling (dashed black) persists for all fit range cut-offs (g).
Supplementary information
Supplementary Information
The Supplementary Information includes Notes 1–13 and Figs. 1–16. In this file, we describe MPS simulations of XY model dynamics, numerical finite-size scaling analysis, alternative correlation fitting schemes and further theoretical analysis of XEB experiments, including computational complexity.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Andersen, T.I., Astrakhantsev, N., Karamlou, A.H. et al. Thermalization and criticality on an analogue–digital quantum simulator. Nature 638, 79–85 (2025). https://doi.org/10.1038/s41586-024-08460-3
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1038/s41586-024-08460-3