Research ArticlePHYSICS

Probing dynamical phase transitions with a superconducting quantum simulator

See allHide authors and affiliations

Science Advances  17 Jun 2020:
Vol. 6, no. 25, eaba4935
DOI: 10.1126/sciadv.aba4935


Nonequilibrium quantum many-body systems, which are difficult to study via classical computation, have attracted wide interest. Quantum simulation can provide insights into these problems. Here, using a programmable quantum simulator with 16 all-to-all connected superconducting qubits, we investigate the dynamical phase transition in the Lipkin-Meshkov-Glick model with a quenched transverse field. Clear signatures of dynamical phase transitions, merging different concepts of dynamical criticality, are observed by measuring the nonequilibrium order parameter, nonlocal correlations, and the Loschmidt echo. Moreover, near the dynamical critical point, we obtain a spin squeezing of −7.0 ± 0.8 dB, showing multipartite entanglement, useful for measurements with precision fivefold beyond the standard quantum limit. On the basis of the capability of entangling qubits simultaneously and the accurate single-shot readout of multiqubit states, this superconducting quantum simulator can be used to study other problems in nonequilibrium quantum many-body systems, such as thermalization, many-body localization, and emergent phenomena in periodically driven systems.


Quantum simulation uses a controllable quantum system to mimic complex systems or solve intractable problems (1, 2). Emergent phenomena in out-of-equilibrium quantum many-body systems (3), e.g., thermalization (4) versus localization (5), and time crystals (6), have all been recently studied using quantum simulation. Recently, the dynamical phase transition (DPT) and the nonequilibrium phase transition in transient time scales have been theoretically studied in the transverse-field Ising model with all-to-all interactions (79). These two transitions can be characterized by the nonequilibrium order parameter (710) and the Loschmidt echo associated with the Lee-Yang-Fisher zeros in statistical mechanics (11), respectively. Moreover, recent experimental progress has allowed the controllable simulation of these exotic phenomena with cold atoms (12, 13) and trapped ions (14, 15). Yet, experimental explorations for the dynamics of entanglement, as a valuable resource in quantum information processing, remain limited in the presence of a DPT.

In our experiments, applying a sudden change of the transverse field with a controllable strength, we drive the system, initially in its ground state, out of equilibrium. Accurate single-shot readout techniques enable us to synchronously record the dynamics of all qubits and to observe essential signatures of DPTs and spin squeezing from the dynamical criticality in the Lipkin-Meshkov-Glick (LMG) model. It is worth mentioning that our experimental system is a 16-qubit device featuring all-to-all connectivity, which complements the type of superconducting circuits used in other simulations of many-body physics (1620), where neighboring couplings dominate. The presence of long-range interactions is essential for realizing the LMG model.

This work presents a systematic quantum simulation of DPTs with two different concepts, providing evidence of the relation between the nonequilibrium order parameter and the Loschmidt echo. We verify entanglement in spin-squeezed states generated from dynamical criticality, directly observing squeezing of −7.0 ± 0.8 dB for 16 qubits.


