Photon blockade and the quantum-to-classical transition in the driven-dissipative Josephson pendulum coupled to a resonator

TL;DR

This paper explores the quantum-to-classical transition in a driven Josephson pendulum coupled to a resonator, revealing nonlinear photon dynamics, multi-photon blockade effects, and the impact of dissipation models on the system's behavior.

Contribution

It provides a quantum mechanical analysis of the Josephson-Rabi model, highlighting the limitations of the Jaynes-Cummings approximation and comparing dissipation formalisms.

Findings

Transition occurs at higher drive powers quantum mechanically.

Photon number dependence on drive power is nonlinear.

Multi-photon blockade is sensitive to state truncation.

Abstract

We investigate the driven quantum phase transition between the oscillating motion and the classical nearly free rotations of the Josephson pendulum coupled to a harmonic oscillator in the presence of dissipation. We refer to this as the Josephson-Rabi model. This model describes the standard setup of circuit quantum electrodynamics, where typically a transmon device is embedded in a superconducting cavity. We find that by treating the system quantum mechanically this transition occurs at higher drive powers than expected from an all-classical treatment, which is a consequence of the quasiperiodicity originating in the discrete energy spectrum of the bound states. We calculate the photon number in the resonator and show that its dependence on the drive power is nonlinear. In addition, the resulting multi-photon blockade phenomenon is sensitive to the truncation of the number of states in the transmon, which limits the applicability of the standard Jaynes-Cummings model as an approximation for the pendulum-oscillator system. We calculate the nth order correlation functions of the blockaded microwave photons and observe the differences between the rotating-wave approximation and the full multilevel Josephson-Rabi Hamiltonian with the counter-rotating terms included. Finally, we compare two different approaches to dissipation, namely the Floquet-Born-Markov and the Lindblad formalisms.

Tables

Table 1. Table 1: Parameters of the driven and dissipative oscillator-pendulum system. The numerical values of the angular frequencies and energies used in the numerical simulations are given in units of ω r subscript 𝜔 r \omega_{\rm r} and ℏ ​ ω r Planck-constant-over-2-pi subscript 𝜔 r \hbar\omega_{\rm r} , respectively. We note that ω q subscript 𝜔 q \omega_{\rm q} is determined by E C subscript 𝐸 C E_{\rm C} and η 𝜂 \eta , see the text.

Symbol	Parameter	Value
${\omega }_q$	qubit frequency	1.0
${\omega }_d$	drive frequency	0.98
${\omega }_p$	plasma oscillation frequency	1.08
$g$	coupling frequency	0.04
$\kappa$	resonator dissipation rate	0.002
$k_{B}T$	thermal energy	0.13
$E_C$	charging energy	0.07
$\eta$	energy ratio $E_J/E_C$	30

Equations

$${\hat{H}}_0={\hat{H}}_{\mathrm{r}}+{\hat{H}}_{\mathrm{t}}+{\hat{H}}_{\mathrm{c}},$$

$${\hat{H}}_{\mathrm{r}}$$

$${\hat{H}}_{\mathrm{t}}$$

$${\hat{H}}_{\mathrm{c}}$$

$${\hat{H}}_0=\hslash {\omega }_{\mathrm{r}}{\hat{a}}^{†}\hat{a}+\sum _{k=0}^{K-1}\hslash {\omega }_k|k⟩⟨k|+\hslash g({\hat{a}}^{†}+\hat{a})\sum _{k,\ell =0}^{K-1}{\hat{\Pi }}_{k\ell }.$$

$$\begin{array}{rl}{\hat{H}}_0\approx & \hslash {\omega }_{\mathrm{r}}{\hat{a}}^{†}\hat{a}+\sum _{k=0}^{K-1}\hslash {\omega }_k|k⟩⟨k|\\ & +\hslash g\sum _{k=0}^{K-2}\left(\hat{a}{\hat{\Pi }}_{k,k+1}^{†}+{\hat{a}}^{†}{\hat{\Pi }}_{k,k+1}\right),\end{array}$$

$$\hat{N}={\hat{a}}^{†}\hat{a}+\sum _{k=0}^{K-1}k|k⟩⟨k|,$$

$${\hat{H}}_{\mathrm{R}}=\hslash {\omega }_{\mathrm{r}}{\hat{a}}^{†}\hat{a}+\hslash {\omega }_{\mathrm{q}}{\hat{\sigma }}_{+}{\hat{\sigma }}_{-}+\hslash g_{01}({\hat{a}}^{†}+\hat{a}){\hat{\sigma }}_{\mathrm{x}}.$$

$${\hat{H}}_{\mathrm{JC}}=\hslash {\omega }_{\mathrm{r}}{\hat{a}}^{†}\hat{a}+\hslash {\omega }_{\mathrm{q}}{\hat{\sigma }}_{+}{\hat{\sigma }}_{-}+\hslash g_{01}({\hat{a}}^{†}{\hat{\sigma }}_{-}+\hat{a}{\hat{\sigma }}_{+}),$$

$$E_{n_r,\pm }$$

$$|n_r,\pm ⟩$$

$$\hat{N}={\hat{a}}^{†}\hat{a}+{\hat{\sigma }}_{+}{\hat{\sigma }}_{-}.$$

$${\hat{H}}_{\mathrm{d}}=\hslash A\cos ({\omega }_{\mathrm{d}}t)[{\hat{a}}^{†}+\hat{a}],$$

$$⟨n_r,\pm |({\hat{a}}^{†}+\hat{a})|0,0⟩$$

$$⟨n_r,\pm |({\hat{a}}^{†}+\hat{a})|{\ell }_r,\pm ⟩$$

$$\hat{H}={\hat{H}}_{\mathrm{S}}+{\hat{H}}_{\mathrm{B}}+{\hat{H}}_{\mathrm{int}},$$

$${\hat{H}}_{\mathrm{B}}$$

$${\hat{H}}_{\mathrm{int}}$$

$${\hat{H}}_{\mathrm{int}}^{\mathrm{t}}=\hslash \hat{n}\sum _k{f}_k({\hat{c}}_k^{†}+{\hat{c}}_k),$$

$$\begin{array}{rl}\frac{\text{d}\hat{\rho }}{\text{d}t}=& -\frac{i}{\hslash }[{\hat{H}}_{\mathrm{S}},\hat{\rho }]+\kappa [n_{\mathrm{th}}({\omega }_{\mathrm{r}})+1]\mathcal{L}[\hat{a}]\hat{\rho }\\ & +\kappa n_{\mathrm{th}}({\omega }_{\mathrm{r}})\mathcal{L}[{\hat{a}}^{†}]\hat{\rho }\\ & +\sum _{k\ell }{\Gamma }_{k\ell }\mathcal{L}[|\ell ⟩⟨k|]\hat{\rho },\end{array}$$

$${\Gamma }_{k\ell }=|⟨\ell |\hat{n}|k⟩{|}^{2}S({\omega }_{k\ell }).$$

$$i\hslash \frac{\text{d}missing}{\text{d}t}|\Psi (t)⟩={\hat{H}}_{\mathrm{S}}(t)|\Psi (t)⟩,$$

$$|\Psi (t)⟩=e^{-i\epsilon t/\hslash }|\Phi (t)⟩,$$

$$\hat{U}(t_2,t_1)|\Psi (t_1)⟩=|\Psi (t_2)⟩,$$

$$i\hslash \frac{\text{d}missing}{\text{d}t}\hat{U}(t,0)={\hat{H}}_{\mathrm{S}}(t)\hat{U}(t,0).$$

$$\hat{U}(\tau ,0)|\Phi (0)⟩$$

$$\hat{U}(t,0)|{\Phi }_{\alpha }(0)⟩=e^{-i{\epsilon }_{\alpha }t/\hslash }|{\Phi }_{\alpha }(t)⟩.$$

$$\begin{array}{rl}{\dot{\rho }}_{\alpha \alpha }(t)& =\sum _{\nu }\left[{\Gamma }_{\nu \alpha }{\rho }_{\nu \nu }(t)-{\Gamma }_{\alpha \nu }{\rho }_{\alpha \alpha }(t)\right],\\ {\dot{\rho }}_{\alpha \beta }(t)& =-\frac{1}{2}\sum _{\nu }\left[{\Gamma }_{\alpha \nu }+{\Gamma }_{\beta \nu }\right]{\rho }_{\alpha \beta }(t),\alpha \ne \beta ,\end{array}$$

$$\begin{array}{rl}{\Gamma }_{\alpha \beta }& =\sum _{\ell =-\infty }^{\infty }\left[{\gamma }_{\alpha \beta \ell }+n_{\mathrm{th}}(|{\Delta }_{\alpha \beta \ell }|)\left({\gamma }_{\alpha \beta \ell }+{\gamma }_{\beta \alpha -\ell }\right)\right],\\ {\gamma }_{\alpha \beta \ell }& =\frac{\pi }{2}\kappa \theta ({\Delta }_{\alpha \beta \ell })\frac{{\Delta }_{\alpha \beta \ell }}{{\omega }_{\mathrm{r}}}|X_{\alpha \beta \ell }{|}^2.\end{array}$$

$$X_{\alpha \beta \ell }=\frac{1}{\tau }{\int }_{t_0}^{t_0+\tau }\text{d}t{e}^{-i\ell {\omega }_{d}t}⟨{\Phi }_{\alpha }(t)|({\hat{a}}^{†}+\hat{a})|{\Phi }_{\beta }(t)⟩,$$

$$p_{\alpha }=\frac{{\sum }_{\nu \ne \alpha }{\Gamma }_{\nu \alpha }{p}_{\nu }}{{\sum }_{\nu \ne \alpha }{\Gamma }_{\alpha \nu }},$$

