Sideband transitions in a two-mode Josephson circuit driven beyond the rotating wave approximation

Driving quantum systems periodically in time plays an essential role in the coherent control of quantum states. The rotating wave approximation (RWA) is a good approximation technique for weak and nearly-resonance driven fields. However, these experiments sometimes require large detuning and strong driving fields, for which the RWA may not hold. In this work, we experimentally, numerically, and analytically explore strongly driven two-mode Josephson circuits in the regime of strong driving and large detuning. Specifically, we investigate beam-splitter and two-mode squeezing interaction between the two modes induced by driving a two-photon sideband transition. Using numerical simulations, we observe that the RWA is unable to correctly capture the amplitude of the sideband transition rates. We verify this finding using an analytical model that is based on perturbative corrections. We find that the breakdown of the RWA in the regime studied does not lead to qualitatively different dynamics, but gives the same results as the RWA theory at higher drive strengths, enhancing the coupling rates compared to what one would predict. This is an interesting consequence compared to the carrier transition case, where the breakdown of the RWA results in qualitatively different time evolution of the quantum state. Our work provides an insight into the behavior of time-periodically driven systems beyond the RWA. We also provide a robust theoretical framework for including these findings in the calculation and calibration of quantum protocols in circuit quantum electrodynamics.


I. INTRODUCTION
Time-periodic driving is a prominent technique for the in − situ coherent control of quantum dynamical processes. However, exactly solving the quantum dynamics of the system with time-varying Hamiltonian is particularly difficult [1,2]. As long as the drive is weak enough and nearly resonant with the quantum state transition of a target observable, then the rotating wave approximation (RWA) may provide a good estimate of the dynamics [3,4]. It is, however, necessary to understand the physics of quantum systems beyond the RWA from both from a fundamental and a practical perspective. In the realm of faithful quantum information processing (QIP), the need for fast gates to suppress quantum errors requires drive strengths that could exceed those that are valid for the RWA. Motivated by this problem, theoretical effort has been directed to understand quantum driven systems beyond the RWA [1,2], along with many experimental works that have used strong driving across many different physical systems. [5][6][7][8][9][10][11][12]. Although the previous studies have explored driven systems beyond the RWA, they focus on a single driven mode and do not address the coupling of different degrees of freedom (e.g., the use of driving fields to induce sideband transitions between modes).
Sideband transitions ubiquitously appear in a variety of physical systems [15, 16, 18-22, 24, 25]. Driving systems with appropriately chosen frequencies can yield engineered interactions among different degrees of freedoms. To engineer a specific interaction, it is important to accurately estimate the transition rates. In many cases the driving parameters for sideband transitions typically satisfy the requirements of the RWA, such as trapped ions, cavity-optomechanics, and Raman transitions [20][21][22]24]. However, this is not always the case for the circuit quantum elecrodynamics (cQED) platform, one of the most promising QIP platforms in recent years, where a strong and far off resonant driving beyond the RWA is sometimes required [26]. Nonetheless, current approaches to quantitative analysis still rely on the application of the RWA. When the sideband driving frequencies are far off-resonant from the transition frequencies of the system, in which the conditions for the application of the RWA should not hold, we may not currently be possible to make reliable predictions of the transition rates.
In this paper, we study the sideband transition rates in a two-mode Josephson circuit that is induced by strong external time-periodical driving. The circuit comprises a transmon [13] that is dispersively coupled to a resonator mode. Specifically, we study beam splitter (BS) and twomode squeezing (TMS) interactions between each mode, which are the simplest forms of sideband transitions in these two-mode systems. For our device, the required driving parameters are close to (TMS coupling), or far beyond (BS coupling) the RWA regime. We confirm a simple relationship between the transition rates and frequency shifts, which explains the data in both regimes.
We perform numerical simulations to support our findings. We also derive an analytical perturbation expansion that goes beyond the RWA, which is validated by our numerical results. Our findings indicate that although the RWA is clearly violated, and significantly underestimates the mode frequency shifts and the sideband transition rates for a known driving strength, the breakdown of the RWA does not result in qualitatively different behaviour but instead its effects in our measurements can be re-produced by the RWA theory using a larger drive field. Although the confirmation of a breakdown of the RWA is only possible to observe experimentally in an accurate independent calibration of the drive field, our results show the importance of including counter-rotating terms for accurate calculations of the sideband transition rates.

