[figure]style=plain,subcapbesideposition=top

††thanks<EMAIL_ADDRESS>
# Superconducting coupler with exponentially large on-off ratio

Catherine Leroux<EMAIL_ADDRESS>Institut quantique &
Département de Physique, Université de Sherbrooke, Sherbrooke J1K 2R1 QC,
Canada  Agustin Di Paolo  Institut quantique & Département de Physique,
Université de Sherbrooke, Sherbrooke J1K 2R1 QC, Canada  Alexandre Blais
Institut quantique & Département de Physique, Université de Sherbrooke,
Sherbrooke J1K 2R1 QC, Canada  Canadian Institute for Advanced Research,
Toronto, ON, Canada

###### Abstract

Tunable two-qubit couplers offer an avenue to mitigate errors in multiqubit
superconducting quantum processors. However, most couplers operate in a narrow
frequency band and target specific couplings, such as the spurious $ZZ$
interaction. We introduce a superconducting coupler that alleviates these
limitations by suppressing all two-qubit interactions with an exponentially
large on-off ratio and without the need for fine-tuning. Our approach is based
on a bus mode supplemented by an ancillary nonlinear resonator mode. Driving
the ancillary mode leads to a coupler-state-dependent field displacement in
the resonator which, in turn, results in an exponential suppression of real
and virtual two-qubit interactions with respect to the drive power. A
superconducting circuit implementation supporting the proposed mechanism is
presented.

##  I Introduction

Two-qubit couplers are useful components for quantum information processing as
they enable fast and high-fidelity operations between qubits while reducing
crosstalk during idle times. Several superconducting coupler designs have been
theoretically proposed and experimentally implemented Niskanen  _et al._
(2007); Majer  _et al._ (2007); Frattini  _et al._ (2017); Yan  _et al._
(2018a); Sete  _et al._ (2021a); Finck  _et al._ (2021); Casparis  _et al._
(2019); Christiane and Homann (2013); Chen  _et al._ (2014); Arute  _et al._
(2019); McKay  _et al._ (2016); Roth  _et al._ (2017); kun Xie  _et al._
(2018); Kounalakis  _et al._ (2018); Chu and Yan (2021); Collodo  _et al._
(2020); Sung  _et al._ (2021); Sete  _et al._ (2021b). These devices can, in
principle, offer precise control of two-qubit interactions while helping to
mitigate frequency crowding effects in multiqubit processors such as to
improve gate speed and fidelity. There are, however, limitations to the
performance of current couplers. For instance, while couplers are designed to
activate interactions between qubits on-demand, spurious interactions can
remain active when the coupler is tuned to its ‘off’ state. A common example
is the ubiquitous always-on cross-Kerr, or $ZZ$, coupling Groszkowski  _et
al._ (2011); Yan  _et al._ (2018b); Mundada  _et al._ (2019); Frattini  _et
al._ (2017); Ku  _et al._ (2020); Zhao  _et al._ (2020). A second difficulty
is that the coupler’s on-off ratio is often sensitive to first order in a
control parameter, such as an external magnetic flux, thus requiring fine-
tuning and frequent calibration. Couplers which do not rely on frequency
tuning do not suffer from this, but the lack of tunability comes with its own
set of challenges such as large crosstalk errors during idle times. Finally,
the impact of these effects could be exacerbated in multiqubit devices where
frequency shifts from spectator qubits can counteract fine-tuning.

Figure 1: Two-qubit coupling strength vs control parameter. For standard
coupler designs, the on-off ratio depends linearly with respect to a control
parameter $\lambda$ (orange line). This results in linear sensitivity to noise
in the control parameter. Our coupler design implements a two-qubit coupling
that can be exponentially suppressed with respect to the control parameter
(green line). Therefore, sensitivity to noise in the ‘off’ state of the
coupler is exponentially small, and the need for fine-tuning is relaxed.

Here, we alleviate these issues by introducing a tunable coupler with an
exponentially large on-off ratio and that does not require fine-tuning of the
coupler or qubit parameters. This is realized by adapting some of the ideas of
protected qubits to coupler designs, as illustrated in Fig. 1. Broadly
speaking, the large on-off ratio is achieved by connecting a bus mode that
wires the two qubits to an ancillary driven nonlinear resonator (NLR), in such
a way that the bus transition matrix elements that control two-qubit
interactions vanish exponentially with respect to the amplitude of the drive
on the ancillary system. This key feature renders the coupler, which includes
the bus and the NLR modes, exponentially insensitive to noise and relaxes the
need for fine-tuning.

This paper is organized as follows. In Sect. II we describe the physical
mechanism enabling the exponential suppression of two-qubit interactions and
introduce a model Hamiltonian realizing this mechanism. In Sect. III, we then
report numerical results demonstrating the exponential suppression of qubit-
qubit interactions mediated by the coupler and discuss implications in the
context of large-scale processors as well as some of the limitations of the
proposed design. Finally, we discuss a superconducting circuit implementation
in Sect. V.

Figure 2: Illustration of the proposed superconducting coupler. a) In the ‘on’
state of the coupler the NLR does not participate in the two-qubit
interactions mediated by the bus mode. Local drives on the qubits or bus
activate two-qubit gates. b) In the ‘off’ state, the NLR is subject to a
microwave drive which strongly suppresses two-qubit interactions that are
mediated by the bus. c) Metapotential of the NLR for the bus states
$\left|0\right\rangle_{b}$ (orange) and $\left|1\right\rangle_{b}$ (purple).
Here $\delta/2\pi=-5.0$ MHz, $\chi/2\pi=-20.0$ MHz and $K_{r}=0$. To help
visualization, the metapotential $\mathcal{E}(I,Q)$ is renormalized as
$\mathcal{E}\cdot\left[(\delta+n\chi)\left|\alpha_{0}/4\right|^{2}\right]^{-1}$
for the $n$th bus state and white corresponds to unity. The global minima of
the metapotentials for $n\geq 2$ are close to the global minimum of the $n=1$
metapotential but not shown for simplicity.

##  II Working principle and Hamiltonian model

Figure 2 a-b) schematically illustrates the proposed device consisting of two
qubits, labeled $Q_{1}$ and $Q_{2}$ coupled by a bus mode $B$ also connected
to a NLR mode $R$. In the absence of a drive on the NLR, the system reduces to
a standard circuit QED setup where the bus mediates energy-exchange
interactions between the qubits Blais  _et al._ (2004, 2021); Yan  _et al._
(2018b); Mundada  _et al._ (2019); Frattini  _et al._ (2017); Ku  _et al._
(2020); Zhao  _et al._ (2020). We assume that the qubit-bus interactions (full
lines) can be modeled by a Jaynes-Cummings-type Hamiltonian. Instead, the bus-
NLR interaction (dashed line) is engineered such that, upon driving the second
mode, the resonator field undergoes a bus-state-dependent displacement
characteristic of a longitudinal interaction. As a result, a distinct
resonator coherent states $\left|\alpha_{n}\right\rangle_{r}$ is associated
with each bus eigenstates $\left|n\right\rangle_{b}$, such that the states
$\left|\psi_{n}\right\rangle=\left|n\right\rangle_{b}\left|\alpha_{n}\right\rangle_{r}$
are stabilized. Then, transitions between the $m$th and $n$th low-energy
eigenstates of the bus are suppressed in the coherent-state amplitude

$\displaystyle\left|n\right\rangle\left\langle m\right|_{b}\otimes
1_{r}\xrightarrow[]{}e^{-\left|\alpha_{n}-\alpha_{m}\right|^{2}/2}\
\left|\psi_{n}\right\rangle\left\langle\psi_{m}\right|.$ (1)

Because all two-qubit interactions are mediated by real or virtual transitions
amongst the bus eigenstates, suppressing these transitions robustly switches
off all interactions mediated by the coupler. As discussed below, if the bus
mode is constrained to its ground state, only the virtual transitions of the
form $0\rightarrow n$ need to be suppressed for all $n$. We also note that
this mechanism is reminiscent of the strategy used to protect cat qubits from
spurious bit flips Cochrane  _et al._ (1999); Puri  _et al._ (2017a);
Mirrahimi  _et al._ (2014).

An effective Hamiltonian realizing this decoupling mechanism can be put in the
form

$\displaystyle\hat{H}=\sum_{j=1}^{2}\hat{H}_{j}+\hat{H}_{b}+\hat{H}_{r}+\sum_{j=1}^{2}\hat{H}_{jb}+\hat{H}_{br}+\hat{H}_{bnlr},$
(2)

where

$\displaystyle\hat{H}_{j}/\hbar=\omega_{j}\hat{q}_{j}^{\dagger}\hat{q}_{j}+\frac{K_{j}}{2}\hat{q}_{j}^{\dagger
2}\hat{q}_{j}^{2},$ (3)
$\displaystyle\hat{H}_{b}/\hbar=\omega_{b}\hat{b}^{\dagger}\hat{b}+\frac{K_{b}}{2}\hat{b}^{\dagger
2}\hat{b}^{2},$ (4)

are the qubits ($j=1,2$) and bus Hamiltonians modeled as Kerr-nonlinear
oscillators, and

$\displaystyle\hat{H}_{r}/\hbar=\omega_{r}\hat{r}^{\dagger}\hat{r}-\varepsilon(t)e^{-i\omega_{d}t}\hat{r}^{\dagger}-\varepsilon^{*}(t)e^{i\omega_{d}t}\hat{r},$
(5)

is the NLR Hamiltonian subject to a drive of amplitude $\varepsilon(t)$ and
frequency $\omega_{d}$. In these expressions, $\hat{q}_{1}$, $\hat{q}_{2}$,
$\hat{b}$, $\hat{r}$ are the annihilation operators of $Q_{1}$, $Q_{2}$, $B$,
and $R$ with mode frequencies $\omega_{1}$, $\omega_{2}$, $\omega_{b}$,
$\omega_{r}$, and anharmonicities $K_{1}$, $K_{2}$, $K_{b}$, respectively.
Although our Hamiltonian model is formulated for the case of transmon qubits
Koch  _et al._ (2007), which can be described as Kerr nonlinear oscillators,
our coupling scheme is in principle applicable to other qubit modalities. The
qubits interact with the bus mode through a Jaynes-Cummings-type Hamiltonian
of the form

$\displaystyle\hat{H}_{jb}/\hbar=g_{j}\left(\hat{q}_{j}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{q}_{j}\right),$
(6)

where $g_{j}$ is the coupling strength, while the bus mode interacts with the
NLR via the cross-Kerr coupling Hamiltonian

$\displaystyle\hat{H}_{br}/\hbar=\chi\hat{b}^{\dagger}\hat{b}\hat{r}^{\dagger}\hat{r}$
(7)

which makes the NLR’s frequency conditional on the bus state via the
dispersive shift $\chi$. For the last term of Eq. 2, we assume the form

$\displaystyle\hat{H}_{bnlr}/\hbar=\sum_{n}\frac{K_{r}}{2}\left|n\right\rangle\left\langle
n\right|_{b}\otimes\hat{r}_{n}^{\dagger 2}(t)\hat{r}_{n}^{2}(t),$ (8)

where $n$ runs over all bus states, and
$\hat{r}_{n}(t)=\hat{r}-\alpha_{n}(t)e^{-i\omega_{d}t}$. This interaction
corresponds to a displaced self-Kerr nonlinearity of the NLR and will be shown
to constrain the system dynamics to a low-energy manifold. Moreover, we show
in Sect. IV how to trade this nonlinear interaction for two additional drives
on the NLR.

Momentarily ignoring the effect of $K_{r}$, the drive on the NLR grows a
coherent state of amplitude $\alpha_{n}$ satisfying

$\displaystyle
i\dot{\alpha}_{n}(t)=(\delta+n\chi-i\kappa/2)\alpha_{n}(t)-\varepsilon(t),$
(9)

which, because of the interaction $\hat{H}_{br}$, is conditional on the bus
state $|n\rangle$. Here, $\delta=\omega_{r}-\omega_{d}$ is the frequency
detuning between the NLR and the drive, and $\kappa$ is the single-photon loss
rate of the NLR. Omitting the qubits, the Hamiltonian of Eq. 2 together with
single-photon loss stabilizes join bus-NLR states of the form
$\left|\psi_{n,k}(t)\right\rangle=\left|n\right\rangle_{b}\left|\alpha_{n}(t)e^{-i\omega_{d}t};k\right\rangle_{r}$,
where
$\left|\alpha;k\right\rangle_{r}=e^{\alpha\hat{r}^{\dagger}-\alpha^{*}\hat{r}}\left|k\right\rangle_{r}$
is the $k$th Fock state displaced by an amplitude $\alpha$. This can be more
clearly seen by plotting the metapotential associated to the Hamiltonian
$\hat{H}$ with the qubit modes traced out and for $K_{r}=0$, obtained by
replacing the operator $\hat{r}$ ($\hat{r}^{\dagger}$) with the complex
variable $I+iQ$ ($I-iQ$). As illustrated in Fig. 2c), this metapotential has a
single well corresponding to a stable point of the system and whose position
in the $I$-$Q$ plane is distinct for each bus state
$\left|n\right\rangle_{b}$. Moreover, because the latter states are associated
with coherent states $\left|\alpha_{n}\right\rangle_{r}$ that are disjoint in
phase space, bus transitions are effectively suppressed. If the system is
energetically constrained to first state $\left|\psi_{n,0}(t)\right\rangle$ of
the metapotential wells, we recover Eq. 1 where the matrix elements of the bus
mode are exponentially suppressed with respect to the drive amplitude. The
nonlinear interaction $\hat{H}_{nlbr}$ of amplitude $K_{r}$ plays the role of
a self-Kerr nonlinearity within each well of the NLR metapotential. As a
result, similarly as in the Kerr-cat qubit Puri  _et al._ (2017b), this Kerr
nonlinearity helps constrain the system’s dynamics to the low-energy states of
each of the metapotential wells.

Rapid switching between the ‘on’ and ‘off’ states of the coupler is realized
by taking advantage of the transitionless-quantum-driving (TQD) method to
rapidly displace the NLR coherent state starting from vacuum Berry (2009). In
the numerical simulations that are discussed below, we use the pulse shape

$\displaystyle\varepsilon(t)=\left(\varepsilon_{0}(t)-\frac{i\dot{\varepsilon}_{0}(t)}{\delta-i\kappa/2}\right)\Theta(\tau-t)+\varepsilon_{0}(\tau)\Theta(t-\tau),$
(10)

where $\varepsilon_{0}(t)$ is a smooth drive amplitude,
$\varepsilon_{0}(0)=0$, and $\Theta(x)$ is the Heaviside step function. With
this choice of drive envelope, the steady-state reached at time $\tau$ takes
the form

$\displaystyle\bar{\alpha}_{n}\simeq\frac{\varepsilon_{0}(\tau)}{\delta+n\chi-i\kappa/2}$
(11)

for each bus state $\left|n\right\rangle_{b}$.

To avoid overlapping metapotential wells and strongly suppress the bus
transitions $0\leftrightarrow n$, the system parameters are chosen such that
$\bar{\alpha}_{0}$ is large with respect to any other $\bar{\alpha}_{n}$. This
last requirement ensures that the bus ground state is well separated in energy
from the NLR excitations, maximizing the exponential suppression of the two-
qubit interaction. More precisely, this is achieved for
$\left|\delta/\chi\right|\ll 1$ and $\left|\kappa/\chi\right|\ll 1$. We note
that choosing the drive such as to make $\bar{\alpha}_{n\neq 0}$ large is also
a valid strategy. However, we numerically find that increasing
$\bar{\alpha}_{0}$ performs better for small to moderate values of $K_{r}$.
Importantly, the TQD protocol is reversible and can be used to bring the
coupler back to the ‘on’ state by emptying the NLR in a time much faster than
$1/\kappa$. The details of this analysis are provided in appendix A.

##  III Numerical experiments

###  III.1 Suppression of bus transitions

Figure 3: Renormalized bus a) Rabi frequency $\tilde{\Omega}/2\pi$ and b)
dephasing time $\tilde{T}_{\varphi}$ under a resonant $\Omega/2\pi=1$ MHz
microwave bus drive as a function of the photon number in the NLR,
$\left|\bar{\alpha}_{0}\right|^{2}$. The results are obtained using a 5-ns-
long TQD scheme with $\kappa/2\pi=100$ KHz, $\delta/2\pi=-5$ MHz,
$K_{b}/2\pi=-300$ MHz. The cross-Kerr $\chi/2\pi$ is varied in the range -5.0
to -30.0 MHz. Fits are done by comparing the time evolution to that of an
effective two-level system, including both $T_{1}$ and $T_{2}$. Dashed lines
correspond to Eqs. 15 and 16. c) Example of the time-evolution traces that are
fitted for $\left|\bar{\alpha}_{0}\right|^{2}=11$.

We now turn to numerical simulations of the concepts presented in the previous
section. To illustrate the working principle—the suppression of bus-state
transitions in the presence of a drive on the NLR—we first simplify the setup
by omitting the qubits. In lieu of the qubits, we add a drive term of the form

$\hat{H}_{\mathrm{Rabi}}/\hbar=\Omega\left(e^{-i\tilde{\omega}_{b}t}\hat{b}^{\dagger}+e^{i\tilde{\omega}_{b}t}\hat{b}\right),$
(12)

where $\Omega$ is the drive amplitude and the drive frequency
$\tilde{\omega}_{b}$ is set to ac-Stark shifted $0-1$ transition frequency of
the bus

$\displaystyle\tilde{\omega}_{b}=\omega_{b}+\frac{\delta(\delta+\chi)-(\kappa/2)^{2}}{\chi}\left|\bar{\alpha}\right|^{2},$
(13)

with

$\bar{\alpha}=\bar{\alpha}_{1}-\bar{\alpha}_{0}=-\frac{\chi}{\delta+\chi-i\kappa/2}\bar{\alpha}_{0}$
(14)

being the distance between the metapotential wells associated with the ground
and first excited states of the bus. In the ‘on’ state of the coupler, the
resonant drive on the bus will result in Rabi oscillations between
$\left|0\right\rangle_{b}$ and $\left|1\right\rangle_{b}$. In the ‘off’ state,
bus transitions, and therefore Rabi oscillations, are exponentially suppressed
with the coherent state amplitude $\bar{\alpha}_{0}$.

Indeed, according to Eq. 1, we expect the Rabi frequency in the presence of a
drive on the NLR to take the form

$\displaystyle\tilde{\Omega}\approx\Omega\exp\left(-\left|\bar{\alpha}\right|^{2}/2\right).$
(15)

Following Eq. 14, $\left|\bar{\alpha}\right|$ is bounded by
$\left|\bar{\alpha}_{0}\right|$. Indeed,
$0\leq\left|\bar{\alpha}\right|<\left|\bar{\alpha}_{0}\right|$ where the lower
bound corresponds to $\chi=0$ or $\bar{\alpha}_{0}=0$, and the upper bound is
reached for $\left|\chi\right|\rightarrow\infty$. As a result, increasing
$\left|\chi\right|$ results in a stronger suppression of $\tilde{\Omega}$.

