Quantum Zeno effect in the strong measurement regime of circuit quantum electrodynamics

The backaction of measurement is a peculiarly quantum mechanical phenomenon which gives rise to striking outcomes, such as the quantum Zeno effect (QZE). In the QZE, the act of measurement inhibits transitions between eigenstates of the measured observable, slowing the state evolution of a 'watched' quantum system. The QZE was described in its modern form in 1977 by Misra and Sudarshan [1], although some related questions were tackled in prior papers [2, 3], and the potential for repeated measurements to influence the state evolution of a quantum system was already known to von Neumann in 1932 [4]. The slowing of state evolution due to the QZE disappears in the classical limit ${\rm{\hslash }}\to 0$, making the QZE a useful test for quantum behavior in a system [5, 6].

The QZE was first observed experimentally in an ensemble of trapped ions [7], and has since been seen in a variety of other systems, including the electronic, nuclear, or motional states of atoms and molecules [8–11], optical photons [12–14], microwave photons [15], and NV centers [16]. In driven superconducting qubits, the QZE has been indirectly inferred from the transition between coherent Rabi oscillations and incoherent exponential population decay with increasing measurement strength [17], and by studying the dependence of this exponential decay on the time between discrete qubit projection pulses [18]. However, theoretical proposals also exist to observe the QZE in the quantum trajectory of a continuously monitored superconducting qubit [19] or in the suppression by measurement of low-frequency superconducting flux qubit dephasing [20].

The QZE and related phenomena can have practical applications in quantum control and the engineering of decoherence. Carefully designed measurements can be used to divide a larger Hilbert space into separate 'Zeno subspaces' [21, 22], where state evolution between subspaces is inhibited by the measurements. The resulting 'quantum Zeno dynamics' have recently been demonstrated experimentally [23–25]. This technique can even be used to generate multiparticle entanglement directly [26–29].

In this work, we report the direct observation of the QZE in a superconducting qubit undergoing continuous strong measurement with simultaneous qubit driving. The measurement record is analyzed to extract quantum jumps indicating individual qubit state transitions, and to determine the rates at which they occur in the presence of simultaneous qubit drive and measurement. The extracted transition rates show inhibition of qubit state transitions due to the measurement, in agreement with both analytical Zeno theory and numerical simulations. Additionally, we examine the excited state decay of the qubit during measurement, finding qualitative agreement with predictions for the Purcell decay of a multi-level qubit [30, 31]. The method for extracting transition rates can also be used to analyze general noisy random telegraph signals.

In the absence of measurement, a resonantly driven qubit will undergo sinusoidal Rabi oscillations between states at frequency ${\rm{\Omega }}/2\pi$, where Ω depends on the strength of the resonant drive. Repeated projective measurements made on this system at time intervals $\tau \ll 1/{\rm{\Omega }}$ will tend to pin the qubit in one state, with occasional state transitions occurring as sudden quantum jumps. The probability per unit time of a quantum jump out of the current state is given by (see e.g. [32]):

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{{{\rm{\Omega }}}^{2}}{4{f}_{{\rm{m}}}},\end{eqnarray} \tag{ 1 }$$

where ${f}_{{\rm{m}}}=1/\tau$ is the frequency of the measurements. A similar result can be derived in the case of a resonantly driven qubit undergoing continuous measurement, parameterized by the measurement rate ${{\rm{\Gamma }}}_{{\rm{m}}}$ at which information is extracted from the qubit [33]:

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{{{\rm{\Omega }}}^{2}}{{{\rm{\Gamma }}}_{{\rm{m}}}}.\end{eqnarray} \tag{ 2 }$$

These expressions show that the rate of quantum jumps between states for a qubit being measured continuously at rate ${{\rm{\Gamma }}}_{{\rm{m}}}$ is the same as that for a qubit undergoing discrete measurements at time intervals τ when ${{\rm{\Gamma }}}_{{\rm{m}}}\tau =4$. This prediction, due to Schulman [33], was first demonstrated experimentally by Streed et al [10].

The two major hallmarks of the QZE can be seen in the expressions (1) and (2). First, Γ is independent of time, meaning that the state evolution under repeated measurement is exponential (linear in time at short times), whereas in the absence of measurement the qubit would exhibit sinusoidal state evolution (quadratic in time at short times). Qubit evolution processes which were already exponential in time, for example T₁ decay, remain exponential and are unaffected by the presence of repeated measurements. Second, the transition rate varies inversely with the measurement frequency, tending to zero in the limit of infinitely frequent measurements. However, the complete 'freezing' of qubit evolution does not occur in real physical systems. This is due to the energy-time uncertainty relationship, which causes the qubit to couple to an arbitrarily large energy spectrum (and thus arbitrarily many decay channels) as the time between measurements goes to zero [34, 35].

Our experiment uses a superconducting qubit in the circuit QED architecture [36]. We bias the qubit in the dispersive regime, where the qubit–cavity detuning Δ is much larger than the qubit–cavity coupling g. In this limit, the cavity resonance frequency depends on the qubit state, enabling a quantum non-demolition measurement of the qubit state by driving the cavity near resonance and observing the response of its steady-state field. The measurement strength is directly related to the distinguishability of the different cavity field amplitudes conditioned on the qubit states. For a linear cavity, the distinguishability is proportional to the amplitude of the measurement drive and could in principle be increased ad infinitum. In real systems, however, cavity nonlinearities (both intrinsic and induced by coupling to the anharmonic qubit circuit) limit the obtainable field separations for readout, even at moderate drive amplitudes [37].

Following [19], and as shown in appendix A, we derive the transition rate of a two-level qubit coupled to a cavity in the Zeno regime starting with the Jaynes–Cummings Hamiltonian in the dispersive regime with an additional qubit drive:

$$\begin{eqnarray}&&\hat{H}={\omega }_{{\rm{r}}}{a}^{\dagger }a+\displaystyle \frac{1}{2}({\omega }_{{\rm{q}}}+\chi ){\sigma }_{z}+\chi {a}^{\dagger }a{\sigma }_{z}+\displaystyle \frac{1}{2}({\epsilon }_{\mathrm{ro}}{a}^{\dagger }{{\rm{e}}}^{{\rm{i}}{\omega }_{\mathrm{ro}}t}+{\epsilon }_{\mathrm{ro}}^{*}a{{\rm{e}}}^{-{\rm{i}}{\omega }_{\mathrm{ro}}t})\\&&+\displaystyle \frac{1}{2}({\rm{\Omega }}{\sigma }_{-}{{\rm{e}}}^{-{\rm{i}}{\omega }_{{\rm{d}}}t}+{{\rm{\Omega }}}^{*}{\sigma }_{+}{{\rm{e}}}^{{\rm{i}}{\omega }_{{\rm{d}}}t}).\end{eqnarray} \tag{ 3 }$$

Here, $a\ ({a}^{\dagger })$ is the annihilation (creation) operator for photons in the cavity, the qubit is described by the Pauli matrices σ, and we have made a rotating wave approximation in the two driving terms. The parameters are the cavity resonance frequency ${\omega }_{{\rm{r}}}$, the qubit frequency ${\omega }_{{\rm{q}}}$, the qubit–cavity coupling g, and the dispersive shift $\chi ={g}^{2}/{\rm{\Delta }}$, where ${\rm{\Delta }}={\omega }_{{\rm{r}}}-{\omega }_{{\rm{q}}}$. We apply a measurement drive to the cavity with frequency ${\omega }_{\mathrm{ro}}$ and strength ${\epsilon }_{\mathrm{ro}}$ and a Rabi drive to the qubit with frequency ${\omega }_{{\rm{d}}}$ and amplitude Ω. We then perform a polaron transformation on the system, which effectively describes the cavity in a semi-classical picture with coherent state amplitudes that are conditioned on the qubit state. This transformation enables us to decouple the dynamical equations of the qubit and cavity and obtain a qubit-only reduced master equation, where the qubit dynamics depend mainly on qubit drive and measurement strength. Because of the presence of the Rabi drive, this treatment is only valid in the case of small measurement strength, when the cavity field amplitudes corresponding to the two qubit states are nearly indistinguishable. Solving this master equation yields a qubit transition rate from the ground state to excited state in the presence of continuous circuit QED measurement [19]:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\uparrow ,\mathrm{drive}}=\displaystyle \frac{{{\rm{\Omega }}}^{2}}{2({\gamma }_{2}+{{\rm{\Gamma }}}_{{\rm{d}}})},\end{eqnarray} \tag{ 4 }$$