Our quantum simulator is a superconducting circuit with 20 fully controllable transmon qubits capacitively coupled to a resonator bus ℛ (Fig. 1A). Sixteen qubits (Q1 to Q16), with XY-control lines, are selected to perform experiments (see Materials and Methods). The resonant frequency of ℛ is fixed at about 5.51 GHz, while the qubit frequencies are individually tunable via their Z-control lines, enabling us to engineer the qubit-qubit interactions induced by ℛ. We detune all 16 qubits from ℛ by, e.g., Δ/2π≃ −450 MHz to switch on the resonator-mediated interactions between two arbitrary qubits (21). Simultaneously, identical resonant microwave drives, with a magnitude of hx, are imposed on all qubits to generate the local transverse fields for the control of a DPT (Fig. 1C). To ensure the uniformity of the local fields, the cross-talk effects of microwave pulses have been precisely corrected, and the microwave phase has been calibrated (see the Supplementary Materials). The effective Hamiltonian of the quenched system isH1/=i<jNλij(σi+σj+σiσj+)+hxj=1Nσjx(1)where N = 16, λijgigj/Δ+λijc is the qubit-qubit coupling strength, gj represents the coupling strength between ℛ and Qj, and gigj/Δ is the resonator-induced virtual coupling strength, which acts as a dominant part (parameters are shown in the Supplementary Materials). Because the values of λij are nearly the same for most pairs of qubits and do not decay over distance ∣ij∣ (Fig. 1B), the quenched system can be reasonably approximated by the LMG model, whose Hamiltonian isHLMG=(J/N)(Sz)2+μSxwith Sx,zjσjx,z/2 and μ = 2hx (see Materials and Methods). The LMG model was first introduced in nuclear physics (22) and then used to describe two-mode Bose-Einstein condensates (23). Recent studies (710, 14) have shown that HLMG has a dynamical critical point separating the dynamical paramagnetic phase (DPP) and the dynamical ferromagnetic phase (DFP) with and without a global ℤ2 symmetry, respectively.

Fig. 1 Quantum simulator and experimental pulse sequences.

(A) False-color optical micrograph of the device highlighting various circuit elements such as qubits (red), the resonator bus (black), qubit XY-control lines (blue), and Z-control lines (green). (B) Connectivity graph of the 16-qubit system when all qubits are equally detuned from the resonator bus by Δ/2π ≃ −450 MHz, with the colored straight lines representing the magnitude of the qubit-qubit couplings. Four pairs of qubits (Q3 and Q14), (Q4 and Q13), (Q5 and Q12), and (Q6 and Q11) have relatively small couplings because of their noticeable cross-talk couplings λijc that neutralize the resonator-induced parts. (C) Experimental pulse sequences for simulating the DPT. First, the qubits are initialized in the ∣00…0〉 state at their corresponding idle frequencies. Then, the rectangular pulses and resonant microwave pulses are applied almost simultaneously to realize the quantum quench. Last, the 16-qubit joint readout is executed, yielding the probabilities {P00…0, P00…1, …, P11…1}, from which σjz can be calculated. When necessary, single-qubit rotation pulses Rj(θj,ϕj)=exp[iθj(cosϕjσjx+sinϕjσjy)/2] (in black dotted box) are applied in advance to bring the axis defined by (θj, ϕj + π/2) in the Bloch sphere of Qj to the σz direction before the readout.

First, we show that our programmable superconducting qubits can simulate and verify the DPT by measuring the magnetization and the spin correlation. The system is initialized at the eigenstate ∣00…0〉 of H1 with hx = 0, where ∣0〉 denotes the ground state of a qubit. Then, we quench the system by suddenly adding a transverse field and monitor its dynamics from t = 0 to 600 ns. With the precise full control and the high-fidelity single-shot readout of each qubit, we are able to omnidirectionally track the evolutions of the average magnetizationσα(t)1Nj=1Nσjα(t)along the x, y, z axes for different strengths of the quenched transverse fields, with α ∈ {x, y, z}. By depicting the trajectory of the Bloch vector σ=[σx,σy,σz], the dynamics of our quantum simulator with two distinct transverse fields is visualized in Fig. 2A. For a small transverse field, e.g., hx/2π ⋍ 2 MHz, 〈σz(t)〉 exhibits a slow relaxation (Fig. 2B). However, given a strong transverse field, e.g., hx/2π ⋍ 8 MHz, 〈σz(t)〉 exhibits a large oscillation at an early time and approaches zero in the long-time limit (Fig. 2B). In Fig. 2C, we show the behavior of the time-averaged magnetizationσz¯(1/tf)0tf dtσz(t)defined as the nonequilibrium order parameter. Figure 2C demonstrates that σz¯0 and σz¯=0 in the DFP and the DPP, respectively. The experimental data of σz¯ for qubits with different detunings Δ are presented in the Supplementary Materials. In addition, the Bloch vector length σ also depends on the strength of the transverse field hx. For large hx, σ decays rapidly to a small value, indicating strong quantum fluctuations in the DPP (Fig. 2D).