Figure 3 (a) shows the Rabi frequency $\tilde{\Omega}$ obtained from numerical
integration of the coupler master equation based on Eq. 2. The result includes
damping in the NLR but excludes decoherence in the bus, and it is computed for
different equilibrium values of $\left|\bar{\alpha}_{0}\right|^{2}$ and cross-
Kerr interactions $\chi$. The data points are extracted from fits to the bus
population $\langle\hat{b}^{\dagger}\hat{b}\rangle(t)$ with the bus and NLR
initialized to the vacuum state, see panel (c). The numerical result (symbols)
is in excellent agreement with Eq. 15 (dashed lines) and display the expected
exponential suppression of the bus Rabi oscillations. This suppression becomes
more significant for increasing cross-Kerr coupling $\left|\chi\right|/2\pi$
which is shown here ranging from 5 to 20 MHz. We note that these results are
obtained for $K_{r}=0$. In the absence of the qubits ($g_{j}=0$), choosing
small $\left|\Omega/\delta\right|$ guarantees that the dynamics is mainly
generated by states $\left|\psi_{n,0}\right\rangle$ for which the exponential
suppression of the Rabi frequency is maximized. Indeed, the states
$\left|\psi_{0,k}\right\rangle$ are separated in energy by $\delta$ and thus,
to prevent transitions to $k\neq 0$ states during a $1\rightarrow 0$ bus
transition, we ideally require the matrix elements of Rabi drive Hamiltonian
in the state basis $\left|\psi_{n,k}\right\rangle$ to be small relative to
$\delta$, i.e.
$\left|\Omega\left\langle\psi_{0,k}|\hat{b}|\psi_{1,0}\right\rangle/k\delta\right|=\left|\tilde{\Omega}\bar{\alpha}^{k}/k\delta\sqrt{k!}\right|\ll
1$ for $k\neq 0$.

In the presence of single-photon loss in the NLR, the distinct coherent states
associated to the different bus states lead to bus dephasing. This originates
from the ‘which-bus-state’ information that is carried by the lost photons,
something that is akin to measurement induced-dephasing in the dispersive
readout of circuit QED Gambetta  _et al._ (2006). With $T_{\varphi}$ denoting
the bare bus dephasing time, the coherence time in the presence of the NLR
drive takes the form

$\displaystyle\tilde{T}_{\varphi}\approx\left(\frac{1}{T_{\varphi}}+\frac{\kappa}{2}\left|\bar{\alpha}\right|^{2}\right)^{-1}.$
(16)

Fig. 3 (b) shows this dephasing time extracted from the numerical simulations
including $\kappa\neq 0$ [symbols]. Similarly to the previous case, we find
excellent agreement with the analytical expression (dashed lines). To isolate
the effects of NLR dissipation on the system, we have omitted intrinsic
relaxation and dephasing of the bus mode. A key observation is that, while
transitions between the bus states are suppressed exponentially with
$\bar{\alpha}$, dephasing only increases polynomially with this quantity.
Moreover, we demonstrate below that the dephasing induced on the bus mode does
not percolate to the qubits.

###  III.2 Suppression of two-qubit interactions

Having numerically confirmed that suppressing the bus transitions by driving
the NLR is possible, we now reincorporate the qubits to the model and explore
the two-qubit decoupling. In particular, we characterize the hybridization
between the qubit and bus modes as a function of the NLR drive parameters and
demonstrate how spurious two-qubit couplings, such as the $ZZ$ interaction,
are exponentially suppressed.

####  III.2.1 Polaron transformation

Analyzing the underlying physics of the model Hamiltonian is made easier by
applying a rotating frame transformation and polaron-like unitaries that
displace the NLR mode conditionally on the state of the bus. Acting with Eq.
62 on Eq. 2, the transformed Hamiltonian takes the form (see appendix A)

$\hat{H}^{P}=\sum_{j=1}^{2}\hat{H}_{j}^{P}+\hat{H}_{br}^{P}+\hat{H}_{\kappa}^{P}+\hat{H}_{g}^{P},$
(17)

with

$\displaystyle\hat{H}_{j}^{P}/\hbar=\tilde{\Delta}_{j}\hat{q}_{j}^{\dagger}\hat{q}_{j}+\frac{K_{j}}{2}\hat{q}_{j}^{\dagger
2}\hat{q}_{j}^{2},$ (18)
$\displaystyle\hat{H}_{br}^{P}/\hbar=\left(\delta+\chi\hat{b}^{\dagger}\hat{b}\right)\hat{r}^{\dagger}\hat{r}+\frac{\tilde{K}_{b}}{2}\hat{b}^{\dagger
2}\hat{b}^{2}+\frac{K_{r}}{2}\hat{r}^{\dagger 2}\hat{r}^{2},$ (19)
$\displaystyle\hat{H}_{\kappa}^{P}/\hbar=\frac{i\kappa}{2}\sum_{n}\left(\alpha_{n}\hat{r}^{\dagger}-\alpha_{n}^{*}\hat{r}\right)\left|n\right\rangle\left\langle
n\right|_{b},$ (20)
$\displaystyle\hat{H}_{g}^{P}/\hbar=\sum_{j,n}g_{j}\hat{q}_{j}^{\dagger}e^{i\phi_{n}}\hat{D}_{n,r}\sqrt{n+1}\left|n\right\rangle\left\langle
n+1\right|_{b}+\mathrm{h.c.},$ (21)

where the NLR decay rate $\kappa$ appears in the displacement transformation
according to Eq. 9. In $\hat{H}_{g}^{P}$ we have defined the bus-state-
dependent phases

$\displaystyle\phi_{n}=\frac{\alpha_{n+1}^{*}\alpha_{n}-\alpha_{n}^{*}\alpha_{n+1}}{2i},$
(22)

the NLR displacement operators

$\displaystyle\hat{D}_{n,r}=e^{\left(\alpha_{n+1}-\alpha_{n}\right)\hat{r}^{\dagger}-\left(\alpha_{n+1}^{*}-\alpha_{n}^{*}\right)\hat{r}}.$
(23)

and the ac-Stark shifted qubit-bus detunings and bus anharmonicity

$\displaystyle\tilde{\Delta}_{j}=\omega_{j}-\omega_{b}-\delta\left|\alpha_{0}\right|^{2}+(\delta+\chi)\left|\alpha_{1}\right|^{2},$
(24)
$\displaystyle\tilde{K}_{b}=K_{b}-\delta\left|\alpha_{0}\right|^{2}+2(\delta+\chi)\left|\alpha_{1}\right|^{2}-(\delta+2\chi)\left|\alpha_{2}\right|^{2},$
(25)

respectively. In what follows, we neglect the transients to focus on times
where the polaronic states are fully grown with $\alpha_{n}=\bar{\alpha}_{n}$
as defined in Eq. 11.

In the polaron frame Hamiltonian $\hat{H}^{P}$, all modes are described by a
Kerr nonlinear oscillator Hamiltonian. Moreover, the interaction between the
qubits and the bus, $\hat{H}_{g}^{P}$, reflects the fact that transitions in
the bus are accompanied by displacements of the NLR field. Damping is
described by the usual Lindblad master equation in the polaron frame, derived
in Sect. A.2. In particular, the NLR photon-loss Lindblad operator
$\hat{L}_{r}=\sqrt{\kappa}\hat{r}$ in the laboratory frame transforms to
$\hat{L}_{r}^{P}=\sqrt{\kappa}\left(\hat{r}+\sum_{n}\bar{\alpha}_{n}\left|n\right\rangle\left\langle
n\right|_{b}\right)$ in the polaron frame. The action of both
$\hat{H}_{\kappa}^{P}$ and $\hat{L}_{r}^{P}$ ensures the stabilization of the
vacuum state in the NLR for all $\left|n\right\rangle_{b}$. Under the
rotating-wave approximation, the Lindblad dynamics of the system can be
further reduced to the effective Hamiltonian
$\hat{H}^{P}-\hat{H}_{\kappa}^{P}$ together with the two Lindblad operators
$\sqrt{\kappa}\hat{r}$ and
$\sqrt{\kappa}\sum_{n}\bar{\alpha}_{n}\left|n\right\rangle\left\langle
n\right|_{b}$. By assuming the dynamics to be restrained to the ground and
first-excited states of the bus, the latter operator takes the simpler form
(c.f. appendix A)

$\displaystyle\hat{L}_{b}^{P}=\sqrt{\frac{\kappa\left|\bar{\alpha}\right|^{2}}{4}}\left(\left|1\right\rangle\left\langle
1\right|_{b}-\left|0\right\rangle\left\langle 0\right|_{b}\right).$ (26)

This result is in agreement with the expression for the bus dephasing time
derived in Eq. 16.

Figure 4: Suppression of $1-\mathrm{IPR}_{1000}$ with respect to the NLR
photon number $\left|\bar{\alpha}_{0}\right|^{2}$. The data points are
obtained from numerical diagonalization with a fixed drive amplitude in the
NLR. Analogous plots for $1-IPR_{0100}$ can be found in Sect. B.2. Dashed
lines correspond to the same system but with the NLR undriven and the bus
frequency tuned to the ac-Stark shifted frequency found in the driven system.
The gray regions are bounded by the analytical estimates in Eq. 28 for
$K_{r}\rightarrow\infty$ and $K_{r}=0$. Here
$(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz, $(\omega_{2}-\omega_{b})/2\pi=14.0$
MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz, $\chi/2\pi=-20.0$ MHz, and
$g/2\pi=2.0$ MHz.

####  III.2.2 Inverse participation ratio

A useful quantity to further characterize the exponential suppression of the
two-qubit interactions is the mode hybridization between the qubits and the
bus in the ‘off’ state. To quantify this effect we make use of the inverse
participation ratio (IPR) Berke  _et al._ (2020); Evers and Mirlin (2008),
which here takes the form

$\displaystyle\mathrm{IPR}_{\mu}:=\frac{\sum_{\nu}\left|\left\langle\psi_{h,\mu}|\psi_{b,\nu}\right\rangle\right|^{4}}{\left|\left\langle\psi_{h,\mu}|\psi_{h,\mu}\right\rangle\right|^{2}},$
(27)

where $\left|\psi_{b,\nu}\right\rangle$ and $\left|\psi_{h,\nu}\right\rangle$
the bare ($g_{j}=0$) and hybridized ($g_{j}\neq 0$) eigenstates of the full
system. The IPR is a measure of how localized the wavefunctions of the circuit
modes are with respect to the bare modes, and it ranges from $1/4$ (maximally
delocalized states) to $1$ (maximally localized states) for our system with
four modes.

In the dispersive regime where the qubits are largely detuned from the bus
mode, analytical expressions for the IPR can be obtained by eliminating
$\hat{H}_{g}^{P}$ in $\hat{H}^{P}$ using a Schrieffer-Wolff transformation.
For the bare states $\left|\psi_{b,\nu}\right\rangle$ with
$\nu\in\\{1000,0100,1100\\}$ and where the state indexing corresponds to
$(Q_{1},Q_{2},B,R)$, we obtain in Sect. B.2

$\displaystyle\mathrm{IPR}_{1000}\approx
1-2e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{1}}{\tilde{\Delta}_{1}}\right)^{2}{}_{4}F_{4}\left(\bm{p}_{1};1+\bm{p}_{1};\left|\bar{\alpha}\right|^{2}\right),$
(28) $\displaystyle\mathrm{IPR}_{0100}\approx
1-2e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{2}}{\tilde{\Delta}_{2}}\right)^{2}{}_{4}F_{4}\left(\bm{p}_{2};1+\bm{p}_{2};\left|\bar{\alpha}\right|^{2}\right),$
(29)

and $\mathrm{IPR}_{1100}\approx\mathrm{IPR}_{1000}+\mathrm{IPR}_{0100}-1$,
where $\bar{\alpha}=\bar{\alpha}_{1}-\bar{\alpha}_{0}$, ${}_{p}F_{q}$ is the
generalized hypergeometric function,
$\bm{p}_{j}=\begin{pmatrix}p_{j-}&p_{j-}&p_{j+}&p_{j+}\end{pmatrix}$ with
$p_{j\pm}=\beta[1\pm\sqrt{1+\smash[b]{2\tilde{\Delta}_{j}\beta^{-2}/K_{r}}}]$,
and $\beta=(\delta+\chi)/K_{r}-1/2$. The expressions for the IPR show that the
degree of mode hybridization decreases with increasing detuning. More
importantly, we also see that hybridization is exponentially suppressed with
increasing photon population of the NLR mode. At large photon number, virtual
transitions to higher-energy states within the NLR’s metapotential wells can
impact the level of exponential suppression, something that is represented by
the contribution from the hypergeometric function. These higher energy
transitions can, however, be prevented by increasing
$\left|\delta+\chi\right|$ and $\left|K_{r}\right|$.

To verify these observations, we investigate the quantity
$1-\mathrm{IPR}_{1000}$ as a function of the photon number
$|\bar{\alpha}_{0}|^{2}$ in the NLR by exact diagonalization of the
Hamiltonian of Eq. 17 for $\kappa=0$ (see Fig. 4). Different colors correspond
to different values of the nonlinearity $K_{r}$. Focusing first on panel (a),
obtained for $\delta/2\pi=-1.5$ MHz, the anticipated suppression of the
hybridization with increasing $|\bar{\alpha}_{0}|^{2}$ is clearly observed,
together with the slowdown of that trend for larger
$\left|\bar{\alpha}_{0}\right|\propto\left|\bar{\alpha}\right|$.

The shaded region is plotted using the analytical expression in Eq. 28 for
$K_{r}$ in the range $\left|K_{r}\right|\rightarrow\infty$ to $K_{r}=0$. In
the former limit, the dynamics is constrained to the low-lying polaronic
states, i.e. the NLR is constrained to the displaced Fock states
$\left|0\right\rangle$ and $\left|1\right\rangle$ in the laboratory frame, and
the exponential suppression persists for large
$\left|\bar{\alpha}_{0}\right|$. As a comparison, the dashed line is obtained
from the usual dispersive factor $(g_{i}/\tilde{\Delta}_{i})^{2}$ taking into
account the change in qubits-bus detuning due to the ac-Stark shift and which
corresponds to usual qubit-bus-qubit couplers without the driven NLR. The very
strong suppression of $1-\mathrm{IPR}$ observed in Fig. 4 (a) for our coupler
design has an important consequence: because of the very small hybridization
of the qubit eigenstates, all real and virtual qubit-qubit interactions
mediated by the coupler are exponentially suppressed in amplitude; see Sect.
B.3 for details. Analogous plots for $1-\mathrm{IPR}_{0100}$ can be found in
Sect. B.2.

Figure 4(b) also shows $1-\mathrm{IPR}_{1000}$ as a function of the number of
photons $|\bar{\alpha}_{0}|^{2}$ in the NLR but now for a positive detuning of
$\delta/2\pi=1$ MHz. In this situation, we observe a divergence in the IPR
associated with a resonance in the ac-Stark shifted detunings
$\tilde{\Delta}_{1}$. As discussed in Sect. B.5, this resonance can be
understood from the poles of the generalized hypergeometric function appearing
in Eqs. 28 and 29, which correspond to frequency collisions with higher energy
levels of the NLR. For negative detunings $\delta$, these collisions are
avoided and the suppression is monotonic with photon number. On the other
hand, choosing $\delta>0$ results in a nonmonotonic IPR but this can lead to a
stronger suppression of unwanted interactions. See Sect. B.5 for a detailed
discussion of these frequency collisions and how to take advantage of them.

We conclude this section with a discussion about the optimal sweep rate across
resonances to minimize their impact on the qubits’ dynamics. Frequency
collisions in Figs. 4 and 9 correspond to pairs of states in the uncoupled
system that are brought into resonance during a sweep of $\bar{\alpha}_{0}$.
To understand the impact of these collisions, let us consider the toy-model
Hamiltonian $\hat{H}_{\lambda}(t)=\lambda
e^{i\int_{0}^{t}dt^{\prime}\Delta_{\lambda}(t^{\prime})}\hat{q}_{1}^{\dagger}\hat{b}+\mathrm{h.c.}$
between $Q_{1}$ and the bus in the interaction picture. Here,
$\Delta_{\lambda}(t)$ corresponds to the instantaneous detuning between these
systems including ac-Stark shifts, and $\lambda$ is an effective coupling
strength. In particular, we focus on the simplest case where
$\Delta_{\lambda}=\Delta_{1}-\delta\nu t$, with
$\nu=\left|\bar{\alpha}_{0}(t)\right|^{2}/t$ constant. The approximate
Bogoliubov angle $\theta$ between $Q_{1}$ and the bus at the end of the sweep
can be estimated using a first-order Magnus expansion of the time-evolution
operator, according to $\theta=\left|-i\int_{0}^{t_{\rm
crit}+\tau}dt^{\prime}H_{\lambda}(t^{\prime})\right|$, where $t_{\rm
crit}=\Delta_{1}/\delta\nu$ is the time at which the resonance occurs for
$\delta/\Delta_{1}>0$, and $t_{\rm crit}+\tau$ is the total time of the sweep.
More precisely, we arrive at

$\theta=\left|\frac{\lambda}{\Delta_{1}}\frac{\Delta_{1}}{\sqrt{\delta\nu/2}}\frac{\sqrt{\pi}}{2}\mathrm{Erf}\left(\frac{e^{i\pi/4}\Delta_{\lambda}(t^{\prime})}{\sqrt{\delta\nu/2}}\right)\Big{|}_{t^{\prime}=0}^{t_{\rm
crit}+\tau}\right|,$ (30)

where $\lambda/\Delta_{1}$ is a measure of the dispersive coupling between the
qubit and the bus away from the resonance. Two limiting cases can be studied
for $\Delta_{\lambda}(t_{\rm crit}+\tau)/\Delta_{\lambda}(0)<0$ and therefore
$\delta/\Delta_{1}>0$. For a slow sweep
($\nu\ll\left|\Delta_{1}^{2}/\delta\right|$) and
$\left|\Delta_{\lambda}(t)\right|\gg\left|\lambda\right|$ at $t=0$ and
$t=t_{\rm crit}+\tau$, we derive the asymptotic expression

$\theta\approx\left|\frac{\sqrt{\pi}\lambda}{\sqrt{\delta\nu/2}}\right|,$ (31)

which is minimized for $\nu\gg\left|2\pi\lambda^{2}/\delta\right|$. If
$\delta\sim\lambda$, we recover the intuitive result
$\nu\gg\left|\lambda\right|$, where the sweep rate must exceed the gap set by
the coupling strength $\lambda$. For faster sweeps
($\nu\gg\left|\Delta_{1}^{2}/\delta\right|$), which are the ideal case,
$\theta\approx 0$. We stress the the sweeping rate $\nu$ is limited by the
cross-Kerr interaction strength $\chi$ [c.f. appendix A] that impacts the TQD
protocol. In the toy model $\hat{H}_{\lambda}$, the envelope is
$\varepsilon_{0}(t)=\delta\sqrt{\nu t}$. It follows that
$d^{k}\varepsilon_{0}(t)/dt^{k}=(1/2)^{(k)}\nu^{k}\varepsilon_{0}(t)/\bar{\alpha}_{0}^{2k}(t)$
must be much smaller in magnitude than $|\chi^{k+1}|$. In other words,
$\left|2\chi^{2}\bar{\alpha}_{0}(t_{\rm
crit}+\tau)/\delta\right|\gg\nu\gg\left|\lambda\right|$. This implies that
sweeping rate across the resonances should ideally exceed the coupling
strengths $g_{j}$ and be smaller than $2\chi^{2}/\delta$, where the upper
limit avoids non-adiabatic errors from an imperfect TDQ protocol.