where ${{\rm{\Gamma }}}_{{\rm{d}}}$ is the measurement-induced dephasing rate and ${\gamma }_{2}$ is the intrinsic qubit dephasing rate in the absence of measurement. The measurement rate is defined as ${{\rm{\Gamma }}}_{{\rm{m}}}=2{{\rm{\Gamma }}}_{{\rm{d}}}$, as in [19]. In the strong measurement limit where ${{\rm{\Gamma }}}_{{\rm{d}}}\gg {\gamma }_{2}$, this expression is the same as equation (2). Although equation (4) is derived in the limit where the cavity linewidth $\kappa \gg {{\rm{\Gamma }}}_{{\rm{m}}}$, the similarity to equation (2)—which has no such constraint in its derivation—suggests that the result may still be valid when $\kappa \lt {{\rm{\Gamma }}}_{{\rm{m}}}$, which is the regime of our experimental data.

The experimental setup is shown schematically in figure 1(a). Our qubit is a planar two-junction transmon qubit [38], tuned with a dc magnetic flux to a Lamb-shifted qubit frequency of ${\tilde{\omega }}_{{\rm{q}}}/2\pi =5.3556\;{\rm{GHz}}$ and an anharmonicity of $\alpha /2\pi =-258\;{\rm{MHz}}$. We denote the two lowest transmon energy levels as $| g\rangle$ and $| e\rangle$, respectively. The qubit is capacitively coupled ($g/2\pi =105.3$ MHz) to a lumped-element readout cavity with bare frequency of ${\omega }_{{\rm{r}}}/2\pi =6.2724\;{\rm{GHz}}$ and linewidth $\kappa /2\pi =7\;{\rm{MHz}}$. Qubit measurement was performed by applying a readout drive to the cavity at ${\omega }_{\mathrm{ro}}/2\pi =6.282\;{\rm{GHz}}$. The readout signal from the cavity was amplified with a near-quantum-limited superconducting parametric amplifier [39] followed by additional cryogenic and room temperature amplification stages, detected with homodyne mixing, and digitized at 10 ns intervals. Both the qubit/cavity system and the parametric amplifier were anchored to the mixing chamber of a dilution refrigerator at 50 mK.

Zoom In Zoom Out Reset image size

The average cavity photon occupation $\bar{n}$ and the measurement-induced dephasing rate ${{\rm{\Gamma }}}_{{\rm{d}}}$, which characterize the measurement strength, were determined by qubit spectroscopy with a simultaneous readout tone. The value of $\bar{n}$ was determined from the ac Stark shift of the qubit frequency, accounting for qubit-induced cavity nonlinearities [30]. Because of the choice of readout frequency, $\bar{n}$ is different for the two qubit states $| g\rangle$ and $| e\rangle ;$ the ac Stark shift measurements give the average photon occupation with the qubit in the ground state, ${\bar{n}}_{{\rm{g}}}$, since the qubit spectroscopy tone was much weaker than the saturation amplitude of the qubit transition. We use numerical simulations to infer the value of ${\bar{n}}_{{\rm{e}}}$ as a function of the measured ${\bar{n}}_{{\rm{g}}}$. The value of ${{\rm{\Gamma }}}_{{\rm{d}}}$ is given by the half-width at half-maximum (HWHM) of the ac-Stark-shifted qubit line [40], which was determined by fitting to the measured line shape for several qubit drive powers and extrapolating to find the HWHM at zero qubit drive power. This method allows us to calibrate out the effects of power broadening due to the qubit drive tone on the measured qubit linewidth. With the readout tone off during qubit evolution, we measured a qubit relaxation time of ${T}_{1}=575$ ns and a pure dephasing time of ${T}_{\varphi }=7.9\;\mu {\rm{s}}$, corresponding to an intrinsic dephasing rate ${\gamma }_{2}=1/2{T}_{1}+1/{T}_{\varphi }=1\;\mu {{\rm{s}}}^{-1}$.

The qubit drive frequency ${\omega }_{{\rm{d}}}$ was chosen to be the ac-Stark-shifted qubit frequency ${\tilde{\omega }}_{{\rm{q}}}({\bar{n}}_{{\rm{g}}})$. Since the qubit drive tone reaches the qubit via the cavity, its amplitude at the qubit (and thus the no-measurement Rabi frequency Ω) depends on both its frequency and its amplitude at the cavity input. We calibrated Ω by measuring the frequency of Rabi oscillations (performed with measurement off during the qubit drive) as a function of qubit drive amplitude for a given qubit frequency. We performed this calibration at several different qubit frequencies spanning the range of the ${\omega }_{{\rm{d}}}$ values used in the Zeno experiments by tuning the qubit frequency with an external flux bias. We interpolated between these data points to determine the value of Ω for a given ${\omega }_{{\rm{d}}}$ and qubit drive amplitude.

Each iteration of the experiment consisted of a single measurement pulse lasting 17.5 μs, with simultaneous qubit drive applied for 10 μs beginning 3 μs after the measurement was turned on, as shown in figure 1(b). These times are chosen to be much longer than the relevant timescales $1/\kappa ,1/{{\rm{\Gamma }}}_{{\rm{m}}},1/{\rm{\Omega }}$, and T₁, so that the experimental record captures the steady-state dynamics of the system. We recorded 10⁴ such iterations, spaced 100 μs apart, for each combination of measurement strength and qubit drive strength. The measurement record was digitized at 10⁸ samples s⁻¹.

The measurement records described above were analyzed to extract the transition rates from $| g\rangle \to | e\rangle$, denoted ${{\rm{\Gamma }}}_{\uparrow }$, and from $| e\rangle \to | g\rangle$, denoted ${{\rm{\Gamma }}}_{\downarrow }$. With multiple gigabytes of raw quantum jump data to analyze, the data processing algorithm must be able to operate with minimal user input and provide reliable output over a broad range of signal-to-noise ratio (SNR) and qubit transition rates, including low-SNR, high-rate scenarios.

We first perform filtering of the raw data to reduce the noise bandwidth and increase the SNR. We use a zero-delay Gaussian finite-impulse-response filter to ensure that state transition edges will remain smooth and will not be shifted in time by the filtering process. The order of the filter must be chosen to reduce the noise without degrading the desired signal, and so will depend both on the initial SNR and the transition rate between states. For each experimental bias point, we histogram all 10⁷ data points taken during simultaneous qubit driving and readout (10⁴ traces of 10³ points each). Example histograms of raw data are shown as the red curves in figures 2(a) and (c). We then perform filtering and histogram the filtered data points, shown as black curves in figures 2(a) and (c). The histogram of the filtered data forms a bimodal distribution⁸ , with a minimum between the two peaks corresponding to the two qubit states. The height of this minimum is dictated both by the amount of noise on the signal and the number of transitions in the measurement record. Filtering will reduce the noise amplitude, and thus the width and overlap of the bimodal peaks, leading to a reduction in the height of the minimum. However, filtering will also slow the rise and fall times of state transition edges, which will increase the height of the minimum between histogram peaks. Based on these considerations, we designate the 'optimal' filter order as the one which gives the lowest minimum between the bimodal peaks of the resulting voltage histogram, and determine it by applying filters of increasing order to the raw data until the height of the minimum between the peaks is as small as possible. This provides a simple, robust method for choosing filter order. The 'optimal' filter order was determined separately for each combination of qubit drive strength and measurement strength.

Zoom In Zoom Out Reset image size

After the data have been filtered, we determine the qubit state at each time point using a thresholding algorithm. For increased robustness to noise, we use two hysteretic thresholds with 'Schmitt trigger' behavior [41]; state transitions are registered when the voltage crosses the higher threshold going upwards (if the state was low), or when it crosses the lower threshold going downwards (if the state was high). The voltage thresholds are shown as dotted horizontal lines in figure 2. Figures 2(b) and (d) show representative individual data traces from the data sets histogrammed in figures 2(a) and (c), respectively, corresponding to two different values of Ω with identical ${{\rm{\Gamma }}}_{{\rm{m}}}$. The red curves show the raw data, while the black curves are after filtering. The blue and white background colors show the extracted qubit state ($| e\rangle$ and $| g\rangle$, respectively) as determined by the thresholding algorithm. The data in (c) and (d) have fewer state transitions, and thus heavier filtering can be employed to increase the SNR.

