key: cord-0979664-yfe38bei
authors: Yin, Zelong; Li, Chunzhen; Allcock, Jonathan; Zheng, Yicong; Gu, Xiu; Dai, Maochun; Zhang, Shengyu; An, Shuoming
title: Shortcuts to Adiabaticity for Open Systems in Circuit Quantum Electrodynamics
date: 2021-07-18
journal: Nature communications
DOI: 10.1038/s41467-021-27900-6
sha: 39938013aef2cd660be415929998de8d12d00d5d
doc_id: 979664
cord_uid: yfe38bei

Shortcuts to adiabaticity (STA) are powerful quantum control methods, allowing quick evolution into target states of otherwise slow adiabatic dynamics. Such methods have widespread applications in quantum technologies, and various STA protocols have been demonstrated in closed systems. However, realizing STA for open quantum systems has presented a greater challenge, due to complex controls required in existing proposals. Here we present the first experimental demonstration of STA for open quantum systems, using a superconducting circuit QED system consisting of two coupled bosonic oscillators and a transmon qubit. By applying a counterdiabatic driving pulse, we reduce the adiabatic evolution time of a single lossy mode from 800 ns to 100 ns. In addition, we propose and implement an optimal control protocol to achieve fast and qubit-unconditional equilibrium of multiple lossy modes. Our results pave the way for accelerating dynamics of open quantum systems and have potential applications in designing fast open-system protocols of physical and interdisciplinary interest, such as accelerating bioengineering and chemical reaction dynamics.

Adiabatic processes -which preserve the nondegenerate instantaneous eigenstates of time-dependent Hamiltonians -have important applications in quantum technologies, including quantum simulation [1, 2] and computation [3, 4] . Though adiabatic evolution is, in principle, relatively robust against parameter fluctuations and environmental noise [5] , in the noisy intermediate-scale quantum era, decoherence is still an obstacle preventing its widespread application. Shortcut to adiabaticity (STA) addresses this issue by finding fast trajectories that connect the initial and final states of slow-paced adiabatic protocols. Since STA was first proposed [6] , it has found many applications, including in atom cooling, trapped atom and ion transportation [7, 8] , spin population transfer [9] [10] [11] [12] , implementing quantum logic gates [13, 14] , and quantum thermodynamics [15] . Due to freedom in choosing intermediate trajectories, time-dependent control parameters of a system can be adjusted in different ways, resulting in various STA protocols [16] . In particular, the method of counterdiabatic (CD) driving [17, 18] adds an auxiliary control H CD to the reference Hamiltonian to suppress unwanted diabatic transitions. This adiabatic-following feature makes CD driving robust to parameter errors [8] and suitable for fast holonomic gates [19] and efficient quantum heat engines [20] .

While STA finds widespread application in closed quantum systems, its generalization to open quantum systems is of fundamental interest. There are two strategies for this generalization: first, one can stick with STA designed for closed systems and mitigate environmental effects by utilizing redundant degrees of freedom [12] ; second, one can directly attempt to accelerate open system adiabatic dynamics. For open classical systems, a swiftequilibration protocol similar to STA was used to accel-erate the equilibration of a Brownian particle [21] . Recently, this idea was extended to the field of biology to guide the probability distribution of genotypes in a population along a specified path and time interval [22] . For open quantum systems, CD driving can be designed theoretically based on non-Hermitian Hamiltonians [23, 24] or Lindblad dynamics [25] . However, it remains challenging to conduct such experiments, as they often require complex controls such as engineered system-bath interactions [26] .

Here, we generalize STA to an open circuit QED (cQED) system consisting of multiple dissipative bosonic modes. When the time-dependent controls are varied sufficiently slowly [27] , the system evolves adiabatically within its time-dependent decoherence-free subspace (DFS) [28] . Analogously to the CD driving for closed quantum systems, we deduce the diabatic part of the Liouvillian which causes non-equilibrium transitions, and add a unitary control to counteract it. Consequently, we can enforce fast adiabatic evolution of the time-dependent DFS, i.e., a system initialized in the DFS remains so at all times [29] . However, when multiple lossy modes are present -as is common in many experimental scenarios -under the one-port driving of our setup (see below), only one hybrid mode at a time can be under perfect CD control.

To realise STA in a multi-mode setting, we further develop an analytical, open-loop control protocol, which we term multi-mode optimal control (MMOC). In MMOC, we make use of non-adiabatic dynamics during ringup(ringdown) to achieve the desired final equilibrium in a duration much less than that required by a slow varying adiabatic reference process (Fig. 2) . Utilising redundant degrees of freedom, we can also minimise a userdefined merit function, which we choose here to be the maximum driving amplitude to avoid undesired qubit excitations [30] . See Fig. 1a for a schematic overview of the system dynamics under the CD and MMOC procedures.

Our setup is illustrated in Fig. 1b and consists of two coupled resonator modes [31] , a and b, with coupling strength J/2π ≈ 10.5 MHz. Mode b is dispersively coupled to a transmon qubit [32] and mode a is coupled to a feedline at an effective temperature of 75 mK with strength κ a /2π ≈ 11.4 MHz, which serves as a cold bath. Given the dispersive coupling to the transmon qubit, the coupled modes a and b can be decoupled into qubitdependent hybrid modes (normal modes) a 0,1 and b 0,1 (see Methods). A coherent drive with amplitude (t) at frequency ω d /2π = 6.44025 GHz is generated by an arbitrary waveform generator and applied through the feedline. The system dynamics is inferred via time-traced output homodyne detection using the input-output theory as detailed in Supplementary Note (SN) 1.

We first implement CD driving designed for lossy hybrid mode b to achieve a fast ringup in a target duration t f , with the qubit kept in the ground state. For a chosen reference drive (t), the required additional CD control can be derived from the Lindblad dynamics or timedependent DFS of the system (see SN 2 and 3). The CD driving in the rotating frame of driving frequency ω d , including the reference (t) and the auxiliary control, takes the form:

where ∆ 0 r ≡ ω 0 r − ω d and κ 0 r are the frequency detuning and decay constant respectively, for the hybrid resonator mode b 0 and qubit state in |0 . For convenience, the superscripts are omitted in what follows. Note that in the dissipation-free limit κ r → 0, the closed system CD solution is readily recovered [8] . We characterise the performance of the CD driving in terms of the quantum speed limit for open quantum systems [33] and conclude that the CD driving has optimal quantum efficiency within the space of all available pulses (see SN 4).

