Simultaneous bistability of qubit and resonator in circuit quantum electrodynamics

We explore the joint activated dynamics exhibited by two quantum degrees of freedom: a cavity mode oscillator which is strongly coupled to a superconducting qubit in the strongly coherently driven dispersive regime. Dynamical simulations and complementary measurements show a range of parameters where both the cavity and the qubit exhibit sudden simultaneous switching between two metastable states. This manifests in ensemble averaged amplitudes of both the cavity and qubit exhibiting a partial coherent cancellation. Transmission measurements of driven microwave cavities coupled to transmon qubits show detailed features which agree with the theory in the regime of simultaneous switching.


I. INTRODUCTION
The generalized Jaynes-Cummings (GJC) model provides a simple basis for describing the interactions between a quantized electromagnetic field and multilevel atoms. Its nonlinearity lies at the heart of cavity quantum electrodynamics (cavity QED), where natural atoms are coupled to cavity photons [1], and circuit quantum electrodynamics (circuit QED), where artificial atoms are coupled to resonators of various dimensionalities [2][3][4][5][6]. The Jaynes-Cummings (JC) interaction also emerges in areas of current interest such as optomechanics in the linearized regime [7] and in the Bose-Hubbard model [8]. There is a large body of work on the resonant and strong-coupling regime of the driven-dissipative JC oscillator [9,10], where driving induces a dynamical Rabi splitting [11,12]. The high excitation strongdispersive regime is also of great interest, for example, in the context of amplifiers [13], squeezing associated with the parametric oscillator [14] and the implementation of qubit readout schemes [15][16][17][18]. In this context, fluctuation-induced switching between metastable states in the driven-dissipative GJC system, which involves two quantum degrees of freedom, has not been directly studied but the theory of quantum activation motivates the interest in this scenario [19][20][21][22][23][24][25].
In this paper, we demonstrate that in a nonlinear intermediate driving regime of circuit QED, the system dynamics exhibits a simultaneous bistability of the qubit and resonator. We support this claim both theoretically, with analytical and numerical results, as well as with experimental measurements obtained from a circuit QED device, consisting of a transmon qubit coupled to a 3D * t.mavrogordatos@ucl.ac.uk FIG. 1. Maxwell-Bloch steady-state bistability of the JC oscillator in the dispersive regime [6] for the coupling strength to detuning ratio g/δ = 0.14, the drive amplitude to photon loss ratio ε d /(2κ) = 25/3 and the photon loss to spontaneous emission ratio 2κ/γ = 12. Here, B(D) denote the bright (dim) semiclassical states. microwave cavity [26]. In particular we report the following findings. (1) The switching process occurs simultaneously for the two coupled quantum oscillators. (2) The ensemble averaged amplitudes of both cavity and qubit exhibit a coherent partial cancellation. Such cancellation, predicted theoretically for a single nonlinear mode by Drummond and Walls in [27], is here verified experimentally and shown to occur for both coupled oscillators. (3) The Duffing oscillator model with one quantum degree of freedom is not sufficient to account for the observed cavity nonlinearity when the two coupled quantum degrees of freedom are involved in the switching process. We also note that, the JC bistable semiclassical intracavity amplitude |α| 2 ss (with α = a ) in the steady state, plotted in Fig. 1(a), does not show any coherent cancellation feature, nor does the bistable average atomic inversion ζ ss = σ z ss shown in Fig. 1(b). Instead, both curves are skewed Lorentzians with the position of their peaks approaching the bare cavity frequency for increasing drive. By increasing the drive strength beyond what is shown in Fig. 1(a) one finds a critical point in the phase space between bistable and linear behavior, lying on the line where the frequency of the drive equals the bare cavity resonance frequency. Beyond that point, the system behaves as a linear oscillator, in contrast to the corresponding Duffing oscillator [17], and exhibits no bistability.