Figure 5: Suppression of the $ZZ$ interaction $\chi_{12}$ between the qubits
as a function of the NLR photon number $\left|\bar{\alpha}_{0}\right|^{2}$.
Data points correspond to numerical diagonalization of the system Hamiltonian
with a fixed drive amplitude in the NLR. Dashed lines correspond to Eq. 33,
i.e. the same system but with the NLR undriven and the bus frequency tuned to
the ac-Stark shifted frequency found in the driven system. Here
$(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz, $(\omega_{2}-\omega_{b})/2\pi=14.0$
MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz, $\chi/2\pi=-20.0$ MHz, and
$g/2\pi=2.0$ MHz.

####  III.2.3 Suppression of spurious interactions

We now analyze how the proposed coupler help to suppress the spurious cross-
Kerr coupling between the qubits, which is given by

$\displaystyle\chi_{12}=\omega_{1100}-\omega_{1000}-\omega_{0100}+\omega_{0000},$
(32)

where
$\omega_{\mu}=\langle\psi_{h,\mu_{\star}}|\hat{H}^{P}-\hat{H}_{\kappa}^{P}|\psi_{h,\mu_{\star}}\rangle$
is the energy associated with the two-qubit eigenstate
$\left|\psi_{h,\mu_{\star}}\right\rangle$ with maximal overlap with the bare
state $\left|\psi_{h,\mu}\right\rangle$, i.e.
$\mu_{\star}=\mathrm{argmax}_{\nu}\left|\left\langle\psi_{h,\nu}|\psi_{b,\mu}\right\rangle\right|^{2}$.

Fig. 5 shows $|\chi_{12}|$ obtained from numerical diagonalization of
$\hat{H}^{P}-\hat{H}_{\kappa}^{P}$ as a function of
$\left|\bar{\alpha}_{0}\right|^{2}$, for different values of $K_{r}$
(symbols). As a comparison, the dashed line shows $|\chi_{12}|$ resulting only
from the change in detuning between the qubits and the bus due to the ac-Stark
shift, and computed using the perturbative expression

$\displaystyle\chi_{12}^{\rm
ac}=\frac{1}{6}\frac{g_{1}^{2}}{\tilde{\Delta}_{1}}\frac{g_{2}^{2}}{\tilde{\Delta}_{2}}\left(\frac{1}{\tilde{\Delta}_{1}}+\frac{1}{\tilde{\Delta}_{2}}\right),$
(33)

valid for $|\tilde{\Delta}_{j}/K_{j}|\ll 1$ and
$|\tilde{\Delta}_{j}/\tilde{K}_{b}|\ll 1$. The latter is obtained from a
Magnus expansion to fourth order in the coupling strengths [c.f. Sect. B.3].
The two resonances observed in the dashed line and the numerical data in panel
(b) correspond to $\tilde{\Delta}_{j}=0$.

As first noticed for $1-\mathrm{IPR}$ in Fig. 4, the suppression of
$\chi_{12}$ is monotonic with photon number for negative detunings $\delta$
[panel (a)], while some nonmonotonic features appear at positive detuning
where the suppression is also stronger [panel (b)]. See Sect. B.5 for a
discussion of the origin of these features.

We also note the presence of more features in Fig. 5b) for $\chi_{12}$ than in
Fig. 4b) for the IPR. The first two dominant peaks in Fig. 5b) result from
accidental resonances between each qubit and the bus, i.e.
$\tilde{\Delta}_{j}=0$, in agreement with the peaks observed for the IPR in
Fig. 4 and Fig. 9. Additional features in Fig. 5 not present in Figs. 4 and 9
result from frequency collisions with higher energy levels in the system,
activated by the ac-Stark shifts in the bus.

Importantly, the suppression of $1-$IPR and the resulting reduction of the
spurious cross-Kerr coupling does not require fine-tuning of the circuit or
drive parameters. Indeed, as illustrated in Figs. 4 and 5, strong suppression
is observed for different choices of circuit parameters including $\delta$ and
$K_{r}$. It is also worth emphasizing that all real and virtual interactions
are suppressed by this scheme. This fact is in stark contrast to other
approaches where cancellation of two-qubit interactions is realized only for a
precise value of a control parameter and where residual virtual interactions
such as $\chi_{12}$ remain present Ku  _et al._ (2020); Zhao  _et al._ (2020);
Sete  _et al._ (2021a).

Finally, we note that it is possible to combine our coupler with other
approaches for suppressing spurious interactions, for instance by using qubits
with opposite sign anharmonicities.

###  III.3 Bus-induced qubit dephasing

At the origin of the suppression of unwanted interaction are the disjoint bus-
state dependent coherent states of the driven NLR. A photon lost from the NLR
carries the ‘which-bus-state’ information and leads to dephasing of the bus
state. Because there exists hybridization between the bus and qubit modes,
this mechanism can introduce additional qubit dephasing. However, as shown in
more details in appendix A, we find that this is not an important contribution
to qubit dephasing. Indeed, by expressing Eq. 26 in the hybridized eigenbasis,
the dephasing rate of the first qubit is given by

$\begin{split}\gamma_{\varphi,1}&=\frac{\kappa\left|\bar{\alpha}\right|^{2}}{2}\left(\sum_{k=0}^{\infty}\left|\left\langle\psi_{h,1000}|\psi_{b,001k}\right\rangle\right|^{2}\right)^{2}\\\
&\approx\frac{\kappa\left|\bar{\alpha}\right|^{2}}{2}\frac{1-\mathrm{IPR}_{1000}}{2},\end{split}$
(34)

where the second line follows from a Schrieffer-Wolff transformation [c.f.
Sect. B.4]. The expression above was obtained with a rotating-wave
approximation, which is valid for $|\gamma_{\varphi,1}/\tilde{\Delta}_{1}|\ll
1$. An expression for the second qubit is obtained by simply replacing the
subscript $1000$ by $0100$.

Similarly to measurement-induced dephasing Gambetta  _et al._ (2006), the
prefactor of Eq. 34 scales with the photon number $|\alpha|^{2}$ in the NLR.
However, because $1-\mathrm{IPR}$ is exponentially suppressed with increasing
$|\bar{\alpha}|^{2}$, the qubit dephasing rate can be made negligible in the
‘off’ state of the coupler. Section B.4 also compares Eq. 34 versus photon
number against the result obtained from numerical diagonalization of
$\hat{H}^{P}-\hat{H}_{\kappa}^{P}$. As with the suppression of unwanted $ZZ$
interactions, the reason for this negligibly small dephasing rate is the very
low hybridization of the qubits’ eigenstates with the bus and NLR modes.

##  IV Effective parametric modulation

In the previous sections, we have seen that large nonlinear interaction
amplitudes $K_{r}$ help in the suppression of the unwanted interactions in the
‘off’ state of the coupler. Here, we explore an alternative strategy that
relies on a two-tone drive on the NLR. Moreover, because the nonlinearity is
not needed in this case, the NLR can be taken to be a linear resonator (LR).
This might also simplify the experimental realization of these ideas.

Our starting point is again the Hamiltonian of Eq. 2 where we now take
$K_{r}=0$ and introduce the following additional drive on the LR

$\displaystyle\hat{H}_{\rm
DD}=-\frac{i\lambda\omega_{m}}{2\bar{\alpha}^{*}}\left[e^{-i\left(\omega_{r}-\omega_{m}\right)t}-e^{-i\left(\omega_{r}+\omega_{m}\right)t}\right]\hat{r}^{\dagger}+\mathrm{h.c.},$
(35)

where $\lambda$ is a real-valued amplitude and the frequency $\omega_{m}$ is
assumed here to be much larger in magnitude than the cross-Kerr interaction
$\chi$. With this additional two-tone drive on the LR, the steady-state bus-
dependent coherent state Eq. 11 becomes

$\displaystyle\bar{\alpha}_{n}\rightarrow\bar{\alpha}_{n}-i\lambda\cos(\omega_{m}t)/\bar{\alpha}^{*},$
(36)

where $\lambda$ plays the role of the amplitude of a modulation around the
steady-state value $\bar{\alpha}_{n}$. Crucially, this modulation changes the
phase $\phi_{n}$ that specifies the bus-state-dependent displacements
Hamiltonian $\hat{H}_{g}^{P}$ in Eq. 21, which can now be written as

$\displaystyle\phi_{n}(t)=\bar{\phi}_{n}-\lambda\mathrm{Re}\left[\frac{\bar{\alpha}_{n+1}-\bar{\alpha}_{n}}{\bar{\alpha}}\right]\cos(\omega_{m}t),$
(37)
$\displaystyle\bar{\phi}_{n}=\frac{\bar{\alpha}_{n+1}^{*}\bar{\alpha}_{n}-\bar{\alpha}_{n}^{*}\bar{\alpha}_{n+1}}{2i}.$
(38)

Moreover, the qubit-bus detunings transform to
$\tilde{\Delta}_{j}=\tilde{\bar{\Delta}}_{j}+\tilde{\Delta}_{j}^{t}$ where

$\displaystyle\tilde{\bar{\Delta}}_{j}=\omega_{j}-\omega_{b}-\delta\left|\bar{\alpha}_{0}\right|^{2}+(\delta+\chi)\left|\bar{\alpha}_{1}\right|^{2}+\frac{\chi\lambda^{2}}{2\left|\bar{\alpha}\right|^{2}},$
(39)
$\displaystyle\tilde{\Delta}_{j}^{t}=-2\chi\lambda\mathrm{Im}\left[\frac{\bar{\alpha}_{1}}{\bar{\alpha}}\right]\cos(\omega_{m}t)+\frac{\chi\lambda^{2}}{2\left|\bar{\alpha}\right|^{2}}\cos(2\omega_{m}t).$
(40)

An additional rotating frame transformation such as to remove the time-
dependence of the qubit-bus detunings, leads to the following approximation
for $\hat{H}_{g}^{P}$

$\displaystyle\hat{H}_{g}^{P\prime}=\sum_{j,n}g_{j}\hat{q}_{j}^{\dagger}e^{i\phi_{n}^{\prime}(t)}\hat{D}_{n,r}\sqrt{n+1}\left|n\right\rangle\left\langle
n+1\right|_{b}+\mathrm{h.c.},$ (41)

where we have introduced
$\phi_{n}^{\prime}(t)=\phi_{n}(t)+\int_{0}^{t}dt\tilde{\Delta}_{j}^{t}$.
Except for the now time-dependent phase $\phi_{n}^{\prime}(t)$, Eq. 41 as the
same form as Eq. 21.

The role of the time-dependent phase $\phi_{n}^{\prime}(t)$ can be understood
by using the Jacobi-Anger expansion

$\displaystyle
e^{i\phi_{n}^{\prime}}=e^{i\bar{\phi}_{n}}\sum_{s=-\infty}^{+\infty}(-i)^{s}J_{s}\left(\lambda\,\mathrm{Re}\left[\frac{\bar{\alpha}_{n+1}-\bar{\alpha}_{n}}{\bar{\alpha}}\right]\right)e^{is\omega_{m}t},$
(42)

where $J_{s}(z)$ is the $s$th Bessel function of the first kind, and where we
considered negligible $\chi/\omega_{m}$ for simplicity. Because the bus mode
is ideally only virtually excited at all times, our goal is to dominantly
suppress the $0\leftrightarrow 1$ transition in the bus. To achieve this, we
adjust the amplitude $\lambda$ to reach a zero of $J_{0}$ in Eq. 42, noticing
that $\bar{\alpha}=\bar{\alpha}_{1}-\bar{\alpha}_{0}$. Higher harmonics of Eq.
42 oscillate rapidly for
$\omega_{m}\gg\left|\chi\right|,|\tilde{\bar{\Delta}}_{j}|$, and result in a
lower bound on the suppression of the two-qubit interactions (c.f. appendix
C).