As shown in Fig. 2 , we apply the reference driving (t) with a sin 2 -shape ringup, i.e. (t) = 0 sin 2 (πt/2t f ) before t f and (t) = 0 afterwards. This sin 2 -shape reference waveform is chosen to give a smooth, hardwarefriendly CD driving pulse. For both t f = 30 ns and 100 ns, we compare the performance of the reference driving (t) and CD driving cd (t). To demonstrate a speedup compared with the adiabatic timescale [34], estimated to be on the order of κ −1 r (see SN 3), we also apply an adiabatic sin 2 driving. Until t f = 800 ns κ −1 r , we observe a relatively good sin 2 -shape response, indicating adiabatic evolution. Our results show an equilibrium time of 100 ns for the CD driving and 350 ns ∼ 5κ −1 r for a quench driving. For t f = 30 ns, the large and rapidly varying CD pulse induces out-of-equilibrium excitation of the untargeted mode a . As a consequence, a requires an extra relaxation time of 75 ns ∼ 5κ −1 f to return to its equilibrium. For t f = 100 ns, equilibrium is achieved almost immediately after t f , and the non-equilibrium excitation of a is negligible. Since the sampled output signal is a coherent superposition of the input driving field and the leakage of the system (see SN 1 for details), the spiked IQ trajectory during the CD driving in Fig. 2c does not imply non-equilibrium dynamics. More data, corresponding to different t f , can be found in SN 9.

The CD drive Eq. 1 can only accelerate adiabatic dynamics for a single, qubit-state-dependent hybrid mode at a time and may excite other modes out of their instantaneous eigenstates. This issue can be mitigated with a relatively small driving-detuning ratio for the untargeted modes (Fig. 1c ). However, due to the time-energy uncertainty, when the protocol duration is reduced, the correspondingly larger drive amplitudes mean such unwanted excitations cannot be avoided entirely.

To eliminate the excitation of the untargeted mode a and remove the qubit-state dependence of the driving, we propose the following MMOC protocol. The MMOC pulse is based on a multi-section ansatz: the protocol duration is divided into m equal-spaced sections, each with a constant complex amplitude, i.e. (t j−1 < t < t j ) = j , j ∈ {1, 2, 3, . . . , m}. During the ringup stage, the control pulse is subject to the boundary conditions (0) = 0 and (t f ) = 0 . During the reset stage, the boundary conditions are reversed. To equilibrate four hybrid modes in time t f , we utilize the underlying Langevin dynamics and obtain four linear equilibrium-transfer equations connecting the pulse vector and the equilibrium difference y with a transfer matrix G as y = G · . If the complex hybrid detuning∆ ≡ ∆ − iκ/2 is found for each hybrid mode, G and y can be analytically derived. These equilibrium-transfer equations can be solved via singular value decomposition of G as long as m ≥ 4. The resulting takes the form = e + m−4 i=1 x i i , where the essential vector e can be analytically solved for, and i are m − 4 vectors orthogonal to e . See SN 5 for a detailed derivation. The extra degrees of freedom x i are chosen numerically to minimize the maximum amplitude component of . As a result, we can obtain a 5 dB reduction in the maximum amplitude compared with the essential vector. See SN 6 for more discussion on the performance (DFS) , which is precise only in the infinite time limit. To accelerate this adiabatic trajectory, a counterdiabatic (CD) driving can be applied to cancel non-adiabatic transitions, making the adiabatic approximation exact. Due to the effect of dissipation, CD driving for closed quantum systems has limited effectiveness, while its open system extension achieves arbitrary speed up for a single mode conditioned on the qubit state. To realize STA for multiple modes and independent of qubit states, we develop the multi-mode optimal control (MMOC) protocol, which includes non-adiabatic trajectories at intermediate times and reaches both final states ρ 0,1 f at a given target time. (b) Setup of the superconducting circuit. The resonator mode b is dispersively coupled to a transmon qubit with strength g/2π ≈ 80 MHz. The Purcell filter mode a is coupled to the resonator mode with strength J/2π = 10.5 MHz and to the feedline with strength κa/2π = 11.4 MHz. From residue qubit excitation, the temperature of the environment is measured to be 75 mK, justifying the cold-bath approximation of the feedline. Driving pulses are applied through the input field ci, and the output field r0 is amplified by an impedance-matched parametric amplifier (IMPA) and homodyne detected to infer the system dynamics. (c) Measured transmission spectrum S21 and hybrid-mode frequencies for both qubit states. The driving frequency ω d = 2π × 6.44025 GHz is allocated roughly in the middle of ω 0 r = 2π × 6.4427 GHz and ω 1 r = 2π × 6.4378 GHz, where the 0 or 1 in the upper right corner denotes the qubit state. The significant shift of resonant dips is partly due to the nonlinearity of the cavity when it is driven with large power. See Supplementary Note 1 for more details on the asymmetry and frequency shifts. and speed limit of this protocol. Fig. 3 shows our results for the fast equilibration of all four lossy hybrid normal modes of the coupled oscillator system. Choosing m = 10, we first apply a t f = 60 ns ringup pulse for fast system equilibration, and then transfer the system to the vacuum state at the end with another 60 ns reset pulse. According to the time-traced IQ trajectories (Fig. 3b , c and insets), different modes undergo different non-equilibrium dynamics and end up in their respective equilibrium states after the MMOC pulse. Our results show that fast unconditional multimode ringup and depletion are almost achieved at the desired time and reduce the duration required to achieve equilibrium compared with natural relaxation.

We have experimentally extended STA to an open quantum system. Our CD method accelerates adiabatic dynamics of the DFS for a single driven-dissipative bosonic mode to occur within 100 ns, compared with the 800 ns of the regular slow varying scheme. Furthermore, by utilizing possible non-adiabatic dynamics, our MMOC protocol -based on methods from optimal control theory -can achieve unconditional adiabatic dynamics for multiple lossy bosonic modes simultaneously in a similar duration. It is worth pointing out that the ringdown pulse cannot be obtained by simply reversing the ringup pulse as can be done in the closed system. Time-reversal symmetry is broken under open-system dynamics.