The location of the thresholds is determined using histograms of the filtered voltages, and depends on the voltage V_m of the minimum between the histogram peaks and the voltages V_h (V_l) and the HWHM w_h (w_l) of the high (low) histogram peaks. The value of w_h (w_l) is derived from the half-height of the histogram above (below) the corresponding peak V_h (V_l), a method that works even when the histogram peaks are not well-separated. We define the SNR of the filtered data to be $\sqrt{2\;\mathrm{ln}\;2}({V}_{{\rm{h}}}-{V}_{{\rm{l}}})/({w}_{{\rm{h}}}+{w}_{{\rm{l}}})$. The thresholds are chosen to be ${V}_{{\rm{m}}}+\tfrac{{w}_{{\rm{h}}}^{2}}{2\;\mathrm{ln}\;2({V}_{{\rm{h}}}-{V}_{{\rm{l}}})}$ and ${V}_{{\rm{m}}}-\tfrac{{w}_{{\rm{l}}}^{2}}{2\mathrm{ln}2({V}_{{\rm{h}}}-{V}_{{\rm{l}}})}$, as suggested in [41], although unlike that work we do not perform iterative refinement of the thresholds. The distance between the thresholds increases as the SNR decreases to reduce the likelihood of spurious state transitions being registered due to noise. For each measurement trace, we determine the initial state by comparing the first data point with V_m. Given a low (high) state, we look for the first upward crossing of the higher threshold (downward crossing of the lower threshold) and note the dwell time before the state change. We continue in this manner, alternating thresholds, until the end of the trace is reached, noting down all dwell times in both the high and low states. For each measurement trace, exactly one dwell time is cut short by the end of the trace and is marked as 'right-censored' [42].

We use maximum likelihood estimation to determine the transition rates between states given the set of observed dwell times in each state. However, the filtering can cause some fast voltage excursions in the raw data, like the one in figure 2(d) near 8.25 μs, not to register as state transitions. In general, the finite bandwidth of the measurement chain, combined with filtering described above, skews the distribution of observed dwell times toward longer times. We compensate for this in our analysis by assuming a probability distribution for the dwell times which takes these effects into account [43]. For right-censored dwell times, the observed time represents a lower bound on the dwell time, rather than an exact value. However, we can include this partial information in the likelihood function as well, which is particularly useful for data sets with relatively few state transitions and thus a higher fraction of censored dwell times. We emphasize that this maximum likelihood method can be applied to dwell times extracted using any technique; for example, it can be used with dwell times determined from wavelet analysis, which is more robust than thresholding for signals with substantial low-frequency drifts [44]. Details on the maximum likelihood estimation and functional forms of the likelihood function are provided in appendix D.

The extraction algorithm was tested on simulated noisy random telegraph data with a variety of experimentally relevant transition rates and SNRs. We note that these simulated data are not the same as the numerical simulations of the qubit dynamics described in the next section. The 3 dB noise bandwidth (14 MHz) and sample rate (10⁸ sample s⁻¹) were the same as in the raw experimental data. The error in the extracted mean transition rate $({{\rm{\Gamma }}}_{\uparrow }+{{\rm{\Gamma }}}_{\downarrow })/2$, relative to the true mean transition rate, is plotted in figure 3 for four values of the initial SNR before filtering. For simulated data at the highest rates shown, the mean dwell time in a given state was 130 ns, while the shortest oscillation timescale of the simulated noise was around 70 ns, making the task of distinguishing signal from noise challenging. Despite this, the mean systematic bias in the extracted rates is below 12%, even at an SNR of only 1.5. This performance is also notable because the algorithm operates autonomously, without the need for externally provided guesses or input, across the entire range of transition rates and SNR shown.

Zoom In Zoom Out Reset image size

The transition rates extracted from the simulated noisy random telegraph data were used to calibrate systematic bias and systematic uncertainty in the rates extracted from the experimental data, by matching experimental traces to simulated data sets with similar extracted transition rates and filtered SNR. The median systematic bias and median systematic uncertainty for the experimental data were −3.5% and 2% of the extracted rates, respectively. The magnitude of the systematic bias (uncertainty) was below 20% (10%) of the extracted rate for all experimental data points, and below 9% (5%) of the extracted rate for 90% of the experimental data points.

We also performed numerical simulations of the qubit/cavity system over a range of bias points, using the generalized Jaynes–Cummings Model for a multi-level qubit coupled to a cavity. We simulated the dynamics numerically, including the measurement, with a stochastic master equation [19, 45, 46], using experimentally measured parameters for the qubit and cavity. The master equation included qubit relaxation and dephasing as well as cavity photon loss. The measurement and qubit drive were simulated as two independent coherent microwave drives acting on the cavity. The simulations yield measurement records for direct comparison with experiment, as well as records of qubit populations, complex cavity field amplitudes, and cavity photon occupation numbers. Details are given in appendix B. Qubit transition rates were then determined from the simulated qubit state population record using the same rate extraction algorithm used to process the experimental data.

The experimental data, and thus the rate extraction algorithm, do not distinguish between the qubit state $| e\rangle$ and higher excited states. This occurs because we operate the parametric amplifier in phase-sensitive mode to achieve the lowest noise performance, and thus only one quadrature of the readout signal from the cavity is amplified [39]. The phase of the amplified quadrature was chosen to maximize the ability to discriminate between states $| g\rangle$ and $| e\rangle ;$ however, because of the parameters of the qubit/cavity system, the projection of the complex cavity amplitudes corresponding to higher excited states onto this choice of amplified quadrature was essentially the same as that for $| e\rangle$, rendering them indistinguishable from $| e\rangle$ in the measurement record [47]. The populations of higher excited states were measured to be $\sim 1\%-5\%$ using a different data set which distinguished $| e\rangle$ from higher qubit states by choosing a different amplified quadrature, at the expense of substantially decreased SNR. These values are corroborated by the numerical simulations. However, the presence of qubit population outside the $\{| g\rangle ,| e\rangle \}$ manifold is not expected to have an effect on the estimates of ${{\rm{\Gamma }}}_{\uparrow }$.

We first examine the qubit excitation rate during simultaneous measurement and qubit driving to look for the QZE. Since ${\bar{n}}_{{\rm{g}}}\ne {\bar{n}}_{{\rm{e}}}$ for our experimental parameters, the qubit drive tone at the ac-Stark-shifted qubit frequency ${\tilde{\omega }}_{{\rm{q}}}({\bar{n}}_{{\rm{g}}})$ was resonant when the qubit was in the ground state, but not when it was in the excited state. For this reason, our analysis of driven transition rates is restricted to the $| g\rangle \to | e\rangle$ transition rate, denoted ${{\rm{\Gamma }}}_{\uparrow }$, which can be expressed as the sum of three rates:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\uparrow }={{\rm{\Gamma }}}_{\uparrow ,\mathrm{drive}}+{{\rm{\Gamma }}}_{\uparrow ,\mathrm{DD}}+{{\rm{\Gamma }}}_{\uparrow ,\mathrm{th}}.\end{eqnarray} \tag{ 5 }$$

Here ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{drive}}$ is the transition rate due to the qubit drive, ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{DD}}$ is the contribution from dressed dephasing [37, 48], and ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{th}}$ represents thermal excitation of the qubit. To isolate ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{drive}}$ for comparison to the predictions in equation (4), we determine ${{\rm{\Gamma }}}_{\uparrow }$ as a function of measurement strength in the absence of qubit drive (${\rm{\Omega }}=0$), which is equal to ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{DD}}+{{\rm{\Gamma }}}_{\uparrow ,\mathrm{th}}$. We can then subtract this contribution from the total rate ${{\rm{\Gamma }}}_{\uparrow }$ in the presence of qubit drive to find ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{drive}}$. This calibration was performed for both experimental data and numerically simulated data. The value of ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{DD}}+{{\rm{\Gamma }}}_{\uparrow ,\mathrm{th}}$ was measured to be between 0.018 and 0.022 $\mu {{\rm{s}}}^{-1}$ for all bias points, which is considerably smaller than ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{drive}}$ for all but the lowest values of Ω.