II. THEORETICAL DESCRIPTION
We derive an analytical expression taking a similar approach in [17,18] but breaking the RWA. The total Hamiltonian on the lab frame is given by, Here, ω (0) t and ω (0) r are the resonant frequencies of each mode.α andβ are the mode destruction operators of the transmon and the resonator modes, respectively. χ t is a Duffing nonlinearity of the transmon mode. g is a transverse coupling between the transmon mode and the resonator mode. In addition toĤ (0) sys , there is the driving HamiltonianĤ where Ω d and ω d are the driving amplitude and frequency, respectively. The total HamiltonianĤ (0) tot is then given byĤ It is often useful to rewrite this Hamiltonian in the normal mode basis (the normal mode annihilation operators areâ andb): With typical circuit QED parameters, χ t is approximately the same as χ (0) t . χ r is the inherited Duffing nonlinearity to the resonator mode by the coupling g. In the dispersive coupling regime (|ω t − ω r | g),α inĤ d can be approximated byâ [18,27]. Then, the driving Hamiltonian can be approximated to beĤ The total Hamiltonian in the normal mode basis is then given byĤ (1) tot =Ĥ (1) sys +Ĥ (1) d . This can be perturbatively diagonalized by taking Schrieffer-Wolff (S-W) transformation [29]Û (t) = eŜ with an appropriate generator In this work, we treat beam splitter (âb † +â †b ) and two-mode squeezing (âb +â †b † ) interactions induced by two-photon driving. These appear with frequency matching conditions 2ω d ≈ |ω After taking Schrieffer-Wolff transformation, collect-ing only the original and relevant derived terms yields, sb is the interaction rate for both the BS and TMS interactions. δω (1) t , δω (1) r and Ω (1) sb can be expressed by, , In the low excitation limit, the total Hamiltonian can be reduced to, (5) ω t,r and A t,r,tr (A tr ≈ √ A t A r ) correspond to the transition frequencies and anharmonicities that we observe in the experiments. We can obtain A t,r,tr by numerically diagonalizing Eq. 2 [23]. The difference between A t,r,tr and χ t,t,tr is due to the off diagonal elements in Eq. 3.
The discrepancy between A t,r,tr and χ t,t,tr suggests that the off diagonal elements also affect the derived quantities after taking the (S-W) transform. We hereby invoke an assumption that the effects of the off diagonal terms can be captured by replacing χ t,t,tr with A t,r,tr . This assumption leads to a conclusion that δω (1) t,r and Ω (1) sb in Eq. 4 should be renormalized to δω t,r and Ω sb , , We provide the supporting information for this finding in Appendix D. When applying the RWA, they are given by, It is also interesting to investigate the case where only the counter-rotating terms inĤ d affect the system. In this case, the frequency shifts and sideband transition rates are given by, The detailed derivation is provided in Appendix A. If δω t is known, then replacing ∆ and Σ with ∆ + δω t and Σ + δω t will provide a more accurate estimate. It is worth pointing out here that many of the previous studies do not seriously distinguish between χ t,r,tr and A t,r,tr . However, the discrepancies between χ t,r,tr and A t,r,tr are sometimes significant, depending on the system's parameters. Renormalization of δω t,r and Ω sb is therefore of great importance for the accurate prediction of the frequency shifts and sideband transition rates. Eq. 6 and Eq. 7 suggests that the RWA significantly underestimates δω t,r and Ω sb when ∆ ∼ Σ but the ratios among them are identical, regardless of whether or not we use the RWA. It is also interesting to note that there is a correlation between the co-and counter-rotating terms in Eq. 6, which makes a significant contribution to the frequency shifts and sideband transition rates.