One issue preventing our protocols from further acceleration is the Kerr nonlinearity correction term 1 2 K(b † ) 2 b 2 to the resonator mode dispersive Hamiltonian [35] . While our methods are designed to work in the linear regime, increasing the protocol speed requires fast-growing driving power, and accounting for nonlinearity becomes essential. However, in the weakly nonlinear regime, i.e., where driving power is well below the first-order dissipative phase transition point [36], we empirically, in both simulation and experiment, find the effect of nonlinearity can be largely mitigated by includ- are used, along with their corresponding CD drivings (blue and orange). The drivings are kept constant at 0 after t f . A slow varying sin 2 driving with t f = 800 ns (purple) is used to illustrate the adiabatic ringup timescale. (b) Output amplitude (in mV) measured by room-temperature FPGA. End times of the CD driving protocols are marked for t f = 30 ns (red) and 100 ns (yellow). The non-equilibrium excitation after t f = 30 ns is caused by excitation of the Purcell filter mode due to the relatively small resonator-filter detuning. The error bar is the standard deviation of the points in the equilibrium states at 990 ns. (c,d) IQ trajectories for the t f = 100 ns sin 2 driving (d) and its CD driving (c). Arrows alongside the simulation indicate the direction of the time evolution. The simulated trajectories are calculated based on the input-output formalism detailed in Supplementary Note 1. To avoid the difficulty of direct tomography, we can infer the bosonic system mean field state from the output signal. The good fit between simulation and experiment data suggests that the CD driving does indeed realize adiabatic following. The mismatch mainly comes from the weak nonlinearity of the resonator. All experimental points are averaged over 3 × 10 4 experiments. A Savitzky-Golay filter with window length 21 and polynomial order 3 is applied to improve the signal-noise ratio.

ing a mean resonator frequency shift ∼ Kn [37] in our protocols, where n is the equilibrium photon number in the resonator mode. If the drive exceeds this dissipative phase transition point, not only does the microwave output suddenly increase, we also observe the transmon qubit becoming excited in a qualitatively similar way to a previously reported result [30] . This result is shown in Fig. 4 , where we stimulate the resonator mode with different amplitudes and durations and measure the remaining ground state population. Our results support the theory [38] that the resonator phase transition and the qubit excitation coincide.

The methods we have presented here are a first step to accelerating open system dynamics in a broader range of settings. The open system CD suppresses diabatic transitions out of the DFS and opens the door for STA to more applications, such as quantum thermodynamics, reservoir engineering [39] , and holonomic quantum computation based on bosonic codes [40] . The openloop MMOC protocol, being independent of qubit-state, is suitable for fast resonator ringup(ringdown) [41, 42] and can be used to reduce measurement cycle times in a b c Supplementary Figure 3 . Multi-mode optimal control (MMOC) by single-port driving. (a) MMOC pulses. A multi-section pulse with target equilibrium time t f = 60 ns is applied to shortcut thermal equilibrium of both resonator and filter modes for different qubit states simultaneously. A different pulse is used to reset to the vacuum state in the end. (b,c) Measured output amplitudes (in mV) for the MMOC pulses with and without reset. The steady output signal is achieved for qubit state |0 (b) and |1 (c) about 30 ns later than the target time, likely due to the high driving amplitude and low-Q energy-storing components in the feedline, such as the impedance-matched Josephson parametric amplifier. This explanation is reinforced by the observation of a similar 30 ns tail with a far-detuned and high-amplitude driving which, in principle, will not excite the multi-mode cQED system (SN 9). In (c), the steady output decays over time due to the T1 decay of the qubit. The error bar is the standard deviation of the points in the equilibrium states at 930 ns. The good final reset performance for the mixed qubit state demonstrates that MMOC works for both qubit states simultaneously. The corresponding IQ trajectories are plotted inset, from which we can infer the system undergoes highly non-equilibrium dynamics during the MMOC pulse. All data are processed in the same way as in Fig. 2. quantum error correction protocols. Both methods can also potentially be applied to other dissipative bosonic systems, such as optical resonant cavities and optomechanics. As recent theoretical developments have shown, open-system STA also has interdisciplinary applications such as accelerating biophysics [43] , bioengineering [22] , and chemical reaction dynamics [44] .

Experimental Setup and Calibration. Our devices consists of a tunable transmon qubit coupled to a λ/4 microwave resonator with coupling g/2π ≈ 80 MHz. The resonator is coupled to an individual λ/4 Purcell filter with coupling J/2π ≈ 10 MHz. The hybrid mode frequencies are ω 0 r /2π = 6.4427 GHz, ω 1 r /2π = 6.4378 GHz and ω 0 f /2π = 6.4634 GHz, ω 1 f /2π = 6.4673 GHz and linewidths κ 0 r = 1/62.88 ns, κ 1 r = 1/77.93 ns and κ 0 f = 1/17.86 ns, κ 1 f = 1/15.64 ns. The frequencies ω r , ω f can a b

Supplementary Figure 4 . Qubit excitation by feedline driving. (a) The pulse sequence used to check the impact of the microwave driving from the feedline on the transmon qubit. We apply a stimulating microwave pulse with fixed frequency ω d and different amplitudes and durations. We wait for 500 ns to clear the resonator-filter population before performing the qubit dispersive measurement. (b) Population of |0 as a function of the amplitude and duration of the stimulating pulse. Once a critical pulse amplitude is exceeded, the qubit is excited. Our results are qualitatively similar to previous findings [30] .

To avoid unwanted excitation, we optimize the driving pulse amplitude and the speed of the protocol. For other linear lossy bosonic mode systems without a dispersively coupled qubit (e.g. optomechanics), STA will not be limited in this way.

be estimated by fitting the transmission S 21 shown in Fig. 1c with the input-output theory of SN 1, or by fitting the measured time-traced IQ plots with numerical trajectories obtained from quantum Langevin dynamics. The hybrid linewidths κ r , κ f are obtained by measuring the decay time constant of the quenched ringdown. All spectrum parameters are characterized to a precision of < 2π × 0.1 MHz. Driving pulses (t) are generated by an arbitrary waveform generator with a sampling rate of 2 GS/s and upconverted to the driving frequency by IQ mixing with the LO microwave signal. The readout signal is first amplified with an impedance-matched Josephson parametric amplifier, then by high electron mobility transistors, and finally by room temperature amplifiers. It is further homodyne demodulated, digitized by an analogue-to-digital converter at 1 GS/s and analyzed by a room-temperature DAQ FPGA. Each experimental point in the figures is an average of 3 × 10 4 experiments. Hybrid-Mode Dynamics. Here we show how the interacting dissipative modes are decoupled so that CD driving is realized for a single lossy hybrid mode. The derivation is based on a simplified version of the circuit model shown in Fig. 1b , where the input gate capacitor is ignored, and the qubit's effect is implicitly accounted for through a state-dependent shift on the filter and resonator modes. The exact input-output relations based on the entire circuit are derived in SN 1, which are equivalent to the results of this simplified model through a parameter renormalization.