Figures 4 and 5 show the extracted ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{drive}}$ as a function of Ω for ten different measurement strengths ranging from ${{\rm{\Gamma }}}_{{\rm{m}}}=134$ $\mu {{\rm{s}}}^{-1}$ to ${{\rm{\Gamma }}}_{{\rm{m}}}=393$ $\mu {{\rm{s}}}^{-1}$, corresponding to values of ${\bar{n}}_{{\rm{g}}}$ between 3.4 and 37. The rates from the experimental data and the numerical simulations are plotted as red squares and black circles, respectively, while the rates predicted by the analytical theory in equation (4), with no adjustable parameters, are plotted as solid blue lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The experimental data and numerics are in good agreement, and both coincide with the prediction of the analytical theory at most points, showing the presence of the QZE across a broad range of measurement strengths and qubit drive amplitudes. Some deviations can be seen above ${\rm{\Omega }}/2\pi \approx 6\;{\rm{MHz}}$, where both the experiment and numerics give consistently lower values of ${{\rm{\Gamma }}}_{\uparrow ,\mathrm{drive}}$ than predicted by the analytical theory. Because the numerical simulations exhibit the same behavior as the experiment, we believe this is a real effect and not an experimental artifact. We postulate that this additional slowing of the transition rate may be due to higher order terms neglected in the derivation of the analytical theory under the approximation ${{\rm{\Gamma }}}_{{\rm{m}}}\ll \kappa$, since this approximation is no longer valid here. The effect is somewhat more pronounced for larger ${{\rm{\Gamma }}}_{{\rm{m}}}$. This could also be associated with the breakdown of the dispersive approximation for ${\bar{n}}_{{\rm{g}}}\gt {n}_{\mathrm{crit}}={{\rm{\Delta }}}^{2}/4{g}^{2}$, which corresponds in our system to ${{\rm{\Gamma }}}_{{\rm{m}}}\approx 340$ $\mu {{\rm{s}}}^{-1}$.

We additionally examine the qubit T₁— the time constant for transitions from $| e\rangle$ to $| g\rangle$—with measurement on but no qubit drive. In the simple two-level picture used to derive equation (4), transition rates from exponential processes such as T₁ decay are unchanged by the QZE. However, in our driven multilevel system other effects are expected to give some dependence of T₁ on measurement strength. In circuit QED, the T₁ decay rate can in general be parameterized as the sum of the Purcell decay rate ${{\rm{\Gamma }}}_{{\rm{P}}}$ and a decay rate ${{\rm{\Gamma }}}_{\mathrm{NR}}$ from nonradiative loss channels such as dielectric loss [49, 50]. The Purcell effect can be thought of as the decay of the photonic part of the qubit eigenstate in the coupled qubit–cavity system, as seen in the expression for the Purcell rate [30, 31]

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{P}}}=\kappa {| \bar{\langle g,n| }a\bar{| e,n\rangle }| }^{2},\end{eqnarray} \tag{ 6 }$$

where $\bar{| g,n\rangle }$ and $\bar{| e,n\rangle }$ are eigenstates of the generalized Jaynes–Cummings Hamiltonian for a multilevel qubit corresponding most closely to n photons in the bare cavity and the bare qubit in $| g\rangle$ or $| e\rangle$, respectively. The matrix elements in equation (6) can be calculated numerically, and become smaller for increasing n. As a result, the Purcell decay of the qubit is suppressed by the presence of photons in the readout cavity [30, 31]. Note that in our case we replace n with ${\bar{n}}_{{\rm{e}}}$, because Purcell decay from $| e\rangle$ occurs with ${\bar{n}}_{{\rm{e}}}$ photons in the cavity initially. Details of the calculation are given in appendix C.

Figure 6 shows the qubit decay rate as a function of ${\bar{n}}_{{\rm{e}}}$. The expected decay rate ${{\rm{\Gamma }}}_{{\rm{P}}}$ due to the Purcell effect, as calculated from equation (6), is plotted as a red solid line. The measured T₁ decay of the qubit (red squares), which shows some suppression with increasing ${\bar{n}}_{{\rm{e}}}$, is larger than ${{\rm{\Gamma }}}_{{\rm{P}}}$, indicating that non-radiative loss channels are contributing to the qubit decay (i.e. ${{\rm{\Gamma }}}_{\mathrm{NR}}\gt 0$). We derive a value for ${{\rm{\Gamma }}}_{\mathrm{NR}}$, which we assume to be independent of ${\bar{n}}_{{\rm{e}}}$, by subtracting the zero-photon Purcell decay rate $\kappa {g}^{2}/{{\rm{\Delta }}}^{2}$ from the experimentally determined value of ${\gamma }_{1}=1/{T}_{1}$ for ${\bar{n}}_{e}=0$. The theoretically predicted total decay rate ${{\rm{\Gamma }}}_{{\rm{P}}}+{{\rm{\Gamma }}}_{\mathrm{NR}}$ is plotted as a red dotted line. This prediction is reasonably close to the experimental data for small ${\bar{n}}_{{\rm{e}}}$, but deviates substantially with increasing ${\bar{n}}_{{\rm{e}}}$. To determine whether this is an experimental artifact (due to readout-power-dependent dielectric loss in the transmon capacitor [51], for example, which would give an ${\bar{n}}_{{\rm{e}}}$ dependence to ${{\rm{\Gamma }}}_{\mathrm{NR}}$), we also performed numerical simulations of the qubit decay with no Rabi drive. The resulting rates (black circles and black triangles) are also consistently lower than the theoretical prediction for ${{\rm{\Gamma }}}_{{\rm{P}}}+{{\rm{\Gamma }}}_{\mathrm{NR}}$, and are in reasonable agreement with the experimental data. Because the numerical calculations assume that ${{\rm{\Gamma }}}_{\mathrm{NR}}$ is independent of ${\bar{n}}_{{\rm{e}}}$, this provides confirmation that the suppression in T₁ decay is not an experimental artifact. The numerics also allow us to determine the effect of qubit population in states higher than $| e\rangle$ on the observed decay rate; the black circles show the observed T₁ decay rate when $| e\rangle$ is not distinguished from higher qubit states, as is the case with our experimental data, while the black triangles show the true T₁ decay rate from $| e\rangle \to | g\rangle$ when $| e\rangle$ is distinguished from higher qubit states. The difference between these traces becomes significant for ${\bar{n}}_{{\rm{e}}}\gtrsim 8$, and shows that the inability to distinguish higher qubit states from $| e\rangle$ in this system causes the qubit decay rate to appear lower than its true value, but the effect does not appear to be large enough to fully account for the discrepancy between the experimental data and the analytical theory. We believe that other effects, such as dressed dephasing [37, 48] and qubit-induced nonlinearities, are responsible for this deviation from the simple theory in equation (6). The blue circles show how the qubit T₁ decay rate would depend on ${\bar{n}}_{{\rm{e}}}$ if it were subject to the QZE and scaled as $1/{{\rm{\Gamma }}}_{{\rm{m}}}$, as in equation (4), relative to the measured value for the lowest ${\bar{n}}_{{\rm{e}}}$. The large disparity between these points and the experimental and numerical data indicate that the variation of ${{\rm{\Gamma }}}_{\downarrow }$ with ${\bar{n}}_{{\rm{e}}}$ is not due to the QZE; this agrees with the expected result that the QZE does not affect exponential processes such as T₁ decay.

Zoom In Zoom Out Reset image size

We have observed the QZE in a continuously measured superconducting qubit and demonstrated quantitative agreement with both analytical theory and numerical simulations. Interestingly, the agreement holds even for measurement strengths much larger than are allowed by the assumptions of the analytical derivation. We have shown that T₁ decay is not subject to the QZE, as predicted. We do observe an gradual increase in the qubit T₁ with increasing measurement strength, in qualitative agreement with the simple picture from the Purcell effect in circuit QED. We have also demonstrated a robust algorithm for determining transition rates from our measurement record, which can be readily applied to the analysis of general noisy random telegraph signals.