Figure 6: Time-averaged inverse participation ratio against photon number
using a two-tone drive on the LR ($K_{r}=0$) with frequencies
$\omega_{r}\pm\omega_{m}$ and amplitudes $\omega_{m}\lambda/\bar{\alpha}^{*}$.
The shaded regions are bounded by $\lambda=\lambda_{0}$ where
$J_{0}(\lambda_{0})=0$ and a 10% error on $\lambda_{0}$. The other parameters
are $\delta/2\pi=-1.0$ MHz, $\chi/2\pi=-20.0$ MHz,
$(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz, $(\omega_{2}-\omega_{b})/2\pi=14.0$
MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz, and $g/2\pi=2.0$ MHz.
$\left<1-\mathrm{IPR}_{0100}\right>_{t}$ can be found in appendix C.

To understand how the two proposed implementation mechanisms compare to each
other, we compute in Sect. B.2 the IPR using a time-dependent Schrieffer-Wolff
transformation to find

$\displaystyle 1-\mathrm{IPR}_{1000}\approx
2g_{1}^{2}e^{-\left|\bar{\alpha}\right|^{2}}\sum_{s_{1},s_{2}=-\infty}^{+\infty}i^{s_{2}-s_{1}}e^{i(s_{2}-s_{1})\omega_{m}t}$
$\displaystyle\times
J_{s_{1}}(\lambda)J_{s_{2}}(\lambda)\frac{{}_{4}F_{4}\left(\bm{p}_{1s_{1}s_{2}};1+\bm{p}_{1s_{1}s_{2}};\left|\bar{\alpha}\right|^{2}\right)}{\left(\tilde{\Delta}_{1}+s_{1}\omega_{m}\right)\left(\tilde{\Delta}_{1}+s_{2}\omega_{m}\right)},$
(43)

where
$\bm{p}_{js_{1}s_{2}}=\begin{pmatrix}p_{js_{1}-}&p_{js_{2}-}&p_{js_{1}+}&p_{js_{2}+}\end{pmatrix}$
with
$p_{js\pm}=\beta[1\pm\sqrt{1+\smash[b]{2(\tilde{\Delta}_{j}+s\omega_{m})\beta^{-2}/K_{r}}}]$
and $\beta=(\delta+\chi)/K_{r}-1/2$.We observe that Eq. 43 is reminiscent of
Eq. 28, and a similar expression for $\mathrm{IPR}_{0100}$ can be derived. In
the large $\omega_{m}\gg|\tilde{\Delta}_{1}|$ limit, the dominant contribution
to Eq. 43 arises from the term with $s_{1}=s_{2}=0$, which is canceled by
adjusting $\lambda$ to reach a zero of $J_{0}$. Importantly, because of the
already suppressed interactions, there is no need for a very fine adjustment
of $\lambda$. The time-averaged IPR according to Eq. 43 is illustrated in Fig.
6, where we take $\omega_{m}=\omega_{0}\left|\bar{\alpha}\right|$ such that
the drive amplitude in Eq. 35 is independent of $\left|\bar{\alpha}\right|$.
The dashed line corresponds to the absence of dynamical decoupling. The shaded
regions correspond to $\pm$ 10% error bounds on the drive amplitude. We
observe a strong suppression of 1-IPR1000 and of the drive amplitude
sensitivity. As discussed further in appendix C, we note that the asymptotic
behavior of the suppression is polynomial in $\bar{\alpha}$. With the very
large suppression of the 1-IPR that is observed in Fig. 6, this is a small
price to pay when trading the nonlinearity $K_{r}$ for an additional drive. We
finally note that the suppression can be further enhanced with the help of a
longitudinal drive in the LR (c.f. appendix C).

##  V Superconducting circuit implementation

Figure 7: Superconducting circuit implementation. a) Circuit design. The
qubits and bus modes are implemented using transmon qubits;
$\hat{\varphi}_{1}$, $\hat{\varphi}_{2}$ and $\hat{\varphi}_{b}$ are the phase
operators of $Q_{1}$, $Q_{2}$ and the bus modes respectively. Here
$2\hat{\varphi}_{\pm}=-\hat{\varphi}_{b}\pm\left(\hat{\varphi}_{r}-\varphi_{\ell}\right)$
where $\varphi_{\ell}$ is a real-valued scalar to be defined. A SNAIL-like
element, representing the NLR mode with phase operator $\hat{\varphi}_{r}$, is
linearly driven by a voltage source $V_{g}$. $\Phi_{0}\varphi_{s2}/2\pi$ and
$\Phi_{0}\left(\varphi_{sN}-\varphi_{s2}\right)/2\pi$ are two external fluxes
that control the bus-NLR interaction. b) Metapotential of the NLR for each of
the bus states $\left|0\right\rangle_{b}$ (orange), $\left|1\right\rangle_{b}$
(purple) and $\left|2\right\rangle_{b}$ (blue). Here $\delta/2\pi=-5.0$ MHz,
$\chi/2\pi=-5.0$ MHz and $K_{r}/2\pi=-10.0$ MHz. c) $1-\mathrm{IPR}_{1000}$
(and $1-\mathrm{IPR}_{0100}$ can be found in appendix D) estimated by
numerical diagonalization of the full system using the effective Hamiltonian
in Eq. 55 with $\delta/2\pi=1$ KHz, $\chi/2\pi=-5.0$ MHz,
$(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz, $(\omega_{2}-\omega_{b})/2\pi=14.0$
MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz, and $g/2\pi=2.0$ MHz.

In this section, we introduce a superconducting quantum circuit realizing our
coupler. To approach the model Hamiltonian of Eq. 2, we draw inspiration from
the Kerr-cat qubit, which exploits the bifurcation physics of driven
Josephson-based devices Puri  _et al._ (2017a); Grimm  _et al._ (2020). A
simplification based on the idea of dynamical decoupling presented in Sect. IV
is also discussed.

###  V.1 Kerr-cat-based circuit model

Figure 7(a) shows a possible circuit realization of our coupler with two
transmon qubits interacting through a transmon-like device playing the role of
bus mode. The latter mode is connected to a driven nonlinear circuit
representing the NLR and consisting of a loop formed by two symmetrical
Josephson junction and a SNAIL-like element which incorporates an array of
$N\sim 3$ junctions Frattini  _et al._ (2017). Omitting the qubits, the
Hamiltonian of the circuit reads

$\displaystyle\hat{H}=\hat{H}_{b}+\hat{H}_{r}+\hat{H}_{br},$ (44)

where

$\begin{split}\hat{H}_{b}&=4E_{C_{b}}\hat{n}_{b}^{2}-E_{J_{b}}\cos\left(\hat{\varphi}_{b}\right)\approx\omega_{b}\hat{b}^{\dagger}\hat{b}+\frac{K_{b}}{2}\hat{b}^{\dagger
2}\hat{b}^{2},\end{split}$ (45)

is the bus Hamiltonian, which we treat as a weakly nonlinear oscillator of
frequency $\omega_{b}=\sqrt{8E_{C_{b}}E_{J_{b}}}-E_{C_{b}}$ and anharmonicity
$K_{b}=-E_{C_{b}}$. We define the the phase operators of the two modes across
the SNAIL-like element as
$2\hat{\varphi}_{\pm}=-\hat{\varphi}_{b}\pm\left(\hat{\varphi}_{r}-\varphi_{\ell}\right)$
where $\hat{\varphi}_{b}$ ($\hat{\varphi}_{r}$) is the phase operator of the
bus (NLR) and $\varphi_{\ell}$ is a real-valued scalar determined from the
minimization of the potential energy of the circuit. We consider two external
flux biases: $\varphi_{s2}$ in the three-node loop and
$\varphi_{sN}-\varphi_{s2}$ in the SNAIL-like circuits. Here, $\varphi_{s2}$
($\varphi_{sN}$) shifts the cosine potential of the two junctions (N
junctions) in the SNAIL-like circuit. Moreover, $\hat{\varphi}_{1}$
($\hat{\varphi}_{2}$) is the phase operator of $Q_{1}$ ($Q_{2}$). The NLR and
bus-NLR Hamiltonians take the form

$\displaystyle\hat{H}_{r}=4E_{C_{r}}\hat{n}_{r}^{2}-NE_{J_{N}}\cos\left(\frac{\hat{\varphi}_{r}-\varphi_{sN}}{N}\right)+2\epsilon(t)\hat{n}_{r}$
$\displaystyle\quad-2E_{J_{\ell}}\cos\left(\frac{\hat{\varphi}_{r}-\varphi_{\ell}}{2}\right)-2E_{J_{2}}\cos\left(\frac{\hat{\varphi}_{r}-\varphi_{s2}}{2}\right),$
(46)
$\displaystyle\hat{H}_{br}=-2E_{J_{\ell}}\left[\cos\left(\frac{3\hat{\varphi}_{b}}{2}\right)-1\right]\cos\left(\frac{\hat{\varphi}_{r}-\varphi_{\ell}}{2}\right).$
(47)

In these expressions, $E_{C_{r}}$ is the NLR charging energy, $E_{J_{\ell}}$
the Josephson energy of the NLR’s symmetrical junctions, $E_{J_{2}}$ the
Josephson energy of each of the two small NLR’s junctions, and $E_{J_{N}}$ the
Josephson energy of each of the $N$ large junctions in the array. Moreover,
$\epsilon(t)$ is the amplitude of the drive of frequency
$2(\omega_{r}-\delta)$ on the NLR, where $\omega_{r}$ is the frequency of the
undriven NLR.

The idea is to stabilize cat states in the NLR with amplitudes that depend on
the bus photon number. Just as in the simplified model discussed in the
previous section, transitions between bus states are associated to
displacements in the NLR. An advantage of this propoposed realization is that
the large anharmonicity in the NLR is now determined by the size of the cat
state. To this end, we follow Frattini  _et al._ (2017) by choosing the
external fluxes and Josehpson energy such as to obtain a cubic nonlinearity of
the form $\hat{b}^{\dagger}\hat{b}\left(\hat{r}^{\dagger}+\hat{r}\right)^{3}$
in $\hat{H}_{br}$. In the presence of a linear drive on the NLR, the cubic
nonlinearity leads to a nearly resonant, bus-photon-number-dependent two-
photon drive in the NLR. The Kerr nonlinearity in the NLR can then stabilize
bus-photon-number-dependent cat-states.

More precisely, we take $E_{J_{2}}=E_{J_{\ell}}$, $\varphi_{\ell}=\pi-2\zeta$
and $\varphi_{s2}=-\pi-2\zeta$ where $\zeta$ is a parameter to be defined.
With these choices, we have

$\displaystyle\hat{H}_{r}\approx
4E_{C_{r}}\hat{n}_{r}^{2}-NE_{J_{N}}\cos\left(\frac{\hat{\varphi}_{r}-\varphi_{sN}}{N}\right)+2\epsilon(t)\hat{n}_{r},$
(48)

and

$\displaystyle\hat{H}_{br}\approx$ $\displaystyle\frac{9\pi
z_{b}E_{J_{\ell}}}{4}\left(2\hat{b}^{\dagger}\hat{b}+1\right)$
$\displaystyle\times\left[\cos\zeta\sin\left(\frac{\hat{\varphi}_{r}}{2}\right)+\sin\zeta\cos\left(\frac{\hat{\varphi}_{r}}{2}\right)\right],$
(49)

where $z_{b(r)}=Z_{b(r)}/R_{Q}$ is the reduced impedance of the bus (NLR) mode
with $R_{Q}\simeq 6.5\,\mathrm{k}\Omega$ the resistance quantum. In
$\hat{H}_{br}$, the sine and cosine terms that depend on $\hat{\varphi}_{r}$
are key for implementing the bus-photon-number-dependent cubic nonlinearity
$\hat{b}^{\dagger}\hat{b}\left(\hat{r}^{\dagger}+\hat{r}\right)^{3}$ and the
cross-Kerr interaction $\hat{b}^{\dagger}\hat{b}\hat{r}^{\dagger}\hat{r}$. In
presence of the drive $\epsilon(t)=\epsilon_{0}\sin[2(\omega_{r}-\delta)t]$
where, for simplicity, we take a constant drive amplitude $\epsilon_{0}$, our
next step is to apply a displacement transformation $\hat{D}[\xi(t)]$ on the
NLR mode to eliminate the drive term. To achieve this, we take
$-i\dot{\xi}+\omega_{r}\xi+i\epsilon(t)/\sqrt{\pi z_{r}}=0$ or, equivalently,
$\sqrt{\pi
z_{r}}\xi(t)\approx-\Omega\cos[2(\omega_{r}-\delta)t]+(i\Omega/2)\sin[2(\omega_{r}-\delta)t]$
with the displacement amplitude $\Omega=2\epsilon_{0}/3\omega_{r}$. By doing
so, we obtain the displaced Hamiltonian
$\hat{H}^{D}=\hat{H}_{b}+\hat{H}_{r}^{D}+\hat{H}_{br}^{D}$ with

$\displaystyle\hat{H}_{r}^{D}\approx\omega_{r}\hat{r}^{\dagger}\hat{r}+\frac{K_{r}}{2}\hat{r}^{\dagger
2}\hat{r}^{2}+\frac{\lambda_{r}}{2}\hat{r}^{\dagger
2}e^{-i2(\omega_{r}-\delta)t}+\text{h.c.}$ (50)
$\displaystyle\hat{H}_{br}^{D}\approx\chi\left(\hat{b}^{\dagger}\hat{b}+\frac{1}{2}\right)\hat{r}^{\dagger}\hat{r}$
$\displaystyle\phantom{\hat{H}_{br}^{D}=}+\frac{\lambda_{\ell}}{2}\left(2\hat{b}^{\dagger}\hat{b}+1\right)\hat{r}^{\dagger
2}e^{-i2(\omega_{r}-\delta)t}+\text{h.c.}$ (51)

for small $\varphi_{sN}/N$. The above expressions are valid for small reduced
mode impedance $\pi z_{r}\approx\sqrt{2NE_{C_{r}}/E_{J_{N}}}$ assume the
rotating-wave approximation. We have introduced the NLR frequency
$\omega_{r}=\sqrt{8E_{C_{r}}E_{J_{N}}/N}-E_{C_{r}}$, the self-Kerr
anharmonicity $K_{r}=-E_{C_{r}}/N^{2}$, the two-photon drive amplitude

$\displaystyle\lambda_{r}=-\frac{(\pi z_{r})^{3/2}\Omega
E_{J_{N}}\varphi_{sN}}{2N^{3}},$ (52)

the cross-Kerr interaction amplitude

$\displaystyle\chi=-\frac{9\pi z_{b}\pi z_{r}E_{J_{\ell}}\sin(\zeta)}{8},$
(53)

and the bus-number dependent two-photon drive amplitude

$\displaystyle\lambda_{\ell}=\frac{9\pi z_{b}(\pi z_{r})^{3/2}\Omega
E_{J_{\ell}}\cos(\zeta)}{64}.$ (54)

We moreover set $\lambda_{r}=-2\lambda_{\ell}$ by a proper choice of the flux
biases.

Finally, in a doubly rotating frame at $\omega_{r}-\delta$ for $\hat{r}$ and
at $\omega_{b}$ for $\hat{b}$, the displaced Hamiltonian reads

$\begin{split}&\hat{H}^{\prime D}\approx\frac{K_{b}}{2}\hat{b}^{\dagger
2}\hat{b}^{2}-\frac{K_{r}\alpha^{4}}{2}\left(2\hat{b}^{\dagger}\hat{b}-1\right)^{2}\\\
&\phantom{\hat{H}^{D}}+\left(\delta+\chi/2+\chi\hat{b}^{\dagger}\hat{b}\right)\hat{r}^{\dagger}\hat{r}\\\
&\phantom{\hat{H}^{D}}+\frac{K_{r}}{2}\left[\hat{r}^{\dagger
2}+\alpha^{2}\left(2\hat{b}^{\dagger}\hat{b}-1\right)\right]\left[\hat{r}^{2}+\alpha^{2}\left(2\hat{b}^{\dagger}\hat{b}-1\right)\right],\end{split}$
(55)

where $\alpha^{2}=\lambda_{\ell}/K_{r}$.

The metapotential associated to this Hamiltonian is illustrated in Fig. 7(b).
As in Fig. 2 for the simplified system, the different Fock states of the bus
mode lead to displaced wells in the $I$-$Q$ plane. However, because of the
combination of the Kerr nonlinearity and the engineered two-photon drive, each
Fock state is associated two metapotential wells Puri  _et al._ (2017a). The
central idea of blocking bus-state transition by entangling those state to
coherent states in the NLR is, however, unchanged. This is confirmed in Fig.
7(c) which shows 1-IPR as a function of photon number. Apart from additional
resonance which can easily be avoided, the overall behavior is the one
expected: we see a an exponential reduction of the bus-state hybridization
with NLR photon number, as originally predicted by the model of Eq. 2.

Furthermore, it is useful to note that the ac-Stark shift on the bus frequency
vanishes in this model, i.e. the ground- and first-excited states of the bus
are both shifted in energy by $K_{r}\alpha^{4}/2$. The bus-NRL entanglement is
therefore entirely responsible for the exponential suppression observed in
Fig. 7c).

###  V.2 Harmonic model with parametric modulation

In Sect. IV we have seen how it is possible to trade the large nonlinear
interaction between the bus and the NLR by additional drives. Here, we show
how this idea can be realized without modifications to the circuit of Fig. 7.
For this second approach, we take $\varphi_{\ell}=0$, $\varphi_{s2}=2\pi$,
$E_{J_{2}}=\lambda E_{J_{\ell}}$ with $\lambda=1-(3/2)^{2}\pi z_{b}/2$, and
$\varphi_{sN}\ \mathrm{mod}\ 2\pi N=0$. With these parameter choices, the
circuit Hamiltonian can now be written as

$\displaystyle\hat{H}_{b}=4E_{C_{b}}\hat{n}_{b}^{2}-E_{J_{b}}\cos\left(\hat{\varphi}_{b}\right)-2E_{J_{\ell}}\cos\left(\frac{3\hat{\varphi}_{b}}{2}\right),$
(56)
$\displaystyle\hat{H}_{r}=4E_{C_{r}}\hat{n}_{r}^{2}-NE_{J_{N}}\cos\left(\frac{\hat{\varphi}_{r}}{N}\right)+2\epsilon(t)\hat{n}_{r}.$
(57)
$\displaystyle\hat{H}_{br}=-2E_{J_{\ell}}\left[\cos\left(\frac{3\hat{\varphi}_{b}}{2}\right)-\lambda\right]\left[\cos\left(\frac{\hat{\varphi}_{r}}{2}\right)-1\right],$
(58)

We note that the reduced mode impedance of the NLR, $\pi
z_{r}\approx\sqrt{2NE_{C_{r}}/E_{J_{N}}}$, needs to be made small to prevent
the drive on the NLR from resulting in appreciable nonlinear terms due to the
cosine potentials and the nearly resonant two-photon and cubic terms. More
precisely, we take $1/\sqrt{\pi z_{r}}$ to be much larger than any
displacement in the NLR field associated with the bus Fock states $n\neq 0$,
and small compared to $N/\sqrt{\pi z_{r}}$ for $n=0$.

As above, we treat the bus and NLR as weakly nonlinear oscillators.
Importantly, for $n=0$, we find that $\hat{H}_{br}\approx 0$ and $\hat{H}_{r}$
is approximately harmonic despite a large displacement in the NLR field.
Moreover, for $n\neq 0$, $\hat{H}_{br}$ effectively implements a large cross-
Kerr interaction that strongly reduces displacements of the NLR field by
rendering the linear drive of the NLR off-resonant. Particularly, we find that

$\begin{split}&\hat{H}\approx\sum_{\nu=b,r}\left(\omega_{\nu}\hat{\nu}^{\dagger}\hat{\nu}+\frac{K_{\nu}}{2}\hat{\nu}^{\dagger
2}\hat{\nu}^{2}\right)+\chi\hat{b}^{\dagger}\hat{b}\hat{r}^{\dagger}\hat{r}\\\
&\phantom{\hat{H}\approx}+i\Omega(t)\left(\hat{r}^{\dagger}-\hat{r}\right).\end{split}$
(59)

In this expression, we have defined the frequencies
$\omega_{\nu}\approx\sqrt{8E_{C_{\nu}}E_{L_{\nu}}}-E_{C_{\nu}}$, the reduced
mode impedances $\pi z_{\nu}\approx\sqrt{2E_{C_{\nu}}/E_{L_{\nu}}}$, the
inductive energies $E_{L_{b}}=E_{J_{b}}+2(3/2)^{2}E_{J_{\ell}}$ and
$E_{L_{r}}=E_{J_{N}}/N$, the anharmonicities
$K_{b}=-[E_{J_{b}}+2(3/2)^{4}E_{J_{\ell}}]E_{C_{b}}/E_{L_{b}}$ and
$K_{r}=E_{C_{r}}/N^{2}$, and the linear drive amplitude
$\Omega(t)=\epsilon(t)/2\sqrt{\pi z_{r}}$.

We also largely reduce the NLR’s anharmonicity $K_{r}$ by choosing a small
$E_{C_{r}}$ and large $N$. It is possible with this model to observe the
exponential suppression of two-qubit coupling by choosing $\epsilon(t)$ to be
nearly resonant with the NLR.

As discussed in Sect. IV, with a reduced anharmonicity in the NLR an
additional drive, which we choose to be of the form

$\displaystyle\Omega(t)=2\delta\bar{\alpha}_{0}\cos\left[(\omega_{r}-\delta)t\right]+\frac{\omega_{m}\lambda}{\bar{\alpha}}\sum_{\nu=\pm}\nu\sin\left[(\omega_{m}+\nu\omega_{r})t\right]$
(60)

can serve as a complementary mechanism to suppress interactions. In the limit
$\left|\delta/\chi\right|\ll 1$, $\left|\chi/\omega_{m}\right|\ll 1$ and
$\left|\omega_{m}/\omega_{r}\right|\ll 1$, we find that the bus-state
dependent displacements take the form

$\displaystyle\alpha_{n}(t)$
$\displaystyle\approx\frac{\delta\bar{\alpha}_{0}}{\delta+n\chi}-\frac{i\lambda}{\bar{\alpha}}\cos(\omega_{m}t),$
(61)

in agreement with Eq. 36. From this point on, the result of Sect. IV follows.

Finally, we emphasize that even though a cross-Kerr type interaction between
the bus and the NLR could, in principle, be implemented using a dispersive
coupling Blais  _et al._ (2021), the dispersive Hamiltonian is invalid at
large photon numbers and yields virtual qubit-qubit interactions through the
driven NLR. We also note that a discussion of the leading effects of stray
couplings can be found in appendix D.

##  VI Conclusion

We introduced a two-qubit coupler with an exponential on-off ratio, realized
by connecting a pair of qubits to a bus mode complemented by a driven
ancillary resonator. The cross-Kerr interaction between the bus and the driven
resonator results in a displacement of the resonator’s field that is
conditional on the bus state. Because the displaced resonator states have
negligible overlap, bus-state transitions are suppressed exponentially in the
amplitude of the drive. In turn, because two-qubit interactions are mediated
by bus transitions, the two-qubit coupling also results strongly suppressed,
leading to a high on-off ratio. As a clear demonstration of this mechanism, we
have shown how the inverse participation ratio, which is a measure of qubit-
bus hybridization, and the spurious cross-Kerr between the qubits are
exponentially reduced with the number of photons in the resonator mode. We
also proposed two complementary superconducting quantum circuit
implementations of our coupler.

The strong reduction in two-qubit couplings demonstrated here can be
advantageous in multiqubit processors, where spectator qubits and long-range
qubit-qubit interactions can have detrimental effects Berke  _et al._ (2020).
For the same reason, the proposed approach can be particularly useful in all-
microwave frequency-fixed qubits architectures with interactions mediated by
frequency-fixed buses. We hope that the mechanisms explored here will pave the
way to a new generation of coupling schemes for diverse platforms.

Finally, we note that possible improvements to the coupler include squeezing
the resonator mode to further reduce the overlaps between the resonator states
associated with distinct bus states, and extending the ancillary system to
multiple modes such that the exponential suppression is now with respect to
multiple modes.

## Acknowledgments

We thank Ross Shillito and Jens Koch for useful discussions. This work was
undertaken in part thanks to funding from NSERC, the Canada First Research
Excellence Fundm the Ministère de l’économie et de l’innovation du Québec and
the U.S. Army Research Office Grant No. W911NF-18-1-0411

## References

  * Niskanen  _et al._ (2007) A. O. Niskanen, K. Harrabi, F. Yoshihara, Y. Nakamura, S. Lloyd, and J. S. Tsai, Science 316, 723 (2007), https://science.sciencemag.org/content/316/5825/723.full.pdf .
  * Majer  _et al._ (2007) J. Majer, J. M. Chow, J. M. Gambetta, J. Koch, B. R. Johnson, J. A. Schreier, L. Frunzio, D. I. Schuster, A. A. Houck, A. Wallraff, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Nature 449, 443 EP (2007).
  * Frattini  _et al._ (2017) N. E. Frattini, U. Vool, S. Shankar, A. Narla, K. M. Sliwa, and M. H. Devoret, Applied Physics Letters 110, 222603 (2017), https://doi.org/10.1063/1.4984142 .
  * Yan  _et al._ (2018a) F. Yan, P. Krantz, Y. Sung, M. Kjaergaard, D. L. Campbell, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Phys. Rev. Applied 10, 054062 (2018a).
  * Sete  _et al._ (2021a) E. A. Sete, A. Q. Chen, R. Manenti, S. Kulshreshtha, and S. Poletto, Phys. Rev. Applied 15, 064063 (2021a).
  * Finck  _et al._ (2021) A. D. K. Finck, S. Carnevale, D. Klaus, C. Scerbo, J. Blair, T. G. McConkey, C. Kurter, A. Carniol, G. Keefe, M. Kumph, and O. E. Dial, “Suppressed crosstalk between two-junction superconducting qubits with mode-selective exchange coupling,” (2021), arXiv:2105.11495 [quant-ph] .
  * Casparis  _et al._ (2019) L. Casparis, N. J. Pearson, A. Kringhøj, T. W. Larsen, F. Kuemmeth, J. Nygård, P. Krogstrup, K. D. Petersson, and C. M. Marcus, Phys. Rev. B 99, 085434 (2019).
  * Christiane and Homann (2013) E. Christiane and M. Homann (2013).
  * Chen  _et al._ (2014) Y. Chen, C. Neill, P. Roushan, N. Leung, M. Fang, R. Barends, J. Kelly, B. Campbell, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey, A. Megrant, J. Y. Mutus, P. J. J. O’Malley, C. M. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, M. R. Geller, A. N. Cleland, and J. M. Martinis, Phys. Rev. Lett. 113, 220502 (2014).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * McKay  _et al._ (2016) D. C. McKay, S. Filipp, A. Mezzacapo, E. Magesan, J. M. Chow, and J. M. Gambetta, Phys. Rev. Applied 6, 064007 (2016).
  * Roth  _et al._ (2017) M. Roth, M. Ganzhorn, N. Moll, S. Filipp, G. Salis, and S. Schmidt, Phys. Rev. A 96, 062323 (2017).
  * kun Xie  _et al._ (2018) J. kun Xie, S. li Ma, Z. peng Yang, Z. Li, and F. li Li, Physics Letters A 382, 2626 (2018).
  * Kounalakis  _et al._ (2018) M. Kounalakis, C. Dickel, A. Bruno, N. K. Langford, and G. A. Steele, npj Quantum Information 4, 38 (2018).
  * Chu and Yan (2021) J. Chu and F. Yan, “Coupler-assisted controlled-phase gate with enhanced adiabaticity,” (2021), arXiv:2106.00725 [quant-ph] .
  * Collodo  _et al._ (2020) M. C. Collodo, J. Herrmann, N. Lacroix, C. K. Andersen, A. Remm, S. Lazar, J.-C. Besse, T. Walter, A. Wallraff, and C. Eichler, Phys. Rev. Lett. 125, 240502 (2020).
  * Sung  _et al._ (2021) Y. Sung, L. Ding, J. Braumüller, A. Vepsäläinen, B. Kannan, M. Kjaergaard, A. Greene, G. O. Samach, C. McNally, D. Kim, A. Melville, B. M. Niedzielski, M. E. Schwartz, J. L. Yoder, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Phys. Rev. X 11, 021058 (2021).
  * Sete  _et al._ (2021b) E. A. Sete, N. Didier, A. Q. Chen, S. Kulshreshtha, R. Manenti, and S. Poletto, “Parametric-resonance entangling gates with a tunable coupler,” (2021b), arXiv:2104.03511 [quant-ph] .
  * Groszkowski  _et al._ (2011) P. Groszkowski, A. G. Fowler, F. Motzoi, and F. K. Wilhelm, Phys. Rev. B 84, 144516 (2011).
  * Yan  _et al._ (2018b) F. Yan, P. Krantz, Y. Sung, M. Kjaergaard, D. L. Campbell, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Phys. Rev. Applied 10, 054062 (2018b).
  * Mundada  _et al._ (2019) P. Mundada, G. Zhang, T. Hazard, and A. Houck, Phys. Rev. Applied 12, 054023 (2019).
  * Ku  _et al._ (2020) J. Ku, X. Xu, M. Brink, D. C. McKay, J. B. Hertzberg, M. H. Ansari, and B. L. T. Plourde, Phys. Rev. Lett. 125, 200504 (2020).
  * Zhao  _et al._ (2020) P. Zhao, P. Xu, D. Lan, J. Chu, X. Tan, H. Yu, and Y. Yu, Phys. Rev. Lett. 125, 200503 (2020).
  * Blais  _et al._ (2004) A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 69, 062320 (2004).
  * Blais  _et al._ (2021) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * Cochrane  _et al._ (1999) P. T. Cochrane, G. J. Milburn, and W. J. Munro, Phys. Rev. A 59, 2631 (1999).
  * Puri  _et al._ (2017a) S. Puri, S. Boutin, and A. Blais, npj Quantum Information 3, 18 (2017a).
  * Mirrahimi  _et al._ (2014) M. Mirrahimi, Z. Leghtas, V. V. Albert, S. Touzard, R. J. Schoelkopf, L. Jiang, and M. H. Devoret, New Journal of Physics 16, 045014 (2014).
  * Koch  _et al._ (2007) J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 76, 042319 (2007).
  * Puri  _et al._ (2017b) S. Puri, S. Boutin, and A. Blais, npj Quantum Information 3, 18 (2017b).
  * Berry (2009) M. V. Berry, Journal of Physics A: Mathematical and Theoretical 42, 365303 (2009).
  * Gambetta  _et al._ (2006) J. Gambetta, A. Blais, D. I. Schuster, A. Wallraff, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 74, 042318 (2006).
  * Berke  _et al._ (2020) C. Berke, E. Varvelis, S. Trebst, A. Altland, and D. P. DiVincenzo, “Transmon platform for quantum computing challenged by chaotic fluctuations,” (2020), arXiv:2012.05923 [quant-ph] .
  * Evers and Mirlin (2008) F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
  * Grimm  _et al._ (2020) A. Grimm, N. E. Frattini, S. Puri, S. O. Mundhada, S. Touzard, M. Mirrahimi, S. M. Girvin, S. Shankar, and M. H. Devoret, Nature 584, 205 (2020).