In the dispersive regime of cQED [45] , after applying the rotating wave approximation, the bare system Hamiltonian in the driving frame (at frequency ω d ) reads

where a, b are bare filter and resonator modes, the superscripts denote the qubit state, ∆ a(b) ≡ ω a(b) − ω d is the driving detuning of mode a(b), J is coupling strength between the modes, and the qubit-induced dispersive shift

Here we have set = 1. Since our protocol time is much shorter than the qubit lifetime, decoherence is dominated by κ a through the cold bath at an effective temperature of 75 mK. In the Heisenberg picture, the system dynamics according to the input-output

where a i is the time-dependent input field. A linear transformation of a and b decouples these equations, resulting in the hybrid modes a 0,1 , b 0,1 with frequencies

where c f , c r are coefficients from the linear transformation. Eq. 6 is equivalent to the master equation in the Schrodinger picture [47]:

where 190503 (2016) .

Here we derive the exact input-output formula used to simulate the output signal of the system shown in Fig. 5 . Our analysis follows that of Ref.

[48], although we additionally account for the distance l from the input capacitor C in to the filter a, which is necessary for the theory to match experimental observations. For wavenumber k, the phase accumulated after passing through this distance is θ = k × l. We find that this phase has a profound effect on the final output signal. Our goal is to determine how the output field r o responds to the input field c i and its interaction with IMPA Supplementary Figure 5 . Modes in the system. The signal network we analyse here includes the feedline, the filter cavity (mode a), the resonator cavity (mode b) and the qubit. In our feed line, there is an input capacitor to reflect the system leakage to the output port. After the filter, there are three circulators and an impedance modified amplifier. All modes and their directions in the calculation are shown.

the system, including filter mode a, the resonator mode b and the qubit state. To simplify the calculation, instead of directly including the qubit state, we will account for its effect by modifying other system mode frequencies.

To build the mode network, we start from the most left input port and consider the transition and reflection of C in as:

where Γ = Z l −Z0 Z l +Z0 is the reflection coefficient, Z 0 is the impedance of the line and the loaded impedance of C in is Z l = 1 iωCin . As a second step, we consider the effect of the microwave length from C in to the system as:

After this, the microwave reaches the T connection between the filter and the feedline. Its scattering matrix is:

Here we assume the impedance of the line connecting the capacitor is also Z 0 . Part of the wave in the feedline will drive the filter mode a, which satisfies the input-output formula:

where κ a is the leakage rate of the mode a. Assuming no output reflection (r i = 0), the relation between the input field c i , the output field r o and the filter mode a is:

To determine the relation between the input c i and the output r o it suffices to deduce how a depends on c i . Under the rotating wave approximation (RWA), the equations of motion in the drive frequency (ω d ) rotating frame are:

where ∆ a/b = ω a/b − ω d is the detuning of mode a/b frequency ω a/b relative to the drive frequency ω d , and J is the coupling strength between modes a and b. It follows from Eq. 8, 9, 10 and 11 that :

Combining Eq. 13 and 14 gives:ȧ

where the effective detuning, leakage rate, and input field of mode a are ∆ a = ∆ a + Im(Γe i2θ )κ a /4, andκ a = κ a [1 + Re(Γe i2θ )]/2, andã i = √ κ a (1 − Γ)e iθ c i /(2 √κ a ) respectively. Note that Eq. 15 is equivalent, up to redefining various parameters, to Eq. 3 and Eq. 4 of the main text. According to the designed values of l and C in , we estimate θ ∼ 0.05 and Γ ∼ 0.98 − 0.17i. To simulate the dynamics of modes a and b, we use a Lindblad master equation with the Hamiltonian:

and the Lindblad operator √κ a a. The effective driving field (t) follows from Eq. 15 and 16 as

which is averaged in the simulation, given the classical (coherent) input field c i . Substituting this into Eq. 12 gives our final input-output formula:

To account for the weak nonlinearity of the resonator and the uncertainty in the estimated design parameters, in the simulation we multiply the √ κ a a term on the right-hand side of Eq. 18 by a complex coefficient, chosen to fit the simulation to experimental data. Physical interpretation. Because the filter mode a driven by (t) can be solved using the Lindblad master equation, we can determine the output mode r o once we know the driving waveform (t). The physical meaning of Eq. 18 can be interpreted as follows. The factor 2 before (t) means only half of the input mode c is used to drive mode a. The signal that finally reaches the output port is twice the driving. The e i2θ Γ term means half of the leakage of a directly goes to the output port, and the other half will go to the input side and be reflected by C in . Finally, these two branches interfere with each other and contribute a complex factor between the input and the system leakage.

Here we give a simple, short derivation of the counterdiabatic (CD) driving (Eq. 1 in the main text) for a single driven bosonic mode coupled to a cold bath, based on a mean-field approximation. We leave a rigorous derivation to Supplementary Note (SN) 3.

As we will see in SN 3, the bosonic mode under consideration can be well approximated by a coherent state, and thus we can use a mean-field approximation for the Heisenberg picture bosonic field a(t), i.e. α(t) = a(t) . Following SN 1, the dynamics are given by the Langevin equation in the drive frame of frequency ω d :

where ∆ r ≡ ω r − ω d is the cavity-drive detuning, κ is the damping rate due to coupling to the readout line, and (t) is the effective drive field. The instantaneous equilibrium state is obtained by settingα = 0. We denote this instantaneous equilibrium state as:ᾱ

If the drive field is varied slowly enough, the adiabatic theorem guarantees thatᾱ(t) be the solution of Eq. 19. Define the instantaneous diabatic excitation δ(t) = α(t) −ᾱ(t). It follows from Eq. 19, and Eq. 20 and its time derivative, that the dynamics of δ(t) satisfies:

where CD is the (new) CD driving. From the boundary conditions δ(0) = 0 andδ(0) = 0, we obtain the desired CD driving as:

Then, for an arbitrarily drive (t) , the instantaneous equilibrium stateᾱ(t) is always the exact dynamic solution of Eq. 19.

In this section, we give rigorous derivations of CD driving for a single driven-dissipative bosonic mode, based on two approaches: (i) Lindblad dynamics [25] and (ii) an adiabatic shortcut of the decoherence free subspace (DFS) [29] . These results justify the mean-field approximation assumed in SN 2, and give additional insight into the adiabatic dynamics of our system.