We thank J Gambetta, O Naaman, J Aumentado, and Y C Hagar for useful discussions and K Lalumière for help with the stochastic master equation simulations. DHS acknowledges support from a Hertz Foundation Fellowship endowed by Big George Ventures. RV and IS acknowledge support from the ARO QCT program. AB acknowledges support from NSERC. We thank Calcul Québec and Compute Canada for computational resources. DHS and RV designed and carried out the experiments. SJW fabricated the qubit sample. DHS analyzed the data and developed the rate extraction algorithm. CM performed the theoretical calculations and numerical simulations. AB and IS supervised the work. DHS and CM wrote the manuscript with input from all authors.

Here we derive the expression for qubit transition rates given in equation (4). We will focus only on the steps necessary to derive the qubit transition rates; more detail can be found in [19]. We note that this derivation is for the case of a two-level qubit, not a multi-level qubit such as a transmon. We describe the qubit–cavity system in the Jaynes–Cummings model (setting ${\rm{\hslash }}=1$):

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{JC}}={\omega }_{{\rm{r}}}{a}^{\dagger }a+\displaystyle \frac{1}{2}{\omega }_{{\rm{q}}}{\sigma }_{z}+g({a}^{\dagger }{\sigma }_{-}+a{\sigma }_{+})+\displaystyle \frac{1}{2}({\epsilon }_{\mathrm{ro}}{{\rm{e}}}^{-{\rm{i}}{\omega }_{\mathrm{ro}}t}{a}^{\dagger }+{\epsilon }_{\mathrm{ro}}^{*}{{\rm{e}}}^{{\rm{i}}{\omega }_{\mathrm{ro}}t}a),\end{eqnarray} \tag{ A1 }$$

where ${\omega }_{{\rm{r}}}$ is the cavity frequency, ${\omega }_{{\rm{q}}}$ the qubit energy, g is the qubit–cavity coupling, and ${\epsilon }_{\mathrm{ro}}$ is the amplitude of the cavity drive. Here a is the annihilation operator for cavity photons, and the Pauli matrices σ describe the qubit. For large qubit–cavity detuning, such that ${\rm{\Delta }}={\omega }_{{\rm{r}}}-{\omega }_{{\rm{q}}}\gg g$, we move into the dispersive frame [36]

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{disp}}=({\omega }_{{\rm{r}}}-{\omega }_{\mathrm{ro}}){a}^{\dagger }a+\displaystyle \frac{1}{2}{\tilde{\omega }}_{{\rm{q}}}{\sigma }_{z}+\chi {a}^{\dagger }a{\sigma }_{z}+\displaystyle \frac{1}{2}({\epsilon }_{\mathrm{ro}}{a}^{\dagger }+{\epsilon }_{\mathrm{ro}}^{*}a),\end{eqnarray} \tag{ A2 }$$

where $\chi ={g}^{2}/{\rm{\Delta }}$ is the dispersive shift and ${\tilde{\omega }}_{{\rm{q}}}={\omega }_{{\rm{q}}}+\chi$ is the Lamb-shifted qubit frequency. We also move the cavity into a rotating frame at frequency ${\omega }_{\mathrm{ro}}$ and assume that ${\omega }_{\mathrm{ro}}$ is far detuned from the qubit transition ${\omega }_{{\rm{q}}}$. For a cavity drive used for a dispersive measurement, this condition is fulfilled. Since we are concerned with the qubit's dynamics under simultaneous qubit drive and measurement, we add a qubit driving term at frequency ${\omega }_{{\rm{d}}}$ to equation (A2), writing

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{qd}}=\displaystyle \frac{1}{2}({\rm{\Omega }}{\sigma }_{-}{{\rm{e}}}^{-{\rm{i}}{\omega }_{{\rm{d}}}t}+{{\rm{\Omega }}}^{*}{\sigma }_{+}{{\rm{e}}}^{{\rm{i}}{\omega }_{{\rm{d}}}t}),\end{eqnarray} \tag{ A3 }$$

where Ω is the qubit Rabi frequency. Hereafter we examine the case of resonant qubit drive, where ${\omega }_{{\rm{d}}}={\tilde{\omega }}_{{\rm{q}}}$. In a frame rotating at the qubit drive frequency, the uncoupled qubit dynamics are described by ${\hat{H}}_{{\rm{q}}}=\tfrac{1}{2}{\rm{\Omega }}{\sigma }_{x}$, leading to coherent qubit flopping at the Rabi frequency. The complete dynamics of the system will be described by a master equation of the form

$$\begin{eqnarray*}&&\dot{\rho }=-{\rm{i}}[{\hat{H}}_{\mathrm{disp}}+{\hat{H}}_{\mathrm{qd}},\rho ]+\kappa { \mathcal D }[a]\rho +{\gamma }_{1}{ \mathcal D }[{\sigma }_{-}]\rho +{\gamma }_{\varphi }{ \mathcal D }[{\sigma }_{z}]\rho ,\end{eqnarray*}$$

with the cavity loss-rate κ, the qubit decay rate ${\gamma }_{1}$, and the qubit pure dephasing rate ${\gamma }_{\varphi }$. The total qubit dephasing rate is then ${\gamma }_{2}=\tfrac{1}{2}{\gamma }_{1}+{\gamma }_{\varphi }$. The dissipative superoperators are defined by ${ \mathcal D }[\hat{o}]\rho$ $=\hat{o}\rho {\hat{o}}^{\dagger }-\tfrac{1}{2}({\hat{o}}^{\dagger }\hat{o}\rho +\rho {\hat{o}}^{\dagger }\hat{o})$.

Next we apply a polaron transformation to the system, defined by [19]

$$\begin{eqnarray}&&P=\displaystyle \sum _{i}D({\alpha }_{i})| {\rm{i}}\rangle \langle {\rm{i}}| ,\end{eqnarray} \tag{ A4 }$$

where $D(\alpha )=\mathrm{exp}\{\alpha {a}^{\dagger }-{\alpha }^{*}a\}$ is the usual field displacement operator and the sum is over the qubit states. The polaron transformation applies a qubit-state-dependent shift to the cavity field, where the complex-valued displacement amplitudes ${\alpha }_{i}$ are the cavity field amplitudes conditioned on the qubit states and are in general time-dependent. The equations of motion for the ${\alpha }_{i}$ follow from the polaron transformation under the constraint that the resulting Hamiltonian represent an un-driven cavity. This treatment and the resulting equations for the cavity field are similar to a semi-classical treatment of the cavity dynamics, while still keeping track of intrinsic and qubit-induced nonlinearities [52]. The polaron frame best captures the dynamics of the full system when the cavity field can be accurately described as a coherent state whose complex amplitude depends on the instantaneous qubit state. In practice the use of the polaron frame is restricted to weak measurement regime, as becomes evident when considering its effect on the qubit ladder operators:

$$\begin{eqnarray}&&{P}^{\dagger }{\sigma }_{-}P={D}_{\beta }{\sigma }_{-}.\end{eqnarray} \tag{ A5 }$$

Here the generalized displacement operator ${D}_{\beta }$ is defined by

$$\begin{eqnarray}&&{D}_{\beta }={D}^{\dagger }({\alpha }_{0})D({\alpha }_{1})=D(\beta ){{\rm{e}}}^{-{\rm{i}}\varphi },\end{eqnarray} \tag{ A6 }$$

with the measurement separability $\beta ={\alpha }_{1}-{\alpha }_{0}$ and the phase $\varphi ={\rm{Im}}({\alpha }_{0}{\alpha }_{1}^{*})$. This follows naturally from the fact that the polaron transformation connects a given qubit state with a corresponding coherent state in the cavity. Any change in the qubit state must be connected with a simultaneous change in the cavity field. Importantly, the transformation equation (A5) results in an expression that, for non-vanishing values of $| \beta |$, contains all orders of the cavity annihilation and creation operators. This makes it impractical to find a reduced description of the dynamics of the qubit alone, except for the weak measurement case, when $| \beta | \ll 1$ and ${D}_{\beta }\sim {\mathbb{1}}$. In this case the system reduces to a set of equations for the classical field amplitudes ${\alpha }_{i}$ conditioned on the qubit state and a Hamiltonian describing the qubit coupled to the ${\alpha }_{i}$ and an effectively un-driven cavity, describing the quantum fluctuations of the cavity field.

By applying the polaron transformation to the dissipators, we immediately find the measurement-induced dephasing since