$$N_{\mathrm{r}}=\sum _{\alpha }{p}_{\alpha }⟨{\hat{a}}^{†}\hat{a}{⟩}_{\alpha },$$

$$⟨{\hat{a}}^{†}\hat{a}{⟩}_{\alpha }=\frac{1}{\tau }{\int }_{t_0}^{t_0+\tau }\text{d}t⟨{\Phi }_{\alpha }(t)|{\hat{a}}^{†}\hat{a}|{\Phi }_{\alpha }(t)⟩,$$

Full text

I Introduction

The pendulum, which can be seen as a rigid rotor in a gravitational potential Baker and Blackburn (2005), is a quintessential nonlinear system. It has two extreme dynamical regimes: the low-energy regime, where it can be approximated as a weakly anharmonic oscillator, and the high-energy regime, where it behaves as a free rotor. Most notably, the pendulum physics appears in systems governed by the Josephson effect, where the Josephson energy is analogous to the gravitational energy, and the role of the momentum is taken by the imbalance in the number of particles due to tunneling across the weak link. Such a system is typically referred to as a Josephson pendulum.

In ultracold degenerate atomic gases, several realizations of the Josephson pendulum have been studied Leggett (2001); Paraoanu et al. (2001); Smerzi et al. (1997); Marino et al. (1999). While the superfluid-fermion case Paraoanu et al. (2002); Heikkinen et al. (2010) still awaits experimental realization, the bosonic-gas version has been already demonstrated Albiez et al. (2005); Levy et al. (2007). Also in this case two regimes have been identified: small Josephson oscillations, corresponding to the low-energy limit case described here, and the macroscopic self-trapping regime Smerzi et al. (1997); Marino et al. (1999), corresponding to the free-rotor situation. Another example is an oscillating $LC$ electrical circuit with a nonlinear inductance realized as a tunnel barrier between two superconducting leads. This is the case of the transmon circuit Koch et al. (2007), which is currently one of the most promising approaches to quantum processing of information, with high-fidelity operations and good prospects for scalability. Its two lowest eigenstates are close to those of a harmonic oscillator, with only weak perturbations caused by the anharmonicity of the potential. The weak anharmonicity also guarantees that the lowest states of the transmon are immune to charge noise, which is a major source of decoherence in superconducting quantum circuits.

In this paper we consider a paradigmatic model which arises when the Josephson pendulum is interacting with a resonator. Circuit quantum electrodynamics offers a rigorous embodyment of the above model as a transmon device coupled to a superconducting resonator - fabricated either as a three-dimensional cavity or as a coplanar waveguide segment. In this realization, the system is driven by an external field of variable frequency, and dissipation affects both the transmon and the resonator.

We study in detail the onset of nonlinearity in the driven-dissipative phase transition between the quantum and the classical regimes. We further compare the photon number in detail with the corresponding transmon occupation, and demonstrate that the onset of nonlinearities is accompanied by the excitation of all bound states of the transmon and, thus, it is sensitive to the transmon truncation. We also find that the onset of the nonlinearities is sensitive to the energy level structure of the transmon, e.g. on the gate charge which affects the eigenenergies near and outside the edge of the transmon potential. The results also show that the full classical treatment is justified only in the high-amplitude regime, yielding significant discrepancies in the low-amplitude regime - where the phenomenology is governed by photon blockade. This means that the system undergoes a genuine quantum-to-classical phase transition. Calculations of the nth order correlation functions support this conclusion. These correlations are very small in the blockaded regime and increase beyond unity in the classical one. Moreover, the effects of the counter-rotating (Bloch–Siegert) terms in the full multilevel Josephson–Rabi system are clearly observable. Overall, our numerical simulations demonstrate that the multi-photon blockade phenomenon is qualitatively different for a realistic multilevel anharmonic system compared to the Jaynes–Cummings case studied extensively in the literature.

To introduce dissipation we use two models, namely the conventional Lindblad master equation and the Floquet–Born–Markov master equation, which is developed especially to capture the effects of the drive on the dissipation. We show that both yield relatively close results. However, we emphasize that the Floquet–Born–Markov approach should be preferred because its numerical implementation is considerably more efficient than that of the corresponding Lindblad equation.

While our motivation is to elucidate fundamental physics, in the burgeoning field of quantum technologies several applications of our results can be envisioned. For example, the single-photon blockade can be employed to realize single-photon sources, and the two-photon blockade can be utilized to produce transistors controlled by a single photon and filters that yield correlated two-photon outputs from an input of random photons Kubanek et al. (2008). In the field of quantum simulations, the Jaynes–Cummings model can be mapped onto the Dirac electron in an electromagnetic field, with the coupling and the drive amplitude corresponding respectively to the magnetic and electric field: then, the photon blockade regime is associated with a discrete spectrum of the Dirac equation, while the breakdown of the blockade corresponds to a continuous spectrum Gutiérrez-Jáuregui and Carmichael (2018). Finally, the switching behavior of the pendulum in the transition region can be used for designing bifurcation amplifiers for the single-shot nondissipative readout of qubits Vijay et al. (2009).

The paper is organized as follows. In Section II we introduce the electrical circuit which realizes the pendulum-oscillator system, calculate its eigenenergies, identify the two dynamical regimes of the small oscillations and the free rotor. In Section III, we introduce the drive and dissipation. We discuss two formalisms for dissipation, namely the Lindblad equation and the Floquet–Born–Markov approach. Section IV presents the main results for the quantum-to-classical transition and the photon blockade, focusing on the resonant case. Here, we also discuss the gate dependence, the ultra-strong coupling regime, and we obtain the nth order photon correlations. Section V is dedicated to conclusions.

II Circuit-QED implementation of a Josephson pendulum coupled to a resonator

We discuss here the physical realization of the Josephson pendulum-resonator system as an electrical circuit consisting of a transmon device coupled capacitively to an $LC$ oscillator, as depicted in Fig. 1(a). The coupled system is modeled by the Hamiltonian

[TABLE]

where

[TABLE]

describe the resonator, the transmon, and their coupling, respectively. We have defined $\hat{a}$ as the annihilation operator of the harmonic oscillator, and used $\hat{n}=-i\partial/\partial\varphi$ as the conjugate momentum operator of the superconducting phase difference $\hat{\varphi}$. These operators obey the canonical commutation relation $[\hat{\varphi},\hat{n}]=i$. The angular frequency of the resonator is given by $\omega_{\rm r}$. We have also denoted the Josephson energy with $E_{\rm J}$, and the charging energy with $E_{\rm C}=e^{2}/(2C_{\Sigma})$ where the capacitance on the superconducting island of the transmon is given as $C_{\Sigma}=C_{\rm B}+C_{\rm J}+C_{\rm g}$. Using the circuit diagram in Fig. 1(a), we obtain the coupling constant $g=2eC_{\rm g}q_{\rm zp}/(\hbar C_{\Sigma}C_{\rm r})$, where the zero-point fluctuation amplitude of the oscillator charge is denoted with $q_{\rm zp}=\sqrt{C_{\rm r}\hbar\omega_{\rm r}/2}$ Koch et al. (2007).

Let us briefly discuss the two components of this system: the Josephson pendulum and the resonator. The pendulum physics is realized by the superconducting transmon circuit Koch et al. (2007) in Fig. 1(a) and described by the Hamiltonian $\hat{H}_{\rm t}$ in Eq. (3). As discussed in Ref. Koch et al. (2007), the Hamiltonian of the transmon is analogous to that of an electrically charged particle whose motion is restricted to a circular orbit and subjected to homogeneous and perpendicular gravitational and magnetic fields. By fixing the x and z directions as those of the gravity and the magnetic field, respectively, the position of the particle is completely determined by the motion along the polar angle in the xy plane. The polar angle can be identified as the $\varphi$ coordinate of the pendulum. Thus, the kinetic energy part of the Hamiltonian (3) describes a free rotor in a homogeneous magnetic field. In the symmetric gauge, the vector potential of the field imposes an effective constant shift for the $\varphi$ component of the momentum, which is analogous to the offset charge $n_{\textrm{g}}$ on the superconducting island induced either by the environment or by a gate voltage. In the following, the ’plasma’ frequency for the transmon is given by $\omega_{\rm p}=\sqrt{8E_{\rm C}E_{\rm J}}/\hbar$ and describes the classical oscillations of the linearized circuit. The parameter $\eta=E_{\rm J}/E_{\rm C}$ is the ratio between the potential and kinetic energy of the pendulum, and determines, by the condition $\eta\gg 1$, whether the device is in the charge-insensitive regime of the free rotor and the gravitational potential.

The eigenvalues $\{\hbar\omega_{k}\}$ and the corresponding eigenvectors $\{|k\rangle\}$, with $k=0,1,\ldots$, of the Hamiltonian in Eq. (3) can be obtained by solving the Mathieu equation, see Appendix A. In general, the eigenvalues of the coupled system Hamiltonian $\hat{H}_{0}$ in Eq. (1) have to be solved numerically. With a sufficient truncation in the Hilbert spaces of the uncoupled systems in Eq. (2) and Eq. (3), one can represent the Hamiltonian $\hat{H}_{0}$ in a matrix form. The resulting eigenvalues of the truncated Hamiltonian $\hat{H}_{0}$ are shown in Fig. 1. We see that the coupling creates avoided crossings at the locations where the pendulum transition frequencies are equal to positive integer multiples of the resonator quantum. Also, the density of states increases drastically with the energy. The nonlinearity in the system is characterized by the non-equidistant spacings between the energy levels. Their origin is the sinusoidal Josephson potential of the transmon. Here, we are especially interested in the regime where the resonator frequency is (nearly) resonant with the frequency of the lowest transition of the pendulum, i.e. when $\omega_{\rm r}\approx\omega_{\rm q}=\omega_{01}=\omega_{1}-\omega_{0}$.