In what follows, we set = 1. After rotating wave approximation, the Hamiltonian in the driving frame is:

where the cavity-drive detuning ∆ r = ω r − ω d depends on the qubit state in the dispersive regime, and (t) is the effective drive amplitude. The transmission line is viewed as a channel for both driving and dissipation, so the dynamics for the cavity density matrix ρ is described by the master equation:

where the dissipator is D[a]ρ(t) = aρ(t)a † − 1 2 {ρ(t), a † a}. Here, only photon decay is considered, since at the effective temperature T mxc = 75 mK, the average photon population is N ≈ 0.015 1 at readout frequency ω d ≈ 2π ×6.5 GHz. Lindblad dynamics approach. In the adiabatic approximation for open systems [34], in the limit where the Liouvillian L(t) is slowly varying, the density matrix ρ evolves independently in each generalized eigenspace of L(t). In other words, ρ can be decomposed into a direct sum of components, one for each independently evolving Jordan block of L(t). The adiabaticity can be made exact by adding a CD Hamiltonian H CD which suppress the inertial part of L(t) that causes transitions between different Jordan blocks [25] .

To determine H CD , we first we find a superoperatorÔ(t) that transforms L(t) into Jordan canonical form (JCF). That is, with respect to a certain (not necessarily Hermitian) basis for the density matrix B = {ρ 1 , ρ 2 , . . .}, we havê

where J i (t) are the Jordan blocks of size n i × n i . Second, we transfer to the adiabatic frame defined by ρ (t) = O(t) −1 ρ(t) and show that the non-JCF part of the new Lindblad superoperator L (t), i.e.ρ (t) = L (t)ρ (t), can be exactly cancelled by adding a specific CD driving Hamiltonian H CD (t) to the system. In the first step, we chooseÔ(t) to take the form of a displacement superoperatorD(α(t))ρ(t) = D(α(t))ρ(t)D(α(t)) −1 where D(α) ≡ exp (αa † − h.c.) is the displacement operator [49] . Using the fact that D(α(t))a = a − α(t), it is straightforward to show that

Choosing

precisely the instantaneous equilibrium state of SN 2, Eq. 20eliminates the time-dependent driving term in H J (t). Thus, in the adiabatic frame defined by D(ᾱ(t)), the Liouvillian L J (t) = L J is time-independent, so a fixed basis B can be chosen in which L J is in JCF.

In the second step, in the adiabatic frame ρ (t) =D(ᾱ(t)) −1 ρ(t) we havė

i.e. the dynamics are exactly in JCF except for an inertial Hamiltonian H i = iḊ(ᾱ) −1 D(ᾱ) which mixes the Jordan blocks of L J . Adding an additional CD term

consistent with SN 2, Eq. 22. Decoherence free subspace approach. The time-dependent decoherence free subspace (DFS) is a subspace of the full system Hilbert space, in which the open system dynamics is unitary and quasi-steady, i.e. its instaneous motion is generated by an effective Hamiltonian H ef f defined within the DFS. By identifying the time-dependent DFS of our system, we derive the CD driving (Eq. 31) and compare it to the Lindblad dynamics approach.

Following the definition in [29], for system dynamics described by a Lindblad master equationρ = L(t)ρ = −i[H(t), ρ] + k D[a k (t)]ρ where the Lindblad operators a k (t) have possible time dependence, the time-dependent DFS is the space spanned by a set of orthonormal states {|φ j (t) }, satisfying: (i) the basis states |φ j (t) are degenerate eigenstates of any Lindblad operator, i.e. a k (t)|φ j (t) = c k (t)|φ j (t) ∀j, k; (ii) the DFS is closed under the effective

. H ef f acting on a state in the DFS results in a state in the DFS. For a single lossy mode described by Eq. 24, the DFS exists and is spanned by the single state |ᾱ(t) for α(t) = i (t)/(∆ r − iκ/2), and H ef f takes the form of a displaced oscillator H ef f (t) = ∆ r (a † −ᾱ * (t))(a −ᾱ(t)).

Suppose evolution of the DFS is given by the unitary transformation U (t), i.e. U (t)|φ j (0) = |φ j (t) . By direct analogy with closed system CD driving, we can transform to the adiabatic frame defined by U (t) and cancel diabatic excitations out of the DFS by adding a CD Hamiltonian H CD (t) = iU (t)U (t) † . In our example, the natural choice for U (t) is the displacement D(ᾱ(t)), from which we can derive the CD Hamiltonian H CD = iḊ(ᾱ(t))D(ᾱ(t)) −1 , equivalent to Eq. 31 obtained from the Lindblad dynamics approach.

Comments. Steady states and adiabatic timescale. In the adiabatic frame, the Liouvillian L J in Eq. 28 takes the form L J ρ = −i[∆ r a † a, ρ] + κD[a]ρ. Its eigenvalues can be found by observing

where |n are the Fock states, indicating L J is upper triangular in the subspace spanned by {|0 j|, |1 j + 1|, |2 j + 2|, . . .} (or their Hermitian conjugates) for natural numbers j. Hence, L J has non-degenerate eigenvalues e j,k = i∆ r j − κ(j/2 + k), k = 0, 1, 2, . . . (or their complex conjugates) in each subspace and can be exactly diagonalized. We note that the only steady state of L J , i.e. the eigenstate of L J with zero eigenvalue, is given by j = k = 0, which is the vacuum state |0 in the adiabatic frame or the coherent state |α(t) in the lab frame. We also comment on the timescale required for the adiabatic approximation to hold in open systems, following the results of [34] . Analogously to closed quantum systems, a sufficient condition for adiabatic evolution of an open system is

where t f is the total evolution time, u, v ≡ tr(u † v) defines the inner product, ρ j,k (t) are (lab-frame) eigenstates of L(t) with eigenvalues e j,k , andρ j,k (t) are eigenstates of L(t) † . Here the adjoint L(t) † is defined as the superoperator that satisfy

The LHS is hard to evaluate in practice, and a crude estimate is obtained by setting ρ j ,k (s), dρ j,k (s)/ds ∼ 1 for normalized time s ≡ t/t f . The adiabatic condition for the total time t f is then derived as

For ∆ r comparable to κ, which is a usual experimental scenario, the adiabatic condition is t f κ −1 , making STA useful for fast protocols operating within unit lifetimes. This adiabatic condition is verified in Fig. 14, where sin 2 -shaped pulses are applied for different durations t f and sin 2 -shaped output signals are observed only for t f > 10κ −1 ≈ 600 ns.