II. THEORETICAL MODELS
In order to develop a comprehensive understanding of the system response, we will consider several different theoretical models: (i) multilevel transmon-cavity GJC model -the most complete in the context of superconducting devices; (ii) two-level 'atom'-cavity JC model which is universal to many strong light-matter coupling scenarios, and (iii) a simplified dressed-cavity Duffing oscillator approximation. To complement the numerical simulations of the three models, we additionally compare the results for the cavity transmission with those coming from an analytical formula, obtained by modeling the transmon itself as a nonlinear oscillator.
We will now define the system Hamiltonians associated with models (i)-(iii). When one cavity field mode of frequency ω c (with corresponding photon annihilation and creation operators a and a † respectively) is coupled to a multilevel system with unperturbed states |n , the coherently-driven GJC Hamiltonian can be written as (setting = 1) [2] H (i) (1) where ε d is the strength of a monochromatic external field with frequency ω d driving the cavity mode. The sum in the third term describes the interaction and is customarily modified to m g m,m+1 |m m + 1| a † + h.c. in the rotating wave approximation (RWA). The interaction energies in the RWA have the approximate form g mn ≈ g √ m + 1 δ m+1,n , with g being the dipole coupling strength. Depending on the range of n and the form of ω n we can distinguish the two-level atom (n = 1, m = 0, ω n = ω q ) -the JC model -from a transmon (n = 1, 2, . . . , N max , ω n = ω(n)) -the GJC model. The JC model in the RWA reads with σ ± the raising (lowering) pseudospin operators and σ z = 2σ + σ − − 1 the inversion operator. In the presence of dissipation, the cavity mode is damped at a rate 2κ [28] while spontaneous emission is present at a rate γ for a qubit dephased at a rate γ φ . It is possible to approximate the Hamiltonian of Eq. (2) further in the strongly dispersive regime, defined by the strong detuning δ = |ω c − ω q | g between the two coupled oscillators. Under an appropriate decoupling transformation [21,29] H JC can be recast in the form H JC = ω c a † a + (1/2)(ω c − ∆)σ z , involving the operator ∆ = δ 2 + 4g 2 N where N = a † a + σ + σ − is the operator of the total number of excitations N . In the dispersive regime, provided that N/N crit 1 [30], we can expand up to the quartic order in the field variables. After normal ordering we obtain the following dressedcavity Duffing oscillator Hamiltonian where setting σ z = −1 is a justifiable approximation for low enough driving amplitudes, yielding a bistable quantum Duffing oscillator [35]. The third term in the parentheses is the leading-order term χ (0) c = (g 2 /δ)σ z which is the familiar Stark shift [36] and provides a valuable tool for qubit readout (with χ (0) c κ). Here, in the bifurcating dispersive region we are studying, the following hierarchy of scales applies [17]: The intracavity excitation number is of the order of N crit , where this perturbation expansion is not strictly valid. However, as we will see later, it gives qualitatively meaningful results.
Having defined the different model Hamiltonians under consideration, we now evolve the corresponding master equations (MEs) in the finite Hilbert state basis numerically, starting from a Fock state of zero photons and the qubit in the ground state, until it reaches a steady state. Note that the steady state obtained is independent of the choice of the initial conditions.

III. ACTIVATED DYNAMICS IN THE DISPERSIVE REGIME
Driving the system beyond the low power regime has a profound effect on the response. We illustrate this fact in  Fig. 2(b)]. The corresponding photon histogram shows quasi-Poissonian statistics obeyed by the two metastable states, one with mean photon occupation of the order of N crit (called 'bright' state), one with mean occupation of about a photon (called 'dim' state) as well as the distribution of a nonclassical (called 'dark') state (for more details see the Supplementary Information). In Fig. 2(c) we draw a sketch illustrating the qubit distribution in the steady state, as viewed from the north pole of the Bloch sphere, for a single quantum trajectory. Red arrows point to the two semiclassical qubit states, corresponding to the two metastable quasicoherent cavity states [depicted by color contour plots in Fig.  2(d)] between which quantum-activated switching takes place. The Q function plot in Fig 2(d) also shows the position in the phase space of the coherently cancelling states. The equal height of the Q function peaks indicates the boundary of a first-order dissipative quantum phase transition [8,39]. This transition is marked by the switching rates to the bright and to the dim state being of the same order of magnitude [8,19]. In our case, as well as in the exact photon statistics section of [27], switching is induced by quantum fluctuations only, as the thermal bath to which the system is coupled is at zero temperature. The trajectories depicted for illustration in Fig. 2 evidence sudden simultaneous jumps. Bistability and synchronisation were studied for the two-level Rabi model in [40]. Note that, in contrast to the cavity field and | σ − |, the exact ME results for the photon number and σ z show no coherent cancellation in the steadystate response. The mean-field behavior, depicted in Fig.  1, shows further that the coherent cancellation is purely a quantum effect at zero temperature, occurring when forming the ensemble averaged quantities, and is already present in the most approximate dressed-cavity Duffing model, even if the qubit is unmonitored.

