# Shortcuts to Adiabaticity for Open Systems in Circuit Quantum
Electrodynamics

Zelong Yin  Tencent Quantum Laboratory, Tencent, Shenzhen, Guangdong 518057,
China  Chunzhen Li  Tencent Quantum Laboratory, Tencent, Shenzhen, Guangdong
518057, China  Jonathan Allcock  Tencent Quantum Laboratory, Tencent,
Shenzhen, Guangdong 518057, China  Yicong Zheng  Tencent Quantum Laboratory,
Tencent, Shenzhen, Guangdong 518057, China  Xiu Gu  Tencent Quantum
Laboratory, Tencent, Shenzhen, Guangdong 518057, China  Maochun Dai  Tencent
Quantum Laboratory, Tencent, Shenzhen, Guangdong 518057, China  Shengyu Zhang
Tencent Quantum Laboratory, Tencent, Shenzhen, Guangdong 518057, China
Shuoming An<EMAIL_ADDRESS>Tencent Quantum Laboratory, Tencent,
Shenzhen, Guangdong 518057, China

###### Abstract

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\,\mbox{ns}$ to $100\,\mbox{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.

###### pacs:

42.50.Pq, 03.65.Yz, 02.30.Yy

##  I Introduction

Adiabatic processes – which preserve the non-degenerate 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 swift-equilibration
protocol similar to STA was used to accelerate 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 user-defined 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.

##  II Experimental setup

Supplementary Figure 1: Experimental setup and principles of shortcut to
adiabaticity (STA) (a) Schematic of relevant dynamics. An adiabatic protocol
transfers the initial state $\rho_{i}$ to the final state $\rho_{f}^{0,1}$
along the light pink trajectories in the decoherence-free subspace (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 $\rho_{f}^{0,1}$ 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\pi\approx
80\,\mbox{MHz}$. The Purcell filter mode $a$ is coupled to the resonator mode
with strength $J/2\pi=10.5\,\mbox{MHz}$ and to the feedline with strength
$\kappa_{a}/2\pi=11.4\,\mbox{MHz}$. From residue qubit excitation, the
temperature of the environment is measured to be $75\,\mbox{mK}$, justifying
the cold-bath approximation of the feedline. Driving pulses are applied
through the input field $c_{i}$, and the output field $r_{0}$ is amplified by
an impedance-matched parametric amplifier (IMPA) and homodyne detected to
infer the system dynamics. (c) Measured transmission spectrum $S_{21}$ and
hybrid-mode frequencies for both qubit states. The driving frequency
$\omega_{d}=2\pi\times 6.44025\,\mbox{GHz}$ is allocated roughly in the middle
of $\omega_{r}^{0}=2\pi\times 6.4427\,\mbox{GHz}$ and
$\omega_{r}^{1}=2\pi\times 6.4378\,\mbox{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.

Our setup is illustrated in Fig. 1b and consists of two coupled resonator
modes [31], $a$ and $b$, with coupling strength $J/2\pi\approx
10.5\,\mbox{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\,\mbox{mK}$ with strength $\kappa_{a}/2\pi\approx 11.4\,\mbox{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 qubit-dependent hybrid
modes (normal modes) $a^{\prime 0,1}$ and $b^{\prime 0,1}$ (see Methods). A
coherent drive with amplitude $\epsilon(t)$ at frequency
$\omega_{d}/2\pi=6.44025\,\mbox{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.

##  III Open System STA by counterdiabatic driving

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

$\epsilon_{CD}(t)=\epsilon(t)-i\frac{\dot{\epsilon}(t)}{\Delta^{0}_{r}-i\kappa^{0}_{r}/2},$
(1)

where $\Delta_{r}^{0}\equiv\omega_{r}^{0}-\omega_{d}$ and $\kappa_{r}^{0}$ are
the frequency detuning and decay constant respectively, for the hybrid
resonator mode $b^{\prime 0}$ and qubit state in $\left|0\right>$. For
convenience, the superscripts are omitted in what follows. Note that in the
dissipation-free limit $\kappa_{r}\to 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 $\epsilon(t)$ with a
$\sin^{2}$-shape ringup, i.e. $\epsilon(t)=\epsilon_{0}\sin^{2}(\pi t/2t_{f})$
before $t_{f}$ and $\epsilon(t)=\epsilon_{0}$ afterwards. This
$\sin^{2}$-shape reference waveform is chosen to give a smooth, hardware-
friendly CD driving pulse. For both $t_{f}=30\,\mbox{ns}$ and
$100\,\mbox{ns}$, we compare the performance of the reference driving
$\epsilon(t)$ and CD driving $\epsilon_{cd}(t)$. To demonstrate a speedup
compared with the adiabatic timescale [34], estimated to be on the order of
$\kappa_{r}^{-1}$ (see SN 3), we also apply an adiabatic $\sin^{2}$ driving.
Until $t_{f}=800\,\mbox{ns}\gg\kappa_{r}^{-1}$, we observe a relatively good
$\sin^{2}$-shape response, indicating adiabatic evolution. Our results show an
equilibrium time of $100\,\mbox{ns}$ for the CD driving and
$350\,\mbox{ns}\sim 5\kappa_{r}^{-1}$ for a quench driving. For
$t_{f}=30\,\mbox{ns}$, the large and rapidly varying CD pulse induces out-of-
equilibrium excitation of the untargeted mode $a^{\prime}$. As a consequence,
$a^{\prime}$ requires an extra relaxation time of $75\,\mbox{ns}\sim
5\kappa_{f}^{-1}$ to return to its equilibrium. For $t_{f}=100\,\mbox{ns}$,
equilibrium is achieved almost immediately after $t_{f}$, and the non-
equilibrium excitation of $a^{\prime}$ 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.

Supplementary Figure 2: Open system shortcut to adiabaticity by
counterdiabatic (CD) driving and bare drivings. (a) Different pulses used in
the experiment. The $\sin^{2}$ ringups with final time $t_{f}=30\,\mbox{ns}$
(green) and $100\,\mbox{ns}$ (red) are used, along with their corresponding CD
drivings (blue and orange). The drivings are kept constant at $\epsilon_{0}$
after $t_{f}$. A slow varying $\sin^{2}$ driving with $t_{f}=800\,\mbox{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\,\mbox{ns}$ (red) and
$100\,\mbox{ns}$ (yellow). The non-equilibrium excitation after
$t_{f}=30\,\mbox{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\,\mbox{ns}$. (c,d) IQ trajectories for the $t_{f}=100\,\mbox{ns}$
$\text{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\times 10^{4}$ experiments. A Savitzky–Golay filter with window length 21
and polynomial order 3 is applied to improve the signal-noise ratio.

##  IV Multi-mode optimal control

To eliminate the excitation of the untargeted mode $a^{\prime}$ 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. $\epsilon(t_{j-1}<t<t_{j})=\epsilon_{j},\
j\in\\{1,2,3,\dots,m\\}$. During the ringup stage, the control pulse is
subject to the boundary conditions $\epsilon(0)=0$ and
$\epsilon(t_{f})=\epsilon_{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 $\vec{\epsilon}$ and the
equilibrium difference $\vec{y}$ with a transfer matrix $G$ as
$\vec{y}=G\cdot\vec{\epsilon}$. If the complex hybrid detuning
$\tilde{\Delta}\equiv\Delta-i\kappa/2$ is found for each hybrid mode, $G$ and
$\vec{y}$ can be analytically derived. These equilibrium-transfer equations
can be solved via singular value decomposition of $G$ as long as $m\geq 4$.
The resulting $\vec{\epsilon}$ takes the form
$\vec{\epsilon}=\vec{\epsilon_{e}}+\sum_{i=1}^{m-4}x_{i}\vec{\epsilon_{i}}$,
where the essential vector $\vec{\epsilon_{e}}$ can be analytically solved
for, and $\vec{\epsilon_{i}}$ are $m-4$ vectors orthogonal to
$\vec{\epsilon_{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 $\vec{\epsilon}$. As a result, we can obtain a $5\,\mbox{dB}$
reduction in the maximum amplitude compared with the essential vector. See SN
6 for more discussion on the performance 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\,\mbox{ns}$ ringup pulse for fast system equilibration, and then
transfer the system to the vacuum state at the end with another
$60\,\mbox{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 multi-mode ringup and
depletion are almost achieved at the desired time and reduce the duration
required to achieve equilibrium compared with natural relaxation.

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\,\mbox{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 $\left|0\right>$ (b) and
$\left|1\right>$ (c) about $30\,\mbox{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\,\mbox{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 $T_{1}$ decay of the qubit. The error bar is the standard
deviation of the points in the equilibrium states at $930\,\mbox{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.

##  V Discussion

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\,\mbox{ns}$, compared with the
$800\,\mbox{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 $\frac{1}{2}K(b^{\prime\dagger})^{2}b^{\prime 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 including a mean resonator frequency shift $\sim 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.

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 $\omega_{d}$ and different amplitudes and durations. We wait for
$500\,\mbox{ns}$ to clear the resonator-filter population before performing
the qubit dispersive measurement. (b) Population of $\left|0\right>$ 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.

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 open-loop 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
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].

##  VI Methods

Experimental Setup and Calibration. Our devices consists of a tunable transmon
qubit coupled to a $\lambda/4$ microwave resonator with coupling
$g/2\pi\approx 80\,\mbox{MHz}$. The resonator is coupled to an individual
$\lambda/4$ Purcell filter with coupling $J/2\pi\approx 10\,\mbox{MHz}$. The
hybrid mode frequencies are $\omega^{0}_{r}/2\pi=6.4427\,\mbox{GHz},\
\omega^{1}_{r}/2\pi=6.4378\,\mbox{GHz}$ and
$\omega^{0}_{f}/2\pi=6.4634\,\mbox{GHz},\
\omega^{1}_{f}/2\pi=6.4673\,\mbox{GHz}$ and linewidths
$\kappa^{0}_{r}=1/62.88\,\mbox{ns},\ \kappa^{1}_{r}=1/77.93\,\mbox{ns}$ and
$\kappa^{0}_{f}=1/17.86\,\mbox{ns},\ \kappa^{1}_{f}=1/15.64\,\mbox{ns}$. The
frequencies $\omega_{r},\ \omega_{f}$ can 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 $\kappa_{r},\
\kappa_{f}$ are obtained by measuring the decay time constant of the quenched
ringdown. All spectrum parameters are characterized to a precision of
$<2\pi\times 0.1\,\mbox{MHz}$.

Driving pulses $\epsilon(t)$ are generated by an arbitrary waveform generator
with a sampling rate of $2\,\mbox{GS}/s$ and up-converted 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\,\mbox{GS}/s$ and analyzed by a room-temperature DAQ
FPGA. Each experimental point in the figures is an average of $3\times 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
$\omega_{d}$) reads

$H^{0,1}=\Delta_{a}a^{\dagger}a+\Delta_{b}^{0,1}b^{\dagger}b+J(a^{\dagger}b+b^{\dagger}a),$
(2)

where $a,\ b$ are bare filter and resonator modes, the superscripts denote the
qubit state, $\Delta_{a(b)}\equiv\omega_{a(b)}-\omega_{d}$ is the driving
detuning of mode $a$($b$), $J$ is coupling strength between the modes, and the
qubit-induced dispersive shift is $2\chi=\Delta_{b}^{1}-\Delta_{b}^{0}$. Here
we have set $\hbar=1$. Since our protocol time is much shorter than the qubit
lifetime, decoherence is dominated by $\kappa_{a}$ through the cold bath at an
effective temperature of $75\,\mbox{mK}$. In the Heisenberg picture, the
system dynamics according to the input-output theory is [46]

$\displaystyle\dot{a}$
$\displaystyle=-i\Delta_{a}a-\frac{\kappa_{a}}{2}a-iJb-\sqrt{\kappa_{a}}a_{i}$
(3) $\displaystyle\dot{b}$ $\displaystyle=-i\Delta^{0,1}_{b}b-iJa,$ (4)

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^{\prime 0,1},\ b^{\prime 0,1}$ with frequencies $\Delta_{f}^{0,1},\
\Delta_{r}^{0,1}$ and linewidths $\kappa_{f}^{0,1},\ \kappa_{r}^{0,1}$. The
hybrid modes are defined so that $[a^{\prime i},a^{\prime
i\dagger}]=[b^{\prime i},b^{\prime i\dagger}]=1$ for $i=0,1$, whereas
$[a^{\prime i},b^{\prime i}]\neq 0$ due to environmental effects. Omitting the
superscripts for convenience, we have

$\displaystyle\dot{a}^{\prime}$
$\displaystyle=-i\Delta_{f}a^{\prime}-\frac{\kappa_{f}}{2}a^{\prime}-c_{f}\sqrt{\kappa_{a}}a_{i}$
(5) $\displaystyle\dot{b}^{\prime}$
$\displaystyle=-i\Delta_{r}b^{\prime}-\frac{\kappa_{r}}{2}b^{\prime}-c_{r}\sqrt{\kappa_{a}}a_{i},$
(6)

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

$\dot{\rho}_{b^{\prime}}=-i[H_{r}(t),\rho_{b^{\prime}}]+\kappa_{r}\mathcal{D}(b^{\prime})\rho_{b^{\prime}},$
(7)

where
$H_{r}(t)=\Delta_{r}b^{\prime\dagger}b^{\prime}-(ic_{r}\sqrt{\kappa_{a}}a_{i}(t)b^{\prime\dagger}+h.c.)$.
For coherent driving fields, we can approximate the Hamiltonian as
$H_{r}(t)=\Delta_{r}b^{\prime\dagger}b^{\prime}-(i\epsilon_{r}(t)b^{\prime\dagger}+h.c.)$,
where the effective drive is $\epsilon_{r}(t)=c_{r}\sqrt{\kappa_{a}}\langle
a_{i}(t)\rangle$. Based on the master equation 7, we can derive (see SN 3) the
CD driving Eq. (1) for a single hybrid mode.

Data availability.
Raw data for the results reported in the text are available from the authors
upon request.

## References

  * [1] Farhi, E. _et al._ A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem.  _Science_ 292, 472–475 (2001).
  * [2] Aharonov, D. _et al._ Adiabatic quantum computation is equivalent to standard quantum computation.  _SIAM review_ 50, 755–787 (2008).
  * [3] Zanardi, P. & Rasetti, M.  Holonomic quantum computation.  _Physics Letters A_ 264, 94–99 (1999).
  * [4] Ekert, A. _et al._ Geometric quantum computation.  _Journal of Modern Optics_ 47, 2501–2513 (2000).
  * [5] Childs, A. M., Farhi, E. & Preskill, J.  Robustness of adiabatic quantum computation.  _Physical Review A_ 65, 012322 (2001).
  * [6] Chen, X. _et al._ Fast optimal frictionless atom cooling in harmonic traps: Shortcut to adiabaticity.  _Physical review letters_ 104, 063002 (2010).
  * [7] Guéry-Odelin, D. & Muga, J.  Transport in a harmonic trap: Shortcuts to adiabaticity and robust protocols.  _Physical Review A_ 90, 063425 (2014).
  * [8] An, S., Lv, D., Del Campo, A. & Kim, K.  Shortcuts to adiabaticity by counterdiabatic driving for trapped-ion displacement in phase space.  _Nature communications_ 7, 1–5 (2016).
  * [9] Bason, M. G. _et al._ High-fidelity quantum driving.  _Nature Physics_ 8, 147–152 (2012).
  * [10] Zhang, J. _et al._ Experimental implementation of assisted quantum adiabatic passage in a single spin.  _Physical review letters_ 110, 240501 (2013).
  * [11] Du, Y.-X. _et al._ Experimental realization of stimulated raman shortcut-to-adiabatic passage with cold atoms.  _Nature communications_ 7, 1–7 (2016).
  * [12] Zhou, B. B. _et al._ Accelerated quantum control using superadiabatic dynamics in a solid-state lambda system.  _Nature Physics_ 13, 330–334 (2017).
  * [13] Martinis, J. M. & Geller, M. R.  Fast adiabatic qubit gates using only $\sigma$ z control.  _Physical Review A_ 90, 022307 (2014).
  * [14] Theis, L., Motzoi, F., Machnes, S. & Wilhelm, F.  Counteracting systems of diabaticities using drag controls: The status after 10 years (a).  _EPL (Europhysics Letters)_ 123, 60001 (2018).
  * [15] An, S. _et al._ Experimental test of the quantum jarzynski equality with a trapped-ion system.  _Nature Physics_ 11, 193–199 (2015).
  * [16] Guéry-Odelin, D. _et al._ Shortcuts to adiabaticity: Concepts, methods, and applications.  _Reviews of Modern Physics_ 91, 045001 (2019).
  * [17] Demirplak, M. & Rice, S. A.  Adiabatic population transfer with control fields.  _The Journal of Physical Chemistry A_ 107, 9937–9945 (2003).
  * [18] Berry, M. V.  Transitionless quantum driving.  _Journal of Physics A: Mathematical and Theoretical_ 42, 365303 (2009).
  * [19] Song, X.-K., Zhang, H., Ai, Q., Qiu, J. & Deng, F.-G.  Shortcuts to adiabatic holonomic quantum computation in decoherence-free subspace with transitionless quantum driving algorithm.  _New Journal of Physics_ 18, 023001 (2016).
  * [20] Beau, M., Jaramillo, J. & Del Campo, A.  Scaling-up quantum heat engines efficiently via shortcuts to adiabaticity.  _Entropy_ 18, 168 (2016).
  * [21] Martínez, I. A., Petrosyan, A., Guéry-Odelin, D., Trizac, E. & Ciliberto, S.  Engineered swift equilibration of a brownian particle.  _Nature physics_ 12, 843–846 (2016).
  * [22] Iram, S. _et al._ Controlling the speed and trajectory of evolution with counterdiabatic driving.  _Nature Physics_ 17, 135–142 (2021).
  * [23] Ibáñez, S., Martínez-Garaot, S., Chen, X., Torrontegui, E. & Muga, J. G.  Shortcuts to adiabaticity for non-hermitian systems.  _Physical Review A_ 84, 023415 (2011).
  * [24] Chen, Y.-H., Xia, Y., Wu, Q.-C., Huang, B.-H. & Song, J.  Method for constructing shortcuts to adiabaticity by a substitute of counterdiabatic driving terms.  _Physical Review A_ 93, 052109 (2016).
  * [25] Vacanti, G. _et al._ Transitionless quantum driving in open quantum systems.  _New Journal of Physics_ 16, 053017 (2014).
  * [26] Villazon, T., Polkovnikov, A. & Chandran, A.  Swift heat transfer by fast-forward driving in open quantum systems.  _Physical Review A_ 100, 012126 (2019).
  * [27] Sarandy, M. & Lidar, D.  Adiabatic approximation in open quantum systems.  _Physical Review A_ 71, 012331 (2005).
  * [28] Lidar, D. A., Chuang, I. L. & Whaley, K. B.  Decoherence-free subspaces for quantum computation.  _Physical Review Letters_ 81, 2594 (1998).
  * [29] Wu, S., Huang, X., Li, H., Yi, X. _et al._ Adiabatic evolution of decoherence-free subspaces and its shortcuts.  _Physical Review A_ 96, 042104 (2017).
  * [30] Sank, D. _et al._ Measurement-induced state transitions in a superconducting qubit: Beyond the rotating wave approximation.  _Physical review letters_ 117, 190503 (2016).
  * [31] Zhou, Y. _et al._ Rapid and unconditional parametric reset protocol for tunable superconducting qubits.  _arXiv preprint arXiv:2103.11315_ (2021).
  * [32] Koch, J. _et al._ Charge-insensitive qubit design derived from the cooper pair box.  _Physical Review A_ 76, 042319 (2007).
  * [33] Taddei, M. M., Escher, B. M., Davidovich, L. & de Matos Filho, R. L.  Quantum speed limit for physical processes.  _Physical review letters_ 110, 050402 (2013).
  * [34] Sarandy, M. & Lidar, D.  Adiabatic approximation in open quantum systems.  _Physical Review A_ 71, 012331 (2005).
  * [35] Blais, A., Grimsmo, A. L., Girvin, S. M. & Wallraff, A.  Circuit quantum electrodynamics.  _Rev. Mod. Phys._ 93, 025005 (2021).
  * [36] Brookes, P. _et al._ Critical slowing down in circuit quantum electrodynamics.  _Science Advances_ 7, eabe9492 (2021).
  * [37] Leghtas, Z. _et al._ Confining the state of light to a quantum manifold by engineered two-photon loss.  _Science_ 347, 853–857 (2015).
  * [38] Mavrogordatos, T. K. _et al._ Simultaneous bistability of a qubit and resonator in circuit quantum electrodynamics.  _Physical review letters_ 118, 040402 (2017).
  * [39] Albert, V. V., Bradlyn, B., Fraas, M. & Jiang, L.  Geometry and response of lindbladians.  _Physical Review X_ 6, 041031 (2016).
  * [40] Albert, V. V. _et al._ Holonomic quantum control with continuous variable systems.  _Physical review letters_ 116, 140502 (2016).
  * [41] McClure, D. T. _et al._ Rapid driven reset of a qubit readout resonator.  _Physical Review Applied_ 5, 011001 (2016).
  * [42] Bultink, C. C. _et al._ Active resonator reset in the nonlinear dispersive regime of circuit qed.  _Physical Review Applied_ 6, 034008 (2016).
  * [43] Chen, K. _et al._ Dynamics of driven polymer transport through a nanopore.  _Nature Physics_ 1–7 (2021).
  * [44] Koch, C. P.  Controlling open quantum systems: tools, achievements, and limitations.  _Journal of Physics: Condensed Matter_ 28, 213001 (2016).
  * [45] Blais, A., Huang, R.-S., Wallraff, A., Girvin, S. M. & Schoelkopf, R. J.  Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation.  _Physical Review A_ 69, 062320 (2004).
  * [46] Gardiner, C. W. & Collett, M. J.  Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation.  _Physical Review A_ 31, 3761 (1985).
  * [47] Breuer, H.-P., Petruccione, F. _et al._ _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * [48] Heinsoo, J. _et al._ Rapid high-fidelity multiplexed readout of superconducting qubits.  _Physical Review Applied_ 10, 034040 (2018).
  * [49] Walls, D. F. & Milburn, G. J.  _Quantum optics_ (Springer Science & Business Media, 2007).
  * [50] Bures, D.  An extension of kakutani’s theorem on infinite product measures to the tensor product of semifinite w*-algebras.  _Transactions of the American Mathematical Society_ 135, 199–212 (1969).
  * [51] Braunstein, S. L. & Caves, C. M.  Statistical distance and the geometry of quantum states.  _Physical Review Letters_ 72, 3439 (1994).
  * [52] Mandelstam, L. & Tamm, I.  The uncertainty relation between energy and time in non-relativistic quantum mechanics.  In _Selected papers_ , 115–123 (Springer, 1991).
  * [53] Sank, D. _et al._ Measurement-induced state transitions in a superconducting qubit: Beyond the rotating wave approximation.  _Physical review letters_ 117, 190503 (2016).

### Acknowledgments

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.

### Author contributions

Z.L.Y., S.M.A., and C.Z.L. developed the theory. S.M.A. and Z.L.Y. performed
the experiment. All authors contributed to the data analysis and writing of
the manuscript.

### Additional information

The authors declare no competing financial interests. Supplementary
information is available for this paper. Correspondence and requests for
materials should be addressed to S.M.A.(shuomingan@tencent.com).

##  VII Supplementary Note 1. Exact Input-output Theory

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 $\theta=k\times l$. We find that this phase
has a profound effect on the final output signal.

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.

Our goal is to determine how the output field $r_{o}$ responds to the input
field $c_{i}$ and its interaction with 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:

$\begin{split}l_{i}^{\prime}&=(1-\Gamma)c_{i}+\Gamma l_{o}^{\prime}\\\
c_{o}&=\Gamma c_{i}+(1-\Gamma)l_{o}^{\prime},\end{split}$ (8)

where $\Gamma=\frac{Z_{l}-Z_{0}}{Z_{l}+Z_{0}}$ is the reflection coefficient,
$Z_{0}$ is the impedance of the line and the loaded impedance of $C_{in}$ is
$Z_{l}=\frac{1}{i\omega C_{in}}$. As a second step, we consider the effect of
the microwave length from $C_{in}$ to the system as:

$\begin{split}l_{i}&=e^{i\theta}l_{i}^{\prime}\\\
l_{o}&=e^{-i\theta}l_{o}^{\prime}.\end{split}$ (9)

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

$\begin{split}l_{o}&=-\frac{1}{3}l_{i}+\frac{2}{3}r_{i}+\frac{2}{3}a_{o}\\\
r_{o}&=\frac{2}{3}l_{i}-\frac{1}{3}r_{i}+\frac{2}{3}a_{o}\\\
a_{i}&=\frac{2}{3}l_{i}+\frac{2}{3}r_{i}-\frac{1}{3}a_{o}.\end{split}$ (10)

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:

$a_{o}=a_{i}+\sqrt{\kappa_{a}}a,$ (11)

where $\kappa_{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:

$r_{o}=(1-\Gamma)e^{i\theta}c_{i}+\frac{1+e^{2i\theta}\Gamma}{2}\sqrt{\kappa_{a}}a.$
(12)

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
($\omega_{d}$) rotating frame are:

$\begin{split}\dot{a}&=-i\Delta_{a}a-iJb-\frac{\kappa_{a}}{2}a-\sqrt{\kappa_{a}}a_{i}\\\
\dot{b}&=-i\Delta_{b}b-iJa,\end{split}$ (13)

where $\Delta_{a/b}=\omega_{a/b}-\omega_{d}$ is the detuning of mode $a/b$
frequency $\omega_{a/b}$ relative to the drive frequency $\omega_{d}$, and $J$
is the coupling strength between modes $a$ and $b$. It follows from Eq. 8, 9,
10 and 11 that :

$a_{i}=\frac{1-\Gamma}{2}e^{i\theta}c_{i}-\frac{(1-e^{i2\theta}\Gamma)}{4}\sqrt{\kappa_{a}}a.$
(14)

Combining Eq. 13 and 14 gives:

$\begin{split}\dot{a}&=-i\Delta_{a}^{\prime}a-iJb-\frac{\tilde{\kappa}_{a}}{2}a-\sqrt{\tilde{\kappa}_{a}}\tilde{a}_{i}\\\
\dot{b}&=-i\Delta_{b}b-iJa,\end{split}$ (15)

where the effective detuning, leakage rate, and input field of mode $a$ are
$\Delta_{a}^{\prime}=\Delta_{a}+\text{Im}(\Gamma e^{i2\theta})\kappa_{a}/4$,
and $\tilde{\kappa}_{a}=\kappa_{a}[1+\text{Re}(\Gamma e^{i2\theta})]/2$, and
$\tilde{a}_{i}=\sqrt{\kappa_{a}}(1-\Gamma)e^{i\theta}c_{i}/(2\sqrt{\tilde{\kappa}_{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 $\theta\sim 0.05$ and $\Gamma\sim
0.98-0.17i$. To simulate the dynamics of modes $a$ and $b$, we use a Lindblad
master equation with the Hamiltonian:

$H(t)=\Delta_{b}b^{\dagger}b+\Delta_{a}^{\prime}a^{\dagger}a+J(a^{\dagger}b+b^{\dagger}a)+i\epsilon^{*}(t)a-i\epsilon(t)a^{\dagger}$
(16)

and the Lindblad operator $\sqrt{\tilde{\kappa}_{a}}a$. The effective driving
field $\epsilon(t)$ follows from Eq. 15 and 16 as

$\epsilon(t)=\sqrt{\tilde{\kappa}_{a}}\tilde{a}_{i},$ (17)

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:

$r_{o}=\frac{2}{\sqrt{\kappa_{a}}}\epsilon(t)+\frac{1+e^{i2\theta}\Gamma}{2}\sqrt{\kappa_{a}}a.$
(18)

To account for the weak nonlinearity of the resonator and the uncertainty in
the estimated design parameters, in the simulation we multiply the
$\sqrt{\kappa_{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 $\epsilon(t)$
can be solved using the Lindblad master equation, we can determine the output
mode $r_{o}$ once we know the driving waveform $\epsilon(t)$. The physical
meaning of Eq. 18 can be interpreted as follows. The factor 2 before
$\epsilon(t)$ means only half of the input mode $c^{\prime}$ is used to drive
mode $a$. The signal that finally reaches the output port is twice the
driving. The $e^{i2\theta}\Gamma$ 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.

##  VIII Supplementary Note 2. Single-Mode Counterdiabatic Driving: Short
Derivation

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.
$\alpha(t)=\mbox{$\langle a(t)\rangle$}$. Following SN 1, the dynamics are
given by the Langevin equation in the drive frame of frequency $\omega_{d}$:

$\dot{\alpha}=-i\Delta_{r}\alpha-\frac{\kappa}{2}\alpha-\epsilon(t)$ (19)

where $\Delta_{r}\equiv\omega_{r}-\omega_{d}$ is the cavity-drive detuning,
$\kappa$ is the damping rate due to coupling to the readout line, and
$\epsilon(t)$ is the effective drive field. The instantaneous equilibrium
state is obtained by setting $\dot{\alpha}=0$. We denote this instantaneous
equilibrium state as:

$\bar{\alpha}(t)=\frac{i\epsilon(t)}{\Delta_{r}-i\kappa/2},$ (20)

If the drive field is varied slowly enough, the adiabatic theorem guarantees
that $\bar{\alpha}(t)$ be the solution of Eq. 19. Define the instantaneous
diabatic excitation $\delta(t)=\alpha(t)-\bar{\alpha}(t)$. It follows from Eq.
19, and Eq. 20 and its time derivative, that the dynamics of $\delta(t)$
satisfies:

$\dot{\delta}(t)=-i(\Delta_{r}-i\kappa/2)\delta(t)-\epsilon_{CD}(t)+\left(\epsilon(t)-i\frac{\dot{\epsilon}(t)}{\Delta_{r}-i\kappa/2}\right).$
(21)

where $\epsilon_{CD}$ is the (new) CD driving. From the boundary conditions
$\delta(0)=0$ and $\dot{\delta}(0)=0$, we obtain the desired CD driving as:

$\epsilon_{CD}(t)=\epsilon(t)-i\frac{\dot{\epsilon}(t)}{\Delta_{r}-i\kappa/2}.$
(22)

Then, for an arbitrarily drive $\epsilon(t)$ , the instantaneous equilibrium
state $\bar{\alpha}(t)$ is always the exact dynamic solution of Eq. 19.

##  IX Supplementary Note 3. Single-Mode Counterdiabatic Driving: Open Quantum
Dynamics Approach

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 $\hbar=1$. After rotating wave approximation, the
Hamiltonian in the driving frame is:

$H(t)=\Delta_{r}a^{\dagger}a-i\left(\epsilon(t)a^{\dagger}-\epsilon^{*}(t)a\right).$
(23)

where the cavity-drive detuning $\Delta_{r}=\omega_{r}-\omega_{d}$ depends on
the qubit state in the dispersive regime, and $\epsilon(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 $\rho$ is
described by the master equation:

$\begin{split}\dot{\rho}(t)&=\mathcal{L}(t)\rho(t)\\\
&=-i[H(t),\rho(t)]+\kappa\mathcal{D}[a]\rho(t)\end{split}$ (24)

where the dissipator is
$\mathcal{D}[a]\rho(t)=a\rho(t)a^{\dagger}-\frac{1}{2}\\{\rho(t),a^{\dagger}a\\}$.
Here, only photon decay is considered, since at the effective temperature
$T_{mxc}=75\,\mbox{mK}$, the average photon population is $N\approx 0.015\ll
1$ at readout frequency $\omega_{d}\approx 2\pi\times 6.5\,\mbox{GHz}$.

Lindblad dynamics approach. In the adiabatic approximation for open systems
[34], in the limit where the Liouvillian $\mathcal{L}(t)$ is slowly varying,
the density matrix $\rho$ evolves independently in each generalized eigenspace
of $\mathcal{L}(t)$. In other words, $\rho$ can be decomposed into a direct
sum of components, one for each independently evolving Jordan block of
$\mathcal{L}(t)$. The adiabaticity can be made exact by adding a CD
Hamiltonian $H_{CD}$ which suppress the inertial part of $\mathcal{L}(t)$ that
causes transitions between different Jordan blocks [25].

To determine $H_{CD}$, we first we find a superoperator $\hat{\mathcal{O}}(t)$
that transforms $\mathcal{L}(t)$ into Jordan canonical form (JCF). That is,
with respect to a certain (not necessarily Hermitian) basis for the density
matrix $B=\\{\rho_{1},\rho_{2},\ldots\\}$, we have

$\hat{\mathcal{O}}(t)^{-1}\mathcal{L}(t)\hat{\mathcal{O}}(t)=\operatorname{diag}(J_{1}(t),J_{2}(t),\ldots)$
(25)

where $J_{i}(t)$ are the Jordan blocks of size $n_{i}\times n_{i}$. Second, we
transfer to the adiabatic frame defined by
$\rho^{\prime}(t)=\hat{\mathcal{O}}(t)^{-1}\rho(t)$ and show that the non-JCF
part of the new Lindblad superoperator $\mathcal{L^{\prime}}(t)$, i.e.
$\dot{\rho}^{\prime}(t)=\mathcal{L^{\prime}}(t)\rho^{\prime}(t)$, can be
exactly cancelled by adding a specific CD driving Hamiltonian $H_{CD}(t)$ to
the system.

In the first step, we choose $\hat{\mathcal{O}}(t)$ to take the form of a
displacement superoperator
$\hat{\mathcal{D}}(\alpha(t))\rho(t)=D(\alpha(t))\rho(t)D(\alpha(t))^{-1}$
where $D(\alpha)\equiv\exp{(\alpha a^{\dagger}-h.c.)}$ is the displacement
operator [49]. Using the fact that
$\hat{\mathcal{D}}(\alpha(t))a=a-\alpha(t)$, it is straightforward to show
that

$\displaystyle\mathcal{L_{J}}(t)\rho(t)$
$\displaystyle\equiv\hat{\mathcal{D}}(\alpha(t))^{-1}\mathcal{L}(t)\hat{\mathcal{D}}(\alpha(t))\rho(t)$
(26) $\displaystyle=-i[H_{J}(t),\rho(t)]+\kappa\mathcal{D}[a]\rho(t)$ (27)
$\displaystyle H_{J}(t)$
$\displaystyle=\Delta_{r}a^{\dagger}a+[\Delta_{r}\alpha(t)-\frac{i\kappa}{2}\alpha(t)-i\epsilon(t)]a^{\dagger}+h.c.$
(28)

Choosing $\alpha(t)=\bar{\alpha}(t)\equiv i\epsilon(t)/(\Delta_{r}-i\kappa/2)$
– i.e. precisely the instantaneous equilibrium state of SN 2, Eq. 20 –
eliminates the time-dependent driving term in $H_{J}(t)$. Thus, in the
adiabatic frame defined by $D(\bar{\alpha}(t))$, the Liouvillian
$\mathcal{L_{J}}(t)=\mathcal{L_{J}}$ is time-independent, so a fixed basis $B$
can be chosen in which $\mathcal{L_{J}}$ is in JCF.

In the second step, in the adiabatic frame
$\rho^{\prime}(t)=\hat{\mathcal{D}}(\bar{\alpha}(t))^{-1}\rho(t)$ we have

$\begin{split}\dot{\rho}^{\prime}(t)&=-i[i\dot{D}(\bar{\alpha})^{-1}D(\bar{\alpha}),\rho^{\prime}(t)]+\hat{\mathcal{D}}(\bar{\alpha})^{-1}\mathcal{L}(t)\hat{\mathcal{D}}(\bar{\alpha})\rho^{\prime}(t)\\\
&=-i[i\dot{D}(\bar{\alpha})^{-1}D(\bar{\alpha}),\rho^{\prime}(t)]+\mathcal{L_{J}}\rho^{\prime}(t).\end{split}$
(29)

i.e. the dynamics are exactly in JCF except for an inertial Hamiltonian
$H_{i}=i\dot{D}(\bar{\alpha})^{-1}D(\bar{\alpha})$ which mixes the Jordan
blocks of $\mathcal{L_{J}}$. Adding an additional CD term $\mathcal{L}_{CD}\
\rho=-i[H_{CD},\rho]$ to $\mathcal{L}$ exactly cancels $H_{i}$ if
$\hat{\mathcal{D}}(\bar{\alpha})^{-1}H_{CD}=-H_{i}$. That is, if

$\displaystyle H_{CD}(t)$
$\displaystyle=i\dot{D}(\bar{\alpha}(t))D(\bar{\alpha}(t))^{-1}$ (30)
$\displaystyle=-\frac{\dot{\epsilon}(t)}{\Delta_{r}-i\kappa/2}a^{\dagger}+h.c.$
(31)

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_{eff}$ 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
$\dot{\rho}=\mathcal{L}(t)\rho=-i[H(t),\rho]+\sum_{k}\mathcal{D}[a_{k}(t)]\rho$
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
$\\{\mbox{$\left|\phi_{j}(t)\right>$}\\}$, satisfying: (i) the basis states
$\left|\phi_{j}(t)\right>$ are degenerate eigenstates of any Lindblad
operator, i.e.
$a_{k}(t)\mbox{$\left|\phi_{j}(t)\right>$}=c_{k}(t)\mbox{$\left|\phi_{j}(t)\right>$}\
\forall j,k$; (ii) the DFS is closed under the effective Hamiltonian
$H_{eff}(t)\equiv H(t)+i/2\sum_{k}c_{k}^{*}(t)a_{k}(t)+h.c.$, i.e. $H_{eff}$
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
$\left|\bar{\alpha}(t)\right>$ for
$\bar{\alpha}(t)=i\epsilon(t)/(\Delta_{r}-i\kappa/2)$, and $H_{eff}$ takes the
form of a displaced oscillator
$H_{eff}(t)=\Delta_{r}(a^{\dagger}-\bar{\alpha}^{*}(t))(a-\bar{\alpha}(t))$.

Suppose evolution of the DFS is given by the unitary transformation $U(t)$,
i.e.
$U(t)\mbox{$\left|\phi_{j}(0)\right>$}=\mbox{$\left|\phi_{j}(t)\right>$}$. 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)=i\dot{U}(t)U(t)^{\dagger}$. In our
example, the natural choice for $U(t)$ is the displacement
$D(\bar{\alpha}(t))$, from which we can derive the CD Hamiltonian
$H_{CD}=i\dot{D}(\bar{\alpha}(t))D(\bar{\alpha}(t))^{-1}$, equivalent to Eq.
31 obtained from the Lindblad dynamics approach.

Comments.

Steady states and adiabatic timescale. In the adiabatic frame, the Liouvillian
$\mathcal{L_{J}}$ in Eq. 28 takes the form
$\mathcal{L_{J}}\rho=-i[\Delta_{r}a^{\dagger}a,\rho]+\kappa\mathcal{D}[a]\rho$.
Its eigenvalues can be found by observing

$\mathcal{L_{J}}\mbox{$\left|m\right>$}\mbox{$\left<n\right|$}=[-i\Delta_{r}(m-n)-\frac{\kappa}{2}(m+n)]\mbox{$\left|m\right>$}\mbox{$\left<n\right|$}+\kappa\sqrt{mn}\mbox{$\left|m-1\right>$}\mbox{$\left<n-1\right|$}$
(32)

where $\left|n\right>$ are the Fock states, indicating $\mathcal{L_{J}}$ is
upper triangular in the subspace spanned by
$\\{\mbox{$\left|0\right>$}\mbox{$\left<j\right|$},\mbox{$\left|1\right>$}\mbox{$\left<j+1\right|$},\mbox{$\left|2\right>$}\mbox{$\left<j+2\right|$},\ldots\\}$
(or their Hermitian conjugates) for natural numbers $j$. Hence,
$\mathcal{L_{J}}$ has non-degenerate eigenvalues
$e_{j,k}=i\Delta_{r}j-\kappa(j/2+k),\ k=0,1,2,\ldots$ (or their complex
conjugates) in each subspace and can be exactly diagonalized. We note that the
only steady state of $\mathcal{L_{J}}$, i.e. the eigenstate of
$\mathcal{L_{J}}$ with zero eigenvalue, is given by $j=k=0$, which is the
vacuum state $\left|0\right>$ in the adiabatic frame or the coherent state
$\left|\alpha(t)\right>$ 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

$\max_{0<t<t_{f}}|\langle\tilde{\rho}_{j^{\prime},k^{\prime}}(t),\dot{\rho}_{j,k}(t)\rangle|\ll|e_{j,k}-e_{j^{\prime},k^{\prime}}|,\
\forall j,k,j^{\prime},k^{\prime}$ (33)

where $t_{f}$ is the total evolution time, $\langle
u,v\rangle\equiv\operatorname{tr}(u^{\dagger}v)$ defines the inner product,
$\rho_{j,k}(t)$ are (lab-frame) eigenstates of $\mathcal{L}(t)$ with
eigenvalues $e_{j,k}$, and $\tilde{\rho}_{j,k}(t)$ are eigenstates of
$\mathcal{L}(t)^{\dagger}$. Here the adjoint $\mathcal{L}(t)^{\dagger}$ is
defined as the superoperator that satisfy
$\langle\mathcal{L}(t)^{\dagger}u,v\rangle=\langle u,\mathcal{L}(t)v\rangle,\
\forall u,v$. The LHS is hard to evaluate in practice, and a crude estimate is
obtained by setting
$\langle\tilde{\rho}_{j^{\prime},k^{\prime}}(s),d\rho_{j,k}(s)/ds\rangle\sim
1$ for normalized time $s\equiv t/t_{f}$. The adiabatic condition for the
total time $t_{f}$ is then derived as

$t_{f}\gg\frac{1}{\min_{j,k,j^{\prime},k^{\prime}}(|e_{j,k}-e_{j^{\prime},k^{\prime}}|)}=\frac{1}{\min(\sqrt{\Delta_{r}^{2}+\kappa^{2}/4},\kappa)}.$
(34)

For $\Delta_{r}$ comparable to $\kappa$, which is a usual experimental
scenario, the adiabatic condition is $t_{f}\gg\kappa^{-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\kappa^{-1}\approx 600\,\mbox{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(\bar{\alpha}(t))$, i.e.
$\mbox{$\left|\phi^{\prime}\right>$}=D(\bar{\alpha}(t))^{\dagger}\mbox{$\left|\phi\right>$}$.
We note that the DFS of our system (i.e. the coherent state
$\left|\bar{\alpha}(t)\right>$) is the vacuum state $\left|0\right>$ in the
adiabatic frame - the only steady eigenstate (i.e. having an eigenvalue with a
non-negative real part) of $\mathcal{L}_{J}$. As a result, the steady
eigenspace of the Liouvillian $\mathcal{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(\omega)\ll 1$ in the frequency band of interest. For bosonic modes in the
high temperature regime ($N(\omega)\sim 1$ or $N(\omega)\gg 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(\alpha(t))$,
the dynamics in Eq. 29 can be rewritten as

$\dot{\rho}^{\prime}=-i[\Delta_{r}a^{\dagger}a+(-i\dot{\alpha}+\Delta_{r}\alpha-i\frac{\kappa}{2}\alpha-i\epsilon)a^{\dagger}+h.c.\
,\rho^{\prime}]+\kappa\mathcal{D}[a]\rho^{\prime}.$ (35)

Choosing $\alpha(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
$\left|\alpha(t)\right>$ in the lab frame.

##  X Supplementary Note 4. Quantum Speed Limit of the CD Driving Protocol

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_{eff}(t)=i\dot{U}(t)U(t)^{\dagger}$. This reduces the QSL to the
Mandelstam-Tamm (MT) bound [52]:

$\operatorname{arccos}|\langle\phi_{i}|\phi_{f}\rangle|\leq\int_{t_{i}}^{t_{f}}\sqrt{\mbox{$\langle\Delta
H_{eff}^{2}(t)\rangle$}}dt$ (36)

where $\left|\phi_{i(f)}\right>$ is the initial(final) state and $\Delta
H_{eff}=H_{eff}-\mbox{$\langle H_{eff}\rangle$}$. 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=\sqrt{F_{Q}(t)/4}=\sqrt{\mbox{$\langle\Delta
H_{eff}^{2}(t)\rangle$}}$ [33]. Here, $F_{Q}(t)$ is the quantum Fisher
information. We define the quantum efficiency of our protocol to be:

$\eta\equiv\frac{\operatorname{arccos}|\langle\phi_{i}|\phi_{f}\rangle|}{\int_{t_{i}}^{t_{f}}\sqrt{\mbox{$\langle\Delta
H_{eff}^{2}(t)\rangle$}}dt}$ (37)

As shown in Eq. 35, for a general driving $\epsilon(t)$ and the system
intialized in the ground state, the dynamics is described by the displacement
operator, i.e. $U(t)=D(\alpha(t))$ where $\alpha(t)$ is the solution to the
Langevin equation
$\dot{\alpha}(t)=(-i\Delta_{r}-\kappa/2)\alpha(t)-\epsilon(t)$ (SN 3, Eq. 35).
We note that $U(t)$ can generate arbitrary dynamics in the space orthogonal to
the system state
$\mbox{$\left|\phi(t)\right>$}=\mbox{$\left|\alpha(t)\right>$}$, but the extra
freedom can be shown to have no contribution to the uncertainty $\langle\Delta
H_{eff}^{2}(t)\rangle$. With this choice of $U(t)$ it is straightforward to
show that $H_{eff}(t)=i\dot{\alpha}(t)a^{\dagger}+h.c.$ and
$\sqrt{\mbox{$\langle\Delta H_{eff}^{2}(t)\rangle$}}=|\dot{\alpha}(t)|$. For
CD driving, $\sqrt{\mbox{$\langle\Delta H_{eff}^{2}(t)\rangle$}}$ is simply
the added drive $|\epsilon_{CD}(t)-\epsilon(t)|$, which provides the resource
for adiabatic speedup in view of the energy-time uncertainty principle.
Identifying
$\mbox{$\left|\phi_{i,f}\right>$}=\mbox{$\left|\alpha_{i,f}\right>$}$ and
applying the triangle inequality, we obtain

$\eta\leq\frac{\operatorname{arccos}(e^{-|\alpha_{f}-\alpha_{i}|^{2}/2})}{|\alpha_{f}-\alpha_{i}|}.$
(38)

Supplementary Figure 6: Quantum Efficiency of the CD Driving. (a) Mean field
trajectory $\alpha(t)$ for a constant drive preceded by a
$t_{f}=100\,\mbox{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 $\eta$ for
resonator-drive detunings $\Delta_{r}/2\pi=3\,\mbox{MHz},\ 6\,\mbox{MHz}$,
plotted for direct (blue for both detunings) and CD driving (orange for
$\Delta_{r}/2\pi=3\,\mbox{MHz}$, green for $\Delta_{r}/2\pi=6\,\mbox{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).

with equality achieved by straight-line trajectories – made possible by CD
driving – in phase space (see Fig. 6 a). The spiral trajectory $\alpha(t)$ in
Fig. 6 a is calculated from Eq. 19 with parameters
$\Delta_{r}/2\pi=3\,\mbox{MHz}$, $\kappa^{-1}=62.88\,\mbox{ns}$. Fig. 6 b
shows $\eta$ as a function of target displacement
$|\Delta\alpha|=|\alpha_{f}-\alpha_{i}|$ for direct and CD driving with two
resonator-drive detunings $\Delta_{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.
$\eta\rightarrow 1$) in the small driving limit $|\Delta\alpha|\rightarrow 0$.
The inefficiency at large $|\Delta\alpha|$ 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.

##  XI Supplementary Note 5. Derivation of the Multi-Mode Optimal Control
Protocol

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}\\}_{i=1}^{n}$ with frequencies $\omega_{i}$ and linear couplings
$J_{ij}$ between modes $a_{i}$ and $a_{j}$. Each mode is coupled to a feedline
with strength $\kappa_{i}$ and driven by an input field $c_{i}$ at frequency
$\omega_{d}$. In general, the fields $c_{i}$ can be linearly dependent if they
come from the same feedline. In the rotating frame with frequency $\omega_{d}$
and after rotating wave approximation (RWA), the system Hamiltonian has the
form

$H_{S}=\sum_{i=1}^{n}\Delta_{i}a_{i}^{\dagger}a_{i}+\sum_{i=1}^{n}\sum_{j=1,j\neq
i}^{n}(J_{ij}a_{i}^{\dagger}a_{j}+h.c.)$ (39)

where $\Delta_{i}=\omega_{i}-\omega_{d}$ is the $i$th detuning. In the
Heisenberg picture, following the input-output formalism [46], the Langevin
dynamics for the $i$th mode is given by
$\dot{a}_{i}=-i[a_{i},H_{S}]-(\kappa_{i}/2)a_{i}-\sqrt{\kappa_{i}}c_{i}$.
Adopting the mean-field approximation $\alpha_{i}\equiv\mbox{$\langle
a_{i}\rangle$}$ for all bosonic modes and defining the effective drive
$\epsilon_{i}\equiv\sqrt{\kappa_{i}}\mbox{$\langle c_{i}\rangle$}$, we can
rewrite the Langevin dynamics in matrix form:

$\frac{d}{dt}\vec{\alpha}+i\Omega\cdot\vec{\alpha}=-\vec{\epsilon},\
\Omega=\begin{pmatrix}\Delta_{1}-i\kappa_{1}/2&J_{12}&\ldots&J_{1n}\\\
J_{12}^{*}&\Delta_{2}-i\kappa_{2}/2&\ldots&J_{2n}\\\ \ldots&\ldots&&\ldots\\\
J_{1n}^{*}&J_{2n}^{*}&\ldots&\Delta_{n}-i\kappa_{n}/2\\\ \end{pmatrix}$ (40)

where $\Omega$ is the (complex) frequency matrix,
$\vec{\alpha}=(\alpha_{1},\ldots,\alpha_{n})^{T}$ is the column vector of the
mean fields, and
$\vec{\epsilon}(t)=(\epsilon_{1}(t),\ldots,\epsilon_{n}(t))^{T}$ is the column
vector of the effective drives. $\Omega$ can be diagonalized as
$\Omega=O^{-1}\Omega_{D}O$, where
$\Omega_{D}=\operatorname{diag}(\Delta_{1}^{\prime}-i\kappa_{1}^{\prime}/2,\ldots,\Delta_{n}^{\prime}-i\kappa_{n}^{\prime}/2)$
defines the hybrid detunings and linewidths as $\Delta_{i}^{\prime}$,
$\kappa_{i}^{\prime}$. 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 $\epsilon_{i}$:

$\begin{split}\epsilon_{i}(t<t_{0})&=\epsilon_{i0}\\\
\epsilon_{i}(t_{j-1}<t<t_{j})&=\epsilon_{ij},\ \ \ j\in\\{1,2,3,\dots,m\\}\\\
\epsilon_{i}(t>t_{f})&=\epsilon_{if}.\end{split}$ (41)

where the $\epsilon_{i0},\epsilon_{ij},\epsilon_{if}$ are constants. Our goal
is to put $\vec{\alpha}(t)$ into the target equilibrium state
$\vec{\alpha}_{f}=i\Omega^{-1}\vec{\epsilon}_{f}$ at final time $t_{f}$,
starting from initial equilibrium state
$\vec{\alpha}_{0}=i\Omega^{-1}\vec{\epsilon}_{0}$. The propagator and general
solution of differential equation Eq. 40 are

$\displaystyle D(t-t^{\prime})$
$\displaystyle=e^{-i\Omega(t-t^{\prime})}\theta(t-t^{\prime})$ (42)
$\displaystyle\vec{\alpha}(t)$
$\displaystyle=e^{-i\Omega(t-t_{0})}\vec{\alpha}(t_{0})-\int_{t_{0}}^{t}e^{-i\Omega(t-t^{\prime})}\vec{\epsilon}(t^{\prime})dt^{\prime}\
(t>t_{0}),$ (43)

where $\theta$ is the step function. Using $\Omega=O^{-1}\Omega_{D}O$ and Eq.
43, our goal can be achieved by solving the equations

$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}O_{ki}G_{kj}\epsilon_{ij}$
$\displaystyle=\frac{1}{i\tilde{\Delta}^{\prime}_{k}}\sum_{i=1}^{n}O_{ki}\left(\epsilon_{if}-\epsilon_{i0}e^{-i\tilde{\Delta}^{\prime}_{k}(t_{f}-t_{0})})\right)$
(44) $\displaystyle G_{kj}$
$\displaystyle\equiv\int_{t_{j-1}}^{t_{j}}e^{-i\tilde{\Delta}^{\prime}_{k}(t_{f}-t^{\prime})}dt^{\prime},$
(45)

for $\epsilon_{ij}$, where
$\tilde{\Delta}^{\prime}_{k}=\Delta^{\prime}_{k}-i\kappa^{\prime}_{k}/2$ is
the complex hybrid detuning. Treating the piece-wise driving $\epsilon_{ij}$
as a vector $\epsilon_{l}$ of dimension $n\times m$ and defining the $n\times
mn$ matrix $M_{kl}\equiv O_{ki}G_{kj}\ (l=1,2,\ldots,mn)$, Eq. 44 reduces to
the linear equations

$\displaystyle\sum_{l=1}^{mn}M_{kl}\epsilon_{l}=y_{k},\
y_{k}\equiv\frac{1}{i\tilde{\Delta}^{\prime}_{k}}\sum_{i=1}^{n}O_{ki}\left(\epsilon_{if}-\epsilon_{i0}e^{-i\tilde{\Delta}^{\prime}_{k}(t_{f}-t_{0})}\right)$
(46)

If complete information of $O$ or $\Omega$ 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
$\epsilon_{i}(t)$ are linearly dependent, and Eq. 41 reduces to

$\begin{split}\epsilon_{i}(t<t_{0})&=c_{i}\epsilon_{0}\\\
\epsilon_{i}(t_{j-1}<t<t_{j})&=c_{i}\epsilon_{j},\ \ \
j\in\\{1,2,3,\dots,m\\}\\\
\epsilon_{i}(t>t_{f})&=c_{i}\epsilon_{f}.\end{split}$ (47)

for constant coefficients $c_{i}$, a single-port driving vector
$\epsilon_{j}$, and boundary conditions $\epsilon_{0}$, $\epsilon_{f}$. In
this case the $\sum_{i}O_{ki}c_{i}$ terms in Eq. 44 cancel, to give

$\displaystyle G\cdot\vec{\epsilon}=\vec{y},\ \
y_{k}\equiv\frac{1}{i\tilde{\Delta}^{\prime}_{k}}(\epsilon_{f}-\epsilon_{0}e^{-i\tilde{\Delta}^{\prime}_{k}(t_{f}-t_{0})}),$
(48)

which takes the form of a linear constraint on $\vec{\epsilon}$.

Eq. 48 can be similarly solved via SVD of the $n\times m$ matrix G, i.e.
$G=U\cdot D\cdot V$ for unitary matrices $U,\ V$ and diagonal matrix
$D=(\operatorname{diag}(s_{1},\ldots,s_{n}),0_{n\times(m-n)})$ ($m\geq n$, 0
is the zero matrix). This gives the general form of $\vec{\epsilon}$ as

$\vec{\epsilon}=\sum_{i=1}^{n}\frac{(U^{-1}y)_{i}}{s_{i}}V^{-1}_{i}+\sum_{i=n+1}^{m}x_{i}V^{-1}_{i},$
(49)

where $x_{n+1},\ldots,x_{m}$ are free complex parameters and $V^{-1}_{i}$ is
the $i$th 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 $\tilde{\Delta}_{i}^{\prime}$, which is simpler to
measure experimentally than the multi-port driving case where full information
of $\Omega$ is required.

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

$\vec{a}=(a_{0},b_{0},a_{1},b_{1})^{T},\vec{\epsilon}=(\epsilon,0,\epsilon,0)^{T},\Omega=\begin{pmatrix}\Delta_{a}-i\kappa_{a}/2&J&0&0\\\
J&\Delta_{b,0}&0&0\\\ 0&0&\Delta_{a}-i\kappa_{a}/2&J\\\ 0&0&J&\Delta_{b,1}\\\
\end{pmatrix}$ (50)

where $a_{i},b_{i}$ are the Purcell filter and readout resonator field
conditioned on qubit state $i=0,1$, $\Delta_{b,i}$ is the resonator detuning
conditioned on qubit state, and $\epsilon$ 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 $\epsilon_{0}=0$ and
$\epsilon_{f}$ in Eq. 48 to be the constant drive amplitude after $t_{f}$. In
the second case, reverse $\epsilon_{0}$ and $\epsilon_{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.

##  XII Supplementary Note 6. Numerical Optimization and Speed Limit of the
MMOC Protocol

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}$,

$E(\\{x_{i}\\})\equiv\langle\vec{\epsilon}\ (\\{x_{i}\\}),\vec{\epsilon}\
(\\{x_{i}\\})\rangle=\sum_{i=1}^{n}\left|\frac{(U^{-1}y)_{i}}{s_{i}}\right|^{2}+\sum_{i=n+1}^{m}|x_{i}|^{2},$
(51)

where $\langle\vec{u},\vec{v}\rangle$ denotes the inner product. From Eq. 51
we see that the minimum energy solution $E_{min}$ is obtained by setting
$x_{i}=0$.

Supplementary Figure 7: Speed-power relation for the MMOC protocol. Maximum
output power $P_{max}$ (in dB) consumed, in units of the steady output power
$P_{0}$ after $t_{f}$, plotted for different pulses and final times $t_{f}$.
(Blue) The minimum energy pulse, by setting $x_{i}=0$ in Eq. 49. (Orange) The
pulse whose maximum power is numerically minimized, by exploiting the
redundant degrees of freedoms $x_{i}$. (Green) Theoretical lower bound
obtained from Eq. 52. The experimental parameters are the same as in the main
text.

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}\\})\equiv\operatorname{max}_{i}(|\epsilon_{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
$mP_{max}\geq\sum_{i=1}^{m}|\epsilon_{i}|^{2}\geq E_{min}$:

$P_{max}(\\{x_{i}\\})\geq\frac{1}{m}\sum_{i=1}^{n}\left|\frac{(U^{-1}y)_{i}}{s_{i}}\right|^{2}\equiv
P_{max,lb}$ (52)

Given the protocol time $t_{f}=60\,\mbox{ns}\approx 1/\kappa_{r}^{0}$ used in
the main text, we find $P_{max}=14.5\,\mbox{dB}$ after numerical optimization,
which is a $4.1\,\mbox{dB}$ reduction from that of the minimum energy pulse.
For speedup beyond unit resonator lifetime $1/\kappa_{r}^{0}$, $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\times 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\times mn$ matrix
$M$, with corresponding complexity $O(mn^{3})$. In either case, the problem
admits an efficient polynomial time solution.

##  XIII Supplementary Note 7. Influences of the Large Drive

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\times 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\kappa_{a}^{-1}$,
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 $\left|0\right>$ state of different pulse amplitudes and
durations. A significant drop in $P_{0}$ is observed above amplitude
$1\,\mbox{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}\equiv(\Delta/2g)^{2}\approx 18$ [45].

Supplementary Figure 8: Output signal comparison for different pulse
amplitudes. (a) The square-wave drive pulse with large (blue) and small (pink)
amplitudes. (b) The average output signal (in mV) for the two pulses, as a
function of time. For the small amplitude pulse, steady output is reached
after the expected equilibrium time $t_{e}=5\kappa_{r}^{-1}\approx
300\,\mbox{ns}$. For the large amplitude pulse, the output signal continues to
grow well beyond this time. (c,d) The I/Q quadrature trajectories for the (c)
large and (d) small amplitude pulse. Various intermediate times are colour-
coded in the figure. Supplementary 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, $P_{0}$, as a function of stimulation pulse amplitude. The pulse
duration is fixed at $2\ \mu\,\mbox{s}$. A resonance peak (of error) is found
between amplitude $1\,\mbox{a.u.}$ to $1.5\,\mbox{a.u.}$, similar to the
observations of [53]. (c) $P_{0}$ 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\,\mbox{a.u.}$ Supplementary Figure 10: Impacts on qubit population induced
by cavity excitation, the full data. Population $P_{0}$ of the
$\left|0\right>$ state, measured as a function of stimulation amplitude and
duration, following the procedures in Fig. 9. A clear drop in $P_{0}$ is found
for amplitude greater than $1\,\mbox{a.u.}$.

##  XIV Supplementary Note 8. Propagator Corrections from the Low-pass Filter
in the AWG Driving Line

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 $\omega_{d}/(2\pi)$ of $200$ to $250\,\mbox{MHz}$, which passes
through the filter with a cutoff frequency
$\omega_{c}/(2\pi)=750\,\mbox{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 $\epsilon^{\prime}(t)$
is described by a convolution $\mathcal{F}$ of the pre-filter waveform
$\epsilon(t)$ with the filter function

$\epsilon^{\prime}(t)=\mathcal{F}[\epsilon](t)=\int_{-\infty}^{\infty}\epsilon(\tau)g(t-\tau)d\tau=\int_{-\infty}^{\infty}\epsilon(\tau)d\tau\int_{\omega_{0}}^{\omega_{1}}\frac{d\omega}{2\pi}e^{-i\omega(t-\tau)},$
(53)

with frequency cutoffs $\omega_{0}=\omega_{d}-\omega_{c}$ and
$\omega_{1}=\omega_{d}+\omega_{c}$. The calculation is done in the rotating
frame defined by $\omega_{d}$. Replacing $\epsilon$ by $\epsilon^{\prime}$ in
Eq. 48 results in a modified constraint matrix $G^{\prime}$, given by
$G^{\prime}\vec{\epsilon}=G\vec{\epsilon}^{\prime}$, which satisfies

$\displaystyle G^{\prime}_{ij}$
$\displaystyle=\int_{t_{0}}^{t_{1}}\mathcal{F}[\theta_{j}](t)e^{-i\tilde{\Delta}^{\prime}_{i}(t_{f}-t)}dt$
(54)
$\displaystyle=\int_{t_{j-1}}^{t_{j}}e^{-i\tilde{\Delta}^{\prime}_{i}(t_{f}-\tau)}d\tau\int_{\omega_{0}}^{\omega_{1}}\frac{d\omega}{2\pi}\frac{e^{i(\tilde{\Delta}^{\prime}_{i}-\omega)(t_{1}-\tau)}-e^{i(\tilde{\Delta}^{\prime}_{i}-\omega)(t_{0}-\tau)}}{i(\tilde{\Delta}^{\prime}_{i}-\omega)}$
(55)
$\displaystyle=e^{-i\tilde{\Delta}^{\prime}_{i}t_{f}}\int_{\omega_{0}}^{\omega_{1}}\frac{d\omega}{2\pi}\frac{e^{i\omega
t_{j}}-e^{i\omega
t_{j-1}}}{i\omega}\frac{e^{i(\tilde{\Delta}^{\prime}_{i}-\omega)t_{1}}-e^{i(\tilde{\Delta}^{\prime}_{i}-\omega)t_{0}}}{i(\tilde{\Delta}^{\prime}_{i}-\omega)},$
(56)

where $[t_{0},t_{1}]$ is the finite time window in which the corrected pulse
$\epsilon^{\prime}$ is integrated over, and $\theta_{j}(t)$ is the $j$th
square function which equals 1 for $t_{j-1}\leq t\leq 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^{\prime}_{ij}$ and $G_{ij}$ is less than $10^{-6}$ and
is generally independent of the choice of window time $t_{0}$, $t_{1}$.

##  XV Supplementary Note 9. Additional Experimental Data

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 $\left|0\right>$ and $\left|1\right>$ 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\,\mbox{ns}$.
However, we see the output reaches the designed stable value around
$t=65\,\mbox{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.

Supplementary Figure 11: Initial larger pulse. (a) Waveform comparison between
the normal driving pulse and a driving pulse with an amplitude $2\times$
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 $\left|0\right>$ on
$\left|1\right>$ state. The driving pulse for $\left|0\right>$ 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 $\left|1\right>$ state will not accelerate the evolution to
equilibrium.  Supplementary Figure 13: Applying the multi-section MMOC
designed for $\left|0\right>$($\left|1\right>$) when the qubit is in the
$\left|1\right>$($\left|0\right>$) 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 $\approx 100$ for
$\left|1\right>$ is caused by the limited lifetime of the qubit.
Supplementary Figure 14: CD and $\sin^{2}$ ringup drives for different
durations. Output signal strength for CD and $\sin^{2}$ pulses with different
ring-up durations, using the same parameters and experimental setup as in Fig.
2 of the main text.  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\,\mbox{ns}$, the large amplitude will excite the
untargeted filter mode, which is indicated by the spiral collapse to the
equilibrium after $30\,\mbox{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\,\mbox{ns}$, the trajectory will not change in shape, which
indicates the system has reached the quasi-equilibrium dynamics. Compared this
with the $<100\,\mbox{ns}$ in the case of CD driving.  Supplementary Figure
17: Output response of the $-2\pi\times 100\,\mbox{MHz}$ detuned CD and
$\sin^{2}$ large-amplitude drivings. When we detune the drive frequency by
$-2\pi\times 100\,\mbox{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. This is particularly pronounced when
the driving power is large.