Exploring the quantum critical behaviour in a driven Tavis–Cummings circuit

Quantum phase transitions play an important role in many-body systems and have been a research focus in conventional condensed-matter physics over the past few decades. Artificial atoms, such as superconducting qubits that can be individually manipulated, provide a new paradigm of realising and exploring quantum phase transitions by engineering an on-chip quantum simulator. Here we demonstrate experimentally the quantum critical behaviour in a highly controllable superconducting circuit, consisting of four qubits coupled to a common resonator mode. By off-resonantly driving the system to renormalize the critical spin-field coupling strength, we have observed a four-qubit nonequilibrium quantum phase transition in a dynamical manner; that is, we sweep the critical coupling strength over time and monitor the four-qubit scaled moments for a signature of a structural change of the system's eigenstates. Our observation of the nonequilibrium quantum phase transition, which is in good agreement with the driven Tavis–Cummings theory under decoherence, offers new experimental approaches towards exploring quantum phase transition-related science, such as scaling behaviours, parity breaking and long-range quantum correlations.

6.05 6.1 6.15 6.2 6.25 6. Illustration of the cubic potential as a function of the junction phase δ in a Josephson-junction device with, e.g., five energy states. The two-level phase qubit is encoded using the two lowest energy states |0⟩ and |1⟩. Transition frequencies ω n,n−1 /2π are separated from each other by about 100 MHz. (b) The swap spectroscopy measurement showing the qubit decay rate at different frequencies: after exciting the qubit to |1⟩, the qubit is biased to a fixed frequency (the x-axis) and its |1⟩-state probability P 1 (color scale) is measured after a delay time (the y-axis). The Lorentzian-shaped chevron indicates the qubit-resonator interaction and its center position is the resonator frequency. (c) A typical T 1 decay curve (points) obtained along the dashed line in (a) at ω q /2π ≡ ω 10 /2π = 6.129 GHz. Line is a fit with T 1 = 511 ns. (d) Fitted T 1 values as a function of the qubit frequency as obtained from data in (a). T 1 values nearby the resonator frequency are not shown since a simple decay fit is not accurate in this region. (e) A histogram analysis of the T 1 values in (c). Line is a Gaussian fit with a mean value of 500 ns and a standard deviation of 30 ns. It is seen that ⟨n r ⟩ at ∆ r /2π = −10 MHz increases quickly with the sweep time t, which may lead to serious state leakage via the qubit-resonator photon swap process. As such ∆ r /2π = −10 MHz should be avoided experimentally. (b) Illustration showing the procedure of obtaining P L by properly lowering the potential barrier to ∆U ′ (∆U ′′ ) and measuring the corresponding tunneled probability P n≥2 (P n≥1 ). P L ≈ P n≥2 − 0.05 · (P n≥1 − P n≥2 ), where the factor 0.05 corresponds to the stray tunneled probability when the qubit is in |1⟩ and the potential barrier is lowered to ∆U ′ . (c) Experimentally obtained P L for each qubit circuit (legend on the far right) as functions of λ/λ c for data shown in   (a) Pulse sequence for the Ramsey-interference-like measurement with the phase of the second π/2 pulse given in Supplementary Equation 10 (Seq. 1), which is equivalent to the Ramsey interference measurement using two identical rotation pulses (Seq. 2) if the steering of ω q (t) is as expected. (b) P 1 versus t using the two pulse sequences in (a) for exemplary qubits as indicated (data for Q 2 and Q 4 are similar to Q 1 's). Using Seq. 1, ω q (t)/2π of each qubit is tuned from 5.87 GHz to 6.14 GHz in 75 ns. Resemblance of the two curves for the same qubit indicates that we achieve a fairly accurate dynamical control of the qubit frequency over a short period. The dip around 55 ns for Q 3 using Seq. 1 is due to the interference by a TLS at 6.06 GHz. Figure 6: Two-qubit (N = 2) ⟨J z ⟩ behaviour in comparison with numerical simulation.
. Lines are numerical simulation results. Value of Ω is calibrated via the on-resonant interaction as described in Methods of the main text.

Unitary transformation into an interacting representation
The four-qubit system in our experiment can be described by a driven TC model in units of = 1, where ω j q (ω r ) is the resonance frequency of the j-th spin (the single mode of the resonator), the latter with raising (lowering) operators a † (a), σ j ± , σ j z are the j th spin Pauli operators, λ j is the coupling strength between the j-th spin and the resonator mode, Ω (ω d ) is the strength (frequency) of the driving applied to the resonator mode. In accordance with the experimental setup, we have neglected the small term regarding the drive on the qubits, which can also be absorbed into the term of Ω under the homogeneous assumption due to the fact that the driving effects on the resonator and the qubits are unitarily equivalent [1]. In the rotating frame with respect to the driving field, the Hamiltonian is rewritten as with detunings ∆ j q = ω j q − ω d and ∆ r = ω r − ω d . Considering the homogeneous case, i.e., identical qubits and couplings, we have with ∆ q = ∆ j q , λ = λ j and the collective spin operators J z = The eigen-solution of Supplementary Equation 3 can be obtained in a displaced frame by displaced Fock states [1]. In the isolated system, to see quantum phase transition (QPT), we only need to focus on the properties of the ground state in the variation of the coupling λ with respect to the critical point λ c = √ ∆ q ∆ r .

In the thermodynamical limit
To understand the finite-spin QPT we have to approach the thermodynamical limit (i.e., infinite number of spins), which is fortunately solvable analytically (see Supplementary Material in [1]). Using scaled ground-state values of ⟨J z ⟩/(N/2) = 2β/N − 1 with β = 0 for λ < λ c and β = (N/2)(1 − ∆ q ∆ r /λ 2 ) for λ ≥ λ c , we find for the latter the cusp-like behaviour around the critical point with the exponent γ z = 1. Similarly, we have ⟨J x ⟩/(N/2) ∼ |λ/λ c − 1| γx with γ x = 1/2 and the scaled mean number We have to mention that the above results in the thermodynamical limit are independent of the driving strength Ω, due to the fact that lim N →∞ (Ω 2 /N λ 2 ) → 0. This implies that, before reaching the critical point, the rise of the ⟨J z ⟩ curve due to the excitation by the driving Ω happens only in the finite-spin cases. It represents the polariton state in the normal phase, rather than in the superradiant phase. The polariton number decreases with N for λ < λ c , and is identically zero in the thermodynamic limit. Therefore, the QPT in the finite-spin case is identified only after traversing the critical point λ c .

Supplementary Note 2: Discussion about the involvement of A 2 term
With respect to a real light-matter interaction we have omitted an A 2 term (with A 2 = κ(a + a † ) 2 ) in our driven TC Hamiltonian H 1 . Actually, although κ is much smaller than other characteristic frequencies, it has been generally considered that A 2 term can forbid the QPT in the Dicke model due to the restriction from the Thomas-Reiche-Kuhn (TRK) sum rule, which is called a 'no-go' theorem [2][3][4]. We have also noticed some recent arguments in this topic [5]. The complete Hamiltonian should be where the last term is the A 2 term. In our case, this no-go theorem is circumvented in the non-equilibrium situation. In the frame rotating with ω d , Supplementary Equation 5 can be rewritten as where we have assumed identical qubits and identical couplings for each qubit to the resonator, as done above in Supplementary Equation 3. This treatment of homogeneity does not change the physical essence of the problem. From Supplementary Equation 6, we know that the extra small quadratic term in Supplementary Equation 5 can be reduced to a small shift in ∆ r in Supplementary Equation 6, following the rotating-wave approximation. This implies that a QPT is allowed under the drive in the non-equilibrium frame as the TRK sum rule is no longer violated. Our experimental observation of the QPT in ⟨J z ⟩ also confirms this conclusion.

Supplementary Note 3: Numerical simulation using master equation
We now consider the influence due to decoherence, which results from the decay/dephasing effects of the phase qubits and the superconducting resonator. In our case, the time evolution of the system can be numerically simulated by the Markovian master equation in which all the decoherence effects are phenomenologically described by the following Lindblad forṁ where κ 1 (κ 2 ) is the decay (pure dephasing) rate of the resonator, and Γ j 1 (Γ j 2 ) is the spontaneous emission (pure dephasing) rate of the j-th phase qubit. In the numerical simulation we assume identical decay for each qubit, with the energy relaxation time T q 1 ≈ 500 ns and the Ramsey dephasing time T q 2 ≈ 200 ns (see Supplementary Note 5 for choosing these parameters in numerical simulations), such that Γ j 1 = 1/T q 1 = 2 MHz and Γ j 2 = 1/T q 2 − 1/(2T q 1 ) = 4 MHz. For the resonator, we have κ 1 = 0.4 MHz and κ 2 → 0 MHz. Considering that the resonator is initially prepared in the vacuum state and the phase qubits are in ground states, we numerically solve Supplementary Equation 7 by the Quantum-Optics Toolbox, where λ/λ c uniformly varies from 0.5 to 2.5 using pulse durations of 600 ns and 1000 ns, respectively, and we consider the drives Ω/2π = 4 and 6 MHz.

Supplementary Note 4: Quantum vs semi-classical treatments
The QPT in the driven TC model was investigated previously using a semi-classical approach [6], and we have also noted some recently developed approaches solving the QPT semi-classically in the light-matter interaction [7][8][9]. Since noise is definitely involved in our experimental observation one may question whether one can describe the experimental data semi-classically.
We argue that there has been no semi-classical solution so far to the off-resonant driven TC model under our consideration. The above mentioned semi-classical solutions are restricted by some conditions and do not apply to our model. For example, only the resonant variant of the driven TC model can be solved and understood [6], and only the generic (non-driven) TC model can be treated semi-classically [7][8][9] since the parity symmetry is necessary in the solutions. In contrast, for our model with detunings and under the condition of parity breaking, it is not clear whether a semiclassical treatment can be found easily to be an appropriate description.
Nevertheless, our quantum treatment has well explained the experimental data with no fitted parameters, as presented in the main text. In contrast, the semi-classical treatment employing continuous variables which works well in the thermodynamic limit may not be applied to our case with only four qubits.

Sample parameters
Our superconducting circuit consists of four phase qubits interconnected by a central resonator. The phase qubit is encoded using a Josephson-junction device, whose cubic potential as a function of the junction phase δ is illustrated in the Supplementary Figure 1a (more about the phase qubit circuit can be found in [10]). The barrier height ∆U can be tuned to change the number of energy levels in the well. The cubic anharmonicity ensures that all transition frequencies ω n,n−1 = (E n − E n−1 )/ are distinct, and is the key for qubit operation, allowing microwaves at frequency ω q ≡ ω 10 = (E 1 − E 0 )/ to drive transitions between |0⟩ and |1⟩ while minimizing leakage to |2⟩ and higher excited states. In each qubit of our device, transition frequencies ω n,n−1 /2π are separated from each other by about 100 MHz.
Different from most previous experiments for gate demonstrations and state generations where a single operation point in frequency for each qubit is enough, here we gradually sweep the qubit frequency ω q (t)/2π over a span of more than 100 MHz for the case of ∆ r /2π = −30 MHz, such that the coherence performance over the entire frequency range matters. For example, Supplementary Figure 1b shows the two-dimensional (2D) swap spectroscopy measurement of Q 1 near the resonator frequency, from which we can fit the qubit T 1 (as in Supplementary Figure 1c) as a function of the bias frequency over a range relevant to the experiment. It is seen that T 1 fluctuates over frequency (Supplementary Figure 1d. T 1 can be significantly reduced, where a spurious two-level state (TLS) is present), and a histogram analysis suggests that the average qubit T 1 in this frequency range is 500 ns, with a standard deviation of 30 ns (Supplementary Figure 1e). Note that due to the strong qubit-resonator interaction at very small ∆ q , T 1 measured close to the resonator, i.e., nearby the chevron in Supplementary Figure 1b, may not be very accurate. Meanwhile, numerical simulation of the ⟨J z ⟩ versus λ/λ c curves based on varying T 1 values for each qubit reveal results very similar to those using a single averaged value for T 1 . Therefore we quote a single value of T 1 = 500 ns for simplified theoretical analysis throughout the work.
For the same reason, we quote an average Ramsey T 2 = 200 ns for all qubits in the analysis. The dephasing process in the phase qubit is non-Markovian due to the 1/f nature of flux noise. As our pulse sequence is typically a few hundred nanoseconds, a Ramsey T 2 of 200 ns should be appropriate: numerical simulation shows that small changes of T 2 would not significantly affect the rise-up magnitude and manner of the ⟨J z ⟩ curves across the critical point.  Figures 2c and 2d) have also caught the key feature of the off-resonantly driven QPT, i.e., a signature rise of the scaled moment ⟨J z ⟩ as λ/λ c increases above 1.
However, difference between the experiment data and numerical simulation at ∆ r /2π = −20 MHz is more apparent than that at ∆ r /2π = −30 MHz. On one hand, it is related to the fact that state leakage to higher excited states of the phase qubit is more serious at ∆ r /2π = −20 MHz as will be discussed next; on the other hand, for ∆ r /2π = −20 MHz, ∆ q has to sweep through a wider range in frequency to achieve the desired trajectory, more likely leading to errors in calculating λ/λ c (in particular for shorter τ values).

