Main

Periodically driven (Floquet) systems can host far-from-equilibrium phenomena that are absent in thermal equilibrium8. Prominent examples include the discrete-time crystals1,2,3, Floquet topological matter9,10 and dynamical phase transitions11,12. Periodic drives have also been widely used for Floquet engineering of many-body interactions13,14,15 and mitigating environment-induced decoherence16,17, serving as a robust and versatile approach to stabilize and control modern quantum simulators18,19. Explorations of non-periodic driving have surged in recent years, with rich discoveries of non-equilibrium phenomena beyond the Floquet lore20,21,22,23,24. For instance, quasi-periodic and structured random drives can lead to the appearance of discrete-time quasi-crystals25,26,27,28 and time rondeau crystals29, 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 state4,5, where the subsystem entanglement entropy also reaches its maximum—the Page value30. 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 disorder31,32. In clean systems, however, heating can also be exponentially suppressed by using high-frequency drives6,7, leading to the transient but long-lived prethermal regime before the eventual heat death33,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.

Fig. 1: Quantum processor and experimental scheme.
figure 1

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.

Full size image

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 driving20,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 accessible35,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 RMD22,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 state40,41 are challenging to simulate classically with current computational resources.

Experimental set-up

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 Qx,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)42. The effective Hamiltonian reads (ħ = 1):

$$\begin{array}{l}\hat{H}\,=\,\sum _{x,y}{\omega }_{x,y}{\hat{n}}_{x,y}+\sum _{x,y}[\,{J}_{{\rm{h}}}{\hat{a}}_{x,y}^{\dagger }{\hat{a}}_{x+1,y}\\ \,+{J}_{{\rm{v}}}{\hat{a}}_{x,y}^{\dagger }{\hat{a}}_{x,y+1}+{\rm{H.c.}}+{\mathcal{O}}(\,{J}^{2}/\eta )],\end{array}$$
(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 Jh and Jv 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 J2 ≪ η, 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 circuits43,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

$${\hat{H}}_{\pm }={\hat{H}}_{\kappa }+(1\pm \delta h){\hat{H}}_{{\rm{p}}},$$
(2)

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

$${\hat{H}}_{{\rm{p}}}={h}_{0}\sum _{x,y}[1-{(-1)}^{y}]{\hat{a}}_{x,y}^{\dagger }{\hat{a}}_{x,y},$$
(3)

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

$${\hat{H}}_{\kappa }=J\sum _{x,y}({\hat{a}}_{x,y}^{\dagger }{\hat{a}}_{x+1,y}+{\hat{a}}_{x,y}^{\dagger }{\hat{a}}_{x,y+1}+{\rm{H.c.}}).$$
(4)

The parameters are chosen as J/2π = 2 MHz, δh = 1.2 and h0/2π = 10 MHz; and T ranges from 3 ns to 8 ns. The average relaxation time (T1) 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. 46), 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 Hamiltonian47. Therefore, the particle number imbalance serves as a good diagnostics of this heating process

$${\mathcal{I}}=\frac{1}{{N}_{0}}\sum _{j\in \{x,y\}}\langle {\psi }_{0}|{\hat{Z}}_{j}(t){\hat{Z}}_{j}(0)|{\psi }_{0}\rangle ,$$
(5)

with \({\hat{Z}}_{j}=2{\hat{n}}_{j}-1\), and the initial particle number N0. In fact, our experiment suffers from weak particle loss owing to the excitation to high energy levels, which can speed up the driving-induced heating. This effect can be suppressed by using an error mitigation scheme (Supplementary Information section 5C). In addition to these local observables, we also study the time evolution of the entanglement entropy, \(S=-\,\text{Tr}[\,\rho \text{ln}\rho \,]\), where ρ denotes the reduced density matrix of a subsystem, to capture the growth of non-local correlations between a subsystem and its complement. It provides valuable information for estimating the numerical complexity in simulating the many-body dynamics48,49,50. In the experiment, we perform quantum-state tomography (QST)51 on a subsystem to reconstruct ρ at different stages of the heating process. This can be achieved by performing statistically complete measurement sets, typically implemented through projective measurements along the Pauli eigenbases.

Characteristics of prethermalization by RMD