##  Appendix A Polaron transformation

To gain intuition about the underlying physics of the model, it is useful to
move to a frame defined by the time-dependent polaron transformation

$\displaystyle\hat{P}(t)=e^{-i\int_{0}^{t}d\tau[\omega_{b}+\Delta_{\rm
ac}(\tau)]\left(\hat{b}^{\dagger}\hat{b}+\sum_{j}\hat{q}_{j}^{\dagger}\hat{q}_{j}\right)}e^{-i\omega_{d}t\hat{r}^{\dagger}\hat{r}}\sum_{n=0}^{\infty}\hat{D}_{r}[\alpha_{n}(t)]\otimes|n\rangle\langle
n|_{b},$ , (62)

where $|n\rangle\langle n|_{b}$ is the projection operator associated with the
eigenstate $|n\rangle$ of the bus mode, and $\Delta_{\rm ac}$ is an ac-Stark
shift that will be defined below. The displacements $\\{\alpha_{n}(t)\\}$ are
determined from the damped classical equation of the NLR [c.f. Eq. 9] that
result from the bilinear Hamiltonian terms only.

###  A.1 Transformed Hamiltonian

Under Eq. 62 the Hamiltonian Eq. 2 transforms to $\hat{H}^{\rm
P}=\hat{P}^{\dagger}\hat{H}\hat{P}-i\hat{P}^{\dagger}\dot{\hat{P}}$, where

$\displaystyle\hat{H}^{P}=\sum_{j=1}^{2}\hat{H}_{j}^{P}+\hat{H}_{br}^{P}+\hat{H}_{\kappa}^{P}+\hat{H}_{g}^{P},$
(63) $\displaystyle\hat{H}_{j}^{P}=(\omega_{j}-\omega_{b}-\Delta_{\rm
ac})\hat{q}_{j}^{\dagger}\hat{q}_{j}+\frac{K_{j}}{2}\hat{q}_{j}^{\dagger
2}\hat{q}_{j}^{2},$ (64)
$\displaystyle\hat{H}_{br}^{P}=\left(\delta+\chi\hat{b}^{\dagger}\hat{b}\right)\hat{r}^{\dagger}\hat{r}+\frac{K_{b}}{2}\hat{b}^{\dagger
2}\hat{b}^{2}+\frac{K_{r}}{2}\hat{r}^{\dagger
2}\hat{r}^{2}+\sum_{n}\Delta_{\rm ac,n}\left|n\right\rangle\left\langle
n\right|_{b},$ (65)
$\displaystyle\hat{H}_{\kappa}^{P}=\frac{i\kappa}{2}\sum_{n}\left(\alpha_{n}\hat{r}^{\dagger}-\alpha_{n}^{*}\hat{r}\right)\left|n\right\rangle\left\langle
n\right|_{b},$ (66)
$\displaystyle\hat{H}_{g}^{P}=\sum_{j,n}g_{j}\left(\hat{q}_{j}^{\dagger}e^{i\phi_{n}}\hat{D}_{n,r}\sqrt{n+1}\left|n\right\rangle\left\langle
n+1\right|_{b}+\mathrm{h.c.}\right),$ (67)

with $i2\phi_{n}=\alpha_{n+1}^{*}\alpha_{n}-\alpha_{n}^{*}\alpha_{n+1}$,
$\hat{D}_{n,r}=\hat{D}_{r}\left(\alpha_{n+1}-\alpha_{n}\right)$,
$\hat{D}_{r}(\alpha)=e^{\alpha\hat{r}^{\dagger}-\alpha^{*}\hat{r}}$ is the
displacement operator in the NLR, and

$\displaystyle\Delta_{\rm
ac,n}=\delta\left|\alpha_{0}\right|^{2}-(\delta+n\chi)\left|\alpha_{n}\right|^{2}-n\Delta_{\rm
ac},$ (68) $\displaystyle\Delta_{\rm
ac}=\delta\left|\alpha_{0}\right|^{2}-(\delta+\chi)\left|\alpha_{1}\right|^{2}.$
(69)

As we only consider Jaynes-Cummings-type interactions, the transformed
Hamiltonian can be reduced to the form in Eq. 17. In addition, driving the NLR
mode introduces an ac-Stark shift $\Delta_{\rm
ac}^{C}(t)=\Delta_{1}^{C}(t)-\Delta_{0}^{C}(t)$ of the bus frequency given by

$\displaystyle\Delta_{\rm
ac}^{C}(t)=\delta\left|\alpha_{0}(t)\right|^{2}-(\delta+\chi)\left|\alpha_{1}(t)\right|^{2}.$
(70)

###  A.2 Transformed Master equation

The full system dynamics can be described by the Lindblad Master equation
formalism,

$\displaystyle\dot{\hat{\rho}}=-i\left[\hat{H},\hat{\rho}\right]+\sum_{j}\hat{L}_{j}\hat{\rho}\hat{L}_{j}^{\dagger}-\frac{1}{2}\left\\{\hat{L}_{j}^{\dagger}\hat{L}_{j},\hat{\rho}\right\\}$
(71)

where $\hat{\rho}$ is the density matrix of the system, $\hat{H}$ is the
Hamiltonian and $\hat{L}_{j}$ are the collapse operators. Under the
transformation Eq. 62 the density matrix transforms as
$\hat{\rho}=\hat{P}\hat{\rho}^{P}\hat{P}^{\dagger}~{}$. It follows that

$\displaystyle\dot{\hat{\rho}}^{P}=-i\left[\hat{H}^{P},\hat{\rho}^{P}\right]+\sum_{j}\hat{L}_{j}^{P}\hat{\rho}^{P}\hat{L}_{j}^{P\dagger}-\frac{1}{2}\left\\{\hat{L}_{j}^{P\dagger}\hat{L}_{j}^{P},\hat{\rho}^{P}\right\\},$
(72)

where we defined the transformed collapse operators
$\hat{L}_{j}^{P}=\hat{P}^{\dagger}\hat{L}_{j}\hat{P}$. As examples, the
collapse operators can take the form $\hat{L}_{r}=\sqrt{\kappa}\hat{r}$,
$\hat{L}_{\nu,\gamma}=\sqrt{\gamma_{\nu}}\left|0\right\rangle\left\langle
1\right|_{\nu}$ and
$\hat{L}_{\nu,\varphi}=\sqrt{\gamma_{\varphi,\nu}/2}\left(\left|1\right\rangle\left\langle
1\right|_{\nu}-\left|0\right\rangle\left\langle 0\right|_{\nu}\right)$ for
$\nu=\\{b,1,2\\}$. We find that qubit collapse operators as well as
$\hat{L}_{b,\varphi}$ are unchanged under the polaron transformation. However,
we have

$\displaystyle\hat{L}_{b,\gamma}^{P}=\sqrt{\gamma_{b}}e^{-i\phi_{0}}e^{-i\int_{0}^{t}d\tau\left[\omega_{b}+\Delta_{\rm
ac}(\tau)\right]}\hat{D}_{0,r}^{\dagger}\left|0\right\rangle\left\langle
1\right|_{b},$ (73)
$\displaystyle\hat{L}_{r}^{P}=\sqrt{\kappa}\left(\hat{r}+\sum_{n}\alpha_{n}\left|n\right\rangle\left\langle
n\right|_{b}\right).$ (74)

from where it follows that $\hat{L}_{b,\gamma}^{P}$ is exponentially
suppressed because of the displacement operator in the NLR. Irrespective of
this observation, it is worth nothing that $\hat{L}_{b,\gamma}^{P}$ does not
prevent the formation of the polaronic states in the coupler and therefore
does not hinder the proposed protocol. In addition, $\hat{L}_{r}^{P}$
corresponds to measurement-induced dephasing in the coupler. It is possible to
further simplify the master equation within the rotating-wave approximation to

$\displaystyle\dot{\hat{\rho}}^{P}=-i\left[\hat{H}^{P}-\hat{H}_{\kappa}^{P},\hat{\rho}^{P}\right]+\sum_{j}\hat{\tilde{L}}_{j}^{P}\hat{\rho}^{P}\hat{\tilde{L}}_{j}^{P\dagger}-\frac{1}{2}\left\\{\hat{\tilde{L}}_{j}^{P\dagger}\hat{\tilde{L}}_{j}^{P},\hat{\rho}^{P}\right\\},$
(75)

where $\hat{\tilde{L}}_{j}^{P}=\hat{L}_{j}^{P}$, except
$\hat{\tilde{L}}_{r}^{P}=\sqrt{\kappa}\hat{r}$ and we define a new collapse
operator
$\tilde{\hat{L}}_{b,\gamma_{\alpha}}=\sqrt{\kappa}\sum_{n}\alpha_{n}\left|n\right\rangle\left\langle
n\right|_{b}$ which captures measurement-induced dephasing in the bus.

###  A.3 Transitionless driving

Large conditional displacements in the NLR (‘off’ state) can be prepared by
controlling the phase of the envelope in Eq. 10 in time. The same mechanism
makes it possible to empty the NLR quickly (‘on’ state). Imperfections in the
envelope lead to deviations in the intended displacements, which we
characterize in this section. Using Eq. 9 we find the displacements

$\displaystyle\alpha_{n}(t)=i\int_{0}^{t}dz\
\varepsilon(z)e^{i(\delta-i\kappa/2+n\chi)(z-t)}+\alpha_{n}(0)e^{-i(\delta-i\kappa/2+n\chi)t}.$
(76)

Extending the envelope in Eq. 10 to include both switching-off and switching-
on events

$\displaystyle\varepsilon(t)=$
$\displaystyle\left(\varepsilon_{0}(t)-\frac{i\dot{\varepsilon}_{0}(t)}{\delta-i\kappa/2}\right)\Theta(\tau-t)+\varepsilon_{0}(\tau)\Theta(t-\tau)\Theta(T+\tau-t)$
$\displaystyle+\left(\varepsilon_{0}(\tau)-\varepsilon_{0}(t-T-\tau)+\frac{i\dot{\varepsilon}_{0}(t-T-\tau)}{\delta-i\kappa/2}\right)\Theta(t-T-\tau)\Theta(T+2\tau-t),$
(77)

where $\tau$ is the ramping time to switch off/on the device and $T$ is the
time during which the drive is on, we find that

$\displaystyle\begin{split}\alpha_{n}(t)=&\frac{1}{\delta-i\kappa/2+n\chi}\begin{cases}0,&t=0\\\
\varepsilon_{0}(t),&0\leq t\leq\tau\\\ \varepsilon_{0}(\tau),&\tau\leq
t\leq\tau+T\\\ \varepsilon_{0}(\tau)-\varepsilon_{0}(t-\tau-T),&\tau+T\leq
t\leq 2\tau+T\\\ 0,&t\leq 2\tau+T\end{cases}\\\
&-\frac{in\chi}{\delta-i\kappa/2}\left(\int_{0}^{t}dz\Theta(\tau-t)-\int_{T+\tau}^{t}dz\Theta(T+2\tau-t)\right)\frac{d\varepsilon_{0}(z)}{dz}\frac{e^{i(\delta-i\kappa/2+n\chi)(z-t)}}{(-i)(\delta-i\kappa/2+n\chi)}\\\
=&\frac{1}{\delta-i\kappa/2+n\chi}\begin{cases}0,&t=0\\\
\varepsilon_{0}(t),&0\leq t\leq\tau\\\ \varepsilon_{0}(\tau),&\tau\leq
t\leq\tau+T\\\ \varepsilon_{0}(\tau)-\varepsilon_{0}(t-\tau-T),&\tau+T\leq
t\leq 2\tau+T\\\ 0,&t\leq 2\tau+T\end{cases}\\\
&-\frac{n\chi}{\delta-i\kappa/2}\sum_{k=1}^{\infty}\frac{d^{k}\varepsilon_{0}(z)}{dz^{k}}\left.\frac{e^{i(\delta-i\kappa/2+n\chi)(z-x)}}{(-i)^{k}(\delta-i\kappa/2+n\chi)^{k+1}}\right|_{z=0}^{z=x}\left(\delta(t-x)\Theta(\tau-t)-\delta(t-\tau-
T-x)\Theta(2\tau+T-t)\right),\end{split}$ (78)

where we assumed $\alpha_{n}(0)=0$. It follows that, for

$\displaystyle\left|\frac{d^{k}\varepsilon_{0}(z)}{dz^{k}}\Big{|}_{z=0}\right|,\
\left|\frac{d^{k}\varepsilon_{0}(z)}{dz^{k}}\Big{|}_{z=\tau}\right|\ll\left|\delta-i\kappa/2\right|\left|\chi\right|^{k},$
(79)

Eq. 78 simplifies to Eq. 11 for $\tau\leq t\leq\tau+T$, and vanishes for
$t\geq 2\tau+T$. Ultimately, the derivatives of the pulse at the endpoints of
the ramp would contribute the most to deviations in the conditional
displacements $\alpha_{n}$. If one has perfect control over the pulse envelope
the ramping time can be made arbitrarily small but limitations could arise
from pulse imperfections. To see this, we define the perturbed envelope

$\displaystyle\varepsilon(t)\rightarrow\varepsilon(t)+\varepsilon_{\rm
err}(t),$ (80)

where $\varepsilon_{\rm err}(t)$ is a small time-dependent perturbation. Using
integration by parts, we find that

$\displaystyle\alpha_{n}(t)\rightarrow\alpha_{n}(t)+i\sum_{k=0}^{\infty}\left.\frac{d^{k}\varepsilon_{\rm
eff}(z)}{dz^{k}}\frac{e^{i(\delta-i\kappa/2+n\chi)(z-t)}}{(-i)^{k}\left(\delta-i\kappa/2+n\chi\right)^{k}}\right|_{0}^{t}.$
(81)

The effects of non-adiabatic errors are quantified by the ratio between the
time derivatives of the drive envelope at the endpoints and powers of
$\left|\delta-i\kappa/2+n\chi\right|$. These errors result in time-dependent
fluctuations of $\alpha_{n}(t)$, which can change the conditional
displacements $\alpha_{n}$ and the ac-Stark shift of the bus. When switching
on, TQD errors could result in residual photons in the NLR and ac-Stark shifts
in the bus that can affect two-qubit interactions in the ‘on’ state. It is
therefore desirable to have a reset scheme for the NLR.

##  Appendix B Numerical experiments and analytical estimates

In this section, we provide details regarding the numerical simulations, the
derivations for the analytical estimates associated with the inverse
participation ratio and the spurious two-qubit interactions. We also report
additional numerical results for the inverse participation ratio and
measurement-induced dephasing.

###  B.1 Rabi drive experiment

The envelope in Eq. 10 used in Fig. 3 is shown in Fig. 8. The ramping time
$\tau$ was set to 5 ns independently of $\bar{\alpha}_{0}$.

Figure 8: Envelope $\varepsilon(t)$ in Eq. 10 used for Fig. 3 to turn off the
coupler. The drive can be turned off after some arbitrary time in order to
turn the coupler back on with the time-reversed pulse shape shown here [c.f.
Sect. A.3].

###  B.2 Inverse participation ratio (IPR)

In analogy to Fig. 4, the inverse participation ratio for the second qubit,
IPR0100, is computed numerically and reported in Fig. 9. The observations made
in Fig. 4 can be extended to Fig. 9. The only difference here is the emergence
of a second divergence at small photon numbers for $K_{r}=0$. This peak
results from the frequency collisions with higher energy levels in the NLR and
the specific choice of parameters. This effect is however absent in the
presence of anharmonicity in the NLR.

We now provide an analytical estimate for the inverse participation ratio
based on a Schrieffer-Wolff (SW) transformation on Eq. 17, where the
Hamiltonian takes the form

$\displaystyle\hat{H}_{I}=e^{i\int
dt\left(\hat{H}^{P}-\hat{H}_{\kappa}^{P}\right)}\ \hat{H}_{g}^{P}\ e^{-i\int
dt\left(\hat{H}^{P}-\hat{H}_{\kappa}^{P}\right)}.$ (82)

To this end we define the generator

$\displaystyle\hat{S}_{I}=i\int dt\ \hat{H}_{I},$ (83)