Fig. 2 Magnetization and spin correlation.

(A) Experimental (left) and numerical (right) data of the time evolution of the average spin magnetization shown in the Bloch sphere for different strengths of the transverse fields. (B) Time evolution of the magnetization 〈σz(t)〉. (C) Nonequilibrium order parameter σz¯, as a function of hx/2π. (D) Dynamics of the Bloch vector length σ. (E) Averaged spin correlation Czz¯ versus hx/2π. The regions with light red and light blue in (C) and (E) show the DFP and DPP, respectively, separated by a theoretically predicted critical point hcx/2π5.7MHz. The solid curves in (B) to (E) are the numerical results using the Hamiltonian of our experimental system without considering decoherence.

Figure 2E shows the averaged spin correlation functionCzz¯(1/tf)0tfdt ijσiz(t)σjz(t)/N2versus hx with a final time tf = 600 ns, where the DPT is characterized by the local minimum of two-spin correlations. We can observe the critical behaviors of σz¯ and Czz¯ as the signatures of the DPT, when the transverse field strength is set near the theoretical prediction hcx/2π=Nλ/8π5.7 MHz, with λλij¯ (see Materials and Methods). In the Supplementary Materials, we also present the experimental results of σz¯ and Czz¯ for 12 and 8 qubits, clarifying finite-size effects. Experimentally, we still observe the DPT signatures down to 8 qubits.

Another perspective on dynamical criticality is based on the Loschmidt echo, defined as L(t)=∣⟨00…0∣eiH1t/∣00…0⟩∣2, where the time t, satisfying L(t) = 0, is a Lee-Yang-Fisher zero. The zero will cause the non-analytical behavior of the rate function r(t) = − N−1 log [L(t)], regarded as the complex-plane generalization of the free-energy density (11, 24). Recent numerical studies (7, 9, 24, 25) have revealed that the existence of Lee-Yang-Fisher zeros closely relates to the DPT between the DFP and the DPP in long-range interacting systems. In Fig. 3A, we show the distinct behaviors of the Loschmidt echo in different dynamical phases. In the DPP (hx/2π ⋍ 8 MHz), the Loschmidt echo decays rapidly to a near-zero value, related to the occurrence of the non-analyticity of the rate function r(t). A clearer signature can be seen from the first minimum of the Loschmidt echo Lmin(1) as a function of hx (Fig. 3B). In the Supplementary Materials, we clarify the reason of choosing Lmin(1) as a circumstantial probe of the dynamical criticality, and the numerical results of the rate function r(t) using the real parameters of our quantum simulator are also presented as reference. The direct observation of the non-analytical points of the rate function r(t), as a diagnostic signature of the dynamical criticality, deserves further experimental investigations.

Fig. 3 Loschmidt echo.

(A) Time evolution of the Loschmidt echo L(t) for different transverse field strengths. (B) Earliest minimum point of L(t) during its dynamics, Lmin(1), as a function of hx/2π. It is shown that Lmin(1) is close to zero in the DPP, while it becomes relatively large in the DFP. The behavior of Lmin(1) versus hx is similar to that in the LMG model (see the Supplementary Materials). The solid curves in (A) and (B) are the numerical results using Eq. 1 without considering decoherence.