We first benchmark our experiments on eight qubits, which can be efficiently simulated numerically. Then, we gradually increase the system size (as shown in Fig. 3a) and finally use all 78 qubits, allowing for a detailed analysis of heating dynamics and the entanglement entropy growth.

Figure 2a shows the measured entanglement entropy dynamics of a single qubit under 1-RMD, in a small-scale system of 8 qubits and 10 couplers, for T = 4 ns and T = 8 ns respectively. The entropy evolution shows several distinct features. (1) Starting from S = 0, as expected for a product initial state, the entanglement entropy rapidly increases within a short time. (2) A distinct prethermal plateau appears within the time interval 20 ns ≲ t ≲ 100 ns, hallmarking the notable suppression of energy absorption. This plateau can be described by a quasi-conserved effective Hamiltonian, which can be perturbatively constructed using a generalized Floquet–Magnus expansion36. (3) At later times, S deviates from the plateau and continues increasing. Around t ≳ 2,000 ns, it saturates at the maximum entropy SM ≈ 0.7, in accordance with the Page value, SP ≈ Lln2 (ref. 30), where L denotes the subsystem size, indicating that the system completely thermalizes. We also perform exact numerical simulations (light dashed curves in Fig. 2a), which match well with our experimental observations. The lifetime of the prethermal plateau, τS, can be quantitatively determined by numerically fitting the entanglement entropy growth to the ansatz \(S\approx {S}_{{\rm{M}}}(1-{{\rm{e}}}^{-t/{\tau }_{S}})\) (ref. 47) in the intermediate region 100 ns ≲ t ≲ 2,000 ns. By reducing the driving period from T = 8 ns to T = 4 ns, τS significantly increases from 0.25 μs to 1.03 μs, indicating a strong suppression of the heating rate.

Fig. 2: Prethermalization by RMD in an eight-qubit system.
figure 2

a,b, Experimental data showing the dynamics of the subsystem entanglement entropy (a) and imbalance (b), for driving periods T = 4 ns and T = 8 ns, where the prethermal plateau is clearly probed (shown as shaded regions). Both the lifetime of the entanglement plateau, τS, and the imbalance decay timescale, τI, can be prolonged by increasing the driving frequency. We implement the 1-RMD protocol, and the system comprises 8 qubits {Qx,y} (the inset in b), indexed by x = {5, 6, 7, 8} and y = {1, 2}. We use QST to determine the von-Neumann entropy of the qubit {Q5,2}. Error bars indicate 1 s.d. of the experimental results, averaged over 10 independent RMD sequences. The light dashed curves depict numerical results for benchmarking.

Full size image

A similar feature also appears during the decay of the imbalance. As shown in Fig. 2b, we observe that, starting from the initial imbalance value \({\mathcal{I}}=1\), the imbalance decays towards 0 at later times, and a high-frequency drive significantly slows down this decay. The characteristic decay timescale τI can be extracted by fitting the imbalance evolution to an exponentially decaying function \({\mathcal{I}}\propto {{\rm{e}}}^{-t/{\tau }_{I}}\). For T = 8 ns and T = 4 ns, the timscales are 1.2 μs and 7.4 μs, respectively.

Tunable heating rate by RMD

To further study the suppression of the heating rate, we perform a similar RMD protocol and examine the dependence of the prethermal lifetime on the driving frequency for a different multipolar order n. Experimentally, we implement the driving protocol on larger systems with distinct geometries, as shown in Fig. 3a. In particular, in Fig. 3b, we present the time evolution of the imbalance for the 78-qubit system driven by the 2-RMD protocol using T = 4 ns and T = 7 ns (orange and red dots, respectively). In both cases, the imbalance decays exponentially, in agreement with both numerical47 and experimental observations during the heating process35. In Fig. 3c,d, we use a log–log scale and show the prethermal lifetimes τI and τS for different 1/T. Data points approximately follow a straight line, suggesting a power-law decay, τI,S ≈ (1/T)α. This power-law behaviour is in sharp contrast to conventional Floquet systems with local interactions, where the prethermal lifetime scales exponentially with the increase of the driving frequency τI,S ≈ e1/T (refs. 37,39). It is noted that the scaling exponent α is now tunable and shows a strong dependence on the multipolar order in the high-frequency regime. This is a significant feature of RMD systems, which cannot be achieved with either Floquet or other quasi-periodic drives. The reason is that the multipolar structure suppresses the low-frequency components of the drive, thereby constraining heating. A generalized Floquet–Magnus theory and a Fermi’s golden rule argument predicts the relation α = 2n + 1 (refs. 22,36).