The Hamiltonian $\hat{H}_{0}$ in Eq. (1) can be represented in the eigenbasis $\{|k\rangle\}$ of the Josephson pendulum as

[TABLE]

We refer to this as the Josephson–Rabi Hamiltonian. Here, $K$ is the number of transmon states included in the truncation. We have also defined $\hat{\Pi}_{k\ell}\equiv\langle k|\hat{n}|\ell\rangle|k\rangle\langle\ell|$ which is the representation of the Cooper-pair-number operator in the eigenbasis of the transmon.

A useful classification of the eigenstates can be obtained by using the fact that the transmon can be approximated as a weakly anharmonic oscillator Koch et al. (2007), thus $\langle k|\hat{n}|\ell\rangle$ is negligible if $k$ and $\ell$ differ by more than 1. Together with the rotating-wave approximation, this results in

[TABLE]

Here, we introduce the total excitation-number operator as

[TABLE]

which commutes with the Hamiltonian in Eq. (6). Thus, the eigenstates of this Hamiltonian can be labeled by the eigevalues of $\hat{N}$, which is a representation that we will find useful in the discussion of the transitions between these states. The terms neglected in the rotating-wave approximation can be treated as small perturbations except for transitions where the coupling frequency $g_{\ell k}=g\langle k|\hat{n}|\ell\rangle$ becomes a considerable fraction of the corresponding transition frequency $\omega_{\ell k}=\omega_{k}-\omega_{\ell}$ and, thus, enters the ultrastrong coupling regime with $g_{k\ell}\geq 0.1\times\omega_{\ell k}$. In the ultrastrong coupling regime and beyond, the eigenstates are superpositions of states with different excitation numbers and cannot, thus, anymore be labeled with $N$.

Another important approximation for the Hamiltonian in Eq.(5) is the two-state truncation ($K=2$), which reduces it to the Rabi Hamiltonian

[TABLE]

Here $g_{01}=g\langle 1|\hat{n}|0\rangle$, the qubit annihilation operator is $\hat{\sigma}_{-}=|0\rangle\langle 1|$, and the Pauli spin matrix $\hat{\sigma}_{\rm x}=\hat{\sigma}_{-}+\hat{\sigma}_{+}$. The Rabi Hamiltonian is a good approximation to the pendulum-oscillator system as long as the corrections for the low-energy eigenvalues and eigenstates, arising from the higher excited states of the pendulum, are taken properly into account Boissonneault et al. (2009, 2012a, 2012b).

Further, by performing a rotating-wave approximation, we obtain the standard Jaynes–Cummings model

[TABLE]

which also results from a truncation of Eq. (6) to the low-energy subspace spanned by the lowest two eigenstates of the transmon. Apart from the non-degenerate ground state $|0,0\rangle$ with zero energy, the excited-state eigenenergies of the Jaynes–Cummings Hamiltonian in Eq. (9) form a characteristic doublet structure. In the resonant case, the excited-state eigenenergies and the corresponding eigenstates are given by

[TABLE]

Here, $n_{r}=1,2,\ldots$ and we have denoted eigenstates of the uncoupled Jaynes–Cummings Hamiltonian with $\{|n_{r},0\rangle,|n_{r},1\rangle\}$ where $|n_{r}\rangle$ are the eigenstates of the resonator with $n_{r}=0,1,\ldots$.

Due to the rotating-wave approximation, the Jaynes–Cummings Hamiltonian commutes with the excitation-number operator in Eq. (7) truncated to two states and represented as

[TABLE]

Thus, they have joint eigenstates and, in addition, the excitation number $N$ is a conserved quantity. For a doublet with given $n_{r}$, the eigenvalue of the excitation-number operator is $N=n_{r}$, while for the ground state $N=0$. We note that the transition energies between the Jaynes–Cummings eigenstates depend nonlinearly on $N$. Especially, the transition energies from the ground state $|0,0\rangle$ to the eigenstate $|n_{r},\pm\rangle$ are given by $n_{r}\hbar\omega_{\rm r}\pm\sqrt{n_{r}}\hbar g_{01}$.

III Models for the driven-dissipative Josephson pendulum coupled to the harmonic oscillator

Here, we provide a master equation approach that incorporates the effects of the drive and dissipation to the coupled system. Previous studies on this system have typically truncated the transmon to the low-energy subspace spanned by the two lowest energy eigenstates Bishop et al. (2010); Reed et al. (2010), or treated the dissipation in the conventional Lindblad formalism Bishop et al. (2009). Recent studies Pietikäinen et al. (2017, 2018); Verney et al. (2019); Lescanne et al. (2019) have treated the dissipation at the detuned limit using the Floquet–Born–Markov approach. We will apply a similar formalism for the case where the pendulum and resonator are in resonance in the low-energy subspace. Especially, we study the driven-dissipative transition between the low-energy and the free rotor regimes of the pendulum in terms of the dependence of the number $N_{\rm r}$ of quanta in the resonator on the drive power.

III.1 Coupling to the drive

The system shown in Fig. 1 and described by the Hamiltonian in Eq. (1) can be excited by coupling the resonator to a monochromatic driving signal modeled with the Hamiltonian

[TABLE]

where $A$ and $\omega_{\rm d}$ are the amplitude and the angular frequency of the drive, respectively. This results in a total system Hamiltonian $\hat{H}_{\rm S}=\hat{H}_{0}+\hat{H}_{\rm d}$. For low-amplitude drive, only the first two states of the pendulum have a significant occupation and, thus, the Hamiltonian $\hat{H}_{0}$ can be truncated into the form of the well-known Rabi Hamiltonian in Eq. (8), which in turn, under the rotating-wave approximation, yields the standard Jaynes–Cummings Hamiltonian in Eq. (9). The transitions induced by the drive in the Jaynes–Cummings system are subjected to a selection rule – the occupation number can change only by one, i.e. $N\rightarrow N\pm 1$. This follows from the relations

[TABLE]

As a consequence, the system climbs up the Jaynes–Cummings ladder by one step at a time. Particularly, a system in the ground state is coupled directly only to states $|1,\pm\rangle$. Indeed, in such a system the Jaynes–Cummings ladder has been observed Fink et al. (2008), as well as the effect of strong drive in the off-resonant Pietikäinen et al. (2017) and on-resonant case Fink et al. (2017). The Jaynes–Cummings model offers a good starting point for understanding the phenomenon of photon blockade in the pendulum-resonator system, which will be discussed later in detail. Indeed, it is apparent from Eq. (10) that, as the system is driven externally by not too intense fields, the excitation to higher levels in the resonator is suppressed by the higher levels being off-resonant, due to the nonlinearity induced by the coupling. This is referred to as photon blockade. As the drive amplitude increases further, the entire Jaynes–Cummings hierarchy breaks down Carmichael (2015).

However, in weakly anharmonic systems such as the transmon, as the drive amplitude is increased, the higher excited states of the Josephson pendulum become occupied and the two-state approximation becomes insufficient. As a consequence, the system has to be modeled by a multilevel Josephson–Rabi model Pietikäinen et al. (2017, 2018); Lescanne et al. (2019); Verney et al. (2019), see Eq. (5), which includes the higher excitations of the transmon as well as the counter-rotating (Bloch–Siegert) terms. In the resonant case, the need to take into account the second excited state of the transmon has been pointed out already in Ref. Fink et al. (2017). Here, at larger drive amplitudes, the pendulum escapes the low-energy subspace defined by the states localized in a well of the cosine potential and the unbound free rotor states also become occupied Pietikäinen et al. (2017, 2018); Lescanne et al. (2019); Verney et al. (2019) even in the case of strongly detuned drive frequency. In the limit of very high drive power, the pendulum behaves as a free rotor and the nonlinear potential can be neglected. Consequently, the resonance frequency of the system is set by the bare resonator frequency, instead of the normal modes.

III.2 Dissipative coupling

The dissipation is treated by modeling the environment as a thermal bosonic bath which is coupled bilinearly to the resonator. The Hamiltonian of the driven system coupled to the bath can be written as

[TABLE]

where

[TABLE]

Above, $\{\hat{b}_{k}\}$, $\{\Omega_{k}\}$, and $\{g_{k}\}$ are the annihilation operators, the angular frequencies, and the coupling frequencies of the bath oscillators. We use this model in the derivation of a master equation for the reduced density operator of the system. We proceed in the conventional way and assume the factorized initial state $\hat{\rho}(0)=\hat{\rho}_{\rm S}(0)\otimes\hat{\rho}_{\rm B}(0)$, apply the standard Born and Markov approximations, trace over the bath, and perform the secular approximation. As a result, we obtain a master equation in the standard Lindblad form.

III.3 Lindblad master equation

Conventionally, the dissipation in the circuit QED setup has been treated using independent Lindblad dissipators for the resonator and for the pendulum. Formally, this can be achieved by coupling the pendulum to another heat bath formed by an infinite set of harmonic oscillators. This interaction can be described with the Hamiltonian

[TABLE]

where $\{f_{k}\}$ and $\{\hat{c}_{k}\}$ are the coupling frequencies and the annihilation operators of the bath oscillators. The bath is coupled to the transmon through the charge operator $\hat{n}$ which is the typical source of decoherence in the charge-based superconducting qubit realizations. By following the typical Born–Markov derivation of the master equation for the uncoupled subsystems, one obtains a Lindblad equation where the dissipators induce transitions between the eigenstates of the uncoupled ($g=0$) system Breuer and Petruccione (2002); Scala et al. (2007); Beaudoin et al. (2011); Tuorila et al. (2017)