DFS from Lindblad dynamics. The derivation of CD driving from both Lindblad dynamics and the DFS approach relies on switching to the adiabatic frame defined by D(ᾱ(t)), i.e. |φ = D(ᾱ(t)) † |φ . We note that the DFS of our system (i.e. the coherent state |ᾱ(t) ) is the vacuum state |0 in the adiabatic frame -the only steady eigenstate (i.e. having an eigenvalue with a non-negative real part) of L J . As a result, the steady eigenspace of the Liouvillian L(t) is equivalent to the DFS, whereas this is not true in general since purity of these steady states requires a zero-temperature approximation or negligible thermal photon number N (ω) 1 in the frequency band of interest. For bosonic modes in the high temperature regime (N (ω) ∼ 1 or N (ω) 1), CD driving is still possible by the Lindblad dynamics approach even though the pure-state DFS does not exist. In this case, although CD driving does not prevent heating into the steady thermal state in the adiabatic frame, it ensures fast transport of this steady state, which is still of practical interest.

Mean Field Approximation. Here we show that the single driven-dissipative mode remains in a coherent state, which justifies the mean-field approximation used in SN 2. For open quantum systems, the coherent state is known to be the consequence of the zero-temperature approximation of the environment [47] . Specifically, in the frame defined by a general displacement D(α(t)), the dynamics in Eq. 29 can be rewritten aṡ

Choosing α(t) that satisfies the Langevin dynamics (Eq. 19) thus eliminates the driving term. Consequently, the system stays in the vacuum state in the displaced frame, corresponding to the coherent state |α(t) in the lab frame.

In this section, we discuss the Quantum Speed Limit (QSL) for a driven-dissipative bosonic mode, and show that our CD driving protocol reaches optimal quantum efficiency among all possible experimental controls.

For open quantum systems, the QSL can be formulated as a geometric constraint, i.e. the total length of the system's trajectory is bounded below by the geodesic connecting its initial and final states, where the geometry is defined in terms of the Bures metric [50] for density matrices. This metric is interpreted as the statistical distinguishability between neighbouring quantum states, expressed in terms of the quantum generalization of Fisher information, i.e. the Fisher information maximized over all choices of quantum measurements [51]. For our system, the dynamics can be equivalently described by a unitary operator, generated by an effective Hamiltonian H ef f (t) = iU (t)U (t) † . This reduces the QSL to the Mandelstam-Tamm (MT) bound [52]:

where φ i(f ) is the initial(final) state and ∆H ef f = H ef f − H ef f . Geometrically, the LHS of Eq. 36 is the Bures length s Bures of the geodesic joining the initial and final state and the RHS is the integrated total length of the system trajectory whose velocity is given by ds Bures /dt = F Q (t)/4 = ∆H 2 ef f (t) [33] . Here, F Q (t) is the quantum Fisher information. We define the quantum efficiency of our protocol to be:

As shown in Eq. 35, for a general driving (t) and the system intialized in the ground state, the dynamics is described by the displacement operator, i.e. U (t) = D(α(t)) where α(t) is the solution to the Langevin equatioṅ α(t) = (−i∆ r − κ/2)α(t) − (t) (SN 3, Eq. 35). We note that U (t) can generate arbitrary dynamics in the space orthogonal to the system state |φ(t) = |α(t) , but the extra freedom can be shown to have no contribution to the uncertainty ∆H 2 ef f (t) . With this choice of U (t) it is straightforward to show that H ef f (t) = iα(t)a † + h.c. and ∆H 2 ef f (t) = |α(t)|. For CD driving, ∆H 2 ef f (t) is simply the added drive | CD (t) − (t)|, which provides the resource for adiabatic speedup in view of the energy-time uncertainty principle. Identifying |φ i,f = |α i,f and applying the triangle inequality, we obtain

with equality achieved by straight-line trajectories -made possible by CD driving -in phase space (see Fig. 6 a) . The

Supplementary Figure 6 . Quantum Efficiency of the CD Driving. (a) Mean field trajectory α(t) for a constant drive preceded by a t f = 100 ns sin 2 ringup, with direct or CD driving. Markers are plotted every 50 ns. CD driving maximizes quantum efficiency by finding the shortest path towards the final state. (b) Quantum efficiency η for resonator-drive detunings ∆r/2π = 3 MHz, 6 MHz, plotted for direct (blue for both detunings) and CD driving (orange for ∆r/2π = 3 MHz, green for ∆r/2π = 6 MHz). Direct driving with higher detuning is less efficient as it induces spiral trajectories with greater length. The gray line shows the efficiencies for the example in (a). spiral trajectory α(t) in Fig. 6 a is calculated from Eq. 19 with parameters ∆ r /2π = 3 MHz, κ −1 = 62.88 ns. Fig. 6 b shows η as a function of target displacement |∆α| = |α f − α i | for direct and CD driving with two resonator-drive detunings ∆ r . For both detunings, the CD driving reaches the optimal quantum efficiency experimentally (as given by the right hand side of Eq. 38). In particular, it saturates the MT bound (i.e. η → 1) in the small driving limit |∆α| → 0. The inefficiency at large |∆α| can be explained by the inability to create direct driving to higher-level Fock states, which is a general issue in applying the MT bound for systems with large numbers of energy levels like the bosonic system we consider. Nevertheless, CD driving achieves optimal quantum efficiency within the space of all available pulses, making it favourable for experimental realisation.

In this section, we derive the multi-mode optimal control (MMOC) protocol used in the main text, which takes the hybrid frequencies of multiple oscillators as input, and generates a single-port waveform that puts these lossy bosonic modes into thermal equilibrium at a desired final time t f . We first present a general framework which can be applied to multiple port driving, and then analyse the simpler single port case which is analytically and experimentally more tractable, and sufficient for our needs in the main text.

General multiple-port framework. Consider n linear bosonic modes {a i } n i=1 with frequencies ω i and linear couplings J ij between modes a i and a j . Each mode is coupled to a feedline with strength κ i and driven by an input field c i at frequency ω d . In general, the fields c i can be linearly dependent if they come from the same feedline. In the rotating frame with frequency ω d and after rotating wave approximation (RWA), the system Hamiltonian has the form

where ∆ i = ω i − ω d is the ith detuning. In the Heisenberg picture, following the input-output formalism [46] , the Langevin dynamics for the ith mode is given byȧ