Fig. 3: Scaling behaviour of the prethermal lifetime with driving frequency.
figure 3

a, Various system configurations and their corresponding power-law scaling exponents. These scaling exponents closely follow the theoretically predicted value of α(n) = 2n + 1 (black horizontal dashed lines). b, Decay of the imbalance shown in a log scale. The markers represent the experimental data, and the dashed and solid curves are obtained using the GMPS with bond dimension χ = 96 and PEPS, for T = 7 ns and T = 4 ns, respectively, with 2-RMD drives. Numerical methods can capture only the early time evolution. c,d, The power-law scaling of the prethermal lifetime obtained from the dynamics of the imbalance τI (c) and the subsystem entanglement entropy τS (d) in the 78-qubit system. A larger multipolar order n significantly suppresses the heating rates. The subsystem for calculating the entanglement entropy comprises qubits {Q5,3, Q6,3}. The dashed lines are linear fits used to extract the scaling exponent, and the shaded regions indicate the corresponding confidence intervals of the fits. Error bars indicate 1 standard deviation of the experimental results, averaged over 10 independent RMD sequences.

Full size image

Experimentally, for the 78-qubit system and n = {0, 1, 2}, we numerically fit the exponent based on the prethermal lifetime of the entanglement entropy plateau and obtain αS = {0.904, 2.976, 4.182}. On the basis of the imbalance decay, we obtain the exponents αI = {0.871, 2.727, 3.344}. In Fig. 3a, we summarize the scaling exponents αI, which reveal a positive correlation between αI and n across all investigated system sizes and configurations. More experimental data on the corresponding time evolutions can be found in Supplementary Information section 7A,B. In particular, the scaling exponents for n = 0 and n = 1 closely follow the theoretical prediction (black horizontal dashed lines in Fig. 3a). For n = 2, it is particularly challenging to obtain a converging scaling exponent, which requires an extremely fast drive with noiseless pulse control, and long evolution time (see detailed discussions in Supplementary Information section 5B). Here we implement approximately 1,000 driving cycles to clearly distinguish the scaling behaviours for n = 1 and n = 2. This also highlights the exceptional long-time stability and controllability of Chuang-tzu 2.0, which allows us to precisely engineer the temporal correlations and manipulate the quantum thermalization process.

To validate our experimental observations, we use advanced numerical techniques based on tensor-network representations52 to approximately simulate the many-body dynamics. We initially developed the PEPS with a single update for simulations. PEPS is a two-dimensional tensor network and is feasible for our circuit system. However, owing to the lack of canonical form, PEPS is not accurate in long-time simulations. The simulation results are shown in Fig. 3b, where we can see that shortly after a few driving cycles (t ≈ 0.1 μs), the PEPS simulations show significant deviations from the experimental data. To improve the numerical accuracy, we further developed the GMPS method. Compared with the conventional matrix product states (MPS) method, the GMPS method merges certain tensors in the conventional MPS into a giant tensor, thereby avoiding truncations inside the giant tensor and the swapping operations in MPS. This strategy is particularly efficient in simulating quasi-one-dimensional systems such as our system. As shown in Fig. 3b, for T = 4 ns, as the staggered potential strongly constrains the entanglement growth, GMPS correctly captures the experimental data at early times (also see results for T = 3 ns in Supplementary Information section 5C with high accuracy) whereas a clear deviation is visible for t > 1 μs. For slow driving (T = 7 ns), notable deviation already appears at early times, t ≲ 0.5 μs. In Supplementary Information section 1C, we systematically study the performance of GMPS and observe that its estimated fidelity significantly decays as heating occurs (fidelity drops below 0.05 when t ≈ 1.5 μs). The rapid growth of entanglement during this process limits both its accuracy and efficiency. Therefore, our quantum processor allows us to simulate the entire heating dynamics towards the highly entangled infinite-temperature states that is extremely costly to simulate classically53,54,55,56,57,58.

Entanglement dynamics and scaling