In addition to demonstrating a DPT, the LMG model is also useful for generating spin-squeezed states with twist-and-turn dynamics (26, 27). Near the equilibrium critical point, spin squeezing can be achieved, originating from quantum fluctuations, according to the Heisenberg uncertainty principle (28). Similarly, we show that spin-squeezed states can also be generated from dynamical criticality. During the dynamics of the quenched Hamiltonian Eq. 1, we can visualize the spin-squeezed state by measuring the quasidistribution Q-function (29) Q(θ, ϕ) ∝ ⟨θ, ϕ∣ρ(t)∣θ, ϕ⟩, where θ,ϕ=j=1N(cosθ20j+sinθ2eiϕ1j) is the spin coherent state. The measurement is realized by applying a single-qubit rotation to bring the axis defined by (θ, ϕ) in the Bloch sphere to the z axis for each qubit before the joint readout. The experimental and numerical data of Q(θ, ϕ) are compared in Fig. 4A, which show spin squeezing with a large strength of the external field, due to stronger quantum fluctuations in the DPP (see also Fig. 2C).

Fig. 4 Quasidistribution Q-function and spin-squeezing parameter.

(A) Experimental and numerical data of Q(θ, ϕ) in spherical coordinates, when the minimum values of the spin-squeezing parameters are achieved during the time evolutions with the strengths of the transverse fields hx/2π ≃ 3 and 6 MHz, respectively. (B) Time evolution of the spin-squeezing parameters with hx/2π ≃ 3 and 6 MHz, respectively. (C) Minimum spin-squeezing parameter ξmin2 as a function of hx. The solid lines in (B) are the numerical results using the Hamiltonian of our experimental system without considering decoherence. The blue shaded area in (B) is only accessible for entangled states. The dotted line in (C) is the piecewise linear fit, whose minimum point is close to the theoretically predicted critical point hcx/2π5.7MHz (dashed line).

We also measured the time-evolved spin-squeezing parameter (26) (see the Supplementary Materials)ξ2=4minn[Var(Sn)]/N()2where n denotes an axis perpendicular to the mean spin direction, and Var(Sn)=(Sn)2Sn2. In Fig. 4B, we show that ξ2 < 1, as a sufficient condition for particle entanglement (30, 31), occurs in the time interval t ≲ 46 ns when hx/2π ⋍ 3 MHz, and for t ≲ 38 ns when hx/2π ⋍ 6 MHz. The minimum spin-squeezing parameter over time, ξmin2, as a function of hx is shown in Fig. 4C, where the minimum value ξmin20.2(7.0 dB) is attained very close to the critical point of the DPT. Compared with the theoretical limit, about N−2/3, of the squeezing parameter for an N-body one-axis twisting Hamiltonian (30), our 16-qubit system achieves a spin-squeezing parameter satisfying ξmin2Nα, with α ≃ 0.58. This indicates the high-efficiency generation of the spin-squeezed state from dynamical criticality and reveals a potential application of the DPT to quantum metrology.


We have presented clear signatures and entanglement behaviors of the DPT in the LMG model with a superconducting quantum simulator featuring all-to-all connectivity, including the nonequilibrium order parameter, Loschmidt echo, and spin squeezing. On the basis of its high degree of controllability, precise measurement, and long decoherence time, our platform with all-to-all connectivity is powerful for generating multipartite entanglement (29, 32) and investigating nontrivial properties of out-of-equilibrium quantum many-body systems, such as many-body localization (33, 34), quantum chaos in Floquet systems (35), and quantum annealing (36).


Device information and system Hamiltonian

The device used here consists of 20 frequency-tunable superconducting qubits capacitively coupled to a central resonator bus. It is the same circuit presented in (29), where more details about the device, the qubit manipulation, and the readout can be found. In table S1, we present the characteristics for the quantum simulator involving 16 of the 20 qubits, with XY-control lines, which have been relabeled in the experiments.