$$\begin{eqnarray*}&&\kappa { \mathcal D }[a]\rho \to \kappa { \mathcal D }[a]\rho +{{\rm{\Gamma }}}_{{\rm{d}}}{ \mathcal D }[{\sigma }_{z}]\rho \;+\cdots ,\end{eqnarray*}$$

where ${{\rm{\Gamma }}}_{{\rm{d}}}=\tfrac{1}{2}\kappa {| \beta | }^{2}$ is the measurement-induced dephasing rate and the omitted terms can be absorbed into the coherent part of the master equation.

We now find the qubit-only dynamics by tracing the master equation over the cavity degrees of freedom in the Fock basis $| n\rangle$, to find the time evolution of the qubit's reduced density matrix:

$$\begin{eqnarray}&&\tilde{\rho }=\displaystyle \sum _{n}\langle n| \rho | n\rangle .\end{eqnarray} \tag{ A7 }$$

Finally, the jump rate in equation (4) can be found from the qubit's Bloch vector elements $\langle {\sigma }_{i}\rangle$. For resonant qubit drive, where ${\omega }_{{\rm{d}}}={\tilde{\omega }}_{{\rm{q}}}$, one finds

$$\begin{eqnarray}\dot{\langle {\sigma }_{x}\rangle } & = & -({\gamma }_{2}+{{\rm{\Gamma }}}_{{\rm{d}}})\langle {\sigma }_{x}\rangle ,\\ \dot{\langle {\sigma }_{y}\rangle } & = & -({\gamma }_{2}+{{\rm{\Gamma }}}_{{\rm{d}}})\langle {\sigma }_{y}\rangle +{\rm{\Omega }}\langle {\sigma }_{z}\rangle ,\\ \dot{\langle {\sigma }_{z}\rangle } & = & {\rm{\Omega }}\langle {\sigma }_{y}\rangle -{{\rm{\Gamma }}}_{1}\langle {\sigma }_{z}\rangle +({\gamma }_{\uparrow }-{\gamma }_{\downarrow }),\end{eqnarray} \tag{ A8 }$$

where ${{\rm{\Gamma }}}_{1}={\gamma }_{\uparrow }+{\gamma }_{\downarrow }$ includes all dissipative transitions between qubit levels (i.e. thermally induced by the environment as well as dressed dephasing). Assuming that the measurement-induced dephasing dominates, such that $\dot{\langle {\sigma }_{x}\rangle }=\dot{\langle {\sigma }_{y}\rangle }=0$, the equation of motion for $\langle {\sigma }_{z}\rangle$ takes the form:

$$\begin{eqnarray}\dot{\langle {\sigma }_{z}\rangle } & = & -\left(\displaystyle \frac{{{\rm{\Omega }}}^{2}}{{\gamma }_{2}+{{\rm{\Gamma }}}_{{\rm{d}}}}+{{\rm{\Gamma }}}_{1}\right)\langle {\sigma }_{z}\rangle +({\gamma }_{\uparrow }-{\gamma }_{\downarrow }),\\ & = & -({{\rm{\Gamma }}}_{\uparrow }+{{\rm{\Gamma }}}_{\downarrow })\langle {\sigma }_{z}\rangle +({\gamma }_{\uparrow }-{\gamma }_{\downarrow }).\end{eqnarray} \tag{ A9 }$$

We identify ${{\rm{\Gamma }}}_{\uparrow }={{\rm{\Gamma }}}_{\uparrow ,{\rm{drive}}}+{\gamma }_{\uparrow }$ and ${{\rm{\Gamma }}}_{\downarrow }={{\rm{\Gamma }}}_{\downarrow ,{\rm{drive}}}+{\gamma }_{\downarrow }$ and find the drive-induced transition rate as in equation (4):

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\uparrow ,{\rm{drive}}}={{\rm{\Gamma }}}_{\downarrow ,{\rm{drive}}}=\displaystyle \frac{{{\rm{\Omega }}}^{2}}{2({\gamma }_{2}+{{\rm{\Gamma }}}_{{\rm{d}}})}.\end{eqnarray} \tag{ A10 }$$

In contrast to the simple two-level qubit model used in the previous description, the transmon qubit used in our experiments is a genuine multi-level system. In the presence of strong driving and measurement, which is the regime of our experiments, the additional levels lead to measurably different behavior of the coupled system. Our numerical simulations therefore take into account the multi-level nature of the transmon and cavity, using experimentally determined parameters. In the rotating wave approximation, the Hamiltonian is given by:

$$\begin{eqnarray}\hat{H} & = & {{\rm{\Pi }}}_{\omega }+{\omega }_{{\rm{r}}}{a}^{\dagger }a+\displaystyle \sum _{i}\;{g}_{i}({\sigma }_{+}^{(i)}a+{\sigma }_{-}^{(i)}{a}^{\dagger })\\ & & +\displaystyle \frac{1}{2}({\epsilon }_{\mathrm{ro}}{{\rm{e}}}^{-{\rm{i}}{\omega }_{\mathrm{ro}}t}{a}^{\dagger }+{\epsilon }_{\mathrm{ro}}^{*}{{\rm{e}}}^{{\rm{i}}{\omega }_{\mathrm{ro}}t}a)+\displaystyle \frac{1}{2}({\epsilon }_{{\rm{q}}}{{\rm{e}}}^{-{\rm{i}}{\omega }_{{\rm{d}}}t}{a}^{\dagger }+{\epsilon }_{{\rm{q}}}^{*}{{\rm{e}}}^{{\rm{i}}{\omega }_{{\rm{d}}}t}a),\end{eqnarray} \tag{ B1 }$$

where ${{\rm{\Pi }}}_{\omega }={\sum }_{i}{\omega }_{i}| {\rm{i}}\rangle \langle {\rm{i}}|$ describes the transmon qubit in its eigenbasis with eigenstates $| {\rm{i}}\rangle$ and corresponding energies ${\omega }_{i}$. The qubit multi-level ladder operators are defined as ${\sigma }_{+}^{(i)}=| {\rm{i}}+1\rangle \langle {\rm{i}}|$ and ${\sigma }_{-}^{(i)}={{\sigma }_{+}^{(i)}}^{\dagger }$. The cavity has frequency ${\omega }_{{\rm{r}}}$ and is coupled with the coupling strengths g_i to the qubit transitions via a Jaynes–Cummings-type interaction. Each of the g_i can be determined from the dipole moments corresponding to different qubit transitions, which are obtained by diagonalizing the transmon qubit Hamiltonian in the charge basis [38]. The cavity measurement drive is applied at frequency ${\omega }_{\mathrm{ro}}$ and with amplitude ${\epsilon }_{\mathrm{ro}}$. In the experiment, driving of the qubit is achieved as a second order effect via off-resonant driving of the cavity at the desired qubit transition frequency; therefore in our model the qubit is driven indirectly via the cavity, with drive frequency ${\omega }_{{\rm{d}}}$ and drive strength ${\epsilon }_{{\rm{q}}}$. In the dispersive limit, $({\omega }_{{\rm{r}}}-{\omega }_{10})={\rm{\Delta }}\gg g$, this term can be written perturbatively as the directly driven model $\sim {\sum }_{i}{{\rm{\Omega }}}_{{\rm{R}},i}({\sigma }_{+}^{(i)}+{\sigma }_{-}^{(i)})$, with the Rabi frequencies for each qubit transition ${{\rm{\Omega }}}_{{\rm{R}},i}\approx {\epsilon }_{{\rm{q}}}{g}_{i}/{\rm{\Delta }}$. However, we use the exact model equation (B1) for our simulations.

To calculate the time evolution of the system under the influence of the environment, we write the master equation for the density matrix as [30]

$$\begin{eqnarray}&&\dot{\rho }={ \mathcal L }\rho =-{\rm{i}}[\hat{H},\rho ]+\kappa { \mathcal D }[a]\rho +2{\gamma }_{\varphi }{ \mathcal D }[{{\rm{\Pi }}}_{\delta {\rm{\Phi }}}]\rho +\displaystyle \sum _{i}{\gamma }_{i}{ \mathcal D }[{\sigma }_{-}^{(i)}]\rho ,\end{eqnarray} \tag{ B2 }$$