[TABLE]

where $\mathcal{L}[\hat{A}]\hat{\rho}=\frac{1}{2}(2\hat{A}\hat{\rho}\hat{A}^{{\dagger}}-\hat{A}^{{\dagger}}\hat{A}\hat{\rho}-\hat{\rho}\hat{A}^{{\dagger}}\hat{A})$ is the Lindblad superoperator and $n_{\rm th}(\omega)=1/[e^{\hbar\omega/(k_{\rm B}T)}-1]$ is the Bose–Einstein occupation. Note that the treatment of dissipation as superoperators acting separately on the qubit and on the resonator is valid if their coupling strength and the drive amplitude are weak compared to the transition frequencies of the uncoupled system-environment. Above, we have also assumed an ohmic spectrum for the resonator bath.

In the Lindblad master equation (20), we have included the effects arising from the coupling $g$ and the drive into the coherent von Neumann part of the dynamics. The first two incoherent terms cause transitions between the eigenstates of the resonator and arise from the interaction Hamiltonian in Eq. (18). The strength of this interaction is characterized with the spontaneous emission rate $\kappa$. The last term describes the relaxation, excitation, and dephasing of the transmon caused by the interaction Hamiltonian in Eq. (19). The transition rates $\Gamma_{k\ell}$ between the transmon eigenstates follow the Fermi’s golden rule as

[TABLE]

In our numerical implementation, we have assumed that the fluctuations of the transmon bath can also be characterised with an ohmic spectrum $S(\omega)=\frac{\gamma_{0}\omega}{1-\exp[-\hbar\omega/k_{\rm B}T]}$, where $\gamma_{0}$ is a dimensionless factor describing the bath-coupling strength. We have also denoted the transition frequencies of the transmon with $\omega_{k\ell}=\omega_{\ell}-\omega_{k}$.

Here, the magnitude of the transition rate from state $|k\rangle$ to the state $|\ell\rangle$ is given by the corresponding matrix element of the coupling operator $\hat{n}$ and the coupling strength $\gamma_{0}$. We note that in a typical superconducting resonator-transmon realization one has $\gamma=\gamma_{0}\omega_{01}\ll\kappa$. In this so-called bad-cavity limit, the effects of the transmon bath are negligible especially if the coupling frequency $g$ with the resonator is large. Thus, the main contribution of the transmon dissipators in the master equation Eq. (20) is that it results to faster convergence in the numerical implementation of the dynamics.

III.4 Floquet–Born–Markov formalism

The dissipators in the Lindblad model above are derived under the assumption of weak driving and weak coupling between the transmon and the resonator. However, both the driving and the coupling affect the eigenstates of the system and, thus, have to be taken into account in the derivation of the master equation. This can be achieved in the so-called Floquet–Born–Markov approach, where the drive and the transmon-resonator coupling are explicitly included throughout the derivation of the dissipators Tuorila et al. (2013); Pietikäinen et al. (2017, 2018); Lescanne et al. (2019); Verney et al. (2019). For this purpose, we represent the system in terms of the quasienergy states which can be obtained only numerically.

Since the drive in Eq. (13) is $\tau=2\pi/\omega_{\rm d}$-periodic, the solution to the time-dependent Schrödinger equation

[TABLE]

corresponding to the Hamiltonian $\hat{H}_{\rm S}(t)$ in Eq. (16), can be written in the form

[TABLE]

where $\varepsilon$ are the quasienergies and $|\Phi(t)\rangle$ are the corresponding $\tau$-periodic quasienergy states. By defining the unitary time-propagator as

[TABLE]

one can rewrite the Schrödinger equation (22) in the form

[TABLE]

Using Eqs. (23) and (24), we obtain

[TABLE]

from which the quasienergies $\varepsilon_{\alpha}$ and the corresponding quasienergy states $|\Phi_{\alpha}(0)\rangle$ can be solved. Using the propagator $\hat{U}$, one can obtain the quasienergy states for all times from

[TABLE]

Due to the periodicity of $|\Phi_{\alpha}(t)\rangle$, it is sufficient to find the quasienergy states for the time interval $t\in[0,\tau]$. Also, if $\varepsilon_{\alpha}$ is a solution for Eq. (26), then $\varepsilon_{\alpha}+\ell\hbar\omega_{\rm d}$ is also a solution. Indeed, all solutions of Eq. (26) can be obtained from the solutions of a single energy interval of $\hbar\omega_{\rm d}$. These energy intervals are called Brillouin zones, in analogy with the terminology used in solid-state physics for periodic potentials.

The master equation for the density operator in the quasienergy basis can be written as Blümel et al. (1991); Grifoni and Hänggi (1998)

[TABLE]

where

[TABLE]

Above, $\theta(\omega)$ is the Heaviside step-function and $\hbar\Delta_{\alpha\beta\ell}=\varepsilon_{\alpha}-\varepsilon_{\beta}+\ell\hbar\omega_{\rm d}$ is the energy difference between the states $\alpha$ and $\beta$ in Brillouin zones separated by $\ell$. Also,

[TABLE]

where $t_{0}$ is some initial time after the system has reached a steady state.

From Eq. (28), we obtain the occupation probabilities $p_{\alpha}=\rho_{\alpha\alpha}(t\rightarrow\infty)$ in the steady state as

[TABLE]

and the photon number

[TABLE]

where

[TABLE]

is the photon number in a single quasienergy state. The occupation probability for the transmon state $|k\rangle$ is given by

[TABLE]

We emphasize that this method assumes weak coupling to the bath but no such restrictions are made for the drive and pendulum-resonator coupling strengths. As a consequence, the dissipators induce transitions between the quasienergy states of the driven coupled system.

III.5 Parameters

The parameter space is spanned by seven independent parameters which are shown in Table 1.

We fix the values of the energy ratio $\eta=E_{\rm J}/E_{\rm C}$ and the coupling strengths $g$ and $\kappa$. The ratio $\eta$ sets the number $K_{\rm b}$ of bound states in the pendulum, see Appendix A, but does not qualitatively affect the response. We have used a moderate value of $\eta$ in the transmon regime, in order to keep $K_{\rm b}$ low allowing more elaborate discussion of the transient effects between the low-energy oscillator and rotor limits. We use the Born, Markov, and secular approximations in the description of dissipation which means that the value of $\kappa$ has to be smaller than the system frequencies. In addition, we work in the experimentally relevant strong coupling regime where the oscillator-pendulum coupling $g\gg\kappa$. The choice of parameters is similar to the recently realized circuit with the same geometry Pietikäinen et al. (2017).

The transition energies of the transmon are determined by the Josephson energy $E_{\rm J}$ and by the charging energy $E_{\rm C}$, which can be adjusted by the design of the shunting capacitor $C_{\rm B}$, see Fig. 1. The transition energy between the lowest two energy eigenstates is given by $\hbar\omega_{\rm q}\approx\sqrt{8E_{\rm J}E_{\rm C}}-E_{\rm C}=E_{\rm C}(\sqrt{8\eta}-1)$. We will study the onset of the nonlinearities for different drive detunings $\delta_{\rm d}=\omega_{\rm d}-\omega_{\rm r}$ as a function of the drive amplitude $A$. We are especially interested in the resonant case $\delta_{\rm q}=\omega_{\rm q}-\omega_{\rm r}=0$. The detuned case has been previously studied in more detail in Refs. Pietikäinen et al. (2017, 2018); Lescanne et al. (2019); Verney et al. (2019). We have used a temperature value of $k_{\rm B}T/(\hbar\omega_{\rm r})=0.13$ which corresponds to $T\approx 30$ mK for a transmon with $\omega_{\rm q}/(2\pi)=5$ GHz.

IV Numerical results

IV.1 Classical system

Classically, we can understand the behaviour of our system as follows: the pendulum-resonator forms a coupled system, whose normal modes can be obtained. However, because the pendulum is nonlinear, the normal-mode frequencies of the coupled-system depend on the oscillation amplitude of the pendulum. The resonator acts also as a filter for the drive, which is thus applied to the pendulum. As the oscillation amplitude of the pendulum increases, the normal-mode frequency shifts, an effect which is responsible for photon blockade. Eventually the pendulum reaches the free rotor regime, where the Josephson energy becomes negligible. As a consequence, the nonlinearity no longer plays any role, and the resulting eigenmode of the system is that of the bare resonator.

We first solve the classical equation of motion (see Appendix B) for the driven and damped resonator-transmon system. We study the steady-state occupation $N_{\rm r}$ of the resonator as a function of the drive amplitude. Classically, one expects that the coupling to the transmon causes deviations from the bare resonator occupation

[TABLE]

We emphasize that $N_{\rm bare}\leq A^{2}/\kappa^{2}$ where the equality is obtained if the drive is in resonance, i.e. if $\delta_{\rm d}=0$. The numerical data for $\delta_{\rm d}/\omega_{\rm r}=-0.02$ is shown in Fig. 2. We compare the numerical data against the bare-resonator photon number in Eq. (35), and against the photon number of the linearized system, see Appendix B,

[TABLE]