Adopting the mean-field approximation α i ≡ a i for all bosonic modes and defining the effective drive i ≡ √ κ i c i , we can rewrite the Langevin dynamics in matrix form:

where Ω is the (complex) frequency matrix, α = (α 1 , . . . , α n ) T is the column vector of the mean fields, and (t) = ( 1 (t), . . . , n (t)) T is the column vector of the effective drives. Ω can be diagonalized as Ω = O −1 Ω D O, where Ω D = diag(∆ 1 − iκ 1 /2, . . . , ∆ n − iκ n /2) defines the hybrid detunings and linewidths as ∆ i , κ i . We use step function driving in our protocol: the time between initial time t 0 and final time t f = t m is divided into m equal-length intervals, over each of which the drive strength is constant. i.e. for each drive i :

where the i0 , ij , if are constants. Our goal is to put α(t) into the target equilibrium state α f = iΩ −1 f at final time t f , starting from initial equilibrium state α 0 = iΩ −1 0 . The propagator and general solution of differential equation Eq. 40 are

where θ is the step function. Using Ω = O −1 Ω D O and Eq. 43, our goal can be achieved by solving the equations

for ij , where∆ k = ∆ k − iκ k /2 is the complex hybrid detuning. Treating the piece-wise driving ij as a vector l of dimension n × m and defining the n × mn matrix M kl ≡ O ki G kj (l = 1, 2, . . . , mn), Eq. 44 reduces to the linear equations

If complete information of O or Ω are given, the general solution of Eq. 46 can be found by performing a singular value decomposition (SVD) of the matrix M. We concentrate instead on the case of single-port driving, which is considerably simpler. Single port driving. In the special case of single-port driving, all drivings i (t) are linearly dependent, and Eq. 41 reduces to

for constant coefficients c i , a single-port driving vector j , and boundary conditions 0 , f . In this case the i O ki c i terms in Eq. 44 cancel, to give

which takes the form of a linear constraint on . Eq. 48 can be similarly solved via SVD of the n × m matrix G, i.e. G = U · D · V for unitary matrices U, V and diagonal matrix D = (diag(s 1 , . . . , s n ), 0 n×(m−n) ) (m ≥ n, 0 is the zero matrix). This gives the general form of as

where x n+1 , . . . , x m are free complex parameters and V −1 i is the ith column of V −1 , which can be chosen to optimize a user-defined objective function such as the maximum power output of the pulse (see SN 6). We note that in the single-port driving case, the only input to the protocol is the complex detuning∆ i , which is simpler to measure experimentally than the multi-port driving case where full information of Ω is required.

Experimental implementation. Our single-port driving experiment in the main text corresponds to

where a i , b i are the Purcell filter and readout resonator field conditioned on qubit state i = 0, 1, ∆ b,i is the resonator detuning conditioned on qubit state, and is the effective driving on the filter port. The drive constants are c 1 = c 3 = 1, c 2 = c 4 = 0 in Eq. 47, and the solution to Eq. 49 determines the two quadratures of the driving function which, after optimization over the parameters x i (discussed in SN 6), yields the waveform used in Fig. 3 of the main text.

Applications. Two applications of the class of waveforms derived above are fast equilibration of the readout cavity and Purcell filter and the fast reset of them to the vacuum state. In the first case, we set 0 = 0 and f in Eq. 48 to be the constant drive amplitude after t f . In the second case, reverse 0 and f . Unlike the continuous driving pulse in the CD case, the MMOC protocol results in many pulse jumps. In SN 8, we estimate the effect of the distortion induced by the filter in the AWG and confirm we can still use the MMOC pulses safely.

This section covers various numerical aspects of the single-port MMOC protocol of SN 5, including optimization over the maximum power needed, the speed limit of the protocol given limited output power, and the computational complexity of calculating the desired pulse. Energy consumption. The total energy consumption (up to an overall constant) of our pulse in Eq. 49 is, due to unitarity of V −1 ,

where u, v denotes the inner product. From Eq. 51 we see that the minimum energy solution E min is obtained by setting x i = 0.

Minimizing the maximum power output. Given the output power limitations of the microwave devices, it is desirable to minimize the maximum output power P max ({x i }) ≡ max i (| i ({x i })| 2 ) of the pulse. To achieve this, we numerically minimize P max (as a function of free parameters x i from Eq. 49) using a differential evolution algorithm. The resulting optimized MMOC pulses for both the ring-up and reset stage are those used in the main text. Fig. 7 shows the numerical results of P max (in dB) in units of the steady power P 0 after t f , plotted for the ring-up stage with different protocol times t f . For comparison, we also plot a lower bound on P max , which follows from Eq. 51 and the fact that

Given the protocol time t f = 60 ns ≈ 1/κ 0 r used in the main text, we find P max = 14.5 dB after numerical optimization, which is a 4.1 dB reduction from that of the minimum energy pulse. For speedup beyond unit resonator lifetime 1/κ 0 r , P max grows rapidly and may induce unwanted qubit transitions, which sets a speed limit for the MMOC protocol, as discussed in SN 7.

Computational complexity. For single-port MMOC, the number of total qubit-state-conditioned bosonic modes n (in our experiment n = 4) is less than the total number m of pulse sections. In this case, the most time-consuming step in computing Eq. 49 is the singular value decomposition of the m × n matrix G, which has time complexity O(mn min(m, n)) = O(mn 2 ). For the general case of n-port driving, G is replaced by the n × mn matrix M , with corresponding complexity O(mn 3 ). In either case, the problem admits an efficient polynomial time solution. 

We observe that the output signal drifts with a large driving power, which sets a limit on the steady-state driving power of our protocols. This can be explained by the nonlinearity of the resonator [36] . At the same time, according to the previous study [53] , higher transmon levels are excited due to the non-RWA part of the qubit-resonator Hamiltonian, which becomes on-resonant as the photon number in the resonator increases through a Raman-like process. These two observations are shown to be closely related in theoretical simulations [38] . Here, we conduct two different experiments to confirm this point and find limitations of our protocol when applied to the transmon-resonator cQED system.