The unused four qubits in this device, without XY-control lines, are detuned far off resonance from the other 16 qubits to avoid interacting with them during the experiments. Thus, they will not be included in the following descriptions. The system Hamiltonian, without applying external transverse fields, can be written asHS1/=ωRaa+j=116[ωj(t)1j1j+gj(σj+a+σja)]+i<j16λijc(σi+σj+σiσj+)where ω and ωj represent the fixed resonant frequency and the tunable frequency of Qj, respectively, while gj is the coupling strength between the Qj and resonator bus. The magnitude of the cross-talk coupling between Qi and Qj beyond the resonator-induced virtual coupling is denoted as λijc. When equally detuning all the 16 qubits from the resonator bus by about Δ/2π ≃ −450 MHz, and simultaneously applying resonant microwaves to each qubit, the system Hamiltonian can be transformed toHS2/=i<j16(λijc+gigj/Δ)(σi+σj+σiσj+)+j=116hjx(σjeiϕj+σj+eiϕj)with gigj/Δ being the magnitude of the resonator-mediated coupling between Qi and Qj. It acts as a dominant part of the qubit-qubit interaction terms, because the cross-talk coupling λijc is much smaller. In Fig. 1B, we plot the connectivity graph of the total coupling strength λij for all the combinations of pairs of qubits. The individually controllable amplitude and the phase of the microwave drive on each Qj are represented by hjx and ϕj, respectively. In our experiments, we set the uniform amplitude and phase for all qubits, leading to the Hamiltonian in Eq. 1. To ensure this uniformity, the calibration process for the microwave drives is described in the Supplementary Materials.

Relation between the quantum simulator and the LMG model

Our device can be described by the Hamiltonian in Eq. 1. With uniform couplings λλij¯, the first term of Eq. 1 can be written asλi<j16(σi+σj+H.c.)=(J/N)[S2(Sz)2]where JNλ. The second term can be directly rewritten as hxi=116σix=μSx, with μ = 2hx. According to [S2, Sα] = 0 (α ∈ {x, y, z}), and the fact that the initial state ∣00…0〉 is an eigenstate of S2, we haveexp[i(H1/)t]000exp(iHLMGt)000indicating that the dynamical properties of the device H1 can be approximately expressed as the ones of the LMG modelHLMG=(J/N)(Sz)2+μSx

The location of the DPT critical point of the LMG model is μc = ∣J∣/2, leading to hcx=Nλ/4. Note that we only roughly estimate the location of the dynamical critical point of the LMG model. The numerical simulations in the main text are based on the Hamiltonian of the quantum simulator described by Eq. 1.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank C. Song, Q. Guo, Z. Wang, and X. Zhang for technical support. Devices were made at the Nanofabrication Facilities at the Institute of Physics in Beijing and National Centre for Nanoscience and Technology in Beijing. The experiment was performed on the quantum computing platform at Zhejiang University. Funding: This work was supported by the National Basic Research Program of China (grant nos. 2016YFA0302104, 2016YFA0300600, and 2017YFA0304300), the National Natural Science Foundation of China (grants nos. 11934018, 11725419, and 11904393), the Strategic Priority Research Program of Chinese Academy of Sciences (grant no. XDB28000000), the China Postdoctoral Science Foundation (grant no. 2018 M640055), the Zhejiang Province Key Research and Development Program (grant no. 2020C01019), the ARO (grant no. W911NF-18-1-0358), the JST Q-LEAP program, the JST CREST (grant no. JPMJCR1676), the JSPS Kakenhi (grant no. JP20H00134), the JSPS Postdoctoral Fellowship (grant no. P19326), FQXi, and the NTT PHI Lab. Author contributions: D.Z., H.F., and H.W. conceived the research. K.X., Z.-H.S., Y.-R.Z., H.F., and H.W. designed the experiment. H.L. and D.Z. fabricated the device. K.X. and W.L. performed the experiment. K.X. and Z.-H.S. did numerical simulations. K.X., Z.-H.S., Y.-R.Z., and H.F. analyzed the results. K.X., Z.-H.S., Y.-R.Z., F.N., H.F., and H.W. wrote the paper. All authors contributed to the experimental setup, discussions of the results, and development of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article