where $\delta_{\rm p}=\omega_{\rm d}-\omega_{\rm p}$, $\hbar\omega_{\rm p}=\sqrt{8E_{\rm J}E_{\rm C}}$, $g_{\rm eff}=g\sqrt[4]{\eta/32}$, and $\gamma$ is the dissipation rate of the pendulum. The above result is obtained by linearizing the pendulum potential which results to system that is equivalent to two coupled harmonic oscillators. We find in Fig. 2 that for small drive amplitude $A/\kappa=0.005$, the steady state of the resonator photon number is given by that of the linearized system. As a consequence, both degrees of freedom oscillate at the drive frequency and the system is classically stable. The small deviation between the numerical and analytic steady-state values is caused by the rotating-wave approximations that were made for the coupling and the drive in the derivation of Eq. (36).

If the drive amplitude is increased to $A/\kappa=7$, the nonlinearities caused by the cosinusoidal Josephson potential generate chaotic behavior in the pendulum. As a consequence, the photon number does not find a steady state but, instead, displays aperiodic chaotic oscillations around some fixed value between those of the bare resonator and the linearized system. This value can be found by studying the long-time average in the steady state. For very high drive amplitude $A/\kappa=500$, the photon number in the classical system is given by that of the bare resonator in Eq. (35). Physically this means that for strong driving, the pendulum experiences rapid free rotations and, as a consequence, its contribution to the photon dynamics is zero on average.

In Fig. 3(a), we study in more detail how the classical steady-state photon number of the resonator changes as the coupled system goes through the transition between the linearized oscillations in the weak driving regime and the bare-resonator oscillations for strong driving. In the absence of driving the steady-state photon number is zero in accordance to Eq. (36). For low drive amplitudes, the resonator-transmon system can be approximated as a driven and damped Duffing oscillator. We show in Appendix B that the system has one stable solution for drive amplitudes $A<A_{\rm min}$ and $A>A_{\rm crit}$, and two stable solutions for $A_{\rm min}<A<A_{\rm crit}$ where

[TABLE]

where we have defined the renormalized oscillator frequency and transmon dissipation rate as $\tilde{\omega}_{\rm r}^{2}=\omega_{\rm r}^{2}-g^{2}\hbar\omega_{\rm r}/(4E_{\rm C})$ and $\tilde{\gamma}=\gamma+gg_{1}\kappa\omega_{\rm d}^{2}/[(\tilde{\omega}_{\rm r}^{2}-\omega_{\rm d}^{2})^{2}+\kappa^{2}\omega_{\rm d}^{2}]$, respectively, the classical oscillation frequency of the linearized transmon as $\hbar\omega_{\rm p}=\sqrt{8E_{\rm J}E_{\rm C}}$, and the renormalized linearized transmon frequency as $\tilde{\omega}_{\rm p}^{2}=\omega_{\rm p}^{2}-g^{2}\omega_{\rm d}^{2}/(\tilde{\omega}_{\rm r}^{2}-\omega_{\rm d}^{2})[\hbar\omega_{\rm r}/(4E_{\rm C})]$.

For amplitudes $A<A_{\rm min}$, the classical system behaves as a two-oscillator system and the photon number has the typical quadratic dependence on the drive amplitude in Eq. (36). As the drive amplitude becomes larger than $A_{\rm min}$ deviations from the linearized model emerge. In addition, the system becomes bistable. If $A\approx A_{\rm crit}$, the number of stable solutions for the Duffing oscillator is reduced from two to one. This is displayed by the abrupt step in the photon number of the classical solution in Fig. 3 around $A_{\rm crit}/\kappa=1.2$. The remaining high-amplitude stable solution appears as a plateau which reaches up to the drive amplitude $A/\kappa\approx 5.6$. If the drive amplitude is further increased, the higher order terms beyond the Duffing approximation render the motion of the classical system chaotic, as described already in Fig. 2. For large drives, the classical photon number approaches asymptotically the photon number of the bare resonator.

IV.2 Quantum description

The transition between the motion of linearized and bare-resonator oscillations is characteristic to oscillator-pendulum systems. However, we show here that in the quantum mechanical context, the onset of the nonlinear dynamical behaviour turns out to be quantitatively different from that provided by the above classical model. This was also observed in recent experimental realization with superconducting circuits Pietikäinen et al. (2017).

In the quantum-mechanical treatment, we calculate the steady-state photon number in the resonator as a function of the drive amplitude using the Floquet–Born–Markov master equation presented in Sec. III.4. We have confirmed that for the used values of the drive amplitude the simulation has converged for the truncation of seven transmon states and 60 resonator states. We compare the quantum results against those given by the classical equation of motion and study also deviations from the results obtained with the two-state truncation of the transmon. In Fig. 3, we present the results corresponding to gate charge $n_{\rm g}=0$, where the resonator, the transmon, and the drive are nearly resonant at low drive amplitudes. The used parameters are the same as in Fig. 2 and listed in Table 1.

First, we notice in Fig. 3(a) that even in the absence of driving there always exists a finite photon occupation of $N_{\rm r}\approx 10^{-3}$ in the ground state, contrary to the classical solution which approaches zero. At zero temperature, the existence of these ground-state photons Lolli et al. (2015) originates from the terms in the interaction Hamiltonian that do not conserve the number of excitations and are neglected in the rotating-wave approximation resulting in Eq. (6). For the two-state truncation of the transmon, one can derive a simple analytic result for the ground-state photon number by treating these terms as a small perturbation. In the second order in the perturbation parameter $g/\omega_{\rm r}$, one obtains that the number of ground-state photons is given by $N_{\rm r}\approx(g/2\omega_{\rm r})^{2}$. We have confirmed that our simulated photon number at zero driving is in accordance with this analytic result if $T=0$ and $g/\omega_{\rm r}\ll 1$. The photon number at zero driving obtained in Fig. 3(a) is slightly higher due to additional thermal excitation - in the simulations we use a finite value for temperature, see Table 1.

As was discussed in the previous section, the resonator photon number of a classical system increases quadratically with the drive amplitude. For amplitudes $A<A_{\rm crit}$, the classical system can be approximated with a linearized model formed by two coupled harmonic oscillators. However, in the quantum case the energy levels are discrete and, thus, the system responds only to a drive which is close to resonance with one of the transitions. In addition, the energy levels have non-equidistant separations which leads to a reduction of the photon number compared to the corresponding classical case, referred to as the photon blockade. This is also apparent in Fig. 3(a).

We emphasize that the photon-blockade is quantitatively strongly-dependent on the transmon truncation. This can be seen as the deviation between the two and seven state truncation results for $A/\kappa>1$ in Fig. 3(a). We further demonstrate this by showing the transmon occupations $P_{\rm k}$ in Fig. 3(b). For weak drive amplitudes, the transmon stays in its ground state. The excitation of the two-level system is accompanied by excitations of the transmon to several its bound states. If $A/\kappa\geq 30$, the transmon escapes its potential well and also the free rotor states start to gain a finite occupation. This can be interpreted as a transition between the Duffing oscillator and free rotor limits of the transmon, see Appendix A. As a consequence, the response of the quantum system resembles its classical counterpart. We will study the photon blockade in more detail in the following section.

Order parameter

In order to characterize the transition between the quantum and classical regime, we can also study the behaviour of the order parameter $\Xi$ defined as the expectation value of the coupling part of the Hamiltonian in Eq. (5), normalized with $\hbar g$, as

[TABLE]

previously introduced and used for the off-resonant case in Ref. Pietikäinen et al. (2017). To get an understanding of its behavior, let us evaluate it for the resonant Jaynes–Cummings model,

[TABLE]

therefore it correctly estimates the absolute value of the cavity field operator. At the same time, when applied to the full Rabi model, it includes the effect of the terms that do not conserve the excitation number.

In Fig. 3(c), we present $\Xi$ as a function of the drive amplitude $A$. Much like in the off-resonant case, this order parameter displays a marqued increase by one order of magnitude across the transition region.

IV.3 Photon blockade: dependence on the drive frequency

Here, we discuss in more detail the phenomenon of photon blockade in the pendulum-resonator system as a function of the drive detuning $\delta_{\rm d}=\omega_{\rm d}-\omega_{\rm r}$. First, we consider the transition between the ground state and the state $|n_{r},\pm\rangle$ [Eq. (11)] of the resonant Jaynes–Cummings system ($\omega_{\rm q}=\omega_{\rm r}$). We recall that the selection rules Eqs. (14) and (15) allow only direct transitions that change the excitation number by one. However, at higher amplitudes the probability of higher order processes is no longer negligible and excited states can be populated by virtual non-resonant single-photon transitions. As a consequence, one obtains the resonance condition for multi-photon transitions as $n_{r}\omega_{\rm d}=n_{r}\omega_{\rm r}\pm\sqrt{n_{r}}g$. Because the energy-level structure is non-equidistant, the drive couples only weakly to other transitions in the system. In the absence of dissipation, the dynamics of the Jaynes–Cummings system can, thus, be approximated in a subspace spanned by the states $\{|0,0\rangle,|n_{r},\pm\rangle\}$.

Thus, one expects that, due to the driving, the system goes through $n_{r}$-photon Rabi oscillations between the basis states of the subspace. The Rabi frequency $\Omega_{n_{r},\pm}$ of such process is proportional to the corresponding matrix element of the driving Hamiltonian in Eq. (13) and the drive amplitude $A$. Consequently, the time-averaged photon number in the system is $N_{\rm r}=(n_{r}-\frac{1}{2})/2$. The driving does not, however, lead into a further increase of the photon number either because the drive is not resonant with transitions from the state $|n_{r},\pm\rangle$ to higher excited states or the matrix element of the resonant transitions are negligibly small. We are referring to this phenomenon as $n_{r}$-photon blockade.