III. EXPERIMENT
Both the BS and the TMS interaction are schematically described in Fig. 1a and Fig. 1b. Two black-wavy arrows indicate the two-photon drive. Fig. 1c and Fig. 1d denote energy diagram descriptions. In all of the descriptions, the resonator and the transmon mode are colored green and magenta, respectively. In addition to the two-photon drive, we have a weak probe field (green) through the resonator mode to estimate Ω sb through the resonator's response. The decay rates of both modes are κ and γ, respectively. The energy levels of the resonator mode are denoted by |0 , |1 , |2 ,... and those of the transmon mode are denoted by |g , |e ,... Fig. 1e depicts a simplified circuit diagram of the system. We drive the transmon mode through a direct driveline and we probe the resonator mode through another feedline coupled to the resonator. Fig. 1f and Fig. 1g show how the probe transmission through the resonator varies with increasing Ω sb for both BS (f) and TMS (g) interaction. The curves are obtained by solving a numerical model based on Eq. 5 with dissipation operators. The decay rates of the resonator and transmon modes in the calculation are κ/2π ≈ 10.2 MHz and γ/2π ≈ 129 kHz. These parameters are similar to those in the experiment. Ω sb /2π is set by 2, 4, and 6 MHz in both BS and TMS interaction. The detailed information on the experimental setup and device is provided in Appendix C 1.
In the experiment, we deliberately design a large κ to facilitate the detection of the interactions through the resonator's transmission, even with small Ω sb . Our system satisfies the condition for electromagnetically induced transparency (EIT) [28] as long as Ω sb is smaller than |γ − κ|. In this regime, Ω sb and the other parameters independently shape the transparency window in the middle of the transmission spectrum of the resonator. Thereby, we extract Ω sb by fitting the resonator's transmission. The resonator's linewidth is overwhelmingly larger than the linewidth of the qubit and therefore the system is in the EIT condition as long as Ω sb is less than around 10 MHz. The observed ω t,r are 2π× 6.8112 and 4.0755 GHz, respectively.
The observed A t is 2π × 150 MHz and A r can be deduced by A tr ≈ √ A t A r . Since the resonator has a broad linewidth, we cannot simply extract A tr from the photon number by splitting the resonator or the transmon spectrum.
In Fig. 2, we present the procedure used for determining the frequency matching conditions. We define ω mat that satisfies 2ω mat = |ω t ± ω r | for both the BS and TMS interactions. In reality, the resonances undergo shifts, ω t,r → ω t,r = ω t,r + δω t,r and in our system we have δω t δω r ≈ 0. Thus, we have modified matching conditions, 2ω mat = |ω t ± ω r |. We swept the driving frequency ω d and find the condition ω d ≈ ω mat . We obtain the matching conditions when the transparency window is located at the center in the transmission spectrum. More quantitatively, ω mat can be obtained by extracting ω t when fitting the transmission data with numerical model given in Appendix A 3. Roughly, ω mat /2π ≈ 1. while keeping ω mat ≈ ω d . Solid curves are fits to the data based on numerical model. From these fits, we extracted Ω sb and δωt. The probe amplitude Ωp is 2π×130.6 kHz, except that Ωp is 10 times larger for the lowest dataset of (a). (c) Observed Ω sb with respect to corresponding δωt (circles). The solid line indicates a theory based on Eq. 6. Fitting errors in Ω sb are around 1% and are not plotted in the figures. and 5.44 GHz are expected for both the BS and the TMS interaction, respectively. For the BS interaction, ω mat is extremely far off-resonant (∆/Σ ≈ 0.6). This regime of the driving parameter has not been explored. Meanwhile, for TMS interaction, ω mat is relatively closer to the RWA regime (∆/Σ ≈ 0.11). In Fig. 3 a-b, we plot a portion of the transmission spectrum observed in the experiment. We scan the sideband driving power preserving the condition ω d ≈ ω mat . The solid curves are the fits based on the numerical model that we used in Fig. 1f-g. In the fitting process, the free parameters are Ω sb , γ and δω mat , while the other parameters are fixed. As we increase the driving amplitudes, we can readily see that the transparency windows be-have as expected from Fig. 1f-g. In Fig. 3c, we plot Ω sb with respect to the corresponding δω t , both of which are extracted from the fitting. The statistical errors in extracting Ω sb from the fitting are around only 1%, and thus not presented in the figures. We can find a linear correlation between δω t and Ω sb . The slope of the solid line is obtained from Eq. 6, with no free parameters. It is of note that both BS and TMS data lie on the same theoretical plot, although the driving parameters for each lives in distinct regimes.
To directly identify the breakdown of the RWA, we need to calibrate Ω d from an independent method not relying on the transmon frequency shifts. If we know the microwave power at the device (P d ), and the coupling rate between the transmon and drive line (γ ex ), then Ω d is simply given by P d γ ex /hω d . However, the uncertainty in the driveline attenuation sets a challenge. An error of only 1 dB in the attenuation induces a 10% error in Ω d , which is critical to our study. In future research, this challenge can be circumvented by using an additional 'sensor' qubit, as recently demonstrated in [50].