We further investigate the subsystem size dependence of the entanglement entropy. Specifically, we perform QST on various subsystems, each comprising up to 4 sites from the 78-qubit system. As shown in Fig. 4a, these configurations show distinct spatial arrangements: the subsystem \({\mathcal{A}}\) is aligned with the x axis, \({\mathcal{B}}\) is oriented along the y axis, and \({\mathcal{C}}\) forms a 2 × 2 square lattice. Figure 4a depicts the entanglement entropy dynamics for different subsystems, all exhibiting the long-lived prethermal regime. Moreover, the entanglement dynamics strongly depend on the specific choice of subsystems: \({\mathcal{A}}\) shows pronounced oscillatory dynamics during the prethermal regime, whereas \({\mathcal{B}}\) and \({\mathcal{C}}\) quickly saturate at a plateau, suggesting a non-uniform entanglement distribution in two dimensions. These oscillations originate from the coherent particle exchange between qubits with even and odd y index, stabilized by the staggered potential in the effective Hamiltonian \({\hat{H}}_{\mathrm{eff}}\), and can also be validated via GMPS simulations (Supplementary Fig. 4).

Fig. 4: Entanglement dynamics and volume-law scaling.
figure 4

a, Dynamics of the entanglement entropy for different subsystem configurations, involving 4 qubits in the 78-qubit system. All subsystems dynamics enter a long-lived prethermal regime, among which \({\mathcal{A}}\) exhibits pronounced oscillatory dynamics. Inset: the entanglement entropy averaged over the prethermal regime (30 ns ≤ t ≤ 100 ns, denoted between two black vertical dashed lines in a), as a function of volume V for subsystems \({\mathcal{A}}\) and \({\mathcal{B}}\). b, The volume (sV) and the area entanglement entropy (sA) per site, numerically fitted by analysing subsystems of varying volumes and areas at different times. The error bars correspond to 1 s.d. of the fitting parameter. c, Ratio sV/sA at different times, as a quantitative measure to distinguish between the area-law and the volume-law scaling. As time evolves, a clear crossover from the area-law to volume-law scaling is observed around 5–10 ns. Here we use 1-RMD protocol and T = 4 ns to perform the experiment.

Full size image