Dissipation modifies somewhat this picture, as it causes transitions outside the resonantly driven subspace. As a consequence, the average photon number decays with a rate which is proportional to $\kappa$. Thus, the steady state of such system is determined by the competition between the excitation and relaxation processes caused by the drive and the dissipation, respectively. At low temperatures, the occupation in the ground state becomes more pronounced as the dissipation causes mostly downward transitions. Thus, the steady-state photon number is reduced compared to the time-averaged result for Rabi-driven non-dissipative transition. This was visible already in Fig. 3 in which the data was obtained with the two-state truncation and corresponds to the 4-photon blockade of the Jaynes–Cummings system.

The diagram in Fig. 4(a) represents the eigenenergies of the Hamiltonian for Eq. (5) in the two-state truncation for the transmon. The states are classified according to the excitation number $N$ from Eq. (12). We note that, here, we do not make a rotating-wave approximation and strictly speaking $N$ is, therefore, not a good quantum number. However, it still provides a useful classification of the states since the coupling frequency is relatively small, i.e. $g/\omega_{\rm r}=0.04$.

In Fig. 4(c), we show the photon blockade spectrum of the resonator-transmon system as a function of the drive detuning $\delta_{\rm d}$, obtained numerically with the Floquet–Born–Markov master equation. Here, one can clearly identify the one-photon blockade at the locations where the drive frequency is in resonance with the single-photon transition frequency of the resonator-transmon system Bishop et al. (2009), i.e. when $\delta_{\rm d}=\pm g$. Two, three, and higher-order blockades occur at smaller detunings and higher drive amplitudes, similar to Ref. Carmichael (2015). Transitions involving up to five drive photons are denoted in the diagram in Fig. 4(a) and are vertically aligned with the corresponding blockades in Fig. 4(c). At zero detuning, there is no excitation as the coupling to the transmon shifts the energy levels of the resonator so that there is no transition corresponding with the energy $\hbar\omega_{\rm r}$. We also note that the photon-number spectrum is symmetric with respect to the drive detuning $\delta_{\rm d}$. We see this same symmetry also in Eq. (36) for the linearized classical system when the classical linearized frequency of the transmon is in resonance with the resonator frequency, i.e. when $\omega_{\rm p}=\omega_{\rm r}$.

However, in experimentally relevant realisations of such systems the higher excited states have a considerable quantitative influence to the photon-number spectrum. We demonstrate this by showing data for the seven-state transmon truncation in Figs. 4(b) and (d). The eigenenergies shown in Fig. 4(b) are those obtained in Fig. 1 at resonance ($\omega_{\rm r}=\omega_{\rm q}$). We have again confirmed that for our choice of drive amplitudes and other parameters, this truncation is sufficient to obtain converged results with the Floquet–Born–Markov master equation. We observe that the inclusion of the higher excited states changes considerably the observed photon number spectrum. However, the states can again be labeled by the excitation number $N$ which we have confirmed by numerically calculating $N=\langle\hat{N}\rangle$ for all states shown in Fig. 4(c). The relative difference from whole integers is less than one percent for each shown state. Corresponding to each $N$, the energy diagram forms blocks containing $N+1$ eigenstates with (nearly) the same excitation number, similar to the doublet structure of the Jaynes–Cummings model. Contrary to the two-state case, these blocks start to overlap if $N>4$ for our set of parameters, as can be seen in Fig. 4(b).

The number of transitions that are visible for our range of drive frequencies and amplitudes in Fig. 4(d) is, thus, increased from ten observed in the Jaynes–Cummings case to 15 in the seven-state system. However, some of these transitions are not visible for our range amplitudes due to the fact that the corresponding virtual one-photon transitions are not resonant and/or have small transition matrix elements. In addition, the spectrum is asymmetric with respect to the detuning as the multi-photon resonances are shifted towards larger values of $\delta_{\rm d}$. As a consequence, the break-down of the photon blockade at $\delta_{\rm d}=0$ occurs at much lower amplitudes as is observed in the Jaynes–Cummings system Carmichael (2015).

IV.4 Approaching the ultrastrong coupling

For most applications in quantum information processing a relative coupling strength $g/\omega_{\rm r}$ of a few percent is sufficient. However, recent experiments with superconducting circuits have demonstrated that it is possible to increase this coupling into the ultrastrong regime ($g/\omega_{\rm r}\sim 0.1-1$) and even further in the deep strong coupling regime ($g/\omega_{\rm r}\geq 1$) Forn-Díaz et al. (2019); Gu et al. (2017); Kockum et al. (2019). While the highest values have been obtained so far with flux qubits, vacuum-gap transmon devices with a similar electrical circuit as in Fig. 1(a) can reach $g/\omega_{\rm r}=0.07$ Bosman et al. (2017a) and $g/\omega_{\rm r}=0.19$ Bosman et al. (2017b). To study the multilevel Josephson–Rabi model at higher couplings, we present in Fig. 5 results for the average number of photons in the resonator for couplings $g/\omega_{\rm r}=0.04,0.06$, and $0.1$, employing the Floquet–Born–Markov approach to dissipation.

At low drive powers the two-level approximation can be used for the transmon, and the Josephson pendulum-resonator system maps into the quantum Rabi model. From Fig. 5 we see that the average number of photons $N_{\rm r}$ in the resonator is not zero even in the ground state; this number clearly increases as the coupling gets stronger. As noted also before, this is indeed a feature of the quantum Rabi physics: differently from the Jaynes–Cummings model where the ground state contains zero photons, the terms that do not conserve the excitation number in $\hat{H}_{\rm c}$ lead to a ground state which is a superposition of transmon and resonator states with non-zero number of excitations. Similar to Sec. IV.2, the perturbative formula $N_{\rm r}\approx(g/2\omega_{\rm r})^{2}$ approximates very well the average number of photons at zero temperature, while in Fig. 5 we observe slightly higher values due to the finite temperature. As the drive increases, we observe that the photon blockade tends to be more effective for large $g$’s. Interestingly, the transition to a classical state also occurs more abruptly as the coupling gets stronger. We have checked that this coincides with many of the upper levels of the transmon being rapidly populated. Due to this effect, the truncation to seven states (which is the maximum that our code can handle in a reasonable amount of time) becomes less reliable and artefacts such as the sharp resonances at some values start to appear.

IV.5 Correlation functions

Photon correlations provide the essential information about the statistics of the quantum field leaking out of the cavity Glauber (1963); Bozyigit et al. (2011); Hamsen et al. (2017). Especially interesting is the zero delay correlation function Fink et al. (2018); Bin et al. (2018). The $n$th order correlation function is defined as

[TABLE]

where $G^{(1)}=N_{\rm r}$. In Fig. 6 we have calculated numerically the correlation functions for the seven state transmon (a) with the rotating-wave approximation and (b) including the counter-rotating (Bloch–Siegert) terms present in the Josephson–Rabi model.

At zero driving and zero temperature, the system is in its ground state. If the rotating-wave approximation is used, this state is $|0,0\rangle$. This means that all $G^{(n)}=0$. This can be seen in Fig. 6(a). Without the rotating-wave approximation, the analytic form of the ground state has not been found even for the two-level Rabi Hamiltonian. We obtain an approximative ground state of the Rabi Hamiltonian by treating the counter-rotating terms $\propto(\hat{a}^{{\dagger}}\hat{\sigma}_{+}+\hat{a}\hat{\sigma}_{-})$ as a small perturbation. We obtain in the second order of the perturbation parameter $g_{01}/\omega_{\rm r}$ that the non-normalized ground state can be written as

[TABLE]

This gives $G^{(1)}=3.7\times 10^{-4}$ and $G^{(2)}=5.2\times 10^{-7}$ with the parameters used here. The corresponding values with seven state Josephson–Rabi Hamiltonian from Fig. 6(b) are $G^{(1)}=3.6\times 10^{-4}$ and $G^{(2)}=7.9\times 10^{-7}$, demonstrating very good agreement. As seen from Fig. 6, raising the temperature reduces the difference between the two cases, since now the system contains thermally-activated photons even at zero drive amplitude. Overall, as expected, in the quantum regime the $G^{(n)}$’s are very small, but once the system approaches the classical regime as $A/\kappa\gg 1$ the values of $G^{(n)}$ become larger than one. The marked difference between the rotating-wave approximation result shown in Fig. 6(a) and that of the full Josephson-Rabi model in Fig. 6(b) suggest the use of statistics as a detection method for the ultrastrong-coupling regime.

IV.6 Dependence on the gate charge

If the transmon is only weakly nonlinear, i.e. $\eta\gg 1$, its lowest bound eigenstates are insensitive to the gate charge, see Appendix A. As a consequence, one expects that the value of the gate charge should not affect the photon-number response to a weak drive. However, as the amplitude of the drive is increased, the higher excited states of the transmon become occupied, as discussed in the context of Fig. 3. Especially, the transition region between the quantum and classical responses should be dependent on the gate-charge dispersion of the transmon states. We demonstrate this in Fig. 7, where we show the simulation data for the gate-charge values $n_{\rm g}=0$ and $n_{\rm g}=0.5$. Clearly, in the weak driving regime, the responses for the two gate-charge values are nearly equal. The deviation in the photon number is of the order of $10^{-3}$, which is explained by our rather modest value of $\eta=30$.

The deviations between the photon numbers of the two gate-charge values are notable if $A/\kappa=10\ldots 20$. In this regime, the transmon escapes the subspace spanned by the two lowest eigenstates and, thus, the solutions obtained with different gate-charge numbers are expected to differ. At very high amplitudes, the free-rotor states with $k\geq 6$ also begin to contribute the dynamics. These states have a considerable gate-charge dispersion but, however, the superconducting phase is delocalized. Accordingly, the gate-charge dependence is smeared by the free rotations of the phase degree of freedom.