IV. NUMERICAL SIMULATION
We performed a comprehensive numerical analysis with the experimental conditions. We simulated the system's time domain dynamics by solving the dρ/dt = without including any dissipation.ρ is the density matrix of the system. As in the experiment, we swept the driving frequency for a given Ω d and find the frequency where a full oscillation takes place in transitions |e0 ⇐⇒ |g1 (BS) or |g0 ⇐⇒ |e1 (TMS). Ω sb is then given by the frequency of the oscillation. More detailed descriptions on the method of the numerical simulation are given in Appendix B.
In Fig. 4, we present the numerical calculation results (circles, triangles and squares) and corresponding analytical calculation results (solid, dashed, and dotted lines). The circles and solid lines refer to the results with thê H (0) sys +Ĥ (0) d . In the plots, triangles and squares refer to the simulation results dropping the counter-rotating and co-rotating driving terms inĤ (0) d . The analytical calculation is based on the Eq. 6. The dashed and dotted lines are obtained by Eq. 7 and Eq. 8, respectively.
In Fig. 4a-b, we compare the sideband transition rates obtained by the numerical simulations (circles, triangles, and squares) with the analytically calculated values (solid, dashed, and dotted lines). In Fig. 4 c-d, we present the frequency shifts of the transmon mode under the matching conditions for given the driving amplitudes in x-axes. The driving frequencies for each data points are set to satisfy the matching conditions for the given driving amplitudes. Although the RWA significantly distorts the Ω sb and δω t , the breakdown of the RWA is not visible in the Ω sb versus δω t relation, as seen in Fig. 4 e-f. The simulation data with the RWA perfectly lie on the data without the RWA. Therefore, a careful treatment is       required when estimating Ω d through Ω sb or δω t . Relying on the RWA results in significant overestimation of Ω d .

V. CONCLUSION
In summary, we performed the quantitative investigation of two-photon assisted four-wave interactions in a superconducting circuit. Over the entire range of the driving amplitudes in this work, our theoretical, numerical and experimental values agree with each other, which suggests that the faithful quantitative estimation of sideband transition rates is possible. This work expands our understanding in the strongly driven quantum systems. The findings through this work are not restricted to the system that we investigate here.

ACKNOWLEDGMENTS
We thank David Theron and Jochem Baselmans for providing us with NbTiN film. Byoung-moo Ann acknowledges support from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 722923 (OMT). This project also has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 828826 -Quromorphic. The data that support the findings of this study are available in [53]. Tab. I lists and defines the symbols that are used in the main text and supplemental material.

Finally,Ĥ
(1) tot can be expressed by, For given ω d , collecting the non-rotating terms at the transmon and resonator rotating frame in Eq. A2 yields Tab. II. We only list the terms that represent the interactions between different modes or the frequency shifts of each mode.

