Prethermalization by random multipolar driving on a 78-qubit processor

Periodically driven (Floquet) systems can host far-from-equilibrium phenomena that are absent in thermal equilibrium⁸. Prominent examples include the discrete-time crystals^1,2,3, Floquet topological matter^9,10 and dynamical phase transitions^11,12. Periodic drives have also been widely used for Floquet engineering of many-body interactions^13,14,15 and mitigating environment-induced decoherence^16,17, serving as a robust and versatile approach to stabilize and control modern quantum simulators^18,19. Explorations of non-periodic driving have surged in recent years, with rich discoveries of non-equilibrium phenomena beyond the Floquet lore^{20,21,22,23,24}. For instance, quasi-periodic and structured random drives can lead to the appearance of discrete-time quasi-crystals^25,26,27,28 and time rondeau crystals²⁹, notably enriching the possible forms of temporal order in non-equilibrium settings.

Owing to the lack of energy conservation, generic time-dependent many-body systems are inherently susceptible to heating, eventually ending in a featureless infinite-temperature state^4,5, where the subsystem entanglement entropy also reaches its maximum—the Page value³⁰. This heating effect thus poses a fundamental challenge for utilizing large-scale quantum simulators and stabilizing sought-after phases, especially over long timescales. In Floquet systems, heating can be suppressed by many-body localization, induced via efficiently strong spatial disorder^31,32. In clean systems, however, heating can also be exponentially suppressed by using high-frequency drives^6,7, leading to the transient but long-lived prethermal regime before the eventual heat death^33,34,35. In contrast, stabilizing non-periodically driven systems is a notoriously difficult task, especially when the driving protocol involves temporal randomness. This typically opens deleterious energy absorption channels, which even many-body localization cannot prevent, and thus heating occurs swiftly.

Here we experimentally demonstrate the existence of a long-lived, doubly tunable prethermal regime on a quantum simulator driven by random but structured protocols, with a universal degree of tunability in the heating rate. We use a superconducting quantum processor, Chuang-tzu 2.0, as shown in Fig. 1, that involves 78 qubits and 137 couplers. Leveraging the precise control and flexibility of this device, we accurately implement stable long-term drives and perform large-scale analogue quantum simulations of a two-dimensional hard-core Bose–Hubbard system.

a, Optical micrograph of the 78-qubit superconducting processor, Chuang-tzu 2.0. The processor is designed as a 6 × 13 square lattice, comprising 78 qubits interconnected by 137 couplers that link all adjacent qubits. b, Schematic diagram of the RMD protocol, characterized by the temporal multipolar order n. The 0-RMD is constructed by randomly selecting elements from the two elementary operators \(\{{\hat{U}}_{+},{\hat{U}}_{-}\}\), and the n-RMD sequence consists of a random selection of two n-multipoles, which are obtained by anti-aligning two (n − 1)th-order operators. c, Elementary operators are generated as \({\hat{U}}_{\pm }=\exp \{-i{\hat{H}}_{\pm }T\}\), with a driving period T, where \({\hat{H}}_{\pm }\) differs in the staggered potential in the y direction. d, Experimental procedure. First, we initialize the system in a density-wave state, where even sites along the y axis are occupied (represented by red spheres), and odd sites are unoccupied (represented by grey spheres). The symbol X denotes X-gate pulses that excite qubits to their first excited states. Next, we implement the RMD protocol that destabilizes the density-wave order, and the system heats up. Finally, we use multiqubit QST to determine the entanglement entropy and perform particle number measurement to characterize the non-equilibrium time evolution.

Furthermore, we implement a family of structured random protocols, known as random multipolar driving (RMD)^22,36. As illustrated in Fig. 1b, the protocol involves two elementary evolution operators \({\hat{U}}_{+}\) and \({\hat{U}}_{-}\), generated by two Hamiltonians \({\hat{H}}_{\pm }\) that differ in the site potential along the y axis (Fig. 1b,c). Starting from an initial density-wave state, the random modulation of \({\hat{U}}_{\pm }\) destabilizes the system and hence induces heating. The heating rate can be significantly suppressed by imposing a dipolar structure into the random driving sequence, such that the elementary building blocks now read \({\hat{U}}_{+}{\hat{U}}_{-}\) and \({\hat{U}}_{-}{\hat{U}}_{+}\). Similarly, the nth multipole can be recursively constructed by anti-aligning two (n − 1)th-order operators, and in the n → ∞ limit it converges to the self-similar Thue–Morse driving^20,24,36.

We first benchmark our experiments on 8 qubits, then we gradually enlarge the system size up to 78 qubits. To quantify the heating process, we experimentally monitor the decay of the particle imbalance. In addition, we measure the evolution of the entanglement entropy and observe distinctive stages of its growth during the heating process, thus going beyond established experimental results on driven systems where typically only the evolution of local observables is accessible^35,37,38,39. We perform more than 1,000 driving cycles, and such a long timescale allows us to capture the long-lived prethermal plateau in the high-frequency regime. Moreover, we experimentally verify the crucial role of the temporal multipolar correlation in stabilizing the system: the heating rates follow a power-law dependence on the driving frequency, with a universal scaling exponent of approximately 2n + 1, in accordance with the original theoretical analyses of the heating processes active for RMD^22,36.