We also note that the photon number response displays two sharp peaks for $n_{\rm g}=0.5$ at $A/\kappa\approx 13$ and $A/\kappa\approx 25$. The locations of the peaks are very sensitive to the value of the gate charge, i.e. to the energy level structure of the transmon. Similar abrupt changes in the transmon occupation were also observed in recent experiments in Ref. Lescanne et al. (2019). They could be related to quantum chaotic motion of the system recently discussed in Ref. Mourik et al. (2018). In this parameter regime, also the Jaynes–Cummings model displays bistability Vukics et al. (2019).

IV.7 Comparison between different master equations

We have also compared our numerical Floquet–Born–Markov method against the Lindblad master equation which was presented in Sec. III.3 and has been conventionally used in the studies of similar strongly driven systems with weak dissipation. We note that for the case of strong coupling to the bath, the possible treatment is by the method of path integrals, as developed by Feynman and Vernon, which has already been applied to describe the dynamics of the Rabi model Henriet et al. (2014).

We recall that in the Lindblad formalism, the environment induces transitions between the non-driven states of the system, whereas in the Floquet–Born–Markov approach the dissipation couples to the drive-dressed states of the system. Thus, one expects deviations from the Floquet–Born–Markov results in the limit of strong driving. In Fig. 8 we show a comparison between the two models in the two-state truncation approximation for the transmon. We see that the largest differences between the models appear when the transition from the quantum to classical response starts to emerge, see Fig. 3. Based on our numerical calculations, the differences are the largest at resonance and both models give equivalent results whenever one of the three frequencies, $\omega_{\rm r}$, $\omega_{\rm q}$, or $\omega_{\rm d}$, is detuned from the other two.

We emphasize, however, that computationally the Floquet–Born–Markov master equation is by two orders of magnitude more efficient than the corresponding Lindblad equation. Moreover, in the case of Fig. 4(d) the computing time of the Floquet–Born–Markov equation was roughly a week with an ordinary CPU. In such cases, the solution of the Lindblad equation becomes impractical and one should use a parallelized implementation of the Floquet–Born–Markov master equation.

V Conclusions

We have given a comprehensive treatment of the driven-dissipative quantum-to-classical phase transition for a Josephson pendulum coupled to a resonator, going beyond the truncated Rabi form of the Hamiltonian through the full inclusion of the higher energy level of the pendulum. We refered to this as the multilevel Josephson–Rabi model. We modelled the open quantum system with the Floquet–Born–Markov method, in which the dissipative transitions occur between the drive-dressed states of the system. We compared our results also against those given by the conventional Lindblad formalism where the dissipation couples to the eigenstates of the non-driven system.

We found that the quantitative description of the multi-photon blockade phenomenon and of the nonlinearities associated with the phase transition in this system requires a systematic inclusion of the higher energy levels of the transmon and a proper model for dissipation. We also studied approximate classical models for this system, and showed that the discrete energy structure of the quantum system suppresses the classical chaotic motion of the quantum pendulum. Indeed, while the classical solution predicts a sudden change between the low and high amplitude solutions, the quantum solution displays a continuous transition from the normal-mode oscillations to the freely rotating pendulum regime. Finally, we analyzed in detail the two models of dissipation and demonstrated that they produce slightly different predictions for the onset of the photon blockade.

VI Acknowledgments

We thank D. Angelakis, S. Laine, and M. Silveri for useful discusssions. We would like to acknowledge financial support from the Academy of Finland under its Centre of Excellence program (projects 312296, 312057, 312298, 312300) and the Finnish Cultural Foundation. This work uses the facilities of the Low Temperature Laboratory (part of OtaNano) at Aalto University.

Appendix A The eigenvalue problem for the Josephson pendulum

The energy eigenstates of the pendulum can be solved from the Mathieu equation Baker and Blackburn (2005); Cottet (2002); Abramovitz and Stegun (1972) which produces a spectrum with bound and free-particle parts. The high-energy unbound states are given by the doubly-degenerate quantum rotor states, which are also the eigenstates of the angular momentum operator.

In analogy with the elimination of the vector potential by a gauge transformation as usually done for a particle in magnetic field, one can remove the dependence on $n_{\rm g}$ from the transmon Hamiltonian in Eq. (3), i.e.

[TABLE]

with the gauge transformation $\hat{U}\hat{H}_{\rm t}\hat{U}^{{\dagger}}$, where

[TABLE]

As a consequence, the eigenstates $|k\rangle$ of the Hamiltonian are modified into

[TABLE]

The transformed Hamiltonian can be written as

[TABLE]

Here, we represent the (Schrödinger) eigenvalue equation for the transformed Hamiltonian in the eigenbasis of the operator $\hat{\varphi}$. As a result, the energy levels of the transmon can be obtained from the Mathieu equation Baker and Blackburn (2005); Cottet (2002); Abramovitz and Stegun (1972)

[TABLE]

where $z=\varphi/2$, $q=-\eta/2=-E_{\rm J}/(2E_{\rm C})$, and $a=E_{k}/E_{\rm C}$. We have also denoted the transformed eigenstate $|k\rangle$ in the $\varphi$ representation with $\psi_{k}(\varphi)=\langle\varphi|e^{-in_{\rm g}\hat{\varphi}}|k\rangle$. Note that $\Psi_{k}(\varphi)=e^{in_{\rm g}\varphi}\psi_{k}(\varphi)$ is the eigenfunction of the original Hamiltonian in Eq. (3). Due to the periodic boundary conditions, one has that $\Psi_{k}(\varphi+2\pi)=\Psi_{k}(\varphi)$. The solutions to Eq. (47) are generally Mathieu functions which have a power series representation, but cannot be written in terms of elementary functions Abramovitz and Stegun (1972). However, the corresponding energy-level structure can be studied analytically in the high and low-energy limits.

In Fig. 9, we present the eigenenergies $E_{k}$ obtained as solutions of the Mathieu equation (47). The eigenstates that lie within the wells formed by the cosine potential are localized in the coordinate $\varphi$, whereas the states far above are (nearly) evenly distributed, see Fig. 9(a). As a consequence, the high-energy states are localized in the charge basis. The data shows that if plotted as a function of the gate charge, the states inside the cosine potential are nearly flat, see Fig. 9(b). This implies that such levels are immune to gate charge fluctuations, which results in a high coherence of the device. Outside the well, the energy dispersion with respect to the gate charge becomes significant, and leads to the formation of a band structure typical for periodic potentials Marder (2000).

A.1 High-energy limit: Free rotor

If the energy in the system is very high due to, e.g., strong driving, the Josephson energy can be neglected and the transmon behaves as a free particle rotating in a planar circular orbit, which can be described solely by its angular momentum $\hat{L}_{\rm z}=\hat{n}$. Since the angular momentum is a good quantum number, the eigenenergies and the corresponding eigenfunctions are given by

[TABLE]

where $k=0,\pm 1,\pm 2,\ldots$. We note that if the magnetic field is zero ($n_{g}=0$), the nonzero free rotor energies are doubly degenerate. The level spacing is not constant but increases with increasing $k$ as Baker and Blackburn (2005)

[TABLE]

In Fig. 9, we show the eigenenergies calculated with Eq. (48). Clearly, with large energies outside the potential, the energy spectrum of the particle starts to resemble that of the free rotor. Also, the eigenfunctions of the free rotor are plane waves in the $\varphi$ eigenbasis, yielding a flat probability density as a function of $\varphi$. On the other hand, in the momentum eigenbasis, the free rotor states are fully localized.

A.2 Low-energy limit: Duffing oscillator

If the pendulum energy is very low, the superconducting phase of the transmon is localized near $\varphi\approx 0$. Thus, the cosine potential can be approximated with the first terms of its Taylor expansion. Consequently, the transmon Hamiltonian reduces to that of a harmonic oscillator with an additional quartic potential

[TABLE]

This is the Hamiltonian operator of the quantum Duffing model. The Duffing model has received a considerable attention in the recent literature Peano and Thorwart (2006); Serban and Wilhelm (2007); Verso and Ankerhold (2010); Vierheilig and Grifoni (2010); DiVincenzo and Smolin (2012); Everitt et al. (2005) especially in the context of superconducting transmon realizations. It is worthwhile to notice that in this regime the potential is no longer periodic and, thus, we can neglect the periodic boundary condition of the wavefunction. As a consequence, the eigenenergies and eigenfunctions are not dependent on the offset charge $n_{\rm g}$.

If $\eta=E_{J}/E_{C}\gg 1$, the quartic term is small and one can solve the eigenvalues and the corresponding eigenvectors perturbatively up to the first order in $\eta$. This regime in which the Josephson energy dominates over the charging energy is referred to as the transmon limit. One, therefore, obtains the eigenenergies

[TABLE]

where $k=0,1,2,\ldots$. Especially, the transition energy between the two lowest Duffing oscillator states can be written as

[TABLE]

This becomes accurate as $\eta\rightarrow\infty$. The anharmonicity of a nonlinear oscillator is typically characterized in terms of the absolute and relative anharmonicity, which are defined, respectively, as

[TABLE]

where $E_{ij}=E_{j}-E_{i}$ and the latter approximations are valid in the transmon limit. We emphasize that in the low-excitation limit the transmon oscillates coherently with frequency $\omega_{\rm q}\approx\omega_{\rm p}-E_{\rm C}/\hbar$. Thus, in the quantum pendulum the nonlinearity is present even in the zero-point energy, whereas the small-amplitude oscillations in classical pendulum occur at the angular (plasma) frequency $\hbar\omega_{\rm p}=\sqrt{8E_{\rm J}E_{\rm C}}$.