IV. THE NONLINEAR RESONATOR TRANSMISSION LINESHAPE
In Fig. 3(a) we compare theoretical and experimental transmission amplitudes of a 3D cavity with embedded transmon (device D 1 , details in the Supplementary Information) for different driving strengths. We observe that, as the driving power is increased, the experimental cavity lineshape develops nonlinear features and a coherent cancellation dip appears. We find perfect agreement with the GJC model. In Fig. 3(b) we show the cavity transmission for the intermediate drive power of −46 dBm for all the models discussed. We observe that the JC model predicts the split of the main peak at the correct position, as opposed to its Duffing reduction, yet fails to capture the position of the dip emerging at a lower frequency. The GJC model with four transmon levels can resolve all the details necessary for a quantitative comparison and provides indeed the most complete description of the cavity nonlinearity.
The behavior of the GJC oscillator depends strongly on the drive strength and frequency, and their relation to the coupling and dissipation rates [14,17]. We find theoretically that the coherent cancellation dip in transmission, discovered by Drummond and Walls for the Duffing oscillator [27], appears in the dispersive response of both the cavity and the qubit within the full nonlinearity of the JC model, where the departure from the mean-field predictions is appreciable and the Duffing oscillator approximation is no longer valid. For the Duffing oscillator the dip is present in the first moment of the field operator, | a |, calculated using the generalized P representation [27]. It is purely a phase effect as the dip does not appear in the number of intracavity photons in the steady state n p = a † a . The coherent cancellation dip appears as well in the qubit projection | σ − | (see the Supplementary Information for more details). The presence of this dip in the cavity response has also been observed in our experimental measurements, which depict the development of nonlinearity for increasing drive strengths within the region of bistability. Similar cancellation effects appear also in classical dissipative systems out of equilibrium, in the presence of thermal fluctuations [27]. The observed dip, appearing progressively in a measurement of complex amplitude, is due to the phase differences between the two metastable states. In the experimental response we can also discern a split of the main peak, alluding to dynamical Rabi splitting [12]. The position of the dip shifts to the lower frequencies with increasing drive strength, while the split gradually fades away in favour of a Duffing-type profile.
In order to gain a further insight, we undertake an analytical approach by identifying an effective Hamiltonian to produce a (second-order) Fokker-Planck equation (FPE) for the transmon, following the adiabatic elimination of the cavity. FPEs in the generalized P representation can be used to solve exactly for the steady state of quantum systems subject to the so-called 'potential conditions', and have been used to study single nonlinear resonator systems [6]. For our two-oscillator model, these conditions are not satisfied, yet in the limit 2κ γ, γ φ the cavity can be eliminated in a similar fashion to the method of [41]. This process leaves a FPE for an effective one-oscillator system, which resembles a driven, damped quantum Duffing oscillator with anharmonicity χ [27] but with parameters that are nontrivial functions of those of the full system. Full details of this method can be found in [42]. The first moment of the cavity field in the steady state is where 0 F 2 (x, y; z) is a generalized hypergeometric function, and we have defined effective decay constants for the cavityγ c = κ + 2i∆ω c and transmonγ q = γ + 2i∆ω q + 2g 2 /γ c respectively (with ∆ω c(q) = ω c(q) − ω d ), effective drive strengthε = −2igε d /γ c and also c =γ q /(2iχ).
The calculated transmission amplitude via Eq. 4 is plotted in Fig. 3(b) and compared to the exact ME results alongside the experimental data. The effective Fokker-Planck model exaggerates the actual nonlinearity in this regime, yet a lower value of χ allows us to capture the essential features of the full transmon-cavity-driven interaction (more details in the Supplementary Information).

V. DISCUSSION AND CONCLUDING REMARKS
We have examined the dispersive interaction of a single qubit and a microwave cavity mode, tracking nonlinearity with increasing drive power. When the regime of bistability is reached, simultaneous switching events allow for both of the metastable states to participate even at zero temperature. Their different phases cause the 'dip' in coherent transmission, for which we have presented theoretical and experimental evidence. Interestingly, the dim quasicoherent state is preceded by a lower amplitude nonclassical state which is not predicted by the mean-field treatment. This state is characterized by very low photon numbers and intense fluctuations in the qubit inversion, which occupies now the north pole of the Bloch sphere. For high excitations, beyond the Duffing oscillator regime, both the cavity and the qubit participate in the switching, and the quantitative comparison with the experiment necessitates the inclusion of more than two levels of the transmon. The superconducting devices we have considered serve as examples of quantum activation with more than one quantum oscillator.
The data underlying this work is available without restriction [43].