In the first experiment (Fig. 8) , we compare the output signal of a small pulse of strength 1 a.u., and another larger pulse of strength 2.66 a.u.. Each point is averaged over 3 × 10 4 measurements and moving averaged with a Savitzky-Golay filter (width 21, order 3). IQ traces of the output signal in Fig. 8(c) show a clear drift even long after 5κ −1 a , which can be qualitatively explained by the nonlinearity of the cavity mode. In Fig. 9 and 10, another experiment is conducted to test the impact on the transition out of the |0 state of different pulse amplitudes and durations. A significant drop in P 0 is observed above amplitude 1 a.u.. At this drive amplitude we estimate the steady-state cavity photon number (via qubit spectroscopy) to be roughly the critical photon number n c ≡ (∆/2g) 2 ≈ 18 [45].

Here we show that corrections to the MMOC pulses imposed by the low-pass fourth-order Chebyshev filter are negligible. In our experiments, the MMOC pulse (Eq. 49) from the arbitrary wave generator (AWG) has a carrier driving frequency ω d /(2π) of 200 to 250 MHz, which passes through the filter with a cutoff frequency ω c /(2π) = Figure 9 . Impacts on qubit population induced by resonator excitation. (a) Pulse sequence for the readout fidelity measurement. A stimulation pulse is first applied with variable strength and duration. After a 500 ns resonator ring-down, a weak measurement pulse is applied to measure the qubit state. (b) Readout fidelity for qubit ground state, P0, as a function of stimulation pulse amplitude. The pulse duration is fixed at 2 µ s. A resonance peak (of error) is found between amplitude 1 a.u. to 1.5 a.u., similar to the observations of [53] . (c) P0 versus stimulation pulse duration, plotted for various amplitudes marked by coloured lines in (b). The fidelity drops drastically in the first 500 ns when the amplitudes are greater than 1 a.u. 750 MHz. The piece-wise constant pulse causes the Gibbs phenomenon, a potential source of error.

Here we give a qualitative evaluation of this error. To simplify our calculations, we assume that the passband's transfer function g is 1, and is 0 outside the passband. The waveform after the filter (t) is described by a convolution F of the pre-filter waveform (t) with the filter function

with frequency cutoffs ω 0 = ω d − ω c and ω 1 = ω d + ω c . The calculation is done in the rotating frame defined by ω d . Replacing by in Eq. 48 results in a modified constraint matrix G , given by G = G , which satisfies

where [t 0 , t 1 ] is the finite time window in which the corrected pulse is integrated over, and θ j (t) is the jth square function which equals 1 for t j−1 ≤ t ≤ t j and 0 elsewhere. The single integral Eq. 56 is easy to evaluate numerically and should be compared to the original matrix elements G ij . We find that the relative difference between G ij and G ij is less than 10 −6 and is generally independent of the choice of window time t 0 , t 1 .

Supplementary Figure 10 . Impacts on qubit population induced by cavity excitation, the full data. Population P0 of the |0 state, measured as a function of stimulation amplitude and duration, following the procedures in Fig. 9 . A clear drop in P0 is found for amplitude greater than 1 a.u..

In this section, we show additional experimental data we have collected. In Fig. 11 , we test the widely used square wave driving with an initial amplitude twice as large as the remaining waveform, and see little decrease in the equilibrium time. In Fig. 12 , we apply the CD pulse designed according to the set of parameters for the |0 and |1 states, and confirm that in the single-port driving situation, CD is only able to accelerate one mode at one time. The main text shows the 4-mode MMOC protocol controlling all four modes with the same pulse. In Fig. 13 , we design the MMOC for two modes corresponding to only one specific qubit state, then apply it on both qubit states. In this case, the method only works when the qubit is in the correct state. In Fig. 14 , we compare the output amplitudes of the sin 2 and the corresponding CD ringup drives of different durations. Fig. 15 and Fig. 16 show the output signal's IQ trajectories for the sin 2 and the corresponding CD ringup drives. In Fig. 17 , a large-amplitude and far off-resonant sin 2 drive , and corresponding CD drive are applied. The far-detuning guarantees the cQED system will not be excited. According to the input-output theory, the output signal will be a simple rescaling of the input signal. The designed input ringup duration t f is 30 ns. However, we see the output reaches the designed stable value around t = 65 ns. This tail indicates the filtering effect of some low-Q and energy-storing microwave components in the feedline. This effect is more noticeable when the input signal is larger and can be simulated by a convolution with a low-pass filter transfer function.

(c) (b) (a) Supplementary Figure 11 . Initial larger pulse. (a) Waveform comparison between the normal driving pulse and a driving pulse with an amplitude 2× larger during the initial 30 ns. In (b) and (c), we see that the initial larger driving pulse does not make the equilibrium process much faster.

Supplementary Figure 12 . Apply the CD driving designed for |0 on |1 state. The driving pulse for |0 is the same as the one used in Fig. 2 in the main text. The same pulse applied when the qubit is prepared in the |1 state will not accelerate the evolution to equilibrium.

(b) (a) Supplementary Figure 13 . Applying the multi-section MMOC designed for |0 (|1 ) when the qubit is in the |1 (|0 ) state. In addition to the 4-mode MMOC we used in the main text, we also test the 2-mode MMOC, designed for the hybrid resonator and filter modes with a fixed qubit state. In (a) and (b), we show the performance of the MMOC for a specific qubit state on both states and see that, as expected, the performance is worse when applied with the qubit in the other state. In (b), the decay of the output from ≈ 100 for |1 is caused by the limited lifetime of the qubit. Supplementary Figure 15 . I-Q plots of CD driving for different ring-up durations. We see a huge spike during the ringup for a short CD driving, which does not mean the system undergoes highly non-equilibrium dynamics. According to the input-output theory, the monitored IQ is the coherent superposition of the input driving and the leakage of the system. When the ringup is only 30 ns, the large amplitude will excite the untargeted filter mode, which is indicated by the spiral collapse to the equilibrium after 30 ns.

Supplementary Figure 16. I-Q plots of the sin 2 drive for different ringup durations. To see how long it takes to reach quasi-equilibrium dynamics for the corresponding bare sin 2 ringup, we measure the IQ trajectories for different t f . Until t f = 800 ns, the trajectory will not change in shape

Output response of the −2π ×100 MHz detuned CD and sin 2 large-amplitude drivings. When we detune the drive frequency by −2π × 100 MHz, the lossy cQED system will not be excited, and the output signal closely follows the input signal. However, we see a delay in the output signal before reaching equilibrium after ringup, and an extended ringdown after the drive has concluded

We thank J.N. Zhang, C.Y. Hsieh, Z.Q. Yin, Z.B. Yang, G.H. Huang, and X. Chen for helpful discussions and comments. We thank the electronics team of Tencent Quantum Lab for preparing the room-temperature electronics.