Modeling transmission spectrum
The resonator transmission spectrum is proportional to Tr t [ρ ssb ]. Here,ρ ss is a steady state density matrix of the transmon and resonator system, and Tr t indicates trace over the transmon states.ρ ss can be calculated based on the below Eq. A3, H p (t) = Ω p cos (ω p t) is the Hamiltonian of the prove field. D[Ô]ρ is defined by 2ÔρÔ † −Ô †Ôρ −ρÔ †Ô . κ is the decay rate of the resonator mode, and γ is that of the transmon mode. We neglect the pure dephasing rate of the transmon mode. Since we employ a single Josephson junction design, it is expected that the coherence time of the transmon mode is only limited to the decay time. For a steady state, we have dρss dt = 0, then we can calculatê ρ ss from Eq. A3.
Transmission spectrum is a function of a set of variables (ω d , ω p , Ω p , Ω sb , ω t , ω r , A t , A r , A tr κ, and γ).
Here, ω p is the independent variable in the fitting process. We fix κ, ω r , A t , A r , and A tr by the values we obtain from the independent measurement without driving field. These quantities are hardly shifted under the driving. ω d is given by the experiment. The free fitting parameters are ω t , Ω sb , Ω p , and γ. These quantities are extracted from the fitting process.
List of a portion of the non-rotating terms at the transmon and resonator rotating frame for given ω d derived from the fourth power term of Eq. A2.