To quantify the entanglement entropy scaling, we calculate the averaged entropy, \({{\mathcal{S}}}_{{\rm{p}}{\rm{r}}{\rm{e}}}\), in the prethermal regime to reduce the temporal fluctuations. In the inset of Fig. 4a, we show \({{\mathcal{S}}}_{{\rm{p}}{\rm{r}}{\rm{e}}}\) as a function of the linear subsystem size V for \({\mathcal{A}}\) (red) and \({\mathcal{B}}\) (blue), and linear dependence can be observed. For other possible subsystem configurations, in practice, one can approximately distinguish different contributions using the ansatz50, \({\mathcal{S}}({\rho }_{{\rm{X}}})={s}_{A}{A}_{{\rm{X}}}+{s}_{V}{V}_{{\rm{X}}},\) where ρX denotes the reduced density matrix of a subsystem X, AX and VX correspond to the subsystem’s area and volume, respectively. The ratio sV/sA quantifies the degree to which the state adheres to area-law or volume-law entanglement scaling. By analysing 12 non-repetitive subsystems (Supplementary Information section 7D) by varying volumes and areas, we numerically fit sV and sA and show their time evolution in Fig. 4b and the ratio sV/sA in Fig. 4c. At early times (t ≤ 30 ns), the system stays close to the density-wave ordered product state, and sV/sA remains small. Indeed, a crossover of the values sV and sA occurs around 5–10 ns (Supplementary Information section 7C. A notable increase of sV/sA occurs around t ≈ 50 ns, and this ratio becomes substantially larger than 1, confirming that the system now exhibits volume-law entanglement scaling. At later times, entanglement further increases, such that the system’s dynamics become further computationally intractable.

Conclusion and outlook

We present a systematic experimental study of the non-equilibrium dynamics of two-dimensional interacting systems driven by the n-RMD protocol on a 78-qubit superconducting processor. Its precise and stable pulse sequence allows for the demonstration of the long-lived prethermal plateau and suppressed heating rates, which show the characteristic algebraic scaling T2n+1 in the high-frequency regime. These features are sufficiently universal, as they should generally appear in locally interacting systems with generic initial states, and are not limited to our current settings. The possibility of driving closed quantum systems with temporal randomness while avoiding heating paves the way for engineering non-equilibrium phases of matter beyond the conventional Floquet paradigm.

By performing QST, we observe the growth of the subsystem entanglement entropy during the entire heating process. Although the system eventually evolves towards the featureless infinite-temperature state, during the prethermal regime, we observe a non-uniform spatial entanglement distribution. In particular, subsystem \({\mathcal{A}}\) shows coherent oscillatory entanglement dynamics. These observations shed light on the understanding of the microscopic generation of entanglement, especially in higher dimensions, where the subsystem configurations can exhibit a rich geometric structure.

Furthermore, by analysing the entanglement entropy for various subsystems, we demonstrate the crossover from area-law to volume-law entanglement scaling as time evolves. The onset of heating further accelerates the entanglement growth. Therefore, our work highlights superconducting quantum processors as a powerful platform for exploring non-equilibrium phases of matter and heating dynamics in driven systems where classical simulation faces formidable challenges.

It will be interesting to further explore the initial-state dependence and the spatial non-homogeneity in the heating process, as well as the stability of many-body localization and anomalous topological phases59 in randomly driven systems. In addition, identifying different prethermalization mechanisms for much broader classes of non-periodic drives remains an open and interesting direction. Finally, although our experiment is performed on a superconducting processor, the underlying heating control mechanism is readily applicable to different quantum simulator platforms and can stabilize sought-after non-equilibrium phenomena in driven systems.

Methods

Experimental set-up

Our experiments are performed on a superconducting quantum processor, named Chuang-tzu 2.0, which comprises 78 qubits arranged in a square lattice configuration with 6 rows and 13 columns. Each pair of nearest-neighbour qubits is interconnected by an adjustable coupler, resulting in a total of 137 couplers within the processor. With all 78 superconducting qubits initialized at their idle points, we prepared the initial state using X gates, and the X-gate pulses are optimized to minimize the leakage to higher energy levels, achieving an average gate fidelity of 99.4%. Then, the RMD pulses are applied on all qubits to engineer the Hamiltonian, for experiments of different parameters. The states of all qubits can be read out simultaneously through the transmission lines coupled to readout resonators. All qubit probabilities are corrected to eliminate the measurement errors.

It is noteworthy that the qubit connectivity in our processor, while maintaining a square lattice structure, differs slightly from processors such as Sycamore53, which feature a square lattice with zig-zag edges. We use a tunable coupler with a capacitively connecting pad architecture45, which facilitates a 1,200-μm spacing between adjacent qubits, thereby ensuring sufficient wiring space. The qubits and couplers are fabricated on the qubit layer chip, whereas the control lines, readout lines and readout cavities are integrated into the wiring layer chip. These two chips are interconnected using flip-chip bonding technology. The readout cavities for the six qubits in each column are multiplexed onto a single readout line and are capacitively coupled to the qubits through interfacial capacitance. The control lines are similarly coupled to the qubits and couplers, enabling excitation and biasing via interfacial capacitance and mutual inductance.

RMD pulse

To implement the Hamiltonian described in the main text, the qubit frequency ωq needs to be rapidly modulated between two values ωq = ωc + δh × h0 for U+, whereas ωq = ωc − δh × h0 for U−, where ωc is the common qubit frequency. For the frequency-tunable transmon qubit, the relationship between ωq and the amplitude of Z pulse z is

$${\omega }_{{\rm{q}}}=\sqrt{8{E}_{\mathrm{JJ}}\,{E}_{{\rm{C}}}| \cos (kz+b)| }-{E}_{{\rm{C}}},$$
(6)

where EJJ denotes the Josephson energy, EC is the charging energy, and kz + b = Φext/Φ0 with the weak external flux Φext (Φext depends linearly on the amplitude of the Z pulse z, with k and b being the slope and intercept, respectively). The RMD pulse consists of a series of square Z waves with duration of T. Owing to the constrained sampling rate of the DAC (digital-to-analogue converter), both the falling and rising edge durations are inherently limited to a minimum of about 0.5 ns. In our experiments, these square waves are substituted to trapezoidal waves with edges of 0.5 ns.

To further characterize the RMD pulse, we measure the dynamical phase induced during its operation. Initially, the qubit is prepared at its idle point and excited using an X/2 gate. Then, the qubit is biased to the working point using the RMD pulse and the flat pulse following a delay, respectively. After turning off all Z pulses to tune the qubit back to its idle point, we apply another rotation \({R}_{\phi }(\frac{{\rm{\pi }}}{2})\), with ϕ ranging from 0 to 2π. The population of \(| 1\rangle \) state reaches its maximum only when ϕ compensates for the accumulated dynamical phase. In Extended Data Fig. 1a,b, we show the accumulated phase in the RMD case

$${\varphi }_{{\rm{r}}}(t)={\int }_{0}^{t}{\rm{d}}t\,[{\omega }_{{\rm{r}}}(t)-{\omega }_{{\rm{i}}{\rm{d}}{\rm{l}}{\rm{e}}}],$$
(7)

where ωr(t) is the qubit frequency under RMD, and the phase in the flat case

$${\varphi }_{{\rm{c}}}(t)={\int }_{0}^{t}{\rm{d}}t\,({\omega }_{{\rm{c}}}-{\omega }_{{\rm{i}}{\rm{d}}{\rm{l}}{\rm{e}}}).$$
(8)

The RMD sequence is precisely recovered by differentiating (φr − φc) with respect to time

$${\omega }_{{\rm{r}}}(t)-{\omega }_{{\rm{c}}}=\frac{{\rm{d}}}{{\rm{d}}t}[{{\varphi }}_{{\rm{r}}}(t)-{{\varphi }}_{{\rm{c}}}(t)].$$
(9)

The result, as depicted in Extended Data Fig. 1c, is well in accordance with the engineered sequence.

GMPS

The commonly used MPS representation for a two-dimensional lattice involves winding the MPS across the lattice, effectively forcing the system into a one-dimensional configuration. This transformation makes the nearest neighbouring interactions non-local, depending on the chosen ordering. To avoid the long-range interactions, we use the GMPS method, where we group Ly lattice sites in a single column, treating them as a single site with a physical dimension of \({2}^{{L}_{y}}\). This allows us to represent the system as an MPS of length Lx. The advantage of this method is that it preserves locality. However, the trade-off is that it is not scalable in the y direction.

We use the second-order Trotter–Suzuki decomposition \({\hat{H}}_{\pm }={\hat{H}}_{\pm }^{1}+{\hat{H}}^{2}\) to approximate the time-evolution operators as follows

$${\hat{H}}_{\pm }^{1}=\sum _{x,y}{h}_{\pm }{\hat{a}}_{x,y}^{\dagger }{\hat{a}}_{x,y}+J({\hat{a}}_{x,y+1}^{\dagger }{\hat{a}}_{x,y}+{\rm{H.c}}\,.\,),$$
(10)
$${\hat{H}}^{2}=\sum _{x,y}J({\hat{a}}_{x+1,y}^{\dagger }{\hat{a}}_{x,y}+{\rm{H.c}}\,.\,).$$
(11)

The single-site operators \({\hat{U}}_{\pm }^{1}=\exp (-{\rm{i}}{\hat{H}}_{\pm }^{1}t)\) are \({2}^{{L}_{y}}\times {2}^{{L}_{y}}\) matrices and can be directly applied to each tensor. For two-site hopping terms, we take Ly = 4 as an example and summarize update rules below (Extended Data Fig. 2).

  1. 1.

    Perform QR decompositions on each tensor to isolate the active physical bonds (i4 and j4 here) into R tensors. This step is essential for ensuring that the calculations are in the most computationally efficient manner.

  2. 2.

    Apply the two-site time-evolution operator and contract into a single tensor.

  3. 3.

    Execute a singular value decomposition and eliminate the singular values that are smaller than a specified threshold. Subsequently, truncate the results by retaining the largest χ singular values.

  4. 4.

    Absorb the singular values and vectors into the updated tensors.

Owing to the canonical form in MPS, it is possible to monitor fidelity throughout the update process. During each truncation step, the singular values s1 ≥ s2 ≥ … ≥sD are arranged in descending order, and only the top χ(<D) singular values are kept. Consequently, the fidelity is given by the ratio:

$$\frac{{\sum }_{k=1}^{\chi }{s}_{k}^{2}}{{\sum }_{k=1}^{D}{s}_{k}^{2}},$$
(12)

This expression provides a measure of how much of the original state’s information is preserved after truncation.