Designing the pulse sequence
The pulse sequence in Figure 1c of the main text is specifically designed considering various experimental limitations. Here we discuss a few important factors other than those already discussed in the main text.
Minimizing state leakage: the driven Tavis-Cummings theory (Equation 1 of the main text) assumes spin-1/2 particles (two-level qubits) interacting with the resonator field. Our phase qubit, though intentionally operated as a two-level system by restricting it to the two lowest energy eigenstates, contains multiple energy eigenstates in its cubic potential as illustrated in Supplementary Figure 1a. Leakage to the |2⟩ and higher excited states is possible during the pulse sequence due to the phase qubit's weak anharmonicity. It is thus important to minimize state leakage for validating the spin-1/2 assumption of the phase qubit.
Since transition frequencies ω n,n−1 to the |2⟩ and higher excited states of the qubits are smaller than ω q , we carry out the experiment by negative values of ∆ q and ∆ r , implying ω d > {ω q , ω r }, to avoid chances of accidentally driving these transitions. Due to symmetry, both positive and negative values of ∆ q,r lead to identical solutions of Supplementary Equation 3. So we may compare the experimental data (using negative detunings) with the ideal case (employing positive detunings) in Figure 2 of the main text.
State leakage can also be due to the qubit-resonator n r -photon Rabi swap, which happens in the subspace of |1, n r + 1⟩ and |2, n r ⟩ (here the first index refers to energy levels in the phase qubit and the second one is for the resonator). Coupling between the qubit's |1⟩-|2⟩ and the resonator's |n r ⟩-|n r + 1⟩ is approximately √ 2(n r + 1)λ as determined in previous experiment [11], and the photon swap probability is ≈ 2(n r +1)λ 2 / [ 2(n r + 1)λ 2 + (ω 21 − ω r ) 2 ] . Therefore state leakage can be significant if the average photon number ⟨n r ⟩ in the resonator is large and meanwhile the qubits and the resonator are close in frequency (the same condition as ∆ q being tuned small for λ/λ c approaching 2.5).
The resonator's photon population can be high under either a strong drive or a weak but continuous (long-time) drive using experimental ∆ r values. In order to minimize state leakage, we choose to keep the drive low and vary ∆ q with time, such that the system is subject to a weak drive at small detunings of ∆ q only for a limited amount of time. As numerically verified in Figure 2b of the main text, linearly sweeping λ/λ c over time offers a dynamical manner to probe the QPT, without exposing the system to external microwave drive at very small detunings of ∆ q for a long time.
Given the linear sweep of λ/λ c , we have to carefully choose values of ∆ r to further minimize state leakage by keeping the resonator's photon population low. For example, Supplementary Figure 3a shows the numerically-simulated average photon number ⟨n r ⟩ in the resonator as functions of the sweeping time t using the pulse sequence parameters of τ = 1000 ns and Ω/2π = 4 MHz. It is seen that ⟨n r ⟩ quickly accumulates with the sweep time t at ∆ r /2π = −10 MHz, which could incur strong qubit-resonator photon swap and should be avoided experimentally. Therefore we perform the experiment with ∆ r /2π = −20 MHz and −30 MHz.
With the negative ∆ q,r , the linear sweep of λ/λ c and the restriction of the driving strength Ω/2π ≤ 6 MHz, the leakage to the |2⟩ and higher excited states for each qubit can be restricted to, on average, within 3% for ∆ r /2π = −30 MHz, as experimentally monitored (Figures 3b and 3c in the main text. Data shown also reflect the readout uncertainty, ∼ 1%). State leakage is slightly higher for ∆ r /2π = −20 MHz, but in most experimental combinations of τ and Ω it is still less than 5% on average. For the small-percentage state leakage, Supplementary  Equation 3 involving two-level qubits is a qualified model to describe the experimental situation.
By taking the next higher excited state |2⟩ into consideration, we can model the experimental situation more accurately, as in the experiment we do not differentiate between P 1 and P n≥1 when calculating ⟨J z ⟩. By including three energy levels for the qubit circuit, i.e., a qutrit, we compare the numerically-calculated (the thick black line) leakage-state probability P L (= P 2 in the qutrit model) with the experimentally obtained values (thin lines) for the case of ∆ r /2π = −30 MHz, Ω/2π = 4 MHz and τ = 600 ns (Supplementary Figure 4a); in Supplementary Figure 4b, we obtain a fairly good agreement between numerical simulation (the dashed line) and the measured ⟨J z ⟩ curve.
Accuracy in tuning and synchronizing all four qubits' trajectory ω q (t): in the experiment to vary λ/λ c linearly from 0.5 to 2.5 within the period of τ , i.e., we adjust the bias of each qubit such that the detuning of each qubit changes according to Note that we cannot use the spectroscopy method to determine the qubit frequency during the steering of ω q (t) as the experimental process happens quickly in short times from a few to a few hundreds nanoseconds. Instead, we use a Ramsey-interference-like measurement to verify the accuracy in tuning ω q (t) by checking the dynamic phase accumulated during the process. Starting with the qubit frequency at ω q (0), we prepare the initial state (|0⟩−i|1⟩)/ √ 2, and then tune the qubit frequency according to Supplementary Equation 9 for a time t, following which we apply a π/2 pulse with an adjusted phase If the steering of ω q is as expected, θ compensates the accumulated dynamic phase and we end up with measuring the envelope of the Ramsey interference fringe. The process is equivalent to the Ramsey interference measurement using two identical π/2 rotation pulses (Supplementary Figure 5a). We note that this method has its limitations. First, due to the relatively short T 1 and T 2 of phase qubits, the probability of finding the qubit in |1⟩, P 1 , gradually decays as t increases, which renders the signal rather weak at larger t values. Second, this method cannot work when the qubit is tuned close to the resonator frequency, due to the Stark shift. As such, we test the accuracy of tuning the qubit frequency in a frequency range from 5.87 GHz to 6.14 GHz in 75 ns. The results for all four qubits are as expected (Supplementary Figure 5b).
Flux bias crosstalk from one qubit to the other is intended to be corrected during the steering of ω q (t). However, we have no conclusive evidence that our synchronized steering of all four qubits is accurate and adequate, nor can we be absolutely sure about the accuracy of correcting the flux bias crosstalk. We estimate an uncertainty of δ [ω q (t)/2π] ≈ 1 MHz, based on all available experimental data and our knowledge of the control accuracy. In particular, the fact that most experimental data measured at the microwave drive strength Ω/2π = 2 MHz show much lower ⟨J z ⟩ signals than expected is consistent with the estimated uncertainty level of ω q (t). For this reason we choose intermediate microwave drive strengths Ω/2π = 4 MHz and 6 MHz in the experiment.
Fortunately, the experimental impact of uncertainties in ω q (t) can be greatly reduced by using larger negative detunings of ∆ r (for ∆ q (t) to sweep through smaller ranges in frequency, which also lowers the crosstalk impact) and longer sweep durations of τ . Using larger negative detunings of ∆ r also helps with minimizing state leakage. Therefore the experimental data agree quite well with numerical simulation that assumes perfect synchronized steering of ω q (t) for most experimental combinations of ∆ r and τ . A slightly large disagreement is seen in the case of ∆ r /2π = −20 MHz and τ = 600 ns, which is compatible with our hypothesis of errors in synchronized tuning of ω q (t).
Two-qubit behaviour and the off-resonant Ω: as a preparation for the four-qubit ⟨J z ⟩ measurements, we quickly perform the two-qubit measurements with Q 1 -Q 2 and Q 1 -Q 3 , respectively, using pulse sequences that are similar to that shown in Figure 1c of the main text for the resonator and two of the qubits. It is observed that the two-qubit ⟨J z ⟩ curves also show the signature rise when λ/λ c increases above 1. The two-qubit ⟨J z ⟩ measurements also confirm that Ω calibrated on-resonantly is applicable off-resonantly for small detunings of ∆ r (Supplementary Figure 6).