VI. EXPERIMENTAL SETUP
Our devices consist of single-junction transmon qubits embedded in aluminium 3D cavities thermally anchored at 10 mK in a dilution refrigerator. We have measured two different devices, referred to as D 1 and D 2 , whose characteristic parameters can be found in Table I. Fig.  4(a) shows the device D 1 consisting of a two-ports Al microwave cavity which holds a lithographically patterned Al transmon qubit [ Fig. 4(b)] on a sapphire substrate. The Josephson junction was fabricated using a doubleangle shadow evaporation technique. Fig. 4(c) shows the experimental setup used for measuring D 1 , a similar setup was used for D 2 . For the device D 1 (D 2 ) we measured the transmitted (reflected) signal over many experiments using a heterodyne detection scheme and computed the average signal.
Each cavity has a bare resonant frequency of f c = ω c /(2π) supporting a mode coupled to a transmon with the lowest transition frequency f q ≡ f 01 = ω 01 /(2π). The interaction between the uncoupled qubit state |n and the cavity with coupling strength g causes a dispersive shift of the cavity resonance f c , depending on the state of the qubit.
In Fig. 5 we show the experimental reflection amplitude for a different device (D 2 , see Table I for more details), as a function of the driving parameters, next to our theoretical predictions. The experimental Lorentzian curve [ Figs. 5(a, b)] shifts towards lower frequencies with increasing power, followed by the emergence of a dip separation, in close agreement with the theoretical predictions shown in Figs. 5(c, d). The horizontal line at −32 dBm marks the onset of bistability, corresponding to a frequency shift of κ. The analytical formulation already captures the pair of emerging dip lines occurring for increasing drive strength, in close agreement with the experiment. The ME results show more conspicuously the shift of the resonance curve to the left for lower drive strengths, but lead to a weaker separation because of thermal broadening.

A. Formulation of the Master Equation
The master equation (ME) corresponding to the open driven GJC model consisting of a multilevel system (transmon qubit) coupled to a single cavity mode with resonance frequency ω c readṡ where L{A, ρ} = 2AρA † −A † Aρ−ρA † A is the Liouvillian super-operator acting upon the Lindblad operator A, ρ is the reduced density operator for the cavity-atom system, 2κ the cavity relaxation rate, γ (γ φ ) the transmon qubit dissipation (dephasing) rate and n(ω c ) = [e ωc/T − 1] −1 is the mean photon number for an oscillator with frequency ω c at temperature T (k B = 1) (we assume no thermal occupation for the transmon qubit). Expressions for the terms α j and β j can be found in [12].
The Hamiltonian H JC of the main text can also be inserted as the system Hamiltonian in (5), where customarily one uses pseudospin operators (σ ± , σ z ) to account for the two-level system. The dressed-cavity system Hamiltonian in the Duffing approximation of the JC model, on the other hand, deserves special attention because the The experimental setups used to measure the transmitted signal for the device D1. An input microwave signal (RF) is sent to the sample through a highly attenuated microwave line. The output signal is first amplified by a cold amplifier and then by a series of room temperature amplifiers before being demodulated and recorded by an analog to digital converter (ADC). A bandpass filter (BP), a low pass filter (LP) and two circulators are used to prevent thermal and amplifier noise from reaching the sample.
Device  generator of the decoupling transformation required to decouple cavity and qubit does not commute with the drive term and the super-operators transform non trivially. As discussed in the main text, this approximation reads The above operator may be inserted in (5) for σ z = −1 together with a drive term only to the leading order in g/δ and N −1/2 [17].