with photon loss from the cavity at rate κ, qubit decay from state $| {\rm{i}}\rangle$ at rate ${\gamma }_{i}$ and pure dephasing acting on the transmon due to fluctuations in its energy levels at rate ${\gamma }_{\varphi }$. Here ${{\rm{\Pi }}}_{\delta {\rm{\Phi }}}={\sum }_{i}\delta {{\rm{\Phi }}}_{i}| {\rm{i}}\rangle \langle {\rm{i}}|$ and $\delta {{\rm{\Phi }}}_{i}\;=\tfrac{\partial ({\omega }_{i}-{\omega }_{0})}{\partial {\rm{\Phi }}}\times {\left(\tfrac{\partial ({\omega }_{1}-{\omega }_{0})}{\partial {\rm{\Phi }}}\right)}^{-1}$, where Φ is the external magnetic flux applied to the qubit loop. Then ${ \mathcal D }[{{\rm{\Pi }}}_{\delta {\rm{\Phi }}}]\rho$ describes dephasing of the transmon due to slow fluctuations in the magnetic field through the qubit loop. As decay of superconducting qubits is mostly due to electromagnetic coupling to spurious environmental modes [52, 54], the relevant matrix elements are proportional to the respective dipole elements of the state transitions. We therefore define the relaxation rates of higher transmon levels as ${\gamma }_{i}={({g}_{i}/{g}_{0})}^{2}{\gamma }_{1}$, where ${\gamma }_{1}$ is the relaxation rate of the first excited qubit level. All rates used in the simulations were measured in independent experiments.

To simulate the measurement process, we use a stochastic master equation technique [19, 45, 46]. For a homodyne measurement, we write the stochastic master equation as

$$\begin{eqnarray}&&\dot{\rho }={ \mathcal L }\rho +\sqrt{\kappa \eta }\xi (t)(\hat{c}\rho +\rho {\hat{c}}^{\dagger }-\langle \hat{c}+{\hat{c}}^{\dagger }\rangle \rho ),\end{eqnarray} \tag{ B3 }$$

where $\hat{c}={{\rm{e}}}^{{\rm{i}}\varphi }a$ is the measurement operator for our system with the local oscillator phase φ (with respect to the cavity driving signal), κ is the measurement rate (due to leakage out of the cavity) and η is the measurement efficiency. Here $\xi (t)$ is a stochastic noise process defined by $\xi (t)={\rm{d}}W(t)/{\rm{d}}t$, where ${\rm{d}}W(t)$ is the Wiener increment and dt is the timestep between successive evaluations of the master equation. dW represents the measurement noise and has the properties

$$\begin{eqnarray}&&{\langle {\rm{d}}W\rangle }_{E}=0\quad \quad {\langle {\rm{d}}{W}^{2}\rangle }_{E}={\rm{d}}t,\end{eqnarray} \tag{ B4 }$$

where the average is over different realizations of the noise process. We then use the simulated stochastic dynamics to numerically calculate the qubit populations ${\langle | {\rm{i}}\rangle \langle {\rm{i}}| \rangle }_{t}$ and determine state transition rates with the same algorithm used on the experimental data. Here the average is defined as ${\langle \hat{o}\rangle }_{t}={\rm{Tr}}\{\rho (t)\hat{o}\}$. We also calculate the cavity field amplitude ${\langle a\rangle }_{t}$, the photon number ${\langle {a}^{\dagger }a\rangle }_{t}$, and the associated ac Stark shifts to calibrate the measurement strength.

The numerical effort involved in these simulations is quite large, as the strong measurement and qubit drive requires one to take into account up to five levels of the transmon qubit as well as up to 50 photon states to ensure proper convergence. In addition, the time steps used in the stochastic numerical simulations must be chosen much smaller than any of the intrinsic timescales of the problem to ensure that no unphysical solutions will be reached. Due to these restrictions, the number of data traces from these numerical simulations is much smaller than the number of experimental traces to which they are compared, and thus the statistical error bars are generally larger.

The model Hamiltonian for the Purcell decay rate of a multilevel qubit in a strongly driven cavity is:

$$\begin{eqnarray}&&\hat{H}={\omega }_{{\rm{r}}}{a}^{\dagger }a+\displaystyle \sum _{i}\;{\omega }_{i}| {\rm{i}}\rangle \langle {\rm{i}}| +\displaystyle \sum _{i}\;{g}_{i}({\sigma }_{i}^{+}a+{\sigma }_{i}^{-}{a}^{\dagger }),\end{eqnarray} \tag{ C1 }$$

where ${\sigma }_{i}^{+}=| {\rm{i}}+1\rangle \langle {\rm{i}}|$ and the sum goes over all relevant qubit levels. Following [30, 31], we define the qubit's Purcell decay rate with the cavity in a Fock state with photon number n as

$$\begin{eqnarray}&&{\rm{\Gamma }}(n)=\kappa {| \bar{\langle g,n| }a\bar{| e,n\rangle }| }^{2},\end{eqnarray} \tag{ C2 }$$

where the eigenstates $\bar{| g,n\rangle }$ and $\bar{| e,n\rangle }$ belong to the subspaces of the Jaynes–Cummings Hamiltonian with n and $n+1$ excitations, respectively. If the cavity is not in a Fock state with definite photon number n, the effective Purcell rate can be found by averaging over the photon number distribution of the cavity state $P(n)$ as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{P}}}=\displaystyle \sum _{n}P(n){\rm{\Gamma }}(n)/\displaystyle \sum _{n}P(n).\end{eqnarray} \tag{ C3 }$$

For a coherent state with mean photon number $\bar{n}$, the photon number distribution is $P(n)={{\rm{e}}}^{-\bar{n}}{\bar{n}}^{n}/n!$. Since in the absence of dissipation the Jaynes–Cummings Hamiltonian preserves the number of excitations, the task of calculating the Purcell rate reduces to diagonalizing the Hamiltonian for each n-excitation subspace. The general form of the Hamiltonian in the n-excitation subspace can be written as

$$\begin{eqnarray}{\hat{H}}^{(n)}=\left(\begin{array}{ccccc}0 & \sqrt{n}{g}_{0} & 0 & 0 & \ldots \\ \sqrt{n}{g}_{0} & \delta {\omega }_{1} & \sqrt{n-1}{g}_{1} & 0 & \\ 0 & \sqrt{n-1}{g}_{1} & \delta {\omega }_{2} & \sqrt{n-2}{g}_{2} & \\ 0 & 0 & \sqrt{n-2}{g}_{2} & \delta {\omega }_{3} & \\ \ldots & & & & \ldots \end{array}\right),\end{eqnarray} \tag{ C4 }$$

with the $\delta {\omega }_{k}={\omega }_{k}-k{\omega }_{{\rm{r}}}$ giving the detuning between the kth qubit level and a k-photon–Fock state. We can write the eigenstates of the n-excitation subspace in the basis of qubit–cavity product states as:

$$\begin{eqnarray}\bar{| g,n\rangle } & = & {b}_{{\rm{g}},n}| g,n\rangle +{b}_{{\rm{e}},n}| e,n-1\rangle +{b}_{{\rm{f}},n}| f,n-2\rangle \;+\cdots ,\\ \bar{| e,n-1\rangle } & = & {c}_{{\rm{g}},n}| g,n\rangle +{c}_{{\rm{e}},n}| e,n-1\rangle +{c}_{{\rm{f}},n}| f,n-2\rangle \;+\cdots ,\end{eqnarray} \tag{ C5 }$$

and so on for all other states involving higher qubit levels. Here $| g\rangle ,| e\rangle$, and $| f\rangle$ denote the first three qubit eigenstates, and the ${b}_{i,n}$ and ${c}_{i,n}$ are complex numbers calculated by diagonalizing equation (C4). Together with the action of the photon annihilation operator on cavity-Fock states, $a| n\rangle =\sqrt{n}| n-1\rangle$, and noting that a will lower the overall excitation number by one, we can then identify the relevant matrix elements of a as

$$\begin{eqnarray}&&\bar{\langle g,n| }a\bar{| e,n\rangle }=\sqrt{n+1}\ {b}_{g,n}{c}_{g,n+1}+\sqrt{n}\ {b}_{e,n}{c}_{e,n+1}+\sqrt{n-1}\ {b}_{f,n}{c}_{f,n+1}\;+\cdots ,\end{eqnarray} \tag{ C6 }$$