under which the Hamiltonian transforms to
$\hat{H}_{S}=e^{\hat{S}_{I}}\hat{H}_{I}e^{-\hat{S}_{I}}+i\dot{\hat{S}}_{I}=\mathcal{O}(g_{j}^{2})$.
In what follows we describe different cases where $\hat{H}_{g}^{P}$ and
$\hat{H}^{P}-\hat{H}_{\kappa}^{P}$ have a particular time dependence used for
parametric modulations in appendix C. Furthermore, we compute the generator
back in the polaron frame as

$\displaystyle\hat{S}=e^{-i\int
dt\left(\hat{H}^{P}-\hat{H}_{\kappa}^{P}\right)}\ \hat{S}_{I}\ e^{i\int
dt\left(\hat{H}^{P}-\hat{H}_{\kappa}^{P}\right)}.$ (84)

This transformation holds for $||\hat{S}||\ll 1$, i.e. if the transition
amplitudes are much smaller in magnitude than the energy gaps. Then, the
hybridized states are approximately given by
$\left|\psi_{h,\nu}\right\rangle=e^{\hat{S}}\left|\psi_{b,\nu}\right\rangle\approx[1+\hat{S}+\mathcal{O}(g_{j}^{2})]\left|\psi_{b,\nu}\right\rangle$,
with $\left|\psi_{b,\nu}\right\rangle$ the bare eigenstates of the full system
for $g_{j}=0$. An estimation of Eq. 27 using the SW transformation follows as

$\displaystyle\mathrm{IPR}_{\mu}\approx
1-2\left\langle\psi_{b,\mu}|\hat{S}^{\dagger}\hat{S}|\psi_{b,\mu}\right\rangle.$
(85)

However, as $\hat{S}$ can be time-dependent, it is useful to also define the
time-averaged quantity

$\displaystyle\left<\mathrm{IPR}_{\mu}\right>_{t}\approx
1-2\left\langle\psi_{b,\mu}|\hat{S}^{\dagger}\hat{S}|\psi_{b,\mu}\right\rangle_{t}=1-2\lim_{t\rightarrow\infty}\frac{1}{t}\int
dt\
\left\langle\psi_{b,\mu}|\hat{S}^{\dagger}\hat{S}|\psi_{b,\mu}\right\rangle.$
(86)

We will now identify the generator of the SW transformation for the different
cases considered in this manuscript.

####  B.2.1 Static $\hat{H}_{g}^{P}$ and $\hat{H}^{P}-\hat{H}_{\kappa}^{P}$

The generator of the SW transformation in this case is

$\displaystyle\hat{S}=\sum_{j,n,m,k,\ell}\frac{\sqrt{(n+1)(m+1)}\left\langle
k|\hat{D}_{n,r}|\ell\right\rangle
g_{j}e^{i\phi_{n}}/\tilde{\Delta}_{j}}{1+(mK_{j}-n\tilde{K}_{b})/\tilde{\Delta}_{j}+q_{j,n,k}-q_{j,n+1,\ell}}\left|m+1\right\rangle\left\langle
m\right|_{j}\left|n\right\rangle\left\langle
n+1\right|_{b}\left|k\right\rangle\left\langle\ell\right|_{r}-\mathrm{h.c.},$
(87) $\displaystyle
q_{j,n,k}=\frac{\delta+n\chi}{\tilde{\Delta}_{j}}k+\frac{K_{r}}{\tilde{\Delta}_{j}}\frac{k(k-1)}{2}.$
(88)

Given that $|\langle k|\hat{D}_{n,r}|\ell\rangle|$ is more strongly suppressed
for small Fock state numbers $\left\\{k,\ \ell\right\\}$, it is clear that the
parameters $(\delta+n\chi)/\tilde{\Delta}_{j}$ and $K_{r}/\tilde{\Delta}_{j}$,
which control the probability of virtually populating larger Fock states of
the NLR, play an important role in the efficiency of the suppression of two-
qubit interactions. More precisely, we find that

$\displaystyle\left\langle\psi_{1000}|\hat{S}^{\dagger}\hat{S}|\psi_{1000}\right\rangle=e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{j}}{\tilde{\Delta}_{j}}\right)^{2}\sum_{\ell=0}^{\infty}\frac{\left|\bar{\alpha}\right|^{2\ell}/\ell!}{\left(1-q_{j,1,\ell}\right)^{2}}=e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{1}}{\tilde{\Delta}_{1}}\right)^{2}{}_{4}F_{4}\left(\bm{p}_{1};1+\bm{p}_{1};\left|\bar{\alpha}\right|^{2}\right).$
(89)

Throughout this work, we consider the two limiting cases including $K_{r}=0$
and $K_{r}\rightarrow\infty$. In these limits, it is possible to derive
asymptotic expressions for the inverse participation ratio, as

$\displaystyle\mathrm{IPR}_{1000}^{\left|K_{r}\right|\rightarrow
0}=1-2\frac{g_{1}^{2}}{\tilde{\Delta}_{1}^{2}}e^{-\left|\tilde{\alpha}\right|^{2}}{}_{2}F_{2}\left(\bm{q}_{1};1+\bm{q}_{1};\left|\tilde{\alpha}\right|^{2}\right),$
(90)
$\displaystyle\mathrm{IPR}_{1000}^{\left|K_{r}\right|\rightarrow\infty}=1-2\frac{g_{1}^{2}}{\tilde{\Delta}_{1}^{2}}e^{-\left|\tilde{\alpha}\right|^{2}}\left(1+\frac{\left|\tilde{\alpha}\right|^{2}}{\left(1+\zeta_{j}\right)^{2}}\right),$
(91)

where $\bm{q}_{j}=\begin{pmatrix}1/\zeta_{j}&1/\zeta_{j}\end{pmatrix}$ and
$\zeta_{j}=-\left(\delta+\chi\right)/\tilde{\Delta}_{j}$. This can be easily
generalized to the states $0100$ and $1100$.

Figure 9: Suppression of $1-IPR_{0100}$ with respect to the NLR photon number
$\left|\bar{\alpha}_{0}\right|^{2}$ in the stabilized ground state. Each data
point is computed from numerical diagonalization with a fixed drive amplitude
in the NLR. Black lines correspond to the same system but with no drive on the
NLR and the bus frequency tuned to the ac-Stark shifted frequency found in the
driven system. The gray regions are bounded by the analytical estimates in Eq.
29 for $K_{r}\rightarrow\infty$ and $K_{r}=0$. We also note the presence of a
resonance in b) for $K_{r}=0$ only for a small photon number for this specific
choice of system parameters. This results from frequency collisions with
higher Fock states in the NLR which can be otherwise prevented by the addition
of a NLR anharmonicity. Here $(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz,
$(\omega_{2}-\omega_{b})/2\pi=14.0$ MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz,
$\chi/2\pi=-20.0$ MHz, and $g/2\pi=2.0$ MHz.

####  B.2.2 Time-dependent $\hat{H}_{g}^{P}$ and static
$\hat{H}^{P}-\hat{H}_{\kappa}^{P}$

We now consider a time-dependent phase

$\displaystyle
e^{i\phi_{n}}=\sum_{s=-\infty}^{+\infty}\zeta_{n,s}e^{is\omega_{m}t},$ (92)

where $\zeta_{s}$ and $\omega_{m}$ are free time-independent parameters. In
this case, the generator takes the form

$\displaystyle\hat{S}=\sum_{j,n,m,k,\ell,s}\frac{\sqrt{(n+1)(m+1)}\left\langle
k|\hat{D}_{n,r}|\ell\right\rangle
g_{j}\zeta_{n,s}e^{is\omega_{m}t}/\tilde{\Delta}_{j}}{1+s\omega_{m}/\tilde{\Delta}_{j}+(mK_{j}-n\tilde{K}_{b})/\tilde{\Delta}_{j}+q_{j,n,k}-q_{j,n+1,\ell}}\left|m+1\right\rangle\left\langle
m\right|_{j}\left|n\right\rangle\left\langle
n+1\right|_{b}\left|k\right\rangle\left\langle\ell\right|_{r}-\mathrm{h.c.}.$
(93)

It then follows that

$\displaystyle\left\langle\psi_{1000}|\hat{S}^{\dagger}\hat{S}|\psi_{1000}\right\rangle=e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{j}}{\tilde{\Delta}_{j}}\right)^{2}\sum_{\ell,s_{1},s_{2}}\frac{\zeta_{0,s_{1}}\zeta_{0,s_{2}}^{*}e^{i(s_{1}-s_{2})\omega_{m}t}\left|\bar{\alpha}\right|^{2\ell}/\ell!}{\left(1+s_{1}\omega_{m}/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)\left(1+s_{2}\omega_{m}/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)},$
(94)

with time average

$\displaystyle\left\langle\psi_{1000}|\hat{S}^{\dagger}\hat{S}|\psi_{1000}\right\rangle_{t}=e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{j}}{\tilde{\Delta}_{j}}\right)^{2}\sum_{\ell,s}\frac{\left|\zeta_{0,s}\right|^{2}\left|\bar{\alpha}\right|^{2\ell}/\ell!}{\left(1+s\omega_{m}/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)^{2}}.$
(95)

####  B.2.3 Static $\hat{H}_{g}^{P}$ and time-dependent
$\hat{H}^{P}-\hat{H}_{\kappa}^{P}$

Here we consider a time-dependent
$\delta\rightarrow\delta-z\omega_{m}\sin(\omega_{m}t)$. The generator takes
the form

$\displaystyle\hat{S}=\sum_{j,n,m,k,\ell,s}\frac{\sqrt{(n+1)(m+1)}\left\langle
k|\hat{D}_{n,r}|\ell\right\rangle
g_{j}e^{i\phi_{n}}\zeta_{k,\ell,s}e^{is\omega_{m}t-i(k-\ell)z\cos(\omega_{m}t)}/\tilde{\Delta}_{j}}{1+s\omega_{m}/\tilde{\Delta}_{j}+(mK_{j}-n\tilde{K}_{b})/\tilde{\Delta}_{j}+q_{j,n,k}-q_{j,n+1,\ell}}$
$\displaystyle\cdot\left|m+1\right\rangle\left\langle
m\right|_{j}\left|n\right\rangle\left\langle
n+1\right|_{b}\left|k\right\rangle\left\langle\ell\right|_{r}-\mathrm{h.c.}$ ,
(96)

where the parameters $\zeta_{k,\ell,s}$ are defined with

$\displaystyle
e^{i(k-\ell)z\cos(\omega_{m}t)}=\sum_{s=-\infty}^{+\infty}i^{s}J_{s}\left[(k-\ell)z\right]e^{is\omega_{m}t}=\sum_{s=-\infty}^{+\infty}\zeta_{k,\ell,s}e^{is\omega_{m}t},$
(97)

where we have used a Jacobi-Anger expansion. We thus arrive at

$\displaystyle\left\langle\psi_{1000}|\hat{S}^{\dagger}\hat{S}|\psi_{1000}\right\rangle=e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{j}}{\tilde{\Delta}_{j}}\right)^{2}\sum_{\ell,s_{1},s_{2}}\frac{(-1)^{s_{1}+s_{2}}J_{s_{1}}(\ell
z)J_{s_{2}}(\ell
z)i^{s_{1}-s_{2}}e^{i(s_{1}-s_{2})\omega_{m}t}\left|\bar{\alpha}\right|^{2\ell}/\ell!}{\left(1+s_{1}\omega_{m}/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)\left(1+s_{2}\omega_{m}/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)},$
(98)

and the time-averaged version

$\displaystyle\left\langle\psi_{1000}|\hat{S}^{\dagger}\hat{S}|\psi_{1000}\right\rangle_{t}=e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{j}}{\tilde{\Delta}_{j}}\right)^{2}\sum_{\ell,s}\frac{J_{s}^{2}(\ell
z)\left|\bar{\alpha}\right|^{2\ell}/\ell!}{\left(1+s\omega_{m}/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)^{2}}.$
(99)

####  B.2.4 Time-dependent $\hat{H}_{g}^{P}$ and
$\hat{H}^{P}-\hat{H}_{\kappa}^{P}$

Finally, we combine the two previous cases, namely we consider a time-
dependent phase

$\displaystyle
e^{i\phi_{n}}=\sum_{s=-\infty}^{+\infty}\zeta_{n,s}^{\phi}e^{is\omega_{m}^{\phi}t},$
(100)

and a time-dependent
$\delta\rightarrow\delta-z\omega_{m}^{\delta}\sin(\omega_{m}^{\delta}t)$ with

$\displaystyle
e^{i(k-\ell)z\cos(\omega_{m}^{\delta}t)}=\sum_{s=-\infty}^{+\infty}i^{s}J_{s}\left[(k-\ell)z\right]e^{is\omega_{m}^{\delta}t}=\sum_{s=-\infty}^{+\infty}\zeta_{k,\ell,s}^{\delta}e^{is\omega_{m}^{\delta}t}.$
(101)

The generator takes the form

$\displaystyle\hat{S}=\sum_{j,n,m,k,\ell,s,r}\frac{\sqrt{(n+1)(m+1)}\left\langle
k|\hat{D}_{n,r}|\ell\right\rangle
g_{j}\zeta_{n,s}^{\phi}\zeta_{k,\ell,r}^{\delta}e^{i\left(s\omega_{m}^{\phi}+r\omega_{m}^{\delta}\right)t-i(k-\ell)z\cos(\omega_{m}^{\delta}t)}/\tilde{\Delta}_{j}}{1+\left(s\omega_{m}^{\phi}+r\omega_{m}^{\delta}\right)/\tilde{\Delta}_{j}+(mK_{j}-n\tilde{K}_{b})/\tilde{\Delta}_{j}+q_{j,n,k}-q_{j,n+1,\ell}}$
$\displaystyle\cdot\left|m+1\right\rangle\left\langle
m\right|_{j}\left|n\right\rangle\left\langle
n+1\right|_{b}\left|k\right\rangle\left\langle\ell\right|_{r}-\mathrm{h.c.}$ .
(102)

With this we find that

$\displaystyle\left\langle\psi_{1000}|\hat{S}^{\dagger}\hat{S}|\psi_{1000}\right\rangle=$
$\displaystyle
e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{j}}{\tilde{\Delta}_{j}}\right)^{2}\sum_{\ell,s_{1},s_{2},r_{1},r_{2}}i^{s_{1}-s_{2}+r_{1}-r_{2}}e^{i(s_{1}-s_{2})\omega_{m}^{\phi}t+i(r_{1}-r_{2})\omega_{m}^{\delta}t}$
$\displaystyle\cdot\frac{(-1)^{r_{1}+r_{2}}\zeta_{0,s_{1}}\zeta_{0,s_{2}}^{*}J_{r_{1}}(\ell
z)J_{r_{2}}(\ell
z)\left|\bar{\alpha}\right|^{2\ell}/\ell!}{\left(1+\left(s_{1}\omega_{m}^{\phi}+r_{1}\omega_{m}^{\delta}\right)/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)\left(1+\left(s_{2}\omega_{m}^{\phi}+r_{2}\omega_{m}^{\delta}\right)/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)},$
(103)

which under time-averaging reduces to

$\displaystyle\left\langle\psi_{1000}|\hat{S}^{\dagger}\hat{S}|\psi_{1000}\right\rangle_{t}=e^{-\left|\bar{\alpha}\right|^{2}}\left(\frac{g_{j}}{\tilde{\Delta}_{j}}\right)^{2}\sum_{\ell,s,r}\frac{\left|\zeta_{0,s}\right|^{2}J_{r}^{2}(\ell
z)\left|\bar{\alpha}\right|^{2\ell}/\ell!}{\left(1+\left(s\omega_{m}^{\phi}+r\omega_{m}^{\delta}\right)/\tilde{\Delta}_{j}-q_{j,1,\ell}\right)^{2}}.$
(104)

Here we assumed $s\omega_{m}^{\phi}+r\omega_{m}^{\delta}=0$ only for $s=0$ and
$r=0$.

###  B.3 Two-qubit interactions

In this section we demonstrate how the matrix elements of the displacement
operator in $\hat{H}_{g}^{P}$ [c.f. Eq. 17] yields exponentially suppressed
two-qubit interactions by deriving an upper bound based on the inverse
participation ratio. Consider the time-evolution operator in the interaction
picture

$\displaystyle\hat{U}_{I}(t)=\mathcal{T}e^{-i\int_{0}^{t}d\tau\hat{H}_{I}(\tau)}=\sum_{n=0}^{\infty}\frac{(-i)^{n}}{n!}\int_{0}^{t}d\tau_{1}\cdots\int_{0}^{t}d\tau_{n}\mathcal{T}\hat{H}_{I}(\tau_{1})\cdots\hat{H}_{I}(\tau_{n}),$
(105)

where $\mathcal{T}$ is the time-ordering operator and we defined the
interaction picture Hamiltonian in Eq. 82. It is convenient to approximate Eq.
105 using a Magnus expansion,

$\displaystyle\hat{U}_{I}=e^{-i\int_{0}^{t}d\tau\hat{H}^{M}_{I}(\tau)},$ (106)

where $\hat{H}_{I}^{M}$ is an effective Hamiltonian. Up to fourth order in
$g_{j}$ we have that

$\displaystyle\hat{H}_{I}^{M}(t)\approx$
$\displaystyle\hat{H}_{I}(t)-\frac{1}{2}\left[\hat{H}_{I}(t),\hat{S}_{I}(t)\right]+\frac{i}{6}\int_{0}^{t}d\tau_{1}\left(\left[\hat{H}_{I}(t),\left[\hat{H}_{I}(\tau_{1}),\hat{S}_{I}(\tau_{1})\right]\right]+\left[\hat{S}_{I}(\tau_{1}),\left[\hat{H}_{I}(\tau_{1}),\hat{H}_{I}(t)\right]\right]\right)$
$\displaystyle+\frac{1}{12}\int_{0}^{t}\int_{0}^{\tau_{1}}d\tau_{1}d\tau_{2}\left(\left[\left[\left[\hat{H}_{I}(t),\hat{H}_{I}(\tau_{1})\right],\hat{H}_{I}(\tau_{2})\right],\hat{S}_{I}(\tau_{2})\right]+\left[\hat{H}_{I}(t),\left[\left[\hat{H}_{I}(\tau_{1}),\hat{H}_{I}(\tau_{2})\right],\hat{S}_{I}(\tau_{2})\right]\right]\right)$
$\displaystyle+\frac{1}{12}\int_{0}^{t}\int_{0}^{\tau_{1}}d\tau_{1}d\tau_{2}\left(\left[\hat{H}_{I}(t),\left[\hat{H}_{I}(\tau_{1}),\left[\hat{H}_{I}(\tau_{2}),\hat{S}_{I}(\tau_{2})\right]\right]\right]+\left[\hat{H}_{I}(\tau_{1}),\left[\hat{H}_{I}(\tau_{2}),\left[\hat{S}_{I}(\tau_{2}),\hat{H}_{I}(t)\right]\right]\right]\right),$
(107)

where based on Eq. 83 we introduced
$\hat{S}_{I}(t)=i\int_{0}^{t}d\tau\hat{H}_{I}(\tau)$.

To compute effective two-qubit interactions, we project $\hat{H}_{I}^{M}$ in
the ground state of the coupler, valid in the dispersive coupling limit. Due
to the form of $\hat{H}_{I}$ [c.f. Eq. 82], only terms of even order in
$g_{j}$ remain in $\hat{H}_{I}^{M}$. Defining the projection operator
$\hat{P}=\left|0_{b},0_{r}\right\rangle\left\langle 0_{b},0_{r}\right|$, we
compute the effective two-qubit Hamiltonian

$\displaystyle\hat{H}_{I}^{q-q}(t)=Tr_{c}\left(\hat{H}_{I}^{M}(t)\hat{P}\right),$
(108)