In Fig. (9), we compare the eigenenergies (51) of the Duffing model obtained with the perturbation theory against the exact solutions of the Mathieu equation (47). We see that in the low-energy subspace the perturbed-Duffing solution reproduces very well the full Mathieu results. For the higher excited states, the momentum dispersion starts to play a dominant role and deviations arise as expected. This starts to occur close to the boundary of the potential. One can estimate the number $K_{\rm b}$ of bound states by requiring $E_{K_{\rm b}-1}\approx E_{\rm J}$ in Eq. (51). This implies that the number of states within the potential scales with $\eta\gg 1$ as

[TABLE]

For the device with parameters listed in Table 1, one has that $\sqrt{\eta}\approx 5$, and the above estimate gives $K_{\rm b}\approx 5$. This coincides with the number of bound states extracted from the numerically exact spectrum of the eigenenergies depicted in Fig. 9.

Appendix B Driven and damped classical system

The classical behaviour of the uncoupled pendulum has been extensively studied in the literature Dykman and Smelyanski (1990); Dykman (2012); Baker and Blackburn (2005). If the driving force is not too strong, one can approximate the pendulum with a Duffing oscillator with a quartic non-linearity, as shown in the previous appendix. The main feature of such an oscillator is the bistability of its dynamics. Namely, in a certain range of drive amplitudes and frequency detunings between the driving signal and the oscillator, two stable solutions with low and high amplitudes of the oscillations are possible. If one gradually increases the driving, the pendulum suddenly jumps from the low to the high amplitude solution at the critical driving strength, at which the low amplitude solution vanishes. In the bistable region the Duffing oscillator may switch between the two solutions if one includes noise into the model Dykman and Smelyanski (1990). This complicated dynamics has been observed in a classical Josephson junction Siddiqi et al. (2004, 2005).

However, differently from the previous work mentioned above, in our setup the pendulum is coupled to a resonator and driven only indirectly. Here we develop the classical theory of the coupled system and show that the basic physics of bistability is present as well. We first linearize the equations of motion and then introduce systematically the corrections due to the nonlinearity. The system Hamiltonian $\hat{H}_{\rm S}=\hat{H}_{0}+\hat{H}_{\rm d}$, defined by Eq. (1) and Eq. (13), can be written in terms of the circuit variables as

[TABLE]

Above, we have denoted the capacitance and inductance of the $LC$ resonator with $C_{\rm r}$ and $L_{\rm r}$, respectively, the effective coupling with $\tilde{g}=\hbar g/q_{\rm zp}$, the effective drive with $\tilde{A}=\hbar A/q_{\rm zp}$, and the zero-point fluctuations with $\phi_{\rm zp}=\sqrt{\hbar/(2C_{\rm r}\omega_{\rm r})}$ and $q_{\rm zp}=\sqrt{C_{\rm r}\hbar\omega_{\rm r}/2}$. Also, the resonance frequency of the bare resonator is defined as $\omega_{\rm r}=1/\sqrt{L_{\rm r}C_{\rm r}}$.

The corresponding equations of motion for the expectation values of the dimensionless operators $\hat{\phi}_{\rm r}=\hat{\phi}/\phi_{\rm zp}$ and $\hat{q}_{\rm r}=\hat{q}/q_{\rm zp}$ can be written as

[TABLE]

where we have denoted the expectation value of operator $\hat{x}$ as $\langle\hat{x}\rangle\equiv x$, applied the commutation relations $[\hat{\phi}_{\rm r},\hat{q}_{\rm r}]=2i$ and $[\hat{\varphi},\hat{n}]=i$, and defined the phenomenological damping constants $\kappa$ and $\gamma=\gamma_{0}\omega_{\rm q}$ for the oscillator and the pendulum, respectively. The exact solution to these equations of motion is unavoidably numerical and is given in Figs. 2 and 3. The resonator occupation is calculated as $N_{\rm r}=\frac{1}{4}(q_{\rm r}^{2}+\phi_{\rm r}^{2})$.

B.1 Solution of the linearized equation

We study Eqs. (56)-(59) in the limit of weak driving. As a consequence, one can linearize the equations of motion by writing $\sin\varphi\approx\varphi$. In addition, by defining

[TABLE]

we obtain

[TABLE]

where we have introduced an effective coupling as $g_{\rm eff}=g\sqrt[4]{\eta/32}$. The above equations describe two driven and dissipative coupled oscillators.

We assume that both oscillators are excited at the drive frequency, i.e. $\alpha=\alpha_{0}\exp(-i\omega_{\rm d}t)$ and $\beta=\beta_{0}\exp(-i\omega_{\rm d}t)$. By making a rotating-wave approximation for the coupling and the drive, we obtain the resonator occupation $N_{\rm lin}=|\alpha_{0}|^{2}$ in the steady state

[TABLE]

where $\delta_{\rm d}=\omega_{\rm d}-\omega_{\rm r}$ and $\delta_{\rm p}=\omega_{\rm d}-\omega_{\rm p}$, with $\omega_{\rm p}=\sqrt{8E_{\rm J}E_{\rm C}}/\hbar$. This appeared already in Eq. (36).

B.2 Correction due to the pendulum nonlinearity

Here, we study the nonlinear effects neglected in the above linearized calculation. We eliminate the variables $\phi_{\rm r}$ and $n$ from Eqs. (56)-(59) and obtain

[TABLE]

where we have denoted $g_{1}=g\hbar\omega_{\rm r}/(4E_{\rm C})$, and defined the renormalized resonator frequency as $\tilde{\omega}_{\rm r}^{2}=\omega_{\rm r}^{2}-g^{2}\hbar\omega_{\rm r}/(4E_{\rm C})$. In Eq. (64), we have included only the term that is proportional to $g^{2}$ as it provides the major contribution to the frequency renormalization, and neglected the other second order terms in $\kappa$, $\gamma$ and $g$ that lead to similar but considerable smaller effects. We write the solutions formally in terms of Fourier transform as

[TABLE]

where $q_{\rm r}[\Omega]$ and $\varphi[\Omega]$ are the (complex valued) Fourier coefficients of $q_{\rm r}(t)$ and $\varphi(t)$, respectively. As a consequence, one can write the equations of motion as

[TABLE]

We solve $q_{\rm r}[\Omega]$ from the first equation and obtain

[TABLE]

By replacing this result into Eq. (69), we obtain

[TABLE]

where we have neglected a constant phase factor. For weak drive amplitudes, $\varphi[\omega_{\rm d}]$ is the only non-zero Fourier component. Thus, one can evaluate the Fourier transform in the above equation at the drive frequency. Consequently, the Fourier component of the third term in the equation can be evaluated as

[TABLE]

where in the second term we have assumed that the dissipation is weak, i.e. $\kappa\ll\sqrt{\tilde{\omega}_{\rm r}^{2}-\omega_{\rm d}^{2}}$ and we taken into account the dominant terms for the real and imaginary parts. As a result, we obtain

[TABLE]

Here, we have defined the renormalized linear oscillation frequency $\tilde{\omega}_{\rm p}$, dissipation rate $\tilde{\gamma}$, and drive amplitude $B$ as

[TABLE]

where the two first equations are valid if $\kappa\ll\sqrt{\tilde{\omega}_{\rm r}^{2}-\omega_{\rm d}^{2}}$.

Thus, we have shown that in the limit of low dissipation, the classical resonator-transmon system can be modeled as a driven and damped pendulum. In the case of weak driving, we expand the sinusoidal term up to the third order in $\varphi$. We obtain the equation of motion for the driven and damped Duffing oscillator:

[TABLE]

This equation can be solved approximatively by applying a trial solution $\varphi(t)=\varphi_{1}\cos(\omega_{\rm d}t)$ into Eq. (77). By applying harmonic balance, and neglecting super-harmonic terms, we obtain a relation for the amplitude $\varphi_{1}$ in terms of the drive amplitude $B$. By taking a second power of this equation and, again, neglecting the super-harmonic terms, we obtain

[TABLE]

The above equation is cubic in $\varphi_{1}^{2}$. It has one real solution if the discriminant $D$ of the equation is negative, i.e. $D<0$. If $D>0$, the equation has three real solutions, two stable and one unstable. The stable solutions can appear only if $\omega_{\rm d}<\tilde{\omega}_{\rm p}$, which is typical for Duffing oscillators with a soft spring (negative nonlinearity). The bistability can, thus, occur for amplitudes $B_{\rm min}<B<B_{\rm crit}$ where the minimal and critical amplitudes $B_{\rm min}$ and $B_{\rm crit}$, respectively, determine the region of bistability and are obtained from the equation $D=0$. By expanding the resulting $B_{\rm min}$ and $B_{\rm crit}$ in terms of $\tilde{\gamma}$ and by taking into account the dominant terms, we find that

[TABLE]

where the last equality holds if $\delta_{\rm p}=\tilde{\omega}_{\rm p}-\omega_{\rm d}\ll\omega_{\rm p}$. The iterative numerical solution of Eq. (77) indicates that the initial state affects the switching location between the two stable solutions. We note that this approximation neglects all higher harmonics and, thus, cannot reproduce any traces towards chaotic motion inherent to the strongly driven pendulum.

Finally, we are able to write the minimal and critical drive amplitudes of the coupled resonator-transmon system using Eqs. (76), (79), and (80). We obtain [see Eq. (38)]

[TABLE]

Note that these equations are valid for $\kappa\ll\sqrt{\tilde{\omega}_{\rm q}^{2}-\omega_{\rm d}^{2}}$.