B. Excitations in the JC oscillator
We will now focus on the behavior of the qubit during the bistable switching. In Fig 6(a) we depict the qubit response in the driving regime where the Duffing approximation is not applicable. As mentioned in the main text, we can observe that | σ − | exhibits the characteristic coherent cancellation dip. In Fig. 6(b) we present the full scatter-plot for the qubit vector in the Bloch sphere, corresponding to the illustration of Fig. 2(c) (but now including the data for the dark state). In the inset of Fig.  7 we show one of the two components (the cavity field) of the system total excitation, the eigenvalue of the operator N for the JC oscillator. None of the constituents of N (n p and σ + σ − ) show the coherent cancellation dip that characterizes the corresponding complex amplitude. Instead, they both present a distorted Lorentzian profile, which marks the onset of amplitude bistability.

C. Bimodality in the GJC
Bimodality is the underlying feature shaping the homodyne cavity response in the driving regime we have considered in the main text for the JC oscillator and its Duffing reduction. We show here that cavity bimodality is also present for a 4-level transmon, as indicated by the bimodal distribution in Fig. 8(a). The emergence of bimodality can be also traced in the transmon level occupation depicted in Fig. 8(b). The 'dark' state is indicated by the departure from the 'dim' state Poissonian  Fig. 2(b). The inset depicts np = a † a as a function of the drive frequency for the same parameters used in Fig. 2. The point on the curve marks the coordinates corresponding to the response depicted in Fig. 2. distribution at P 0 [inset of Fig. 8(c)].

D. Approaches to single-atom bistability
Single-atom absorptive bistability was first reported at resonance in [32], where two distinct states were identified that are long-lived in terms of the characteristic decay times of the two-level atom and cavity. Adiabatic elimination of the atomic variables has been used in [44] to extract approximate Gaussian distributions able to adequately calculate mean square photon fluctuations for absorptive optical bistability at resonance. However, adiabatic elimination underestimates the dynamical rôle of the qubit in shaping the system response. There has been evidence that the effect of the qubit on the cavity mode is significant in a situation, where the drive is constantly tuned to maintain resonance [15]. The two degrees of freedom describing the system and the inherent square root nonlinearity mark the departure from the Duffing oscillator [17]. Even in the dressed model describing dynamic Rabi splitting for the driven JC oscillator [14], with one degree of freedom, the effective Hamiltonian leads to an equivalent description for the master equation of resonance fluorescence which, as is well known, cannot be reduced to a Fokker-Planck equation (FPE) [45] following the adiabatic elimination used in Ref. [44].

E. The effective Fokker-Planck model
We will now provide more details for the adiabatic elimination of the cavity in the transmon-cavity system. In this model we approximate the level structure of the transmon itself by a Duffing oscillator. We set ω n = ω q n − χn(1 − n)/2 in Eq. (1), which is valid when the transmon occupation is relatively low. This Duffing approximation is a consequence of retaining the lowest order terms in the expansion of the Josephson potential [2]. We then adiabatically eliminate the cavity from the system in the limit 2κ γ, and solve the Fokker-Planck equation analytically for the steady-state distribution. The solution obtained retains much of the structure of the full system despite the elimination process, and also allows us to recover the moments of the cavity field in the steady state. Thus, we model the transmon-cavity system as a nonlinear oscillator, with the first moment given by the result in Eq. (4).
We can also use the intracavity amplitude to plot the reflection of the cavity. In Fig. 5(c) we show the calculated reflection for the device D 2 , compared with experimental results. For a further comparison we solve numerically the corresponding ME for the same parameters, and account for a finite temperature T 200 mK in the presence of dissipation. For this device we see that the features of the analytical reflection spectra agree qualitatively with both experiment and ME at all four plotted powers.
For the device D 1 , this analytical solution does not provide a good fit to the experimental data, since we are no longer operating in the limit 2κ γ. When κ ≈ γ, this solution appears to overestimate the effect of the transmon nonlinearity on the cavity transmission spectrum. In Fig. 9 we plot the cavity amplitude for different values of the anharmonicity coefficient, comparing with the experimental data. We see that for the true device anharmonicity of χ/2π = −150 MHz the response resembles the JC oscillator in Fig. 3(b). If a much lower anharmonicity coefficient χ/2π = −20 MHz is selected, then the response is in closer agreement with the GJC model and experiment. This is consistent with the idea that the nonlinearity is overestimated in this approximation -if the experimental nonlinearity happens to be large it is then additionally magnified by the approximation so that the transmon behaves like a true two-level system. The solution to make the model useful in this regime is to choose an effective anharmonicity parameter χ smaller then the one measured by the experiment, so that it is enhanced to the level of the true experimental value by the approximation.