corresponding to $\hat{H}_{I}^{M}$, where the coupler mode is traced out
assuming that it is stabilized in its ground state. For compactness, we define
$\hat{H}_{I}^{M}=\sum_{n=1}^{4}\hat{H}_{I,n}^{M}$ where $\hat{H}_{I,n}^{M}$
incorporates the couplings $g_{j}$ at $n$th order. In what follows we shall
approximate the qubits and coupler as two-level systems. This approximation is
valid in the limit where the qubit-coupler detunings are small compared to the
anharmonicities of the qubits and the coupler. Our results can, however, be
extended to include the effects of finite anharmonicity. We therefore write

$\displaystyle\hat{H}_{I}(t)\approx\sum_{j,k,\ell}g_{j}e^{i\tilde{\Delta}_{j}\left(1+q_{j,0,k}-q_{j,1,\ell}\right)t}\hat{A}_{k,\ell}\hat{\sigma}_{+,j}+\text{h.c.},$
(109) $\displaystyle\hat{A}_{k,\ell}=e^{i\phi_{0}}\left\langle
k|\hat{D}_{0,r}|\ell\right\rangle\left|k\right\rangle\left\langle\ell\right|_{r}\hat{\sigma}_{-,b},$
(110) $\displaystyle
q_{j,n,k}=\frac{\delta+n\chi}{\tilde{\Delta}_{j}}k+\frac{K_{r}}{\tilde{\Delta}_{j}}\frac{k(k-1)}{2}=\frac{r_{n,k}}{\tilde{\Delta}_{j}},$
(111)

where $\hat{\sigma}_{-,j}$ ($\hat{\sigma}_{-,b}$) are the spin ladder
operators, the two-level approximation of $\hat{q}_{j}$ ($\hat{b}$). The
interaction picture generator $\hat{S}_{I}$ and the polaron frame generator
$\hat{S}$ then take the form

$\displaystyle\hat{S}_{I}(t)\approx\sum_{j,k,\ell}\frac{g_{j}}{\tilde{\Delta}_{j}}e^{i\tilde{\Delta}_{j}\left(1+q_{j,0,k}-q_{j,1,\ell}\right)t}\frac{\hat{A}_{k,\ell}\hat{\sigma}_{+,j}}{1+q_{j,0,k}-q_{j,1,\ell}}-\text{h.c.},$
(112)
$\displaystyle\hat{S}\approx\sum_{j,k,\ell}\frac{g_{j}}{\tilde{\Delta}_{j}}\frac{\hat{A}_{k,\ell}\hat{\sigma}_{+,j}}{1+q_{j,0,k}-q_{j,1,\ell}}-\text{h.c.},$
(113)

in analogy to Eq. 83 and Eq. 84. Our goal is to find an upper bound on the
amplitude of two-qubit interactions in Eq. 108 corresponding to a partial
trace over Eq. 107. This upper bound can be obtained from the Cauchy-Schwartz
inequality,

$\displaystyle\left|Tr\left(\hat{A}\hat{B}\right)\right|\leq\sqrt{Tr\left(\hat{A}^{\dagger}\hat{A}\right)}\sqrt{Tr\left(\hat{B}^{\dagger}\hat{B}\right)}=\left|\left|\hat{A}\right|\right|_{F}\left|\left|\hat{B}\right|\right|_{F},$
(114)

where $\left|\left|\bullet\right|\right|_{F}$ stands for the Frobenius norm.
We will separate the two-qubit interactions by the order in the coupling
strengths $g_{j}$ and stop at fourth order. However, the following analysis
can be extended to higher orders in the coupling strengths.

####  B.3.1 Second-order interactions

We now focus on the two-qubit interactions of second order in the coupling
strengths in Eq. 108. To this end, it is convenient to expand

$\displaystyle
Tr\left(\hat{H}_{I,2}^{M}(t)\hat{P}\right)=-\frac{1}{2}Tr_{c}\left(\hat{H}_{I}(t)\hat{S}_{I}(t)\hat{P}\right)+\text{h.c.},$
(115)

where

$\displaystyle
Tr_{c}\left(\hat{H}_{I}(t)\hat{S}_{I}(t)\hat{P}\right)=\sum_{s_{1},s_{2}=1}^{2}e^{i\left(\tilde{\Delta}_{q_{s_{1}}}-\tilde{\Delta}_{q_{s_{2}}}\right)t}g_{s_{1}}Tr\left(\sum_{\ell=0}^{\infty}\hat{A}_{0,\ell}\frac{g_{s_{2}}}{\tilde{\Delta}_{q_{s_{2}}}}\frac{\hat{A}_{0,\ell}^{\dagger}}{1-q_{s_{2},1,\ell}}\hat{P}\right)\hat{\sigma}_{+,s_{1}}\hat{\sigma}_{-,s_{2}}.$
(116)

Interestingly, applying Eq. 114 yields

$\displaystyle\left|Tr\left(\sum_{\ell=0}^{\infty}\hat{A}_{0,\ell}\frac{g_{s_{2}}}{\tilde{\Delta}_{q_{s_{2}}}}\frac{\hat{A}_{0,\ell}^{\dagger}}{1-q_{s_{2},1,\ell}}\hat{P}\right)\right|\leq\left|\left|\sum_{\ell=0}^{\infty}\hat{A}_{0,\ell}\right|\right|_{F}\left|\left|\sum_{\ell=0}^{\infty}\frac{g_{s_{2}}}{\tilde{\Delta}_{q_{s_{2}}}}\frac{\hat{A}_{0,\ell}^{\dagger}}{1-q_{s_{2},1,\ell}}\hat{P}\right|\right|_{F}$
. (117)

where

$\displaystyle\left|\left|\sum_{\ell=0}^{\infty}\hat{A}_{0,\ell}\right|\right|_{F}=\left|\left|\left|0\right\rangle\left\langle
0\right|_{r}e^{i\phi_{0}}\hat{D}_{0,r}\right|\right|_{F}\cdot\left|\left|\hat{\sigma}_{-,b}\right|\right|_{F}=1.$
(118)

Furthermore,

$\displaystyle\left|\left|\sum_{\ell=0}^{\infty}\frac{g_{s_{2}}}{\tilde{\Delta}_{q_{s_{2}}}}\frac{\hat{A}_{0,\ell}^{\dagger}}{1-q_{s_{2},1,\ell}}\hat{P}\right|\right|_{F}$
$\displaystyle=\sqrt{\left(\frac{g_{s_{2}}}{\tilde{\Delta}_{q_{s_{2}}}}\right)^{2}\sum_{\ell=0}^{\infty}\frac{\left|\left\langle\ell|\hat{D}_{0,r}|0\right\rangle\right|^{2}}{\left(1-q_{s_{2},1,\ell}\right)^{2}}}\approx\sqrt{\frac{1-\mathrm{IPR}_{\Omega(s_{2})}}{2}},$
(119)

with $\Omega(s_{2})=1000$ for $s_{2}=1$ and $0100$ for $s_{2}=2$, and where we
used the dispersive limit expression for $\mathrm{IPR}$ [c.f. Eq. 89]. Since
each term in Eq. 115 is bounded in magnitude by
$\text{max}_{\nu\in\\{1000,0100\\}}\sqrt{1-\mathrm{IPR}_{\nu}}$ it follows
that Eq. 115 is equally bounded in magnitude by
$\text{max}_{\nu\in\\{1000,0100\\}}\sqrt{1-\mathrm{IPR}_{\nu}}$. Intuitively,
this means that the virtual interaction is always smaller in magnitude than
$\left|g_{i}g_{j}/\tilde{\Delta}_{j}\right|$.

####  B.3.2 Fourth-order interactions

We will now demonstrate that the
$\text{max}_{\nu\in\\{1000,0100\\}}\sqrt{1-\mathrm{IPR}_{\nu}}$ bound holds
true for higher-order two-qubit interactions. First, we expand the commutators
such that

$\displaystyle Tr_{c}\left(\hat{H}_{I,4}^{M}\hat{P}\right)$
$\displaystyle=\frac{1}{12}\int_{0}^{t}\int_{0}^{\tau_{1}}d\tau_{1}d\tau_{2}\left[2Tr_{c}\left(\hat{H}_{I}(t)\hat{H}_{I}(\tau_{1})\hat{H}_{I}(\tau_{2})\hat{S}_{I}(\tau_{2})\hat{P}\right)\right.$
$\displaystyle-
Tr_{c}\left(\left[\hat{H}_{I}(\tau_{1}),\hat{H}_{I}(t)\right]\hat{H}_{I}(\tau_{2})\hat{S}_{I}(\tau_{2})\hat{P}\right)-Tr_{c}\left(\left[\hat{H}_{I}(\tau_{2}),\hat{H}_{I}(t)\right]\hat{H}_{I}(\tau_{1})\hat{S}_{I}(\tau_{2})\hat{P}\right)$
$\displaystyle-
Tr_{c}\left(\left[\hat{H}_{I}(\tau_{1}),\hat{H}_{I}(\tau_{2})\right]\hat{H}_{I}(t)\hat{S}_{I}(\tau_{2})\hat{P}\right)-Tr_{c}\left(\left[\hat{H}_{I}(t),\hat{H}_{I}(\tau_{1})\right]\hat{S}_{I}(\tau_{2})\hat{H}_{I}(\tau_{2})\hat{P}\right)$
$\displaystyle-
Tr_{c}\left(\left[\hat{H}_{I}(t),\hat{H}_{I}(\tau_{2})\right]\hat{S}_{I}(\tau_{2})\hat{H}_{I}(\tau_{1})\hat{P}\right)\left.-Tr_{c}\left(\left[\hat{H}_{I}(\tau_{1}),\hat{H}_{I}(\tau_{2})\right]\hat{S}_{I}(\tau_{2})\hat{H}_{I}(t)\hat{P}\right)\right]+\text{h.c.}.$
(120)

Second, we find that

$\displaystyle
Tr_{c}\left(\hat{H}_{I}(t_{1})\hat{H}_{I}(t_{2})\hat{H}_{I}(t_{3})\hat{S}_{I}(t_{4})\hat{P}\right)=\sum_{s_{1},\cdots,s_{4}=1}^{2}e^{i\tilde{\Delta}_{q_{s_{1}}}t_{1}}e^{-i\tilde{\Delta}_{q_{s_{2}}}t_{2}}e^{i\tilde{\Delta}_{q_{s_{3}}}t_{3}}e^{-i\tilde{\Delta}_{q_{s_{4}}}t_{4}}g_{s_{1}}g_{s_{2}}g_{s_{3}}$
$\displaystyle\cdot
Tr\left(\sum_{k,\ell,m=0}^{\infty}e^{-ir_{1,k}(t_{1}-t_{2})-ir_{0,\ell}(t_{2}-t_{3})-ir_{1,m}(t_{3}-t_{4})}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\hat{A}_{\ell,m}\frac{g_{s_{4}}}{\tilde{\Delta}_{q_{s_{4}}}}\frac{\hat{A}_{0,m}^{\dagger}}{1-q_{s_{4},1,m}}\hat{P}\right)\hat{\sigma}_{+,s_{1}}\hat{\sigma}_{-,s_{2}}\hat{\sigma}_{+,s_{3}}\hat{\sigma}_{-,s_{4}},$
(121)

and similarly,

$\displaystyle
Tr_{c}\left(\hat{H}_{I}(t_{1})\hat{H}_{I}(t_{2})\hat{S}_{I}(t_{3})\hat{H}_{I}(t_{4})\hat{P}\right)=-\sum_{s_{1},\cdots,s_{4}=1}^{2}e^{i\tilde{\Delta}_{q_{s_{1}}}t_{1}}e^{-i\tilde{\Delta}_{q_{s_{2}}}t_{2}}e^{i\tilde{\Delta}_{q_{s_{3}}}t_{3}}e^{-i\tilde{\Delta}_{q_{s_{4}}}t_{4}}g_{s_{1}}g_{s_{2}}g_{s_{4}}$
$\displaystyle\cdot
Tr\left(\sum_{k,\ell,m=0}^{\infty}e^{-ir_{1,k}(t_{1}-t_{2})-ir_{0,\ell}(t_{2}-t_{3})-ir_{1,m}(t_{3}-t_{4})}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\frac{g_{s_{3}}}{\tilde{\Delta}_{q_{s_{3}}}}\frac{\hat{A}_{\ell,m}\hat{A}_{0,m}^{\dagger}}{1+q_{s_{3},0,\ell}-q_{s_{3},1,m}}\hat{P}\right)\hat{\sigma}_{+,s_{1}}\hat{\sigma}_{-,s_{2}}\hat{\sigma}_{+,s_{3}}\hat{\sigma}_{-,s_{4}}.$
(122)

As it was done in the previous section we apply Eq. 114 to find an upper bound
on Eq. 121

$\displaystyle\left|Tr\left(\sum_{k,\ell,m=0}^{\infty}e^{-ir_{1,k}(t_{1}-t_{2})-ir_{0,\ell}(t_{2}-t_{3})-ir_{1,m}(t_{3}-t_{4})}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\hat{A}_{\ell,m}\frac{g_{s_{4}}}{\tilde{\Delta}_{q_{s_{4}}}}\frac{\hat{A}_{0,m}^{\dagger}}{1-q_{s_{4},1,m}}\hat{P}\right)\right|$
$\displaystyle\leq\left|\left|\sum_{k,\ell,m=0}^{\infty}e^{-ir_{1,k}(t_{1}-t_{2})-ir_{0,\ell}(t_{2}-t_{3})-ir_{1,m}(t_{3}-t_{4})}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\hat{A}_{\ell,m}\right|\right|_{F}\left|\left|\sum_{m=0}^{\infty}\frac{g_{s_{4}}}{\tilde{\Delta}_{q_{s_{4}}}}\frac{\hat{A}_{0,m}^{\dagger}}{1-q_{s_{4},1,m}}\hat{P}\right|\right|_{F}$
$\displaystyle\leq\left|\left|\sum_{k,\ell,m=0}^{\infty}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\hat{A}_{\ell,m}\right|\right|_{F}\left|\left|\sum_{m=0}^{\infty}\frac{g_{s_{4}}}{\tilde{\Delta}_{q_{s_{4}}}}\frac{\hat{A}_{0,m}^{\dagger}}{1-q_{s_{4},1,m}}\hat{P}\right|\right|_{F},$
(123)

where it can also be verified that
$\left|\left|\sum_{k,\ell,m=0}^{\infty}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\hat{A}_{\ell,m}\right|\right|_{F}=1$
and where

$\displaystyle\left|\left|\sum_{m=0}^{\infty}\frac{g_{s_{4}}}{\tilde{\Delta}_{q_{s_{4}}}}\frac{\hat{A}_{0,m}^{\dagger}}{1-q_{s_{4},1,m}}\hat{P}\right|\right|_{F}\approx\sqrt{\frac{1-\mathrm{IPR}_{\Omega(s_{4}))}}{2}}.$
(124)

Similarly, Eq. 114 can be applied on Eq. 122

$\displaystyle\left|Tr\left(\sum_{k,\ell,m=0}^{\infty}e^{-ir_{1,k}(t_{1}-t_{2})-ir_{0,\ell}(t_{2}-t_{3})-ir_{1,m}(t_{3}-t_{4})}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\frac{g_{s_{3}}}{\tilde{\Delta}_{q_{s_{3}}}}\frac{\hat{A}_{\ell,m}\hat{A}_{0,m}^{\dagger}}{1+q_{s_{3},0,\ell}-q_{s_{3},1,m}}\hat{P}\right)\right|$
$\displaystyle\leq\left|\left|\sum_{k,\ell=0}^{\infty}e^{-ir_{1,k}(t_{1}-t_{2})-ir_{0,\ell}(t_{2}-t_{3})}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\right|\right|_{F}\left|\left|\sum_{\ell,m=0}^{\infty}e^{-ir_{1,m}(t_{3}-t_{4})}\frac{g_{s_{3}}}{\tilde{\Delta}_{q_{s_{3}}}}\frac{\hat{A}_{\ell,m}\hat{A}_{0,m}^{\dagger}}{1+q_{s_{3},0,\ell}-q_{s_{3},1,m}}\hat{P}\right|\right|_{F}$
$\displaystyle\leq\left|\left|\sum_{k,\ell=0}^{\infty}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\right|\right|_{F}\left|\left|\sum_{\ell,m=0}^{\infty}\frac{g_{s_{3}}}{\tilde{\Delta}_{q_{s_{3}}}}\frac{\hat{A}_{\ell,m}\hat{A}_{0,m}^{\dagger}}{1+q_{s_{3},0,\ell}-q_{s_{3},1,m}}\hat{P}\right|\right|_{F},$
(125)

where
$\left|\left|\sum_{k,\ell=0}^{\infty}\hat{A}_{0,k}\hat{A}_{\ell,k}^{\dagger}\right|\right|_{F}=1$
and where

$\displaystyle\left|\left|\sum_{\ell,m=0}^{\infty}\frac{g_{s_{3}}}{\tilde{\Delta}_{q_{s_{3}}}}\frac{\hat{A}_{\ell,m}\hat{A}_{0,m}^{\dagger}}{1+q_{s_{3},0,\ell}-q_{s_{3},1,m}}\hat{P}\right|\right|_{F}=\sqrt{\sum_{\ell,m,m^{\prime}=0}^{\infty}\left(\frac{g_{s_{3}}}{\tilde{\Delta}_{q_{s_{3}}}}\right)^{2}\frac{\left\langle
0|\hat{D}_{0,r}|m^{\prime}\right\rangle\left\langle
m^{\prime}|\hat{D}_{0,r}^{\dagger}|\ell\right\rangle\left\langle\ell|\hat{D}_{0,r}|m\right\rangle\left\langle
m|\hat{D}_{0,r}^{\dagger}|0\right\rangle}{\left(1+q_{s_{3},0,\ell}-q_{s_{3},1,m^{\prime}}\right)\left(1+q_{s_{3},0,\ell}-q_{s_{3},1,m}\right)}}.$
(126)

Considering this norm to be bounded from above by the limiting case $\delta=0$
and $K_{r}=0$, we find that

$\displaystyle\left|\left|\sum_{\ell,m=0}^{\infty}\frac{g_{s_{3}}}{\tilde{\Delta}_{q_{s_{3}}}}\frac{\hat{A}_{\ell,m}\hat{A}_{0,m}^{\dagger}}{1+q_{s_{3},0,\ell}-q_{s_{3},1,m}}\hat{P}\right|\right|_{F}$
$\displaystyle\lessapprox\sqrt{\sum_{\ell,m=0}^{\infty}\left(\frac{g_{s_{3}}}{\tilde{\Delta}_{q_{s_{3}}}}\right)^{2}\frac{\left|\left\langle
m|\hat{D}_{0,r}|0\right\rangle\right|^{2}}{\left(1-q_{s_{3},1,m}\right)^{2}}}\approx\sqrt{\frac{1-\mathrm{IPR}_{\Omega(s_{3})}}{2}}.$
(127)

As previously found for second-order interactions, all contributions are
bounded in magnitude by
$\text{max}_{\nu\in\\{1000,0100\\}}\sqrt{1-\mathrm{IPR}_{\nu}}$. These
observations can be straightforwardly extended to higher orders in a similar
fashion.

This implies that exponentially suppressing $1-\mathrm{IPR}_{\nu}$ will
equally exponentially suppress virtual two-qubit interactions. However, we
stress that this only provides an estimate for the order of magnitude.

###  B.4 Measurement-induced dephasing

In order to quantify the dephasing induced by the coupler drive we express Eq.
26 in the hybridized eigenbasis