where we obtain the coefficients numerically, using the experimental parameters for our transmon qubit. Figure 7 shows a plot of the calculated Purcell decay rates for our experimental parameters with a coherent state in the cavity, accounting for various numbers of qubit levels. The Purcell decay is suppressed by the presence of photons in the cavity, although the effect is less pronounced when higher qubit levels are accounted for. The calculated rates for more than five qubit levels (not plotted) are very similar to those for five qubit levels.

Zoom In Zoom Out Reset image size

Our qubit is a two-state system where state transitions obey Poissonian statistics, but the finite bandwidth of the measurement apparatus—which results from the inherent bandwidth κ of circuit QED measurement, the bandwidth of the superconducting parametric amplifier and other elements in the measurement chain, and the Gaussian filtering employed to increase SNR in post-processing—tends to skew the observed dwell times in each state toward longer times and reduces the number of events seen with dwell times smaller than the measurement bandwidth. A mathematical model to compensate for these effects has been presented by Naaman and Aumentado [43].

Their model allows the states of the qubit (denoted A and B) to be independent of the corresponding states of the readout (denoted ${A}^{*}$ and ${B}^{*}$). A diagram of this model, adapted from [43], is shown in figure 8. For each readout state, the qubit can be in either state A or B. Transitions occur between A and B with rates ${{\rm{\Gamma }}}_{A}$ and ${{\rm{\Gamma }}}_{B}$. However, when the qubit and readout 'disagree', i.e. states $(B,{A}^{*})$ and $(A,{B}^{*})$, the readout can also make a transition (to $(B,{B}^{*})$ or $(A,{A}^{*})$, respectively) with rate ${{\rm{\Gamma }}}_{\mathrm{det}}$. Assuming that the readout does not change states unless the qubit has changed states (no false positives), one can write the probability distribution function $h(t)$ for the observed dwell times in ${A}^{*}$ as a function of the underlying rates ${{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B}$, and ${{\rm{\Gamma }}}_{\mathrm{det}}$ [43]:

$$\begin{eqnarray}&&h(t;{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}})=\displaystyle \frac{2}{{\theta }_{A}}{{\rm{\Gamma }}}_{A}{{\rm{\Gamma }}}_{\mathrm{det}}{{\rm{e}}}^{-\lambda t/2}\mathrm{sinh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right),\end{eqnarray} \tag{ D1 }$$

where $\lambda ={{\rm{\Gamma }}}_{A}+{{\rm{\Gamma }}}_{B}+{{\rm{\Gamma }}}_{\mathrm{det}}$ and ${\theta }_{A}=\sqrt{{\lambda }^{2}-4{{\rm{\Gamma }}}_{A}{{\rm{\Gamma }}}_{\mathrm{det}}}$. The probability distribution for dwell times in state ${B}^{*}$ is identical but with ${{\rm{\Gamma }}}_{A}$ and ${{\rm{\Gamma }}}_{B}$ interchanged (defining ${\theta }_{B}=\sqrt{{\lambda }^{2}-4{{\rm{\Gamma }}}_{B}{{\rm{\Gamma }}}_{\mathrm{det}}}$), as can be seen from the symmetry of the model. We can calculate the complementary cumulative distribution function, often called the survival function, by integrating the expression in (D1):

$$\begin{eqnarray}&&s(t;{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}})={\displaystyle \int }_{t}^{\infty }h({t}^{\prime }){\rm{d}}{t}^{\prime }=\displaystyle \frac{{{\rm{e}}}^{-\lambda t/2}}{{\theta }_{A}}\left[\lambda \mathrm{sinh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right)+{\theta }_{A}\mathrm{cosh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right)\right].\end{eqnarray} \tag{ D2 }$$

Zoom In Zoom Out Reset image size

The survival function is useful for handling right-censored dwell times, where the observed dwell time represents a lower bound on the true dwell time in the state. The survival function $s(\tau )$ gives the probability of a dwell time $t\gt \tau$, which is precisely the probability distribution function needed for right-censored events. Given a data set of ${n}_{{A}^{*}}$ dwell times in state ${A}^{*}$ and ${n}_{{B}^{*}}$ dwell times in state ${B}^{*}$, denoted $\{{t}_{i}\}$ and $\{{t}_{j}\}$, respectively, as well as corresponding censoring variables $\{{\delta }_{i}\}$ and $\{{\delta }_{j}\}$, where $\delta =0(1)$ indicates an uncensored (right-censored) dwell time, we can write the likelihood function as:

$$\begin{eqnarray}{ \mathcal L }({{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}}) & = & \displaystyle \prod _{i=1}^{{n}_{{A}^{*}}}h{({t}_{i};{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}})}^{1-{\delta }_{i}}s{({t}_{i};{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}})}^{{\delta }_{i}}\\ & & \times \displaystyle \prod _{j=1}^{{n}_{{B}^{*}}}h{({t}_{j};{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{\mathrm{det}})}^{1-{\delta }_{j}}s{({t}_{j};{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{\mathrm{det}})}^{{\delta }_{j}}.\end{eqnarray} \tag{ D3 }$$

We then use nonlinear optimization methods to maximize the value of ${ \mathcal L }$ by varying the parameters ${{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B}$, and ${{\rm{\Gamma }}}_{\mathrm{det}}$. Because of the limitations of floating point arithmetic, we actually perform the maximization on the log-likelihood function $L=\mathrm{ln}({ \mathcal L })$, which has the form:

$$\begin{eqnarray}L({{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}}) & = & \displaystyle \sum _{i=1}^{{n}_{{A}^{*}}}[(1-{\delta }_{i})\mathrm{ln}[h({t}_{i};{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}})]+{\delta }_{i}\mathrm{ln}[s({t}_{i};{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}})]]\\ & & +\displaystyle \sum _{j=1}^{{n}_{{B}^{*}}}[(1-{\delta }_{j})\mathrm{ln}[h({t}_{j};{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{\mathrm{det}})]+{\delta }_{j}\mathrm{ln}[s({t}_{j};{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{\mathrm{det}})]].\end{eqnarray} \tag{ D4 }$$

The functional forms of $\mathrm{ln}[h(t)]$ and $\mathrm{ln}[s(t)]$ for dwell times in state ${A}^{*}$ are:

$$\begin{eqnarray}&&\mathrm{ln}[h(t;{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}})]=\mathrm{ln}\left(\displaystyle \frac{2{{\rm{\Gamma }}}_{A}{{\rm{\Gamma }}}_{B}}{{\theta }_{A}}\right)-\displaystyle \frac{\lambda t}{2}+\mathrm{ln}\left(\mathrm{sinh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right)\right),\end{eqnarray} \tag{ D5 }$$

$$\begin{eqnarray}&&\mathrm{ln}[s(t;{{\rm{\Gamma }}}_{A},{{\rm{\Gamma }}}_{B},{{\rm{\Gamma }}}_{\mathrm{det}})]=-\mathrm{ln}({\theta }_{A})-\displaystyle \frac{\lambda t}{2}+\mathrm{ln}\left[\lambda \mathrm{sinh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right)+{\theta }_{A}\mathrm{cosh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right)\right].\end{eqnarray} \tag{ D6 }$$

The corresponding functions for dwell times in state ${B}^{*}$ can be found by interchanging ${{\rm{\Gamma }}}_{A}$ with ${{\rm{\Gamma }}}_{B}$, and ${\theta }_{A}$ with ${\theta }_{B}$. For sufficiently large values of ${\theta }_{A}t$, these expressions may cause overflows in double-precision floating point arithmetic. Therefore, for ${\theta }_{A}t\gt 40$, we make the following approximations, which introduce fractional errors of less than 10⁻¹⁸ for each approximated term:

$$\begin{eqnarray}&&\mathrm{ln}\left(\mathrm{sinh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right)\right)\approx \displaystyle \frac{{\theta }_{A}t}{2}-\mathrm{ln}2,\end{eqnarray} \tag{ D7 }$$

$$\begin{eqnarray}&&\mathrm{ln}\left[\lambda \mathrm{sinh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right)+{\theta }_{A}\mathrm{cosh}\left(\displaystyle \frac{{\theta }_{A}t}{2}\right)\right]\approx \displaystyle \frac{{\theta }_{A}t}{2}+\mathrm{ln}(\lambda +{\theta }_{A})-\mathrm{ln}2.\end{eqnarray} \tag{ D8 }$$