Then by selecting different subsystem configurations, we demonstrate a non-uniform spatial entanglement distribution and observe the crossover from the area-law to the volume-law entanglement scaling within the prethermal regime. The onset of heating further accelerates the entanglement growth, and advanced tensor-network numerical techniques, such as grouped matrix product states (GMPS) and projected entangled pair states (PEPS), struggle to keep pace with the rapid entanglement growth. Therefore, the experimentally observed prethermalization with controllable heating rates and the entire heating dynamics towards the maximally entangled infinite-temperature state^40,41 are challenging to simulate classically with current computational resources.

Our experiments are performed on a flip-chip superconducting processor, arranged in a 6 × 13 square lattice array (Fig. 1a), with 78 qubits and 137 couplers between all neighbour qubits. The qubits are labelled as Q_x,y, with x ranging from 1 to 13 and y ranging from 1 to 6 (Fig. 1d). As the anharmonicity η is designed to be much larger than the hopping strength J, with an average value of −2π × 200 MHz, the system can be described as the non-integrable two-dimensional hard-core Bose–Hubbard model (or equivalently, the two-dimensional XY model)⁴². The effective Hamiltonian reads (ħ = 1):

where ω_x,y is the on-site potential, \({\hat{n}}_{x,y}\) denotes the particle number operator at the site (x, y), \({\hat{a}}_{x,y}^{\dagger }\) and \({\hat{a}}_{x,y}\) are the creation and annihilation operators, respectively, and J_h and J_v are the horizontal and vertical hopping strength between two nearest-neighbour qubits, respectively. H.c. denotes the Hermitian conjugate of the preceding coupling terms. In our set-up J² ≪ η, the term \({\mathcal{O}}(\,{J}^{2}/\eta )\) is negligible. Device information can be found in Supplementary Information section 3. Owing to the significant progress in coherence time, scalability and controllability of superconducting quantum circuits^43,44,45, this platform has emerged as a powerful system for exploring complex quantum phenomena that require precise manipulation.

Our target elementary operators in the n-RMD protocol are generated as \({\hat{U}}_{\pm }=\exp (-i{\hat{H}}_{\pm }T)\), where T denotes the characteristic timescale (noted as the driving period below). In our experiments

with a dimensionless parameter δh characterizing the driving amplitude. Here the on-site term (h₀ denotes the site potential) exhibits a uniform potential along the x axis and a staggered arrangement along the y axis

with the hopping term isotropic in both spatial directions, which reads

The parameters are chosen as J/2π = 2 MHz, δh = 1.2 and h₀/2π = 10 MHz; and T ranges from 3 ns to 8 ns. The average relaxation time (T₁) of our device is 26.4 μs (Supplementary Information section 3B), allowing us to experimentally implement more than 1,000 driving cycles, before any notable decoherence occurs.

For 1-RMD, the protocol involves a random sequence of two dipolar operators \({\hat{U}}_{+}{\hat{U}}_{-}\) and \({\hat{U}}_{-}{\hat{U}}_{+}\); and for 2-RMD, the elementary building blocks become \({\hat{U}}_{+}{\hat{U}}_{-}{\hat{U}}_{-}{\hat{U}}_{+}\) and \({\hat{U}}_{-}{\hat{U}}_{+}{\hat{U}}_{+}{\hat{U}}_{-}\). Experimentally, implementing the driving protocol requires precise temporal modulation of the pulse signal, especially in the high-frequency regime (T ≈ 3 ns): qubits with an odd y index are biased to the working point, whereas the pulse signal following an RMD sequence is implemented on qubits with an even y index. Through careful calibration of the Z-pulse distortion and crosstalk, combined with precise timing alignment via Floquet engineering, we achieve highly accurate RMD control (Supplementary Information section 4 and Methods).

As shown in Fig. 1d, we initialize the system as a density-wave ordered product state \(| {\psi }_{0}\rangle \), where all lattice sites with an even y index are occupied, resulting in the maximal Hilbert space dimension. In the high-frequency limit, T → 0, the early time evolution of the system can be described by an effective Hamiltonian, \({\hat{H}}_{\mathrm{eff}}=({\hat{H}}_{+}+{\hat{H}}_{-})/2\). After a short transient period, the system evolves to the prethermal state captured by the Gibbs ensemble \(\exp (-\beta {\hat{H}}_{\mathrm{eff}})\) (ref. ⁴⁶), with the inverse temperature β determined by the initial-state energy. The expectation value of a given local operator can also be determined accordingly. In our current setting, the staggered potential along the y axis (equation (3)), stabilizes the initial density-wave order during the prethermal regime. However, for any finite T, switching between H₊ and H₋ induces heating and destabilizes the prethermal state. As the system heats up to the infinite-temperature state with β(t) gradually dropping to zero, the expectation value of the local observable (and the entanglement entropy as defined later) generally follows its expected thermal value, with respect to the effective Hamiltonian