Appendix B: Numerical simulations
In this section, we describe the detail procedures of the time-domain numerical simulations. The dynamics of the system are governed by the equation, dρ sys /dt = −i[Ĥ  Here,ρ sys is density matrix of the transmon and resonator. We do not take the dissipation into consideration in the time-domain dynamic simulations. Fig. 5a shows the simulated dynamics (blue line) when the driving frequency satisfies the matching condition for two-mode squeezing (TMS) interaction. The system parameters used in the simulation are the same with the experimental conditions. The sideband drive (Ω d (t), green line) is given as a pulse with 10ns of Gaussian rising and falling. The arrow indicates the length of the pulse. Fig. 5b shows the area enclosed by the dashed square in Fig. 5a. One can identify the qubit and resonator states significantly vary during the rising and falling duration of the sideband pulse. In Fig. 5c, we sweep the length of the sideband pulse and plot the states of the system at the end of the pulse. We obtain    Fig. 5d shows the area enclosed by the dashed square in Fig. 5c.
We sweep the driving frequency for each simulation data point and find the optimal frequency that yields the resonant sideband transitions. This procedure is described in Fig. 6. We chose the w d when the oscillation has a maximum contrast. We present the simulation data with different driving Hamiltonian in Fig. 7. The solid lines refer to the results with a full driving Hamiltonian containing both co-and counter-rotating terms. The dotted lines (dashed lines) are obtained by the simulations with only co-rotating (counter-rotating) terms in the driving Hamiltonian. See the caption for the detail conditions in the simulations.
Appendix C: Experimental methods

Experimental setup
An optical microscope image of the device is given in Fig. 8a. The device is comprised of a transmon and two co-planar waveguide resonators. The design of the device is the same as the one used in our previous work [51]. Only one of the resonators was used in this experiment. In addition, there is a drive line directly coupled to the transmon. The base layer of the circuit is fabricated from 100 nm niobium titanium nitride (NbTiN) film on a Silicon substrate. The detailed procedure to prepare the NbTiN film is described in [52]. The transmon is comprised of a Al-AlOx-Al Josephson junction and a finger capacitor. The transmon is not flux tuneable and there- fore the frequency is insensitive to the external magnetic field noise.
A cryogenic wiring diagram and measurement electronics are given in Fig. 8b. The device is mounted at the mixing chamber plate of a Bluefors LD-400 dilution fridge. The temperature of the plate is around 10 mK during the measurements. The device is enclosed within a cylindrical cooper shield to block the infrared radiation. To block the external magnetic fields, the copper can is enclosed by a Aluminum shield and two Mu-metal shields. The shields are not represented in the figure. We used a vector network analyzer (Keysight N5222A) to measure the resonator transmissions. An additional microwave source (Keysight N5183B) was used for sideband drivings. We used a non-dissipative low pass filter (Minicircuit VFL-3800+) in the drive line (third column).

Device parameter extraction
In this section, we provide the procedure to calibrate the cross-anharmonicity (A tr ) between the transmon and resonator modes in the experiment. We use the fact that the EIT transmission spectrum of the resonator depends on the A tr in the nonlinear response regime. In Fig. 9, we simulate the resonator's transmission spectrum with a beam splitter interaction (Ω sb /2π 1.2 MHz). The model that we used in the simulation is based on Eq. (3) in the main text including dissipation operators. In addition, we set δω mat /2π by -300 kHz. In the simulation, the linewidths of the resonator and transmon modes are the same with those in the experiment. We simulate in both linear response (Fig. 9a) and nonlinear response (Fig. 9b regime. Fig. 10a shows the measured resonator transmission spectrum while sweeping the probe power. P p is the resonator probe power measured at the output port of the vector network analyzer (VNA). Note that the contrast of the transparency window near the center decreases with increasing probe power. We first fit the resonator's transmission data in the linear response regime (solid line), setting Ω sb , ω d , ω t , κ and γ as free parameters. Then, we fit the data in the nonlinear response regime (dashed line) while fixing all the parameters obtained from the first fitting and only Ω p and A tr are free fitting parameters. When fitting the data in the linear response regime, we set A tr = 0 and Ω p /2π = 10 kHz. The choice of A tr can be justified since we already know A tr hardly affects the transmission in the linear response regime. The fitting results in both regimes are given in Fig. 10b and Fig. 10c. We obtain A tr /2π = 497 kHz and Ω p /2π = 4.35 MHz from the data in the nonlinear regime.
We can also obtain A tr from the fact that the resonator's transition frequency depends on the transmon's quantum states [37]. Fig. 11 shows how the resonator's transmission spectrum changes as we populate the transmon's first excited state. We drive the transmon mode with its resonant frequency and increase the power until we cannot see any further shift in the resonator's frequency. With this drive power, we can approximate the transmon's state 50:50 mixed state between the ground and first excited states. We observe a frequency shift of 520 kHz, which can be interpreted as A tr .
A tr extracted from Fig. 11 is slightly larger than the value obtained from Fig. 10. The discrepancy of the expected sideband transition rates based on both is about 2 percent. In the main text, we use A tr /2π = kHz obtained from Fig. 10. This approach is advantageous because we can extract the resonator probe power and A tr simultaneously, and consequently it guarantees more consistency. Comparison between the qubit decay rates (γ) extracted from the EIT spectrum fitting (dots) and the two-tone spectroscopy with a low probe and spectroscopy power (solid line).

Transmon decay rate analysis
In the fitting process to extract the sideband transition rates, the free fitting parameters other than Ω sb are γ and δω mat . We also present the extracted values for γ and δω mat in [53]. In this section, we especially focus on the γ. Fig. 12 shows the fitted γ (dots) with respect to corresponding Ω sb . These values are consistent with the γ from the low power two-tone spectroscopy (dashed line) in general. For BS interaction case, some data points far deviate from the dashed line. We attribute this to the undesired higher order sideband interactions. The matching frequency for BS interaction is close to the matching frequency for single-photon assisted sideband interaction between |e0 and |g2 . Since the resonator mode has a much larger decay rate, this undesired interaction can increase the effective decay rate of the transmon mode. The rightmost two data of TMS interaction case also far deviate from the solid line. We cannot find the systematic reason for the discrepancy. We could attribute this to the fluctuation of the transmon's decay rate with respect to time.

Appendix D: Additional analytical and numerical analysis
In this section, we confirm that Eq. 6 more accurately predicts the δω t and Ω sb than Eq. 4. In Fig. 13, we compare the analytical calculation based on Eq. 4 and numerical simulation results in Fig. 4. We can clearly see the discrepancy between the analytical and numerical results becomes larger than that in Fig. 4.
In Fig 14, we perform the additional simulation with various system parameters and compare the numerically simulated sideband transition rates (Ω sb -Sim) to the theoretical calculations (Ω sb -Th). We compare two different theoretical approaches based on Eq. 4 and Eq. 6, respec- tively. Aside from one case (Fig 14-d), Ω sb -Th based on Eq. 6 are closer to Ω sb -Sim. Even in Fig 14-d, Ω sb -Th based on Eq. 6 is more accurate with low driving amplitudes. * byoungmoo.ann@gmail.com