$\displaystyle\hat{L}_{b}^{P}=\sqrt{\frac{\kappa\left|\bar{\alpha}\right|^{2}}{4}}\sum_{\nu,\nu^{\prime}}c_{\nu;\nu^{\prime}}\left|\psi_{h,\nu}\right\rangle\left\langle\psi_{h,\nu^{\prime}}\right|,$
(128)

with

$c_{\nu;\nu^{\prime}}=\left\langle\psi_{h,\nu}|(\left|1\right\rangle\left\langle
1\right|_{b}-\left|0\right\rangle\left\langle
0\right|_{b})|\psi_{h,\nu^{\prime}}\right\rangle.$ (129)

To study the effects of bus dephasing on $Q_{1}$ alone, we first trace out
$Q_{2}$, $B$ and $R$ in the hybridized basis where all three have zero
excitation

$\displaystyle\hat{L}_{1}^{P}\approx\sum_{i,j}\left\langle\psi_{h,i000}|\hat{L}_{b}^{P}|\psi_{h,j000}\right\rangle\left|\psi_{h,i000}\right\rangle\left\langle\psi_{h,j000}\right|.$
(130)

We then apply a rotating-wave approximation to obtain

$\displaystyle\hat{L}_{1}^{P}\approx\sqrt{\frac{\kappa\left|\bar{\alpha}\right|^{2}}{4}}\sum_{i}c_{i000;i000}\left|\psi_{h,i000}\right\rangle\left\langle\psi_{h,i000}\right|.$
(131)

Because the total excitation number in $Q_{1}$, $Q_{2}$ and $B$ is conserved
under Eq. 17 and the interaction in Eq. 17 dominantly yields hybridization
between $1$ and $B$, we neglect
$\left\langle\psi_{h,1000}|\psi_{b,010k}\right\rangle\approx 0$. We also use
the fact that
$\sum_{k}\left|\left\langle\psi_{h,0000}|\psi_{b,000k}\right\rangle\right|^{2}=1$
due to total excitation number conservation, and finally,
$\sum_{k}\left|\left\langle\psi_{h,1000}|\psi_{b,100k}\right\rangle\right|^{2}\approx
1-\sum_{k}\left|\left\langle\psi_{h,1000}|\psi_{b,001k}\right\rangle\right|^{2}$.
It follows that $c_{1000;1000}\approx
2\sum_{k}\left|\left\langle\psi_{h,1000}|\psi_{b,001k}\right\rangle\right|^{2}-1$
and $c_{0000;0000}\approx-1$, leading to

$\displaystyle\hat{L}_{1}^{P}\approx$
$\displaystyle\sqrt{\frac{\kappa\left|\bar{\alpha}\right|^{2}}{4}}\sum_{k=0}^{\infty}\left|\left\langle\psi_{h,1000}|\psi_{b,001k}\right\rangle\right|^{2}\left(\left|\psi_{h,1000}\right\rangle\left\langle\psi_{h,1000}\right|-\left|\psi_{h,0000}\right\rangle\left\langle\psi_{h,0000}\right|\right).$
(132)

The dephasing rate of $Q_{1}$ can be estimated by computing the dephasing rate
associated with Eq. 132, which takes the form

$\begin{split}\gamma_{\varphi,1}&=\frac{\kappa\left|\bar{\alpha}\right|^{2}}{2}\left(\sum_{k=0}^{\infty}\left|\left\langle\psi_{h,1000}|\psi_{b,001k}\right\rangle\right|^{2}\right)^{2}\approx\frac{\kappa\left|\bar{\alpha}\right|^{2}}{2}\frac{1-\mathrm{IPR}_{1000}}{2},\end{split}$
(133)

where the second line follows from the Schrieffer-Wolff transformation [c.f.
Eq. 85]. The expression above was obtained under a rotating-wave
approximation, which is valid for $|\gamma_{\varphi,1}/\tilde{\Delta}_{1}|\ll
1$. An expression for the second qubit is obtained by simply replacing the
subscript $1000$ by $0100$. Measurement-induced dephasing rates in the qubits
are estimated numerically from diagonalization and are reported in Fig. 10.

Figure 10: Measurement-induced dephasing rates in the qubits $1$ and $2$.
Dephasing rate for a)-b) qubit $1$, and c)-d) qubit $2$. The gray regions
correspond to the analytical estimates in Eq. 34 with Eqs. 90 and 91. Left
column (a), c) and e) with $\delta/2\pi=-1.5$ MHz) and right column (b), d)
and f) for $\delta/2\pi=1.0$ MHz) correspond to two different parameter
regimes set by the sign of $\delta$. Here $(\omega_{1}-\omega_{b})/2\pi=7.0$
MHz, $(\omega_{2}-\omega_{b})/2\pi=14.0$ MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$
MHz, $\chi/2\pi=-20.0$ MHz, $g/2\pi=2.0$ MHz, and $\kappa/2\pi=100.0$ kHz.

###  B.5 Parameter regimes

We have seen in Figs. 4 and 5 that the choice of frequency detuning between
the NLR and the drive $\delta$, which controls the sign of the ac-Stark shift
between the qubits and the bus, has both a quantitative and qualitative impact
on the system. To better understand this behavior, we consider again Eqs. 28
and 29 where we now express the generalized hypergeometric functions as

$\displaystyle{}_{4}F_{4}(\bm{p}_{j},1+\bm{p}_{j};\left|\bar{\alpha}\right|^{2})=\sum_{n=0}^{\infty}\frac{\left|\bar{\alpha}\right|^{2n}}{n!}\prod_{k=\pm}\left(1+\frac{n}{p_{jk}}\right)^{-2}.$
(134)

We recall that in this expression
$p_{j\pm}=\beta[1\pm\sqrt{1+\smash[b]{2\tilde{\Delta}_{j}\beta^{-2}/K_{r}}}]$
with $\beta=(\delta+\chi)/K_{r}-1/2$ and where $\tilde{\Delta}_{j}$ is the ac-
Stark shifted detuning given in Eq. 24. Physically, Eq. 134 corresponds to a
weighted sum over all virtual transitions to higher energy levels in the NLR
for a displacement $\bar{\alpha}$ during a $0\rightarrow 1$ transition in the
bus. The poles $1+n/p_{jk}$ in Eq. 134 correspond to frequency collisions with
these higher energy levels.

Given that we wish to maximize the exponential suppression of IPR [i.e. the
exponential factor in Eqs. 28 and 29], we want
${}_{4}F_{4}(\bm{p}_{j},1+\bm{p}_{j};\left|\bar{\alpha}\right|^{2})$ to be as
close as possible to unity. To have a monotonic suppression with respect to
$\left|\bar{\alpha}\right|^{2}$, this parameter should be chosen such as to
avoid frequency collisions corresponding to the poles of Eq. 134. With these
constraints, we define two key parameter regimes below. In what follows we
take $(\kappa/2)^{2}\ll\left|\delta(\delta+\chi)\right|$, but we note that
$\kappa$ also plays a role in controlling the ac-Stark shift of the bus and
can be seen as an additional knob. To simplify the analysis below, it is
useful to note that for $K_{r}=0$ we find $p_{j+}\rightarrow\infty$ and
$p_{j-}\rightarrow-\tilde{\Delta}_{j}/(\delta+\chi)$, while for
$\left|K_{r}\right|\gg\\{|\tilde{\Delta}_{j}|,\ |\delta+\chi|\\}$ we find
$p_{j+}\rightarrow-1+2(\delta+\chi-\tilde{\Delta}_{j})/K_{r}$ and
$p_{j-}\rightarrow 2\tilde{\Delta}_{j}/K_{r}$.

####  B.5.1 Monotonic suppression

We first focus on the situation illustrated in panels (a) of Figs. 4 and 5
where there are no frequency collisions between the different modes of the
system and the suppression factor is a smooth function of the photon number.
There are two possible sources of frequency collisions: i) the ac-Stark
shifted qubit-bus detuning $\tilde{\Delta}_{j}$ of Eq. 24 and ii) higher
energy levels in the NLR as captured by Eq. 134.

First, to avoid a collision where $\tilde{\Delta}_{j}=0$, the ac-Stark-shifted
qubit-bus detuning should ideally grow in magnitude with respect to
$\bar{\alpha}$. This, in turn implies that $\omega_{qj}-\omega_{b}$ and
$\delta(\delta+\chi)/\chi\approx\delta$ should have opposite signs.
Considering the second source of frequency collisions, $p_{j\pm}>0$ in Eq. 134
suppresses frequency collisions with higher energy levels in the NLR. For
$K_{r}=0$ this condition is achieved with
$p_{j-}=-\tilde{\Delta}_{j}/(\delta+\chi)>0$, meaning that
$\tilde{\Delta}_{j}$ and $\delta+\chi$ have opposite signs. Combining this
with the above finding, it follows that $\delta$ and $\chi$ have the same
sign, and opposite sign to $\omega_{qj}-\omega_{b}$. We stress that this
choice of parameters is equally compatible with the large $\left|K_{r}\right|$
limit because $\left|1/p_{j-}\right|$ is then large.

Combining the above with the fact that $\left|\delta/\chi\right|\ll 1$ for
large conditional displacements, allows us to define the parameter regime for
monotonic suppression of the IPR as

$\displaystyle\left|\delta/\chi\right|\ll
1,\quad\left(\omega_{j}-\omega_{b}\right)\delta<0,\quad\left(\omega_{j}-\omega_{b}\right)\chi<0.$
(135)

Because the nonlinearity of Josephson junctions is negative, in panels (a) of
Figs. 4 and 5 we chose $\chi<0$ and $\omega_{j}-\omega_{b}>0$.

####  B.5.2 Nonmonotonic and strong suppression

As illustrated in panels (b) of Figs. 4 and 5, working in a parameter regime
where the behavior of the IPR with photon number is nonmonotonic can lead to
stronger suppression. A first observation to understanding this effect is
that, taking $p_{j-}\propto\tilde{\Delta}_{j}$ and $p_{j-}\rightarrow 0$,
forces Eq. 134 to be unity. Moreover, we can exploit the fact that symmetric
two-qubit interactions, such as the $ZZ$ interaction, are suppressed when the
bus frequency lies between the qubit frequencies. As a result, the choice
$\tilde{\Delta}_{1}+\tilde{\Delta}_{2}=0$ minimizes the $ZZ$ interaction in
Eq. 33 in the limit of large anharmonicities with respect to the detunings,
i.e. $|\tilde{\Delta}_{j}/K_{j}|\ll 1$ and
$|\tilde{\Delta}_{j}/\tilde{K}_{b}|\ll 1$.

Based on these arguments, we define a second parameter regime according to

$\displaystyle\left|\delta/\chi\right|\ll
1,\quad\left(\omega_{j}-\omega_{b}\right)\delta>0,\quad\left(\omega_{j}-\omega_{b}\right)\chi<0.$
(136)

The essential difference with respect to Eq. 135 is that $\delta$ has now
changed sign. In this regime, the ac-Stark shift of the bus mode changes the
sign of the qubit-bus detunings for some $\bar{\alpha}$, something that can
help suppressing spurious interactions.

It is worth highlighting that the protocol presented here works best for small
qubit-bus detunings $|\tilde{\Delta}_{j}|$ relative to $\left|\chi\right|$.
This is because the drive on the NLR, close to its resonance frequency,
renders the energy gaps between the stabilized states in the NLR smaller in
magnitude than in the bare energy spectrum of the NLR. Without anharmonicity
in the NLR, these energy gaps are predominantly set by $\left|\chi\right|$.

Finally, although we have focused on $\bar{\alpha}_{0}$ with the TQD protocol,
one could have chosen to grow $\bar{\alpha}_{1}$ for example, such as to still
suppress the dominant $0\leftrightarrow 1$ transition in the bus. The main
reason for the focus on $\bar{\alpha}_{0}$ is that, in this case, $\delta$ is
made small but $\delta+\chi$ is large, thus preventing transitions to higher
energy states of the NLR.

##  Appendix C Effective parametric modulation

In this section we provide supporting analytical derivations and numerical
results for the parametric modulation (PAM) scheme presented in Sect. IV and
discuss another decoupling scheme based on a longitudinal drive (LD) in the
NLR to mimic the effects of an anharmonicity. We also show that the two
schemes can be potentially combined to offer stronger suppression of qubit-
qubit interactions.

###  C.1 PAM: additional tones in the NLR drive

We comment on the effects of the fast-oscillating contributions in Eq. 42. As
previously stated, these terms impose a lower bound on the time-averaged
$1-$IPR. For large $\omega_{m}$ as compared to $\delta$ and $\chi$, this lower
bound is approximately $2\left(g/\omega_{m}\right)^{2}$. In the particular
case of $\omega_{m}=\omega_{0}\left|\bar{\alpha}\right|$, such that the
voltage drive amplitude in the NLR is displacement-independent, we observe
that the asymptotic behavior in the suppression is polynomial in
$\left|\bar{\alpha}\right|$. We also stress that an appreciable anharmonicity
$K_{r}$ in the NLR could in principle result in additional frequency
collisions due to negative $s\omega_{m}/\tilde{\Delta}_{1}$ where $s$ is an
integer. Finally and in addition to Fig. 6 we report the inverse participation
of the second qubit in Fig. 11. For both Figs. 6 and 11 we plot the time-
averaged inverse participation ratio obtained in Sect. B.2.2.

Figure 11: Dynamical decoupling in the bus: time-averaged inverse
participation ratio against photon number in the ground state. Here we assume
two additional tones $\omega_{r}\pm\omega_{m}$ with amplitude
$\omega_{m}\lambda/\bar{\alpha}$, and we choose $\omega_{m}$ such that the
voltage-drive amplitude is fixed for all $\bar{\alpha}$s. The regions are
bounded by $\lambda=\lambda_{0}$ where $J_{0}(\lambda_{0})=0$ and $\pm$10%
error on $\lambda_{0}$. Note that we exclude anharmonicity in the NLR
($K_{r}=0$). Parameters: $\delta/2\pi=-1.0$ MHz, $\chi/2\pi=-20.0$ MHz,
$(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz, $(\omega_{2}-\omega_{b})/2\pi=14.0$
MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz, and $g/2\pi=2.0$ MHz.

###  C.2 LD: longitudinal drive in the NLR

We now discuss how a longitudinal drive in the NLR can help recover a strong
suppression in absence of large anharmonicity and without the need for fine-
tuning. Our starting point is Eq. 99. In the $z\rightarrow\infty$ limit, the
asymptotic behaviour of the Bessel functions is

$\displaystyle J_{\nu}(z)\sim\sqrt{\frac{2}{\pi
z}}\cos\left(z-\frac{2\nu+1}{4}\pi\right),$ (137)

and we also have $J_{\nu}(-z)=(-1)^{\nu}J_{\nu}(z)$. It follows that the
matrix elements of Eq. 96 are renormalized by $\left|(k-\ell)z\right|^{-1/2}$
for $k\neq\ell$. This result is appealing given that the suppression of two-
qubit interactions is strongest if the NLR is constrained to Fock states of
small photon number and it is desirable to suppress $0\rightarrow k$
transitions in the NLR, as shown in Fig. 12.

A key advantage of this scheme is that it might be possible to relax parameter
constraints for $\tilde{\Delta}_{q}$ and $\chi$ due to the simulated
anharmonicity. Regarding the physical implementation, it can range from a
modulated detuning in the voltage drive to flux-modulating a superconducting
loop with junctions. The results shown in Fig. 12 correspond to the time-
averaged participation ratio obtained in Sect. B.2.3.

Figure 12: Dynamical decoupling in the NLR: time-averaged inverse
participation ratio against photon number in the ground state. We assume an
additional longitudinal drive
$\delta\rightarrow\delta-z\omega_{m}\sin\left(\omega_{m}t\right)$ in the NLR
with $\omega_{m}=1.0$ GHz. Note that no anharmonicity is considered for the
NLR ($K_{r}=0$). Parameters: $\delta/2\pi=-1.0$ MHz, $\chi/2\pi=-20.0$ MHz,
$(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz, $(\omega_{2}-\omega_{b})/2\pi=14.0$
MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz, and $g/2\pi=2.0$ MHz.

###  C.3 PAM and LD in parallel

A key challenge in combining PAM and LD is to prevent frequency collisions.
Having simultaneously large modulation frequencies for both schemes is
therefore not a good option. Given that the lower bound on $1-$IPR is
dominantly set by the modulation frequency in PAM, we choose this frequency to
be the largest. Since the modulation frequency in LD is small, we can choose
$z$ to be large. We indeed observe in Fig. 13 that it is possible to grow
stronger suppression factors by increasing $z$. Fig. 13 illustrates the time-
averaged inverse participation obtained in Sect. B.2.4.

Figure 13: Dynamical decoupling in the bus and in the NLR: time-averaged
inverse participation ratio against photon number in the ground state. Here we
assume two additional tones $\omega_{r}\pm\omega_{m}^{\phi}$ in the voltage
drive with amplitude $\omega_{m}^{\phi}\lambda/\bar{\alpha}$, and we choose
$\omega_{m}^{\phi}=0.5\bar{\alpha}$ such that the voltage-drive amplitude is
fixed for all $\bar{\alpha}$s with $\lambda=\lambda_{0}$ where
$J_{0}(\lambda_{0})=0$. We assume an additional longitudinal drive
$\delta\rightarrow\delta-z\omega_{m}^{\delta}\sin\left(\omega_{m}^{\delta}t\right)$
in the NLR with $\omega_{m}^{\delta}=10.0$ MHz. Parameters: $\delta/2\pi=-1.0$
MHz, $\chi/2\pi=-20.0$ MHz, $(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz,
$(\omega_{2}-\omega_{b})/2\pi=14.0$ MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz,
and $g/2\pi=2.0$ MHz.

##  Appendix D Superconducting implementation

The inverse participation ratio for the second qubit is illustrated in Fig. 14
and should be contrasted to that of panel c) in Fig. 7.

Figure 14: $1-\mathrm{IPR}_{0100}$ estimated with the numerical
diagonalization of the full system with the effective Hamiltonian in Eq. 55
with $\delta/2\pi=19.9$ MHz, $\chi/2\pi=-20.0$ MHz (parameter regime of c)
where $\left|\alpha_{0}\right|\rightarrow 0$). Here
$(\omega_{1}-\omega_{b})/2\pi=7.0$ MHz, $(\omega_{2}-\omega_{b})/2\pi=14.0$
MHz, $K_{1}/2\pi=K_{2}/2\pi=-300.0$ MHz, and $g/2\pi=2.0$ MHz.

We conclude with a remark on stray couplings, not captured by Eq. 2, but most
likely present in a superconducting circuit implementation. For instance,
direct coupling $g_{1-2}$ between the qubits cannot be suppressed by
manipulating the coupler and, in that case, $\mathrm{IPR}$ includes an
additional term $-2g_{1-2}^{2}/(\omega_{1}-\omega_{2})^{2}$ which bounds
$1-\mathrm{IPR}$. However, this bound can be conveniently lowered by detuning
the qubits and improving circuit design such that $g_{1-2}$ is minimized. It
is also worth noticing that this interaction is typically very small compared
to desired couplings. Another possibility for the presence of stray couplings
is to have spurious qubit-NLR interactions. In this case, virtual two-qubit
transitions mediated by the NLR are not exponentially suppressed. However, if
the NLR is far detuned in frequency with respect to the qubits, these
interactions can be greatly reduced. Finally, stray dispersive coupling
between the bus and the NLR can also exist. However, this type of
nonidealities are not particularly detrimental, as the resulting weak
hybridization between the bus and the NLR does not prevent the suppression of
two-qubit interactions.