Observation of strong and weak thermalization in a superconducting quantum processor
Fusheng Chen, Zheng-Hang Sun, Ming Gong, Qingling Zhu, Yu-Ran Zhang, Yulin Wu, Yangsen Ye, Chen Zha, Shaowei Li, Shaojun Guo, Haoran Qian, He-Liang Huang, Jiale Yu, Hui Deng, Hao Rong, Jin Lin, Yu Xu, Lihua Sun, Cheng Guo, Na Li, Futian Liang, Cheng-Zhi Peng, Heng Fan, Xiaobo Zhu, Jian-Wei Pan
OObservation of strong and weak thermalization in a superconducting quantum processor
Fusheng Chen,
1, 2, 3, ∗ Zheng-Hang Sun,
4, 5, ∗ Ming Gong,
1, 2, 3, ∗ Qingling Zhu,
1, 2, 3
Yu-Ran Zhang, Yulin Wu,
1, 2, 3
Yangsen Ye,
1, 2, 3
Chen Zha,
1, 2, 3
Shaowei Li,
1, 2, 3
Shaojun Guo,
1, 2, 3
Haoran Qian,
1, 2, 3
He-Liang Huang,
1, 2, 3
Jiale Yu,
1, 2, 3
Hui Deng,
1, 2, 3
Hao Rong,
1, 2, 3
Jin Lin,
1, 2, 3
Yu Xu,
1, 2, 3
Lihua Sun,
1, 2, 3
Cheng Guo,
1, 2, 3
Na Li,
1, 2, 3
Futian Liang,
1, 2, 3
Cheng-Zhi Peng,
1, 2, 3
Heng Fan,
4, 5, 7, 8
Xiaobo Zhu,
1, 2, 3 and Jian-Wei Pan
1, 2, 3 Hefei National Laboratory for Physical Sciences at Microscale and Department of Modern Physics,University of Science and Technology of China, Hefei, Anhui 230026, China Shanghai Branch, CAS Center for Excellence and Synergetic Innovation Center in Quantum Information and Quantum Physics,University of Science and Technology of China, Shanghai 201315, China Shanghai Research Center for Quantum Sciences, Shanghai 201315, China Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100190, China Theoretical Quantum Physics Laboratory, RIKEN Cluster for Pioneering Research, Wako-shi, Saitama 351-0198, Japan Songshan Lake Materials Laboratory, Dongguan 523808, Guangdong, China CAS Center for Excellent in Topological Quantum Computation,University of Chinese Academy of Sciences, Beijing 100190, China
We experimentally study the ergodic dynamics of a 1D array of 12 superconducting qubits with a transversefield, and identify the regimes of strong and weak thermalization with different initial states. We observe con-vergence of the local observable to its thermal expectation value in the strong-thermalizaion regime. For weakthermalization, the dynamics of local observable exhibits an oscillation around the thermal value, which canonly be attained by the time average. We also demonstrate that the entanglement entropy and concurrence cancharacterize the regimes of strong and weak thermalization. Our work provides an essential step towards ageneric understanding of thermalization in quantum systems.
Statistical mechanics is developed to describe the ther-modynamics of both classical and quantum systems. Ifa balloon, containing a large number of air molecules,is pierced in an evacuated chamber, the molecules willmove around all possible states, and with the long-timeaverage, they will satisfy the Maxwell velocity distribu-tion, which is independent of the initial condition [1].In quantum cases, thermal equilibrium states, describedby the statistical-mechanical prescription, are expected toemerge in the non-equilibrium dynamics of a closed non-integrable many-body system [2–7]. Different from theclassical case, the time average is not necessary for quan-tum thermalization, and the quenched states could con-verge to their thermal expectations at a time after a shortrelaxation [3]. This phenomenon is regarded as strongthermalization. However, strong thermalization cannotbe certainly achieved by a non-integrable system drivenout of equilibrium, whose occurrence is also relevant tothe choice of initial states. Numerical works [8, 9] haverevealed that in a non-integrable 1D Ising model, strongthermalization happens when the effective inverse tem-perature of initial states is close to 0. In contrast, if the ef-fective inverse temperature of initial states is sufficientlyfar away from 0, the temporal evolution of the local ob-servable shows an obvious oscillation, with the long-time ∗ Those authors contributed equally to this work. average attaining the thermal expectation value. Thisphenomenon is known as weak thermalization. Recently,it has been numerically shown that regimes of strong andweak thermalization exist in the long-range Ising modeldescribing trapped ions [10]. Nevertheless, the experi-mental observation of both strong and weak thermaliza-tion remains absent.On the basis of the high-precision control, long co-herence time, and the accurate readout, a superconduct-ing quantum processor is an excellent platform for gen-erating multipartite entangled states [11–13], character-izing quantum supremacy [14–16], and demonstratingvariational quantum computation [17, 18]. Moreover,by performing analog quantum simulations, the platformis also employed to study the phenomena in quantummany-body systems out of equilibrium, including quan-tum walks [19], many-body localization [20–22], dynam-ical phase transitions [23], and ergodic-localized junc-tions [24].Here, we realize a non-integrable system using a 1Darray of 12 superconducting qubits with a controllabletransverse field. We experimentally observe the signa-tures of strong and weak thermalization via measuringthe local observable with different initial states. Sincethe description of the local observable, using statisticalmechanics, relies on the local entropy created by entan-glement, the dynamics of the entanglement entropy (EE)plays a key role in thermalization [25–29]. Thus, westudy the EE of the single-qubit subsystems and show a r X i v : . [ qu a n t - ph ] F e b Q Q Q | ⟩ | ⟩ | ⟩ Initialization Time evolution Measurement Density of states, ( ) N o r m a li z ed ene r g y , energy, Normalized Q i | ⟩ QST Q Q Q Q Q Q Q Q Q Q Q Q Qubit (MHz)
Initial time, (a)(c)(b) (d)(e)
FIG. 1. (a) Architecture of the 24-qubit superconducting circuit. The qubits Q - Q are employed to realize a 1D non-integrablemodel. (b) The schematic diagram of the pulse sequence used to observe strong and weak thermalization, which consists of threeparts, i.e., initialization, evolution and measurement. For the initialization, we prepare all qubits at ∣ θ , φ ⟩ , i.e., Eq. (4). Startingfrom the ∣ ⟩ state, all qubits are prepared at ∣ θ , φ ′ ⟩ via the gate ˆ R ( θ , φ ′ ) rotating ∣ ⟩ around the axis ˆ n = cos ( φ ′ ) ˆ σ y − sin ( φ ′ ) ˆ σ x by a angle θ . However, due to the dynamical phases accumulated in tunning qubits to the evolution point, φ ′ is not equal to φ . We calibrate the dynamical phases, and then correct all qubits to the same initial state ∣ θ , φ ⟩ . (c) The control waveformscorresponding to the pulse sequence in (b). The single-qubit pulses after ∣ ⟩ refer to the gate ˆ R ( θ , φ ) . To realize the non-equilibrium quantum dynamics, the Z pulses and resonant microwave pulses (the sinusoidal line) are simultaneously applied to eachqubit. After the evolution, quantum state tomography measurements are performed at idle points of the qubits and single-qubitpulses are required. (d) For the Hamiltonian (3), the density of states ρ ( (cid:15) ) as a function of the normalized energy (cid:15) obtained by thenumerical simulation. The normalized energy of two initial states ∣ π, ⟩ and ∣ π / , π / ⟩ are highlighted. (e) The normalized energyof the initial state ∣ θ , φ ⟩ as a function of θ and φ . that the EE can distinguish the strong-thermalizationregime from the weak one. Furthermore, we measurethe concurrence [30] of the reduced density matrices oftwo nearest qubits, employing the tomographic readout,and observe thermal entanglement [31] in the presence ofweak thermalization.The experiment is performed on a chain of 12 super-conducting transmon qubits [see Fig. 1(a)], described bythe Hamiltonian of the Bose-Hubbard model [19, 32–34] ˆ H = J ∑ j = ( ˆ a † j ˆ a j + + H.c. )+ U ∑ j = ˆ n j ( ˆ n j − ) + ∑ j = µ j ˆ n j , (1)where ˆ n j = ˆ a † j ˆ a j is the bosonic number operator, with ˆ a † j ( ˆ a j ) being the creation (annihilation) operator, and J , U , and µ n denote the nearest-neighbor coupling strength,the on-site nonlinear interaction, and tunable on-site po-tential, respectively. More details of our device can befound in Supplementary Materials [35]. With U / J → ∞ and all µ n being tuned to the samefrequency, the Hamiltonian (1) can be rewritten as [36,37] ˆ H = λ ∑ j = ( ˆ σ xj ˆ σ xj + + ˆ σ yj ˆ σ yj + ) , (2)with λ = J / and ˆ σ αj ( α ∈ { x, y, z }) being the Paulimatrices. The Hamiltonian (2) describes an integrablemodel, which can be exactly solved via introducing theJordan-Wigner transformation [38]. To realize a non-integrable system where thermalization occurs, we im-pose resonant microwave drives with a magnitude g ≃ λ on all qubits, generating a local transverse field [23]. Thefinal effective Hamiltonian of the system is ˆ H = λ ∑ j = ( ˆ σ xj ˆ σ xj + + ˆ σ yj ˆ σ yj + ) + g ∑ j = ˆ σ yj . (3)For details regarding the Hamiltonian (3), see Supple-mentary Materials [35].To observe strong and weak thermalization, we initial-ize the system by preparing each qubit in the direction Time (ns) Lo c a l ob s e r v ab l e , Time (ns) Lo c a l ob s e r v a t i on , Time (ns) Lo c a l ob s e r v ab l e , Time (ns)Time (ns) E n t ang l e m en t en t r op y E n t ang l e m en t en t r op y Time (ns) E n t ang l e m en t en t r op y (a) (b) (c)(d) (e) (f) FIG. 2. (a) Experimental data of the time evolution of the local observable ˆ σ z ( t ) with the initial state ∣ π, ⟩ . (b) As in (a), but forthe initial state ∣ π / , π / ⟩ . (c) As in (a), but for the initial state ∣ π / , π / ⟩ . (d) Experimental data of the time evolution of theentanglement entropy with the initial state ∣ π, ⟩ . (e) As in (d), but for the initial state ∣ π / , π / ⟩ . (f) As in (d), but for the initialstate ∣ π / , π / ⟩ . The horizontal lines in (d)-(f) denote the Page value of the EE. The solid lines in (a)-(f) are numerical results withoutconsidering decoherence. The solid line in (f) is numerical results considering decoherence (see Supplementary Materials [35] forthe effects of decoherence). The shaded region shows the errorbars of the numerical results, taking the uncertainties of the localfield into consideration (see Supplementary Materials [35]). The initial states are presented in Bloch spheres, where ∣ X +⟩ , ∣ Y +⟩ ,and ∣ Z +⟩ are the eigenstate of σ x , σ y and σ z with the eigenvalue + , respectively. ( θ , φ ) , which can be described as the spin coherentstate ∣ θ , φ ⟩ = ∏ j = ( cos θ ∣ + Z ⟩ j + e − iφ sin θ ∣ − Z ⟩ j ) , (4)where ∣ + Z ⟩ j ( ∣ − Z ⟩ j ) denotes the eigenstate of ˆ σ zj with the eigenvalue + ( − ). Next, all qubits are bi-ased to the working point to start the quench dynamics, ∣ Ψ t ⟩ = e − i ˆ Ht ∣ θ , φ ⟩ . We then tune the qubits to their idlepoints, and perform the quantum state tomography to re-construct the one- and two-qubit density matrices. Theexperimental pulse sequence and control waveforms areshown in Fig. 1(b) and (c), respectively. There are threeessential experimental requirements to be satisfied: ( i )To realize the time evolution, all qubits should be tunedto the same frequency; ( ii ) The initial states of all qubitsshould be uniform at the start point of the time evolu-tion; ( iii ) The local transverse fields of all qubits shouldbe uniform during the evolution. These requirementsare fulfilled with specific calibrations (see Supplemen-tary Materials for details [35]).The occurrence of strong or weak thermaliza-tion relates closely to the effective inverse tempera-ture β of ∣ θ , φ ⟩ , which can be obtained by solv-ing Tr [ ˆ H (∣ θ , φ ⟩⟨ θ , φ ∣ − ˆ ρ β )] = , with ˆ ρ β = e − β ˆ H / Tr ( e − β ˆ H ) being the thermal state [8]. Moreover,the quasiparticle explanation of weak thermalizationindicates that initial states in the weak-thermalizationregime are near the edge of the energy spectrum [9].Here, we first consider two initial states ∣ θ , φ ⟩ =∣ π / , π / ⟩ and ∣ π, ⟩ whose effective inverse temperatureare numerically estimated as Jβ ≃ − . and , lyingin the weak- and strong-thermalization regime, respec-tively. In addition to the effective inverse temperature, theregimes of strong and weak thermalization can be identi-fied by defining the normalized energy of the initial state ∣ θ , φ ⟩ (cid:15) = ⟨ θ , φ ∣ ˆ H ∣ θ , φ ⟩ − E min E max − E min (5)with E max and E min being the maximum and minimumeigenvalue of ˆ H , respectively. In Fig. 1(d), the (cid:15) ofthe initial state ∣ π, ⟩ corresponds to the maximum den-sity of states (DoS) ρ ( (cid:15) ) , while the (cid:15) of the initial state ∣ π / , π / ⟩ is close to the edge with ρ ( (cid:15) ) ≃ . Moreover,in Fig. 1(e), we plot the normalized energy of differentinitial states ∣ θ , φ ⟩ , i.e., (cid:15) ( θ , φ ) .We start by characterizing strong and weak ther-malization employing the local observable ˆ σ z ( t ) = ∑ j = ⟨ Ψ t ∣ ˆ σ zj ∣ Ψ t ⟩ . Figure 2(a) and (c) present the ex- T i m e - a v e r aged en t ang l e m en t en t r op y (b)(a) FIG. 3. (a) Experimental data of the time-averaged EE with ini-tial states ∣ θ , φ ⟩ plotted in the φ − θ plane. (b) Experimentaldata of the time-averaged EE with initial states ∣ π / , φ ⟩ . Thewhite dashed line in (a) highlights θ = π / . The solid linein (b) shows the numerical result considering decoherence (seeSupplementary Materials [35] for the effects of decoherence).The shaded region shows the errorbars of the numerical results,taking the uncertainties of the local field into consideration (seeSupplementary Materials [35]). perimental results of the time evolution of ˆ σ z ( t ) withinitial states ∣ π, ⟩ and ∣ π / , π / ⟩ , respectively. It isshown that for the initial state ∣ π, ⟩ in the strong-thermalization regime, ˆ σ z ( t ) stably achieve the thermalvalue Tr ( ˆ ρ β ˆ σ zj ) = after t ≃ ns. In contrast, for theinitial state ∣ π / , π / ⟩ in the weak-thermalization regime, ˆ σ z ( t ) strongly oscillates around the thermal value 0. Inaddition, we measure the dynamics of ˆ σ z ( t ) with theinitial state ∣ π / , π / ⟩ , which also lies in the strong-thermalization regime, since its effective inverse tem-perature is Jβ ≃ . The results, depicted in Fig. 2(b),show that even the behavior of short relaxation is differ-ent from that with the initial state ∣ π, ⟩ , the local observ-able also has a stationary value near the thermal value after t ≃ ns, which is a hallmark of strong thermal-ization.Next, we consider the von Neumann entanglement en-tropy (EE), S = − Tr [ ˆ ρ j ln ( ˆ ρ j )] , where ˆ ρ j is the reduceddensity matrix of the j -th qubit. We average the EE overall qubit sites. The dynamics of the EE, with the ini-tial states ∣ π, ⟩ and ∣ π / , π / ⟩ in the regime of strongthermalization, are displayed in Fig. 2(d) and (e), respec- tively. We observe that for strong thermalization, the EErapidly reaches the Page value S Page ≃ . as the max-imum EE of a single-qubit subsystem of a total systemin the random pure state [39]. However, for weak ther-malization, the non-equilibrium dynamics gains the EEsmaller than the Page value [Fig. 2(f)].Furthermore, we study the time-averaged EE between ns and ns with different initial states ∣ θ , φ ⟩ .In Fig. 3(a), we show the experimental data of time-averaged EE with different initial states ∣ θ , φ ⟩ , whichbears a close resemblance of to the normalized energy inFig. 1(e). Specifically, with θ = π / , around φ ≃ π / , (cid:15) ≃ . , and the DoS ρ ( (cid:15) ) becomes the maximum [seeFig. 1(d) and (e)]. Thus, it can be predicted that strongthermalization occurs in this regime. Additionally, ac-cording to the results in Fig. 1(d) and (e), the normalizedenergy is near 1 at φ = π / , where the DoS is close to 0,and the degree of thermalization is the weakest. The ex-perimental data of the time-averaged EE, with θ = π / and different φ , are presented in Fig. 3(b). There is aminimum of the EE around φ = π / , corresponding tothe weakest thermalization. Moreover, the maximum EEreveals a regime of strong thermalization with θ = π / and φ ∈ [ . π, . π ] .The trace distance between the non-equilibrium state ˆ ρ t = ∣ Ψ t ⟩⟨ Ψ t ∣ and the thermal state ˆ ρ β , with the β be-ing the effective inverse temperature of the initial state,i.e., Tr (∣ ˆ ρ t − ˆ ρ β ∣) , can diagnose quantum thermaliza-tion [25]. It has been numerically shown that the dis-tance monotonically decays to 0 for strong thermal-ization. With initial states in the weak-thermalizationregime, the decay of the distance can also be observedbut with a strong fluctuation [8].We measure the reduced density matrix ˆ ρ j,j + t of thesubsystem consisting of the j -th and ( j + ) -th qubit, us-ing the quantum state tomography. For the initial state ∣ π, ⟩ , the effective inverse temperature is Jβ = , andthe corresponding two-qubit thermal state is ˆ I / with ˆ I being an identity matrix. Then, the trace distance be-tween ˆ ρ j,j + t and ˆ I / , averaged over the qubit site j , canbe directly obtained. Similarly, by considering the ther-mal states ˆ ρ j,j + β with Jβ ≃ − . , the dynamics ofthe trace distance, with the initial state ∣ π / , π / ⟩ , canalso be measured. As shown in Fig. 4(a) and (b), thetrace distance decays during the time evolution for bothstrong and weak thermalization, indicating the tendency ˆ ρ j,j + t ≃ ˆ ρ j,j + β . However, for weak thermalization, thetrace distance strongly oscillates [Fig. 4(b)].Finally, we experimentally investigate the concurrenceof the two-qubit reduced density matrix ˆ ρ j,j + t , which isdefined as E ( ˆ ρ ) = max { , √ γ − √ γ − √ γ − √ γ } ,where γ , ..., γ are eigenvalues listed in decreasing orderof the matrix Γ = ˆ ρ ( ˆ σ y ⊗ ˆ σ y ) ˆ ρ ∗ ( ˆ σ y ⊗ ˆ σ y ) [30]. The time Time (ns) D i s t an c e Time (ns) D i s t an c e Time (ns) C on c u rr en c e Time (ns) C on c u rr en c e (a)(c) (b)(d) FIG. 4. (a) Experimental data of the trace distance between thenon-equilibrium and thermal states with the initial state ∣ π, ⟩ .(b) As in (a), but for the initial state ∣ π / , π / ⟩ . (c) Experimen-tal data of the concurrence with the initial state ∣ π, ⟩ . (d) Asin (c), but for the initial state ∣ π / , π / ⟩ . The solid lines in (a)and (c) are numerical results without considering decoherence.The solid lines in (b) and (c) are numerical results consideringdecoherence (see Supplementary Materials [35] for the effectsof decoherence). The shaded region shows the errorbars of thenumerical results, taking the uncertainties of the local field intoconsideration (see Supplementary Materials [35]). The initialstates are presented in Bloch spheres, where ∣ X +⟩ , ∣ Y +⟩ , and ∣ Z +⟩ are the eigenstate of σ x , σ y and σ z with the eigenvalue + , respectively. evolution of the concurrence, with the initial state ∣ π, ⟩ and ∣ π / , π / ⟩ , are presented in Fig. 4(c) and (d), respec-tively. We observe that the concurrence vanishes after t ≃ ns for strong thermalization. Whereas, the concur-rence preserves a finite value with the initial state in theweak-thermalization regime, which can be interpreted asthe thermal entanglement, i.e., the concurrence in thermalstates ˆ ρ β [31], according to ˆ ρ j,j + t ≃ ˆ ρ j,j + β as a result of the ergodic dynamics in the weak-thermalization regime.The numerics of the concurrence in thermal states ˆ ρ β with different β are presented in Supplementary Mate-rials [35].In summary, we have provided a clear experimental ev-idence for characterizing the regimes of strong and weakthermalization. Weak thermalization, with a slow growthof the EE, has the potential for generating states withlong-lived coherence and stabilizing the phases of matterfar away from equilibrium, such as Floquet symmetry-protected topological phases [40], discrete time crys-tals [41], many-body localized phase [6, 42, 43], and dy-namical paramagnetic and ferromagnetic phases [23, 44].Our work also indicates that the 1D array of supercon-ducting qubits can be a promising platform for exploringthe issues at the foundation of quantum thermodynamics. ACKNOWLEDGMENTS
The authors thank the USTC Center for Micro- andNanoscale Research and Fabrication. The authors alsothank QuantumCTek Co., Ltd. for supporting the fab-rication and the maintenance of room temperature elec-tronics. This research was supported by the National KeyR & D Program of China (Grants No. 2018YFA0306703,No. 2017YFA0304300, No. 2016YFA0302104, No.2016YFA0300600), the Chinese Academy of Sciences,and Shanghai Municipal Science and Technology Ma-jor Project (Grant No. 2019SHZDZX01), the StrategicPriority Research Program of Chinese Academy of Sci-ences (Grant No. XDB28000000), Japan Society forthe Promotion of Science (JSPS) Postdoctoral Fellow-ship (Grant No. P19326), JSPS KAKENHI (Grant No.JP19F19326), the National Natural Science Foundationof China (Grants No. 11574380, No. 11905217, No.11934018, No. 11774406), the Key-Area Research andDevelopment Program of Guangdong Province (GrantNo. 2020B0303030001), and Anhui Initiative in Quan-tum Information Technologies. [1] G. Gallavotti,
Statistical mechanics: A Short Treatise (Springer, Berlin, 1999).[2] J. M. Deutsch, Phys. Rev. A , 2046 (1991).[3] M. Srednicki, Phys. Rev. E , 888 (1994).[4] M. Rigol, V. Dunjko, and M. Olshanii, Nature , 854(2008).[5] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat-tore, Rev. Mod. Phys. , 863 (2011).[6] R. Nandkishore and D. A. Huse, Annu. Rev. Condens.Matter Phys. , 15 (2015).[7] J. Eisert, M. Friesdorf, and C. Gogolin, Nature Physics , 124 (2015).[8] M. C. Ba˜nuls, J. I. Cirac, and M. B. Hastings, Phys. Rev.Lett. , 050405 (2011).[9] C.-J. Lin and O. I. Motrunich, Phys. Rev. A , 023621(2017).[10] F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang,C. Monroe, and A. V. Gorshkov, Phys. Rev. Lett. ,150601 (2019).[11] M. Gong et al. , Phys. Rev. Lett. , 110501 (2019).[12] C. Song et al. , Science , 574 (2019).[13] K. X. Wei et al. , Phys. Rev. A , 032343 (2020). [14] F. Arute et al. , Nature , 505 (2019).[15] S. Boixo et al. , Nature Physics , 595 (2018).[16] C. Neill et al. , Science , 195 (2018).[17] F. Arute et al. , Science , 1084 (2020).[18] M.-C. Chen et al. , Phys. Rev. Lett. , 180501 (2020).[19] Z. Yan et al. , Science , 753 (2019).[20] Q. Guo et al. , Nature Physics (2020), 10.1038/s41567-020-1035-1.[21] K. Xu et al. , Phys. Rev. Lett. , 050507 (2018).[22] B. Chiaro et al. , (2019), arXiv:1910.06024 [cond-mat.dis-nn].[23] K. Xu et al. , Science Advances , eaba4935 (2020).[24] C. Zha et al. , Phys. Rev. Lett. , 170503 (2020).[25] A. M. Kaufman et al. , Science , 794 (2016).[26] C. Neill et al., Nature Physics , 1037 (2016).[27] Y. O. Nakagawa, M. Watanabe, H. Fujita, and S. Sugiura,Nature Communications , 1635 (2018).[28] H. Kim and D. A. Huse, Phys. Rev. Lett. , 127205(2013).[29] S. Popescu, A. Short, and A. Winter, Nature Physics ,754 (2006).[30] W. K. Wootters, Phys. Rev. Lett. , 2245 (1998).[31] M. C. Arnesen, S. Bose, and V. Vedral, Phys. Rev. Lett. , 017901 (2001).[32] Q. Zhu et al. , (2021), arXiv:2101.08031 [quant-ph].[33] P. Roushan et al. , Science , 1175 (2017).[34] J. Koch et al. , Phys. Rev. A , 042319 (2007).[35] Supplementary Material.[36] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, andM. Rigol, Rev. Mod. Phys. , 1405 (2011).[37] M. Cramer, A. Flesch, I. P. McCulloch, U. Schollw¨ock,and J. Eisert, Phys. Rev. Lett. , 063001 (2008).[38] E. Lieb, T. Schultz, and D. Mattis, Annals of Physics ,407 (1961).[39] D. N. Page, Phys. Rev. Lett. , 1291 (1993).[40] A. C. Potter, T. Morimoto, and A. Vishwanath, Phys. Rev.X , 041001 (2016).[41] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Annu.Rev. Condens. Matter Phys. , 467 (2020).[42] E. Altman, Nature Physics , 979 (2018).[43] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[44] J. Zhang et al. , Nature , 601 (2017). Supplementary Materials for ‘Observation of strongand weak thermalization in a superconducting quan-tum processor’
I. DEVICE INFORMATION AND SYSTEMHAMILTONIAN
The experiment is performed on a ladder-type su-perconducting circuit consisting of 24 transmon qubits,which is the same circuit presented in Ref. [1]. One cansee more details about the device, the readout and thequbit manipulation in Tab. S1 or Ref. [1]. As shown inthe Fig. 1(a) in the main text, the qubits Q − Q are em-ployed to the quantum simulation of an one-dimensionalsystem.Without the external transverse field, the Hamiltonianof the one-dimensional (1D) array of superconductingqubits can be written as ˆ H = ∑ j = J j,j + ( ˆ a † j ˆ a j + + H.c. )+ ∑ j = U j n j ( ˆ n j − ) + ∑ j = µ j ˆ n j . (6)The values of the hopping interaction between the j − thand ( j + )− th qubit, i.e., J j,j + are presented in the Fig.1(a) in the main text and Tab. S1, with an average value J / π ≃ . MHz. The average value of the onset non-linear interactions U j is U / π ≃ − . MHz. Becauseof ∣ UJ ∣ ≃ , the Bose-Hubbard model (6) is quite close tothe hard-core limit where the bosonic creation and anni-hilation operator are mapped to the spin raising and low-ering operator. When all qubits are tuned to the sameworking point, the Hamiltonian with the hard-core limitcan be written as ˆ H XX = ∑ j = J j,j + ( ˆ σ + j ˆ σ − j + + H.c. )≃ λ ∑ j = ( ˆ σ xj ˆ σ xj + + ˆ σ yj ˆ σ yj + ) (7)with λ = J / . By employing the Jordan-Wigner trans-formation ˆ σ − j = exp (− iπ ∑ m < j ˆ c † m ˆ c m ) ˆ c j and the Fouriertransformation ˆ d k = √ N ∑ Nj = ∑ Nj = exp (− ijk ) ˆ c j , theHamiltonian can be rewritten as ˆ H = ∑ k J cos k ˆ d † k ˆ d k with the momentum k taking N (the number of qubit)discrete values from the Brillouin zone k = πnN with n ∈ {− N / + , ..., N / } . Consequently, a 1D array ofsuperconducting qubits can be regarded as an integrablesystem described by free fermions [3].Numerical works have shown that strong and weakthermalization occur in 1D non-integrable systems. To observe strong and weak thermalization, we apply res-onant microwaves to each qubit, which induces a localtransverse field described by [2] ˆ H drive = ∑ j = g j ( e − iϕ j ˆ σ + j + e iϕ j ˆ σ − j ) (8)with g j and ϕ j as the strength and phase of the localtransverse field imposed on the j − th qubit. In our ex-periment, we adopt that g / π = g j / π = MHz, and ϕ = ϕ j = π / . Finally, the Hamiltonian of the qubit chainwith a local transverse field reads ˆ H = ˆ H XX + ˆ H drive = λ ∑ j = ( ˆ σ xj ˆ σ xj + + ˆ σ yj ˆ σ yj + ) + g ∑ j = ˆ σ yj . (9)It is recognized that when we tune all qubits to theworking point employing Z pluses, there is a disorderin the chemical potential, i.e., ˆ H disorder = ∑ j = µ j ˆ σ + j ˆ σ − j ,due to the crosstalk of the Z pluses. Here, we estimatethat µ j / π ∈ [− , ] MHz. In addition, the crosstalkof the local transverse field should also be considered.The strength of the local transverse field g j would satisfy g j / π = . ± MHz (see below comparison between thenumerics and experimental data).We then discuss the integrability breaking originatedfrom the local transverse field. Without the local trans-verse field, in the system (7), the total number ∑ j = ˆ σ + j ˆ σ − j is conserved. The local transverse field can break theconservation. Moreover, we can study the level spac-ing distribution to characterize the non-integrability ofthe system with the local transverse field. Let E n as theeigenenergy of a Hamiltonian, and the level separationuniformity is defined as r n = min { s n , s n − } max { s n , s n − } , (10)where s n = E n + − E n is the nearest-neighbor spacings.It has been postulated that for a non-integrable system inthe ergodic phase, the statistics of energy levels followsthe Gaussian orthogonal ensemble (GOE) [4], i.e., P GOE ( r ) = r + r ( + r + r ) / . (11)As shown in Fig. S1, the level statistics of the systemHamiltonian (9) follows the GOE and the superconduct-ing qubit chain with the local transverse field is a non-integrable system. II. CALIBRATION
In this experiment, in order to realize a non-integrablesystem where thermalization occurs, there are three re-quirements we have to meet: (1) All qubits should be Q Q Q Q Q Q Q Q Q Q Q Q ω read / π (GHz) 6.688 6.729 6.790 6.832 6.886 6.927 6.701 6.754 6.801 6.855 6.909 6.954 ω maxq / π (GHz) 4.928 5.536 4.962 5.600 4.887 5.600 4.941 5.562 4.904 5.602 4.905 5.587 ω idleq / π (GHz) 4.835 5.31 4.693 5.38 4.82 5.26 4.68 5.32 4.74 5.27 4.8 5.4 T ( µ s) 20.6 24.3 28.1 24.0 26.3 27.3 21.4 25.4 21.4 18.8 22.2 18.7 T ∗ ( µ s) 3.9 2.2 2.0 2.0 4.9 1.8 2.1 2.0 2.3 1.9 4.4 2.8 U / π (MHz) − − − − − − − − − − − − J / π (MHz) 12.4 12.3 12.4 12.4 12.2 13.3 13.6 13.7 13.8 13.7 13.6 f (%) 97.6 98.7 98.2 98.7 98.9 97.9 98.9 98.9 98.6 99.1 97.8 96.7 f (%) 91.2 94.2 94.1 91.6 93.6 88.5 93.8 94.8 92.3 95.6 94.1 92.6 Q RB fidelity (%) 99.86 99.88 99.87 99.88 99.93 99.91 99.94 99.87 99.74 99.94 99.94 99.91
TABLE S1. Qubit performance. Parameters of the device. w read / π is the resonant frequency of the readout resonator; w max q is themaximum frequency of the qubit; w idle q is qubit’s idle frequency; T and T ∗ are the energy relaxation time and Ramsey dephasingtime measured at the idle frequency; U is the anharmonicity of qubit. J / π is the coupling strength of the corresponding qubit-pairmeasured at the working frequency ( . GHz). Q RB fidelity is the average gate fidelity of single-qubit X / gate measured atidle frequency. The gate time is 30 ns. FIG. S1.
Level statistics of the system Hamiltonian.
Thedashed line is the probability distribution of r n following theGOE. The solid line is the numerics of P ( r ) of the system (9). tuned to the same interacting frequency; (2) The initialstate of all qubits should be uniform at the start point ofthe time evolution; (3) The local transverse field of allqubits should be unifrom. We use the multi-qubit exci-tation propagation to calibrate frequency alignment of allqubits [1], meeting the requirement (1). In the followingsection, we explain how to calibrate the initial states [therequirement (2)] and the transverse fields [the require-ment (3)]. A. Calibration of the initial state
In our experiment, the initial state of each qubit is pre-pared using a single qubit gate ˆ R ( θ , φ ) , which is re-alized by resonant microwave pulses. It is required thatall 12 qubits are in the same initial state at the start pointof the time evolution, as shown in the Fig. 1(b) and (c)of the main text. Therefore, for each qubit, both the θ and φ need extra calibrations. The θ is determined bythe amplitude and length of the microwave pulse. We canperform Rabi oscillation measurements on each qubit atthe idle frequencies to ensure the uniformity of the θ .Moreover, the φ is determined by the phase of the mi-crowave pulse. As the qubits will accumulate extra dy-namical phases during the process tuned to the workingfrequency, we still need to calibrate the overall phasesbefore the time evolution.The main idea of the calibration is that when twonearest-neighboring qubits are tuned in resonance, if both θ and φ are uniform, there will not be any swapping be-tween them. The calibration process of the initial state isshown in Fig. S2 a . We take the qubit Q and Q as anexample. We select the phase of the microwave pulse im-posed on Q as the frame of reference. We firstly applyan X π / rotation pulse to Q and a π / pulse whose rota- T i m e ( n s ) Time abc
FIG. S2.
Phase alignment of the rotation pulse employed tothe initialization. a,
Experimental sequence for phase align-ment of the initial state. b, The experimental result for phasealignment of the initial state. We plot the measured excitedprobabilities P Q as a function of the phase φ and the evolu-tion time t . c , The standard deviation of P Q with different φ .The dashed lines in b and c refer to the phase φ with minimumstandard deviation of P Q . tion axis is of phase φ ′ in the x - y plane to Q . Then, webias Q and Q to the working frequency ω I with a rect-angular pulse. After a period of evolution, we bias themto their idle frequencies ω idleq and measure the populationof state ∣ ⟩ , i.e., P Q . When the initial states of Q and Q are the same, there will not be any swapping and themeasured P Q will not oscillate during the dynamics.The experimental results of P Q are displayed inFig. S2 b . To extract the phase φ with the weakest oscilla-tion, we also study the standard deviation of P Q with re-spect to the evolution time t , i.e., σ ( P Q ) . The σ ( P Q ) as a function of φ is displayed in Fig. S2 c . The phase used T i m e ( n s ) -0.2 -0.4 abc FIG. S3.
Phase alignment of the transverse field a,
Exper-imental sequence for phase alignment of the transverse field. b, The experimental result for phase alignment of the trans-verse field. We plot the measured difference of probabilities ∆ P = P Q − P Q as a function of the phase φ and the evolu-tion time t . c , The standard deviation of ∆ P with different φ .The dashed lines in b and c refer to the phase φ with minimumstandard deviation of ∆ P . for the correction is marked by the black dashed verti-cal line. We calibrate all nearest-neighboring qubits pairswith the same method. After that, all qubits are preparedto the same initial state before the time evolution. B. Calibration of the transverse field
The experiment requires the uniformity of the localtransverse field for each qubit, which means the g j and ϕ j in Eq. (8) should be the same for all qubits. The pa-rameters g j is determined by the amplitude of microwave0pulse A j .There are two factors we need to consider in calibra-tion: (1) The difference of the microwave phases causedby the length difference of qubits¡ ¯ control lines. (2)The XY-crosstalk caused by the microwave leakage. Ourgoal is to get 12 pairs of microwave pulse parameters ( A j , ϕ ′ j ) , which ensure the uniformity of all transversefields. The initial value of ϕ ′ j is set to 0. We use thefollowing procedure to calibrate the transverse field pa-rameters.1. By tuning each qubit Q j to the interacting frequencyand performing Rabi oscillation measurements, we canget the initial values of the driving amplitude A j ( j = , , ..., ) with Ω j / π = MHz.2. We use sequence in Fig. S3 a to calibrate the phaseof transverse field. Taking Q and Q as an example.We start with initializing them in the ground state, thenbias them to the interacting frequency, and apply resonantmicrowave drives on them with the same magnitude buta phase difference of ϕ ′ j . Note that when driving Q and Q , we simultaneously drive other qubits to realize theeffect of XY-crosstalk from other qubits. After a periodof evolution, we bias Q and Q to their idle frequenciesand measure their population of ∣ ⟩ . When the phases ofthe transverse fields are aligned, the difference betweenthe population of Q and Q , which is defined as ∆ P = P Q − P Q , will not oscillate. The experimental resultsof ∆ P are displayed in Fig. S3 b . The standard deviationof ∆ P with respect to the evolution time t , i.e., σ ( ∆ P ) ,is plotted in Fig. S3 c . We perform this calibration to allnearest-neighboring qubits pairs. Then, the microwavephases ϕ ′ j are updated.3. Again, we tune each qubit Q j to the interacting fre-quency and run Rabi oscillation measurements. There isa little difference from step 1 to take the XY-crosstalkinto account. When driving Q j , we simultaneously ap-ply microwave pulses to other qubits which are in idlefrequencies. Then we can update driving amplitudes A j .4. We iterate over step (2) and step (3) until the mi-crowave pulse parameters A j and ϕ ′ j converge. Typi-cally, we only need two to three iterations. III. NUMERICS OF THE LOCAL OBSERVABLE ANDENTANGLEMENT ENTROPY
For the local transverse field ˆ H driven = ∑ j = g j ˆ σ yj , weexperimentally set g / π = g j / π = MHz. However,if we numerically simulate the dynamics of the local ob-servable σ z ( t ) using g / π = g j / π = MHz, there isan obvious discrepancy between the numerics and ex-perimental data (see the solid line and circle points inFig. S4 c ). Thus, the crosstalk effect should be consid-ered in detail. Time (ns) -0.4-0.200.20.4 0 300150
Time (ns) Lo c a l ob s e r v ab l e , Lo c a l ob s e r v ab l e , Time (ns) -0.6-0.4-0.200.20.4 a bc
FIG. S4.
Numerics of the local observable in comparisonwith the experimental data with the initial state ∣ π / , π / ⟩ .a, Numerical data of the local observable with different disor-der strength W . b, Numerical data of the local observable withdifferent strength of the local transverse field g / π . c, The com-parison between the experimental and numerical data.
Time (ns) E n t ang l e m en t en t r op y Exp. dataDecoherenceIsolated system
FIG. S5.
The dynamics of the EE with the initial state ∣ π / , π / ⟩ . The solid line is the numerical data obtained fromunitary evolution. The dashed line is the numerical data consid-ering decoherence effects. The diamond points are the experi-mental data.
We assume that the crosstalk can induce a disorder ofthe local transverse field, i.e., g j / π ∈ ( g + [− W, W ])/ π MHz with g / π = MHz and W as the strength of dis-order. The numerical results of the local observable withdifferent W is shown in Fig. S4 a . One can see that theamplitude of σ z ( t ) is suppressed with the increase of W . Additionally, we assume that the strength of the lo-cal transverse field g / π can be enlarged by the crosstalkeffect. As shown in Fig. S4 b , the frequency of the os-cillation of σ z ( t ) is closely related to g / π . We considerthe two parameters g and W , and find that when we chose1 a b c T i m e - a v e r aged en t ang l e m en t en t r op y FIG. S6.
The time-averaged EE. a,
Numerical data of the time-averaged EE for the initial states ∣ π / , φ ⟩ ( φ ∈ [ , π ] ), withoutdecoherence, i.e., considering unitary evolution. b, Numerical data of the time-averaged EE for the initial states ∣ π / , φ ⟩ , takingthe decoherence effects into consideration. c, Experimental data of the time-averaged EE for the initial states ∣ π / , φ ⟩ . g / π = . MHz and W / π = MHz, the numerics is ingood agreement with the experimental data.We also present the numerics of the local observableand entanglement entropy (EE) with the initial states instrong-thermalization regime [see the Fig. 2(a), (b), (d)and (e) in the main text]. The numerical simulations alsoconsider g / π = . MHz and W / π = MHz and theresults are consistent with the experimental data.
IV. THE IMPACT OF DECOHERENCE ON THEENTANGLEMENT ENTROPY
The above numerics mainly consider an isolated sys-tem evolved by the unitary transformation, i.e., ∣ Ψ t ⟩ = e − i ˆ Ht ∣ θ , φ ⟩ . However, the decoherence is inevitablein experiments. The results in Fig. S4 c show a goodagreement between the numerical and experimental data,indicating that the dynamics of local observable in theregimes of both strong and weak thermalization, and thedynamics of the EE for strong thermalization are robustagainst decoherence effects.As shown in Fig. S5, there exists an obvious differ-ence between the numerics without considering decoher-ence and the experimental data. Consequently, we shouldnumerically estimate the impacts of decoherence on thedynamics of the EE in the weak-thermalization regime.The dynamics with decoherence can be modeled by theLindblad master equation [5,6]d ˆ ρ t d t = − i [ ˆ H, ˆ ρ t ]+ ∑ j = ( L j ˆ ρ t ˆ L † j − { ˆ L † j ˆ L j , ˆ ρ t }) (12)with ˆ L j ( ˆ L † j ) as the Lindblad operators and ˆ ρ t = ∣ Ψ t ⟩⟨ Ψ t ∣ .There are two effects of decoherence: the energy relax-ation and dephasing effect, identified by the energy life-time T and dephasing time T , respectively. The corre-sponding Lindblad operators are ˆ L j = ˆ σ − j /√ T (energy Time (ns) C on c u rr en c e Time (ns) C on c u rr en c e , isolated system, isolated system, decoherence, decoherence a b FIG. S7.
Numerics of the concurrence in comparison withthe experimental data. a,
Experimental data of the concur-rence with different initial states. b, Numerical results of theconcurrence with and without considering decoherence effects. relaxation effect) and ˆ L j = ˆ σ zj /√ T (dephasing effect).For our device, the averaged energy lifetime is T ≃ . µ s, and the averaged dephasing time T ≃ . µ s. Thenumerical results of the dynamics of EE with decoher-ence effects displayed in Fig S5 are obtained by solvingEq. (12) with T ≃ . µ s and T ≃ . µ s. It is shownthat the numerics with decoherence effects are consis-tent with the experimental data. The numerics of time-averaged entanglement entropy without and with consid-ering the decoherence effects are presented in Fig S6 a and b , respectively. In comparison with the experimentaldata in Fig S6 c , the numerical data, taking the decoher-ence effects into consideration, are in better agreementwith the experimental data. V. NUMERICS OF THE CONCURRENCE
In Fig. S7, we present the comparison between the nu-merical and experimental results of the concurrence. Weshow that the dynamics of concurrence with the initialstate ∣ π / , π / ⟩ in weak-thermalization regime is also ap-parently influenced by decoherence effects. We find anoverall decrease of the concurrence as a function of timein the presence of decoherence effects.2 -2 -1 0 1 200.10.20.3 C on c u rr en c e FIG. S8.
Numerics of the thermal entanglement quantifiedby the concurrence.
The concurrence of thermal states ˆ ρ j,j + β ,averaged over the qubit site j , as a function of the inverse tem-perature β . In addition, we calculate the concurrence of thethermal state ˆ ρ j,j + β ≡ Tr A ∉{ j,j + } ˆ ρ β with ˆ ρ β ≡ exp (− β ˆ H )/ Tr [ exp (− β ˆ H )] , and average the concur-rence over the qubit site j = , , ..., . The averagedconcurrence with different inverse temperature β is dis-played in Fig. S8. As shown in Fig. 1 in the maintext, the β of the initial state ∣ π / , π / ⟩ in the weak-thermalization regime is estimated as Jβ ≃ − . . InFig. S8, the concurrence of the thermal state ˆ ρ j,j + β with Jβ ≃ − . , while the concurrence vanishes for thethermal state ˆ ρ j,j + β with − . ≤ Jβ ≤ . . Therefore,the finite value of the concurrence with the initial state ∣ π / , π / ⟩ , shown in Fig. S7, can be interpreted as a ther-mal entanglement in the state ˆ ρ j,j + β with Jβ ≃ − . ,based on ˆ ρ j,j + t → ˆ ρ j,j + β (see the main text). Next, we discuss the trace distance between thequenched states and the thermal states. For thequenched states, the reduced density matrices ˆ ρ j,j + t ≡ Tr A ∉{ j,j + } ˆ ρ t can be measured via the two-qubit quantumstate tomography performed on the j − th and ( j + )− thqubit. Whereas, the thermal state ˆ ρ j,j + β can only be nu-merically estimated. We emphasize that the thermal state ˆ ρ β is dependent on the Hamiltonian ˆ H . As we discussedin Section 2, there are random parameters in the Hamilto-nian ˆ H , i.e., ˆ H drive = ∑ j = ( g + δg j ) ˆ σ yj , where g / π = . MHz and δg j ∈ [− W, W ] with W / π = MHz.For the initial state ∣ π, ⟩ in the strong-thermalizationregime, the estimated inverse temperature of the initialstate ∣ π, ⟩ is β = which is independent of the ran-dom parameters δg j . Therefore, the thermal state is ˆ I = diag ( , , , ) . Nevertheless, for the initial state ∣ π / , π / ⟩ in the weak-thermalization regime, the esti-mated inverse temperature β is dependent on the randomparameters. Here, we use 20 samples of δg j and estimatethe β of ∣ π / , π / ⟩ . Then we can calculate the thermalstates ˆ ρ j,j + β for the 20 Hamiltonian and average themover the 20 samples.[1] Q. Zhu, et al. , arXiv: 2101.08031.[2] K. Xu, et al. , Sci. Adv. , eaba4935 (2020).[3] E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. ,407 (1961).[4] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux,Phys. Rev. Lett. , 084101 (2013).[5] R. Johansson, P. Nation, and F. Nori, Comput. Phys.Commun. , 1760 (2011).[6] R. Johansson, P. Nation, and F. Nori, Comput. Phys.Commun.184