2023

[1,2]Yao Lu These authors contributed equally to this work.

[1,2]Aniket Maiti These authors contributed equally to this work.

[1,2]Robert J. Schoelkopf

1]Departments of Applied Physics and Physics, Yale University, New Haven,
06511, CT, USA

2]Yale Quantum Institute, Yale University, New Haven, 06511, CT, USA

# High-fidelity parametric beamsplitting with a parity-protected converter

<EMAIL_ADDRESS><EMAIL_ADDRESS>John W. O. Garmon  Suhas Ganjam
Yaxing Zhang  Jahan Claes  Luigi Frunzio  Steven M. Girvin
<EMAIL_ADDRESS>[  [

###### Abstract

Fast, high-fidelity operations between microwave resonators are an important
tool for bosonic quantum computation and simulation with superconducting
circuits. An attractive approach for implementing these operations is to
couple these resonators via a nonlinear converter and actuate parametric
processes with RF drives. It can be challenging to make these processes
simultaneously fast and high fidelity, since this requires introducing strong
drives without activating parasitic processes or introducing additional
decoherence channels. We show that in addition to a careful management of
drive frequencies and the spectrum of environmental noise, leveraging the
inbuilt symmetries of the converter Hamiltonian can suppress unwanted
nonlinear interactions, preventing converter-induced decoherence. We
demonstrate these principles using a differentially-driven DC-SQUID as our
converter, coupled to two high-Q microwave cavities. Using this architecture,
we engineer a highly-coherent beamsplitter and fast ($\sim$ 100 ns) swaps
between the cavities, limited primarily by their intrinsic single-photon loss.
We characterize this beamsplitter in the cavities’ joint single-photon
subspace, and show that we can detect and post-select photon loss events to
achieve a beamsplitter gate fidelity exceeding 99.98$\%$, which to our
knowledge far surpasses the current state of the art.

###### keywords:

Bosonic, Beamsplitting, dual-rail qubit, multi-mode, parametric converter

## Introduction

The precise manipulation of high-Q bosonic modes is vital to experimentally
explore a wide range of phenomena in many-body physics and quantum
information. A crucial component in this endeavor is a programmable two-mode
interaction, exemplified by the time-dependent beamsplitter Hamiltonian
$\hat{\mathcal{H}}_{\text{BS}}/\hbar=g_{\text{BS}}(t)\left(e^{i\varphi_{\text{BS}}}\hat{a}^{\dagger}\hat{b}+e^{-i\varphi_{\text{BS}}}\hat{a}\hat{b}^{\dagger}\right)$.
This controlled photon-exchange coupling is a requirement for continuous
variable quantum computation Gaussian_QI , and has direct applications in
bosonic simulations of hopping models and lattice gauge theories
YaoMottInsulator ; CreutzLadderExperiment ; JensCQEDLattices ; Roushan2016 .
It is particularly appealing to implement such a Hamiltonian in circuit-QED, a
flexible platform offering readily available nonlinear control of high-Q modes
in superconducting resonators ReagorStubCavity . These resonators primarily
lose coherence through single-photon loss ECDControlPaper ;
SergeMushroomCavity , a ‘noise-bias’ that has been utilized in important
demonstrations including bosonic error-correction ChuangBosonicCodes ;
CochraneCatCode ; Leghtas4CatCode ; JamesDualRailPaper  beyond break-even
BeyondBreakeven2016 ; VladGKPBeyondBreakeven ; BinomialBeyondBreakeven .
However, the experimental implementation of a fast, high-fidelity beamsplitter
that preserves the long lifetime of these resonators and does not introduce
additional decoherence has remained a challenge. If realized, this interaction
would help construct logical entangling gates between qubits encoded in the
resonators Yvonne2019_Nature_eSWAP ; TakaPaper ; Joshi2021 , and enhance long-
distance interactions through microwave quantum buses
PhillipeDeterministicRemoteEntanglementPitchAndCatch ;
AxlineOnDemandStateTransfer .

Within the circuit-QED framework, interactions between linear superconducting
resonators can be implemented by coupling them via a Josephson junction-based
nonlinear ‘converter’ JPC_Paper ; SPA_1 ; JAMPA . Driving this converter with
RF drives can actuate a ‘parametric beamsplitter’ AumentadoDualRail ;
SiroisStorageRetrieval2015 ; Yvonne2018_PRX_ProgrammableInterference ; FPJA ;
ParametricCQEDPaper ; HatridgeRouter  between modes that are widely separated
in frequency space, offering large on-off ratios. The amplitude of the drives
sets the strength of the beamsplitter, so one could ideally improve both the
speed and fidelity of beamsplitting-based gates by simply driving harder.
However, controlling the dynamics of strongly driven, dissipative, nonlinear
systems can be challenging KorotkovBeyondRWAPaper ;
Yaxing_PRA_BilinearModeCoupling . The wide-bandwidth of the Josephson
nonlinearity can activate numerous parasitic processes that cause exchanges
with modes lossier than the resonators, including with the converter itself,
spoiling fidelity. The converter can also be incoherently excited through a
dressing of its natural decoherence in the presence of the drives, which can
directly dephase the beamsplitter interaction and harm the noise-bias of the
resonators. These parasitic processes and drive-induced coupler excitations
have been significant barriers in previous implementations
Yvonne2018_PRX_ProgrammableInterference ; Yvonne2019_Nature_eSWAP  and
suppressing them is crucial to engineering a clean, high-fidelity
beamsplitter.

One approach that has been used to tame this nonlinearity is engineering
multi-junction converters with useful symmetries JPC_Paper ; ATS_Paper , that
prevent a significant fraction of the nonlinear processes allowed in a single-
junction circuit like the transmon JensTransmonPaper . In this work we employ
a similar approach, using the familiar circuit of the symmetric DC
Superconducting QUantum Interference Device (SQUID), with careful engineering
to make it compatible with a high-Q environment. Making full use of this
symmetry requires the SQUID to be driven in a purely differential manner, and
we introduce an architecture for delivering this drive through an auxiliary
‘buffer’ mode.

In this work, we use the differentially-driven SQUID (DDS) to perform a fast,
highly coherent beamsplitter between two high-Q 3D superconducting cavity
resonators. We show that the coupler experiences almost no drive-induced
excitation, independent of the beamsplitter rate. The beamsplitter fidelity is
then limited only by the cavities’ single-photon decay, thus being highly
compatible with existing bosonic encodings ATS_Paper ; Grimm2020 ; Puri2020 .
As a first demonstration, we characterize the pulsed operation of this
beamsplitter with single photons in a microwave implementation of the dual-
rail qubit subspace $\\{|0_{a}1_{b}\rangle$, $|1_{a}0_{b}\rangle\\}$, where
the beamsplitter interaction provides universal control. We achieve an average
gate fidelity of $99.92\%$, which on detecting single-photon loss events is
boosted to $99.98\%$. This paves the way for quantum computing architectures
based on erasure-limited dual-rail qubits DivincenzoDualRail ;
AumentadoDualRail ; JamesDualRailPaper , and for cleaner parametric processes
in general.

Figure 1:  The differentially driven SQUID as a parity protected converter. a,
The symmetric DC-SQUID contains two orthogonal modes SQUIDModesLecocq2011 ;
KamalSQUIDModes , the common mode (coupler) and the differential mode
(actuator). We selectively couple the former to two bosonic modes and the
latter to the drives to take advantage of the natural symmetries of the
Hamiltonian in eq. 1. b, Implementing the purely differential drive through a
3D buffer post-cavity (figure is exaggerated for illustrative purposes). The
natural separation in electric and magnetic fields in the $\lambda/4$ mode is
used to purely drive the actuator, without exciting the coupler. The sensitive
quantum information is stored in two high-Q $\lambda/4$ post-cavities (Alice
and Bob) that participate in the coupler, enabling parametric beamsplitting
between them. The inset shows an optical micrograph of the SQUID device,
displaying the purposely offset antenna pad that counters residual drive-
asymmetry. c, Frequency stack for relevant modes in the system. The difference
of the two drive frequencies ($\Delta_{d}$) is fixed to be equal to the cavity
detuning ($\Delta_{ab}$) for resonant beamsplitting. The drives are placed
symmetrically around the buffer mode resonance, which is engineered to be far-
red detuned from the coupler frequency.

fig]fig1

## Results

### Engineering a differentially-driven SQUID

We now outline the key insight of properly leveraging the symmetry of the
SQUID for a cleaner beamsplitter, and describe an architecture that
accomplishes it. The SQUID has two orthogonal modes, a common ‘coupler’ mode
with symmetric Josephson phases across its two junctions, and a differential
‘actuator’ mode with anti-symmetric phases. Fully utilizing the symmetry of
these modes requires a precisely-engineered selective coupling of the coupler
and actuator to the resonators and the drive respectively (LABEL:fig1a), at
the sweet-spot of zero DC flux. This is different from a conventional flux-
driven SQUID OliverJPA ; FirstOrderSideband ; MckayCoupler ;
YaoUniversalStabilization , since extra care must be taken to ensure that the
drive does not excite the common mode. This drive must also be introduced in a
minimally invasive manner that does not spoil the lifetime of either the
resonators or the coupler.

The Hamiltonian of such a differentially-driven SQUID is (Supplementary Note
1):

$\hat{\mathcal{H}}_{\text{DDS}}=4E_{C}\hat{n}_{c}^{2}-E_{J}\underbrace{\cos(\phi_{d})\cos(\hat{\theta}_{c})}_{\text{Even
Parity}},$ (1)

where $\hat{\theta}_{c}=\left(\hat{\theta}_{1}+\hat{\theta}_{2}\right)/2$ and
$n_{c}$ are the conjugate phase and charge variables describing the coupler
mode, and $\phi_{d}=\left<\hat{\theta}_{1}-\hat{\theta}_{2}\right>/2$ is the
classical response of the driven actuator, for junction phases
$\hat{\theta}_{1,2}$. $E_{J}$ and $E_{C}$ are the total Josephson and charging
energies respectively. The advantage of engineering this Hamiltonian is two-
fold. First, the Hamiltonian is protected from processes that involve an odd
number of mode quanta, while preserving the strength of desired quadratic
process that provides the beamsplitting interaction. This ‘parity protection’
forbids up to half the parasitic coherent processes allowed by a single-
junction converter, including the swapping of sensitive information into the
coupler that limited previous implementations
Yvonne2018_PRX_ProgrammableInterference . Second, driving through the
orthogonal actuator port provides flexibility in choosing the drive
frequencies, which we can utilize to drive far away from non-protected
processes that could otherwise limit fidelity. Combining these two advantages
has the potential to allow the suppression of not only unwanted coherent
processes, but also incoherent coupler excitations
Yaxing_PRA_BilinearModeCoupling ; DykmanHeating  due to coupler decoherence
dressed by the drives (see Supplementary Note 2 for a detailed description).

Implementing the symmetric Hamiltonian in eq. 1 amounts to fabricating a SQUID
with symmetric junctions, calibrating it to zero DC flux, and delivering a
purely differential drive. The first two conditions are possible to optimize
in fabrication or calibrate out respectively, and have been described in
Supplementary Note 7. Engineering the differential drive, however,
necessitates the delivery high-frequency flux in a superconducting package,
while simultaneously eliminating stray electric fields that couple to the
common mode of the SQUID and accounting for spatial gradients in the time-
dependent flux (see Methods “Fine-tuning a differential drive” for more
details).

We simultaneously fulfill these constraints using the $\lambda/4$ mode of a
stub cavity (dubbed the ‘buffer mode’), integrated into the same monolithic
package that hosts the high-Q storage modes (LABEL:fig1b). Once driven, the
buffer mode provides an oscillating magnetic field that serves as the flux
drive penetrating the SQUID loop. We place the SQUID at the base of the buffer
mode, which functions as a virtual ground, minimizing the common mode drive
due to the driven electric field while maximizing the flux penetrating the
SQUID loop. The coupling to any remnant electric field is further suppressed
by orienting the SQUID’s electric dipole moment to be perpendicular to the
driven field. In addition to this minimization of the common-mode drive,
driving through the buffer mode also has the benefit of imposing a finite
bandwidth ($\sim 250$ MHz) on the noise spectrum of the environment seen
through the drive line.

Finally, the full elimination of the common mode drive is achieved by a
deliberate offset of the SQUID’s capacitive pad. This counters the effects of
flux gradients in our system, by engineering an additional flux-dependent
electromotive force, which can be directly optimized using finite-element-
method simulations (see Methods). In experiment, imperfections in fabrication
or placement of the SQUID chip might affect the achievable drive
orthogonality. However, we are able to directly estimate this orthogonality in
our actual device by experimentally comparing the strength of an allowed
process to one that should be forbidden by parity protection , finding a
residual common drive on the order of $|\phi_{c}/\phi_{d}|\sim 1\%$ (see
Methods “Experimentally characterizing residual drive asymmetry”).

### Demonstrating a high-coherence beamsplitter

We now present the full experimental realization (LABEL:fig1b) of a high-
fidelity beamsplitter that strongly suppresses undesirable coupler heating.
Our construction consists of a high-purity aluminum package hosting three
coaxial $\lambda/4$ stub cavity modes: the buffer mode and the two high-Q
storage modes ReagorStubCavity . The storage modes are capacitively coupled to
the Y-shaped antenna of the SQUID Ymon_Paper  and have a negligible mutual
inductance to the SQUID loop, ensuring that they exclusively participate
BBQPaper ; EPRPaper  in the coupler mode:

$\hat{\theta}_{c}\approx\left(\frac{2E_{C}}{E_{J}}\right)^{\frac{1}{4}}\left(\frac{g_{a}}{\Delta_{a}}\;\hat{a}+\frac{g_{b}}{\Delta_{b}}\;\hat{b}+\hat{c}\right)+\text{h.c.}$
(2)

Here $\hat{a}$, $\hat{b}$ are the ladder operators for our dressed storage
modes, named Alice and Bob respectively, while $\hat{c}$ represents the
dressed coupler. The coupling strengths ($g_{a,b}$) and mode detunings
($\Delta_{a,b}$, LABEL:fig1c) between the storage modes and the coupler are
chosen to be in the dispersive regime
$\left(\frac{g_{a}}{\Delta_{a}},\frac{g_{b}}{\Delta_{b}}\sim 0.1\right)$.
Additionally, we can prepare and readout Fock states in Bob through a
dispersively coupled ancilla transmon DispersiveUniversalControl  and a
dedicated stripline readout resonator. The SQUID coupler also includes a
dedicated readout resonator and drive pin, for explicit characterization of
frequency shifts and drive-induced excitation. The measured device parameters
are presented in Supplementary Note 12, and details of device fabrication can
be found in Supplementary Note 11.

We activate and control the amplitude and phase of our beamsplitter with a bi-
chromatic drive on the actuator:
$\phi_{d_{1,2}}(t)=|\phi_{d_{1,2}}|\cos(\omega_{d_{1,2}}t+\varphi_{d_{1,2}})$.
When the difference in our drive frequencies ($\Delta_{d}$) is close to our
cavity detuning ($\Delta_{ab}$), eq. 1 and eq. 2 combine to create a tunable
beamsplitter Hamiltonian (see Methods for full derivation):

$\displaystyle\hat{\mathcal{H}}_{\text{BS}}/\hbar=\;$
$\displaystyle\Delta_{\text{BS}}\;\hat{a}^{\dagger}\hat{a}+g_{\text{BS}}\;(e^{i\varphi_{\text{BS}}}\hat{a}^{\dagger}\hat{b}+e^{-i\varphi_{\text{BS}}}\hat{a}\hat{b}^{\dagger}),$
(3) $\displaystyle\text{with \hskip 10.00002pt}g_{\text{BS}}$
$\displaystyle\approx\;\dfrac{\omega_{c}}{2}\>\dfrac{g_{a}\;g_{b}}{\Delta_{a}\Delta_{b}}\;J_{1}\left(|\phi_{d_{1}}|\right)J_{1}\left(|\phi_{d_{2}}|\right),$
(4) $\displaystyle\Delta_{\text{BS}}$
$\displaystyle=\Delta_{ab}-\Delta_{d}+\Delta_{\text{Z},ab}$ (5)

where $J_{1}(|\phi_{d_{1,2}}|)$ is the first-order Bessel function of the
drive amplitudes, and $\varphi_{\text{BS}}$ is the beamsplitter phase,
controlled by the relative phase of the drives. The drives also induce an
additional frequency offset, the relative AC-Zeeman shift of the cavities,
$\Delta_{\text{Z},ab}$, and we experimentally find the amplitude-dependent
condition $\Delta_{\text{BS}}=0$ to execute a resonant beamsplitter.

Figure 2:  Beamsplitting with the differentially-driven SQUID. a,
Beamsplitting implements an effective driven Rabi evolution in the Bloch
sphere of the dual-rail qubit formed by the single photon subspace Alice and
Bob, where decay can be detected by monitoring the vacuum state. b, Resonant
evolution of a single-photon prepared in Bob. The data is normalized for
readout infidelity, and state preparation fidelity is shown as a dashed grey
line (Supplementary Note 10). The fast coherent oscillations (black dots)
between the cavities are fitted to eq. 6 (green lines show envelope) to obtain
the decay and dephasing time-scales. The evolution for the first $1.5\;\mu$s
is plotted separately to better illustrate the oscillations, and fit to a
sinusoid to extract $g_{\text{BS}}$. c, Sweeping both drive amplitudes
simultaneously and repeating experiment a lets us quantify $g_{\text{BS}}$
(blue crosses), and the decoherence limit on beamsplitter infidelity (red
diamonds) at various drive strengths. We choose a drive strength with
simultaneously low infidelity and high beamsplitter rate as our operating
point (yellow dashed line). d, The coupler’s driven excitation ($P_{c}$) after
evolving for 10 swaps is directly quantified through a dedicated on-chip
readout mode. We observe no monotonic correlation with respect to drive
amplitude, and driven populations mostly remain within the range of the
undriven population (grey region). e, Coupler population as a function of
number of swaps at the operational driving point. The heating rate is nearly
immeasurable, with a fitted (pink line) slope of $(1.2\pm 2.4)\times 10^{-5}$
excitation per swap, which is within expectation for our natural thermal
background ($\gamma_{c,\uparrow}\sim(3.3\textrm{ ms})^{-1}$). The non-zero
offset of the fit arises from preparation and readout infidelities. Error bars
in both d and e represent fit errors from the protocol described in
GeerlingsRPM .

fig]fig2

We characterize our beamsplitter interaction using the joint single-photon
subspace of our storage cavities, which forms a microwave implementation of a
dual-rail qubit (LABEL:fig2a). We initialize a single photon in Bob with a
preparation fidelity of $\sim 94\%$, by displacing Bob to a coherent state
($\alpha_{b}=\sqrt{2}$) and using number-resolved measurements through the
ancilla to post-select the desired Fock state. We then apply the resonant
beamsplitter interaction for a range of times up to 32 $\mu$s (LABEL:fig2b),
to estimate both the beamsplitting strength and the driven decoherence time of
the dual-rail qubit. The evolution of the time-dependent photon population in
Bob follows (see Methods “Decoherence in the dual-rail subspace” for
derivation):

$P_{\text{Bob}}=\dfrac{1}{2}e^{-\kappa_{1}t}\left(1+e^{-\kappa_{\varphi}t}\;\cos{(2g_{\text{BS}}t)}\right)$
(6)

Here, $\kappa_{1}$ is the mean of the driven cavities’ single-photon decay
rates, and is the effective rate of population leakage out of the dual-rail
subspace into vacuum $|0_{a}0_{b}\rangle$. The dual-rail qubit may also
experience dephasing at a rate $\kappa_{\varphi}$, which would drive this
evolution towards an evenly mixed state within the qubit subspace. Combining
these lets us place a lower bound on the expected decoherence limit on the
fidelity of a single beamsplitter operation (see Methods):

$\mathcal{F}\approx
1-\frac{\pi}{4}\frac{\kappa_{\text{BS}}}{g_{\text{BS}}},\hskip
10.00002pt\kappa_{\text{BS}}=\kappa_{1}+\frac{\kappa_{\varphi}}{2},$ (7)

which is a more accurate metric than the one used in
Yvonne2018_PRX_ProgrammableInterference .

We choose our operating point to maximize this expected fidelity $\mathcal{F}$
with respect to drive strength, which we characterize by analyzing sections of
the resulting driven long-time evolution at various drive amplitudes
(LABEL:fig2c). At each amplitude, we extract the $g_{\text{BS}}$ and
$\kappa_{\text{BS}}$ by fitting two short sections of the evolution (Methods)
to eq. 6. We are able to obtain a maximum $g_{\text{BS}}/2\pi$ exceeding 5
MHz, with the effective decoherence-limited fidelity surpassing $99.9\%$ for a
wide range of beamsplitting strengths. We observe that upon increasing the
drive strength, the coupler frequency shifts closer to the cavities (by $\sim
200$ MHz), which may be accompanied by direct sideband interactions between
the coupler and the cavities (Supplementary Note 5), leading to an increased
hybridization between the modes. This effect results in a faster-than-
quadratic dependence of the beamsplitting strength on our drive amplitudes,
but limits the fidelity of the beamsplitter at higher amplitudes due to the
aggravated Purcell loss.

At the operating point, we fit the evolution in LABEL:fig2b to find
$g_{\text{BS}}/2\pi=2.16\pm 0.01\textrm{MHz}$, $\kappa_{1}=(197\pm
8\;\mu\mathrm{s})^{-1}$ and $\kappa_{\varphi}=(313\pm
40\;\mu\mathrm{s})^{-1}$. Our effective $\kappa_{\text{BS}}=(150\pm
25\;\mu\text{s})^{-1}$ places a decoherence-based upper bound on the fidelity
of $\mathcal{F}=99.96\pm 0.01\%$, which is almost two orders of magnitude
better than the previous transmon-based implementation. Crucially, this
fidelity is also limited primarily by photon loss in the cavities, preserving
their advantageous noise bias.

To directly quantify the suppression of drive-induced coupler excitation, we
measure the coupler population as a function of drive amplitude after
resonantly evolving for ten swaps. With our coupler prepared in the ground
state, we apply this pulse and measure the coupler’s population through a
protocol that is robust to readout infidelity GeerlingsRPM . We measure no
correlated increase of its driven population as a function of drive amplitude,
up to our measurement uncertainty of $\sim 0.2\%$ (LABEL:fig2d). At the
operating point, we explicitly quantify the increase in coupler excitation as
a function of number of swaps (LABEL:fig2e). We evolve the system up to 100
swaps and observe a total heating rate below $\sim 4\times 10^{-5}$
excitations per swap, which is consistent with the undriven heating rate of
the coupler, implying no additional drive-induced heating. This substantial
suppression (three orders of magnitude better than transmon-based
implementations Yvonne2018_PRX_ProgrammableInterference ) eliminates
limitations placed by coupler-induced dephasing on the fidelity of the
beamsplitter, allowing us to harness the long lifetimes and even longer
dephasing times of the 3D cavities.

While this device was optimized for performance in the single-photon manifold,
any drive-induced increase in the self-Kerr of the cavities would lead to
coherent errors in higher-photon manifolds. The sideband interaction resulting
from our choice of mode frequencies exacerbates this effect, with up to 128
KHz of inherited Kerr at the operating point. Numerous avenues exist to
minimize this driven non-linearity if desired (see Supplementary Note 6),
including choosing a slightly higher coupler frequency, arraying multiple
SQUIDs, or dynamical Kerr cancellation YaxingKerrCancellation . If minimizing
inherited self-Kerr is vital, such as in schemes utilizing large coherent
states, then implementing these improvements or using an alternative scheme
like Kerr-free three-wave mixing SPA_1 ; ChapmanSNAILBeamsplitterPaper  may be
required.

Figure 3:  Randomized benchmarking with a calibrated beamsplitter pulse. a,
Wigner function of Bob after preparing $|0_{a}1_{b}\rangle$ and implementing
1, 10 and 60 calibrated swaps. Each swap is a combination of two identical
beamplitter pulses. b, Probability of ending in the target state
$|0_{a}1_{b}\rangle$ after executing the RB protocol with randomly selected
gates from $G_{\textrm{DR}}$ (yellow circles). The curve is normalized to
account for state preparation and readout imperfections (Supplementary Note
10), and is in good agreement with a single-exponential, with a decay constant
of $\tau_{\text{RB}}=1271\pm 4$ gates. This ‘Raw RB’ is practically
indistinguishable from the sequences post-selected on the coupler ground state
(black crosses), with the difference of the two curves shown below the main
plot (pink crosses). c, Focusing on the first $2250$ gates, we use
measurements of both cavities to post-select on sequences in which no photon
loss event occurred (green diamonds). We compare these sequences to the raw RB
in b (yellow), showing an improvement in average gate infidelity from
$0.078\pm 0.001\%$ to $0.020\pm 0.001\%$. d, The gate-sets required for the
above protocols are generated from calibrated beampslitter pulses with tanh-
shaped ramps, where different $U_{\text{BS}}(\varphi)$ are obtained by
changing the relative phase of our drives. e, The benchmarking sequences
consist of randomly generated pulses that, under ideal operation, map
$|0_{a}1_{b}\rangle$ back to itself. After each sequence, we measure whether
the coupler is in its ground state, the presence of a photon in Bob, and the
presence of a photon in Alice using an additional swap gate. This lets us
generate the raw, coupler-selected and leakage-detected RB datasets. All
sequences are also conditioned on Bob’s ancilla ending in its ground state, to
discount first-order effects of ancilla heating.

fig]fig3

### Benchmarking fidelity in the single-photon subspace

We precisely characterize the fidelity and noise bias of our beamsplitter by
using it to implement universal control of the dual-rail qubit subspace,
allowing techniques akin to standard randomized benchmarking (RB) protocols
RBKnill ; EmersonRBPrecursor ; RBKristine . First, we identify the amplitude
and frequency required for a fixed-length pulse (LABEL:fig3d) to achieve the
beamsplitter unitary

$U_{\text{BS}}\left(\varphi\mathord{=}0\right)=e^{i\pi/4(\hat{a}^{\dagger}\hat{b}+\hat{a}\hat{b}^{\dagger})},$
(8)

by repeating the pulse to perform swaps between Alice and Bob, and iteratively
checking its performance up to $\sim 1000$ such repetitions (Supplementary
Note 8). This calibrates an effective $X_{\pi/2}$ gate for the dual-rail
qubit, with a relative amplitude precision of less than $3\times 10^{-6}$ and
frequency precision of less than $3\times 10^{-7}$. Through phase control and
repetition, we use this pulse to construct a set of native gates,

$G_{\text{DR}}=\\{X_{\pi/2},Y_{\pi/2},X_{-\pi/2},Y_{-\pi/2},X_{\pi},Y_{\pi}\\},$
(9)

that generate the Clifford group for the dual rail qubit (LABEL:fig3d,
Supplementary Note 9). This allows a form of direct randomized benchmarking
DirectRB , which under uniform sampling should convert both dephasing and
coherent control errors into an effective depolarization channel. The dominant
but detectable error of cavity photon loss appears as a leakage to the
orthogonal state $|0_{a}0_{b}\rangle$, which is not converted to
depolarization under this protocol, but can be separately quantified and
selected out in post-processing.

The RB protocol consists of initializing the system in $|0_{a}1_{b}\rangle$
with the coupler and ancilla prepared in their ground states, and running
sequences of varying lengths of randomly chosen gates from $G_{\text{DR}}$.
Each sequence ends with an additional gate from $G_{\text{DR}}$ that maps the
state back to $|0_{a}1_{b}\rangle$, after which the presence of an excitation
in the coupler, Bob, and Alice are measured (LABEL:fig3e). We explicitly
discount the effects of the ancilla in all results shown by separately
measuring the ancilla and only including sequences where it ends in the ground
state. The sequences range up to 8100 gates, where we choose up to 900 random
gates (limited by FPGA memory), and repeat each gate nine times to fully
capture the fidelity decay timescale (this is roughly equivalent to performing
nine times as many random gates and provides a lower-bound for the fidelity of
a non-repeated sequence, see Supplementary Note 9). We average over $\sim
10^{5}$ such semi-random sequences at each sequence length.

We observe that the average success probability of returning to Bob
(LABEL:fig3b) decays exponentially, with an effective time constant of
$\tau_{\text{RB}}=1271\pm 4$ gates. This curve is normalized to account for
effects from state preparation and readout infidelity. We compare this curve
to post-selected sequences where the coupler explicitly ended in the ground
state, finding remarkable similarity. This comparison strongly suggests that
the coupler-heating induced errors on actual beamsplitter-based gates have
been suppressed to below the measurement accuracy, agreeing with our earlier
measurement of diminished driven coupler excitation.

Finally, we quantify our beamsplitter infidelity and noise bias by comparing
the raw protocol to sequences where the system was observed to have not
experienced a photon loss error. We focus on the span of the first $2250$
gates for this comparison (LABEL:fig3c, plotted on a log scale), where the raw
decay of Bob’s population suggests an un-selected single gate infidelity of
$0.078\pm 0.001\%$. On detecting out leakage events, our evolution changes to
a curve that decays toward a steady-state fidelity of $\sim 0.5$ instead of 0,
as expected (Supplementary Figure 7). For a fair comparison between this curve
and the raw RB, we re-scale this fidelity to also have a steady-state of zero
(Supplementary Note 10). This helps illustrate the clear improvement under
error-detection to an infidelity per gate of $0.020\pm 0.001\%$, which shows
that discarding only one out of every $\sim 1300$ shots per gate can lead to a
$3.9\pm 0.2\times$ increase in gate performance, enabled by the cavity noise-
bias. We reiterate that the fact that the dominant errors remain photon loss,
which is detectable or correctable with various bosonic encoding schemes,
satisfies one of the key goals for a high-performance inter-cavity
interaction.

Notably, because our gate-set is crafted from nearly identical beamsplitter
pulses, with a random gate from $G_{DR}$ containing $4/3$ beamsplitters on
average, we can directly convert these gate infidelities into an effective
single beamsplitter fidelity. Our measurements thus imply an effective un-
selected beamsplitter fidelity of $99.941\pm 0.001\%$, which improves on
leakage detection to $99.985\pm 0.001\%$. The remaining errors after leakage
detection can be due to intrinsic dephasing of the cavities, drifts in our
control electronics, or other effects that are not treated by the post-
selection protocols, like cascaded heating and decay events of the ancilla
JamesDualRailPaper .

## Discussion

In conclusion, we constructed a microwave implementation of a tunable cavity-
cavity beamsplitter and characterized its performance within the cavities’
joint single-photon subspace. We obtained a beamsplitter gate fidelity
exceeding $99.94\%$ in this subspace, and were limited by detectable single-
photon loss in the cavities. This performance was enabled by carefully
engineering the drive frequencies and leveraging fundamental symmetries of the
nonlinear converter to keep the coupling mode in its ground state even when
driving a fast beamsplitter. Such high-fidelity control of a strongly driven
nonlinear element in the presence of decoherence is a significant step forward
for fast parametric operations in circuit-QED.

While this system was optimized to engineer a beamsplitter between 3D
cavities, the generality and versatility of the design framework far exceed
this specific application. One could apply such a beamsplitter to on-chip
resonators, or to phononic modes in hybrid architectures. Other low-order
mixing processes allowed by the converter’s symmetry, such as two-mode or
single-mode squeezing, could also be implemented while still reaping the
benefits of parity protection and suppressed coupler-induced infidelity.
Parity protection is also not exclusive to the SQUID Hamiltonian JPC_Paper ;
ATS_Paper , and future converter designs could leverage more advantageous
forms of this selection rule.

Beyond the context of parametric interactions, we have also demonstrated the
delivery of AC flux in a high-Q 3D environment. This can be used to control
other devices in similar architectures that require fast-flux modulation. The
electromagnetic simulation techniques we have developed are also readily
applicable to other work involving driven circuit-QED systems, for both 3D and
planar devices.

Finally, the demonstration of high-fidelity control in the dual-rail subspace
motivates the hypothesis that this subspace could itself be used as a
computational qubit DivincenzoDualRail ; AumentadoDualRail ;
JamesDualRailPaper . The error-hierarchy of detectable decay over dephasing
makes the dual-rail qubit amenable to erasure conversion
ShrutiErasureConversion , an approach potentially yielding high thresholds in
the surface-code architecture. The single-qubit control demonstrated here can
be extended to realize a high-fidelity gate-set TakaPaper  for multi-qubit
control, charting a path towards a general dual-rail qubit-based architecture
in circuit-QED. The performance in higher-photon manifolds is the subject of
further research, with promising avenues including arraying the SQUID element,
and creating parity-protected couplers with suppressed Kerr nonlinearity.

## Methods

### Deriving the programmable beamsplitter Hamiltonian

We now show how differential driving at the correct resonance condition
generates a beamsplitter between Alice and Bob. For the sake of simplicity in
derivations, we set $\hbar=1$. The differentially-driven SQUID Hamiltonian in
eq. 1 can be expanded in the ladder operators of the bare SQUID coupler
($\hat{\tilde{c}},\hat{\tilde{c}}^{\dagger}$), in a frame rotating at the
coupler frequency ($\omega_{c}$), and gives (up to quartic order):

$\displaystyle\hat{\mathcal{H}}_{tot}^{(4)}=$
$\displaystyle-\frac{E_{C}}{2}\;{\hat{\tilde{c}}^{\dagger}}^{2}{\hat{\tilde{c}}}^{2}$
$\displaystyle+E_{J}\left(\cos\left(\phi_{d}\right)-1\right)\;\theta_{c,\text{zpf}}^{2}\;{\hat{\tilde{c}}^{\dagger}}{\hat{\tilde{c}}}.$
(10)

Here, $\phi_{d}$ is assumed to be driven at frequencies far off-resonant from
the coupler squeezing condition, and
$\theta_{c,\text{zpf}}=\left(\frac{2E_{C}}{E_{J}}\right)^{\frac{1}{4}}$ is the
zero-point fluctuations of the bare coupler phase.

The drive ($\phi_{d}$) modulates the bare coupler frequency, which enacts a
beamsplitting process between the cavities Alice and Bob ($\omega_{a,b}$)
through their linear capacitive couplings to the coupler ($\omega_{c}$):

$\hat{\mathcal{H}}_{\text{coupling}}=\sum_{\hat{\tilde{k}}=\hat{\tilde{a}},\hat{\tilde{b}}}g_{k}(\hat{\tilde{k}}+\hat{\tilde{k}}^{\dagger})(\hat{\tilde{c}}+\hat{\tilde{c}}^{\dagger})$

Here $\hat{\tilde{k}}=\hat{\tilde{a}},\hat{\tilde{b}}$ represent the bare
modes of Alice and Bob respectively. The coupling strengths
$g_{a,b}/2\pi\approx 80$ MHz are designed to place us in the dispersive regime
($\Delta_{a,b}=\omega_{a,b}-\omega_{c}\gg g_{a,b},\;E_{C}$) and therefore we
can follow standard circuit-QED derivations Krantz2019 ; Blais2021  to
diagonalize the linear part of the Hamiltonian, yielding eq. 2 in the main
text. We then move into a frame rotating at the dressed frequencies of Alice,
Bob and the coupler respectively and keep relevant terms:

$\displaystyle\hat{\mathcal{H}}_{\text{tot}}^{(\text{RWA})}=$
$\displaystyle\;\hat{\mathcal{H}}_{\text{static}}^{(\text{RWA})}+\hat{\mathcal{H}}_{\text{driven}}^{(RWA)}$
$\displaystyle\hat{\mathcal{H}}_{\text{static}}^{(\text{RWA})}=$
$\displaystyle\;-\sum_{k=a,b}2E_{C}\beta_{k}^{2}\hat{c}^{\dagger}\hat{c}\hat{k}^{\dagger}\hat{k}-\frac{E_{C}}{2}{\hat{c}^{\dagger^{2}}}\hat{c}^{2}$
$\displaystyle\hat{\mathcal{H}}_{\text{driven}}^{(RWA)}=$
$\displaystyle\;\frac{\omega_{c}}{2}(\cos\left(\phi_{d}\right)-1)\times$
$\displaystyle\sum_{k,k^{\prime}=a,b,c}\beta_{k}\beta_{k^{\prime}}\;\hat{k}^{\dagger}\hat{k^{\prime}}\;e^{i\left(\omega_{k}-\omega_{k^{\prime}}\right)t}$
(11)

where $\hat{a},\hat{b},\hat{c}$ are the ladder operators for the dressed
Alice, Bob and coupler modes respectively,
$\beta_{a,b}=\frac{g_{a,b}}{\Delta_{a,b}}$ are the participations of Alice and
Bob in the coupler mode, and $\beta_{c}\approx 1$.

Next, we introduce a pair of drive tones whose response on the actuator is
given by:

$\phi_{d}=\sum_{1,2}\phi_{d_{1,2}}(t)=\sum_{1,2}|\phi_{d_{1,2}}|\sin\left(\omega_{d_{1,2}}t+\varphi_{d_{1,2}}\right)$
(12)

which are IQ-modulated signals, with Local Oscillators (LOs) that are directly
derived by down-converting Bob and Alice’s LOs, respectively, using a
reference signal $\omega_{\Delta}$ (see Supplementary Note 13). This places
the difference of the two drive frequencies close to the Alice-Bob detuning by
design:

$\displaystyle\omega_{d_{2,1}}=\omega_{\Delta}-\omega_{a,b}+\Delta_{d_{2,1}}^{(\text{ssb})}$
(13) $\displaystyle\implies$
$\displaystyle\Delta_{d}=\omega_{d_{2}}-\omega_{d_{1}}=\Delta_{ab}+\Delta_{d_{2}}^{(\text{ssb})}-\Delta_{d_{1}}^{(\text{ssb})}$
(14)

Figure 4:  Measured driven Zeeman shift and cavity swaps. a, We measure the
coupler frequency through direct spectroscopy under a single drive tone, as a
function drive amplitude. Comparing the Zeeman shift for either drive tone
(pink, purple) to the prediction from Floquet simulation (grey solid line)
allows accurate calibration of the drive strength in terms of the driven
junction phase $\phi_{d1,2}$. b, Measured population in Bob in the presence of
both drive tones as a function of the drive-detuning
($\Delta_{d}-\Delta_{ab}$) and the time of evolution. A single photon is
prepared in Bob, and the drives swap this photon between Alice and Bob under
the detuned beamsplitter interaction. c, Fitting the oscillations as a
function of drive-detuning to the Rabi model allows us to calibrate amplitude-
dependent shift in the resonance condition ($\Delta_{Z,ab}$) and strength
$g_{\text{BS}}$ of the beamsplitting interaction.

fig]fig:Chevron_Zeeman

Here $\omega_{\Delta}\approx\omega_{a}+\omega_{b}-\omega_{\text{buffer}}$ is
chosen to place the the drive tones symmetrically about the buffer mode
resonance, and $\Delta_{d_{1,2}}^{(\text{ssb})}$ are the detunings in the
sideband modulation frequency for each drive tone with respect to the
sidebands used to control Alice and Bob respectively. We set
$\Delta_{d_{2}}^{(\text{ssb})}=0$ and use
$\Delta_{d_{1}}^{(\text{ssb})}\ll\Delta_{d}$ as a finely tunable parameter to
find the true beamsplitter resonance condition. Similarly, the drive phases
($\varphi_{d_{1,2}}$) are defined with respect to the control of Alice and Bob
respectively, with $\varphi_{d_{2}}=0$ and
$\varphi_{d_{1}}=\left\\{0,\pi/2,\pi,-\pi/2\right\\}$ being used to generate
the various pulses required for the randomized benchmarking sequence.

In this limit of
$\Delta_{d}\approx\Delta_{ab}\ll\omega_{d_{1}}+\omega_{d_{2}}$, we can
simplify the drive using the Jacobi-Anger expansion and only keep terms that
are either static or rotating at the detuning frequency (up to second order in
drive strength):

$\displaystyle\cos(\phi_{d})$ $\displaystyle\approx
J_{0}\left(|\phi_{d_{1}}|\right)J_{0}\left(|\phi_{d_{2}}|\right)$
$\displaystyle+2J_{1}\left(|\phi_{d_{1}}|\right)J_{1}\left(|\phi_{d_{2}}|\right)\cos\left(\Delta_{d}t+\varphi_{d_{1}}\right),$
(15)

where $J_{0}$ and $J_{1}$ are the zeroth and first order Bessel functions.
Substituting this into eq. 11, we find, keeping only slow-rotating terms:

$\displaystyle\hat{\mathcal{H}}_{\text{driven}}^{(\text{RWA})}\approx$
$\displaystyle\frac{\omega_{c}}{2}\left[J_{0}\left(|\phi_{d_{1}}|\right)J_{0}\left(|\phi_{d_{2}}|\right)-1\right]\sum_{k=a,b,c}\beta_{k}^{2}\hat{k}^{\dagger}\hat{k}$
$\displaystyle+\frac{\omega_{c}}{2}J_{1}\left(|\phi_{d_{1}}|\right)J_{1}\left(|\phi_{d_{2}}|\right)\beta_{a}\beta_{b}$
$\displaystyle\times\left(e^{i(\Delta_{d}-\Delta_{ab})t+\varphi_{d1}}\hat{a}\hat{b}^{\dagger}+\text{h.c.}\right)$
(16)

Importantly, this approximation only holds when the drive tones don’t activate
any other higher-order resonances. The first part of this driven Hamiltonian
is a driven frequency shift, which we call the AC-Zeeman shift, and is
reminiscent of the Stark shift in a regular charge-driven transmon. This is
easier to see when rewriting it for small amplitude, where
$\left(1-J_{0}(x)\right)\approx x^{2}/4$:

$\Delta_{\text{Z},k}=-\beta_{k}^{2}E_{C}\left(\frac{|\phi_{d_{1}}|^{2}+|\phi_{d_{2}}|^{2}}{\theta_{c,\text{zpf}}^{2}}\right)$
(17)

Specifically, the frequency shift of the coupler ($\Delta_{\text{Z},c}$)
provides a simple way to experimentally calibrate the strength of the drives
(see LABEL:fig:Chevron_Zeemana).

Finally, we derive the beamsplitting Hamiltonian by assuming the coupler to be
in the ground state ($\langle c^{\dagger}c\rangle=0$), which eliminates
contributions from $\hat{\mathcal{H}}_{\text{static}}^{(\text{RWA})}$. We move
into frames rotating at Zeeman-shifted frequencies for both Alice and Bob
($\omega_{a,b}+\Delta_{Z,(a,b)}$). We carry out an additional frame
transformation for Alice, to a frame rotating at
$\Delta_{\text{BS}}=\Delta_{ab}-\Delta_{d}+\Delta_{\text{Z},ab}$, which gives
us the beamsplitter Hamiltonian:

$\hat{\mathcal{H}}_{\text{BS}}/\hbar=\Delta_{\text{BS}}\hat{a}^{\dagger}\hat{a}+g_{\text{BS}}\left(e^{i\varphi_{\text{BS}}}\hat{a}^{\dagger}\hat{b}+\text{h.c.}\right).$
(18)

Here, $\varphi_{\text{BS}}=\varphi_{d_{1}}$,
$\Delta_{Z,ab}=\Delta_{Z,a}-\Delta_{Z,b}$ is the relative Zeeman shift of the
cavities, and $g_{\text{BS}}$ is the value described in the main text in eq.
4.

We find $\Delta_{Z,ab}$ to be on the order of $g_{\text{BS}}$, which needs to
be taken into account when finding the resonance condition. We approximately
calibrate for it by preparing a single photon in Bob and examining the
system’s resonant evolution as a function of $\omega_{d_{1}}$. This produces a
Chevron-like pattern (LABEL:fig:Chevron_Zeemanb), whose oscillations follow a
detuned Rabi model:
$\omega_{\text{osc}}=\sqrt{4g_{\text{BS}}^{2}+\Delta_{\text{BS}}^{2}}$.
Fitting to this model (LABEL:fig:Chevron_Zeemanc) lets us find both the
$g_{\text{BS}}$ at any amplitude and the drive detuning required for
$\Delta_{\text{BS}}=0$.

Figure 5:  Minimizing residual drive asymmetry in the SQUID via fine-tuning
the circuit geometry. a, The SQUID device in the buffer cavity package driven
by an oscillating B-field. The cylindrical geometry of the buffer cavity
(outer wall radius greater than the inner wall) dictates that the B-field has
a non-uniform distribution along the radial direction (from left to right).
The dashed lines with capacitor symbol represents the capacitance between the
antenna pads of the SQUID and the wall of the package. b, The lumped-element
circuit model of the SQUID device, taking into account the geometry of SQUID
and the spatial distribution of the B field. c, HFSS simulation of the quality
factor of the common mode (red) and drive asymmetry (blue, defined in eq. 21)
as functions of the top pad displacement from center to right, $\delta$. As a
result of introducing this asymmetry in the device geometry, the optimization
of the drive asymmetry and the quality factor are simultaneously achieved at
$\delta\approx 350\mu$m.

fig]HFSS_asymmetry

### Fine-tuning a differential drive

We realize a purely differential drive by utilizing the separation of the
coupler and actuator in electromagnetic space — the coupler has a large
electric dipole moment, while the actuator appears primarily as a magnetic
dipole. This provides a convenient way to address these two modes in a 3D
architecture, which is not immediately present in other parity-protected
devices like the asymmetrically threaded SQUID (ATS) ATS_Paper  or the
Josephson parametric converter (JPC) JPC_Paper . We require our drive to
imitate a spatially uniform AC flux in the SQUID loop, with minimal stray
capacitive coupling to its electric dipole. The introduction of such a high-
frequency flux-bias line into a superconducting package is challenging with
existing techniques, since it must both preserve the lifetime of the cavities
and also effectively screen parasitic driven electric fields.

We are able to conveniently address both of these issues by utilizing the
natural geometric separation of the electric and magnetic fields in a
$\lambda/4$ co-axial post-cavity. We engineer our drive antenna as a 3D cavity
that we call the ’buffer mode’, that is incorporated into the same monolithic
package that contains our high-Q storage cavities (LABEL:fig1b), thus
introducing minimal additional seam or dielectric loss. The base of this
buffer cavity forms a virtual electric ground where we simultaneously have an
E-field node and a B-field anti-node. Placing the SQUID loop in this region,
with its electric dipole moment oriented perpendicular to the direction of any
residual E-field, automatically smothers any stray capacitive coupling between
the buffer mode and the coupler.

Any additional effects due to a spatial non-uniformity of the magnetic field
threading the loop can be addressed through a purposely-engineered asymmetric
capacitance in the on-chip SQUID device (LABEL:HFSS_asymmetrya). This
asymmetry in the capacitance serves as a control knob for tuning the coupling
between the electromotive force and the coupler mode in the presence of an
alternating magnetic field, allowing to compensate for residual drive
asymmetry KochACFluxQuantization . To illustrate this effect, we analyze the
circuit model in LABEL:HFSS_asymmetryb that captures this effect to represent
the actual experimental device. We take into account a non-uniform AC magnetic
field that is distributed across the device, not only inducing a flux in the
SQUID loop, but also inducing an electromotive force on the shunting
capacitors. Assuming the geometric inductance of the loop is much smaller than
the Josephson inductance (Supplementary Note 1), we arrive at the effective
Hamiltonian of

$\hat{\mathcal{H}}=\frac{\left(\hat{Q}_{c}-C_{\Sigma}V_{\text{emf}}\right)^{2}}{2C_{\Sigma}}-E_{J}\cos{\frac{\Phi_{\text{ext}}^{\left(1,\Sigma\right)}}{2\phi_{0}}}\cos{\frac{\hat{\Phi}_{c}}{\phi_{0}}}.$
(19)

Here, $\hat{Q}_{c}$ is the charge operator of the common mode,
$C_{\Sigma}=C_{1}+C_{2}+C_{3}+C_{4}+C_{5}$ is its total shunting capacitance,
and
$\Phi_{\text{ext}}^{\left(1,\Sigma\right)}=\Phi_{\text{ext}}^{\left(1,l\right)}+\Phi_{\text{ext}}^{\left(1,r\right)}$
is the total external flux penetrating the SQUID loop. The charge operator
$\hat{Q}_{c}$ and the flux operator $\hat{\Phi}_{c}$ are related to the
Cooper-pair number operator $\hat{n}_{c}$ and the superconducting phase
operator $\hat{\theta}_{c}$ through $\hat{Q}_{c}=2e\hat{n}_{c}$,
$\hat{\Phi}_{c}=\phi_{0}\hat{\theta}_{c}$, where $e$ and $\phi_{0}$ are the
electron charge and the reduced flux quantum, respectively. The electromotive
voltage, $V_{\text{emf}}$, is given by

$\displaystyle V_{\text{emf}}$
$\displaystyle=\frac{C_{1}}{2C_{\Sigma}}\dot{\Phi}_{\text{ext}}^{\left(1,\delta\right)}+\frac{C_{l}}{C_{\Sigma}}\left(\frac{1}{2}\dot{\Phi}_{\text{ext}}^{\left(1,\Sigma\right)}+\dot{\Phi}_{\text{ext}}^{\left(2\right)}\right)$
$\displaystyle-\frac{C_{r}}{C_{\Sigma}}\left(\frac{1}{2}\dot{\Phi}_{\text{ext}}^{\left(1,\Sigma\right)}+\dot{\Phi}_{\text{ext}}^{\left(3\right)}\right),$
(20)

with
$\Phi_{\text{ext}}^{\left(1,\delta\right)}=\Phi_{\text{ext}}^{\left(1,l\right)}-\Phi_{\text{ext}}^{\left(1,r\right)}$,
$C_{l}=C_{2}+C_{4}$, $C_{r}=C_{3}+C_{5}$.

Shifting the upper antenna pad of the SQUID from the center to the right
simultaneously increases $C_{3}$ and decreases $C_{1}$ and $C_{2}$ in Eq. 20.
Fine-tuning this offset ($\delta$) can result in the minimization of
$V_{\text{emf}}$, as demonstrated in LABEL:HFSS_asymmetryc, where we simulate
the drive asymmetry $l$ (Supplementary Note 2) at various offsets in Ansys
finite-element High-Frequency Simulation Software (HFSS). We obtain this
asymmetry from the currents across the two junctions,

$l=\frac{\left|\phi_{c}\right|}{\sqrt{\left|\phi_{c}\right|^{2}+\left|\phi_{d}\right|^{2}}}=\frac{\left|I_{J1}+I_{J2}\right|}{\sqrt{2\left[\left|I_{J1}\right|^{2}+\left|I_{J2}\right|^{2}\right]}},$
(21)

with $I_{J1,2}=\phi_{0}\phi_{1,2}/L_{J1,2}$. When the pad is shifted towards
the right by $\delta\approx 350\mu$m, the asymmetry in the drive is almost
completely eliminated, resulting in the desired purely differential drive. At
the same time, the quality factor of the common mode (coupler) reaches a
maximum, since the Purcell effect from the coupler’s linear coupling to the
buffer mode is highly suppressed. In the experiment, this lets us strongly
couple the buffer mode to a coaxial drive line and achieve a quality factor
($\mathsf{Q}$) as low as 12 while maintaining a $\mathsf{Q}_{\text{c}}\sim
2\times 10^{6}$ for the coupler, and $\mathsf{Q}_{\text{a}}\sim 1.4\times
10^{7}$, $\mathsf{Q}_{\text{b}}\sim 1.2\times 10^{7}$ for Alice and Bob
respectively. We also make sure that other types of back-actions of the buffer
mode, such as the coherent- or thermal-photon-induced dephasing, is not
limiting the coherence of the coupler and the high-Q storage modes
(Supplementary Note 3).

### Experimentally characterizing residual drive asymmetry

Non-idealities in the experimental implementation of our differential drive
inevitably result in a small but finite common-mode drive strength. One way to
experimentally characterize this residual drive asymmetry is to probe the
strength of interactions that are uniquely generated by the common-mode drive.
A candidate for this is the ‘ge/3’ interaction, where a single drive tone at
$\approx\omega_{c}/3$ drives Rabi oscillations of the coupler between its
ground state and its excited state through three drive photons. The rate of
this Rabi-like process is given by (Supplementary Note 1):

$\Omega_{\text{ge/3}}=\frac{E_{J}}{\hbar}\;\theta_{c,\text{zpf}}\phi_{c}\left(|\phi_{c}|^{2}/3+|\phi_{d}|^{2}\right).$
(22)

Figure 6:  Experimentally characterizing drive asymmetry. We directly probe
the coupler ‘ge/3’ transition due to any residual common-mode drive. The
coupler is prepared in the ground state at $t=0\mu s$, followed by a single-
tone drive through the buffer-mode near $\omega_{c}/3$. The rate of coherent
oscillation of the coupler population (color-scale) compared to the effective
Zeeman shift (x-axis) of this resonance frequency bounds the drive asymmetry
to $<1\%$

fig]fig:Meth_ge3

To obtain the relative drive asymmetry, we compare the strength of the ge/3
process to the Zeeman-shift induced by the same single-tone drive, which
depends on both $\phi_{c}$ and $\phi_{d}$ :

$\displaystyle\Delta_{\text{Z},c}$
$\displaystyle\approx\left[J_{0}\left(|\phi_{c}|\right)J_{0}\left(|\phi_{d}|\right)-1\right]\omega_{c}/2$
$\displaystyle\approx-\frac{\omega_{c}}{8}\left(\lvert\phi_{c}\rvert^{2}+\lvert\phi_{d}\rvert^{2}\right).$
(23)

Note the difference between this and eq. 17 which assumes a purely
differential drive with two tones. Measuring the strength of both these
processes (eqs. 22 and 23) is sufficient to extract $|\phi_{c}|$ and
$|\phi_{d}|$. As illustrated in LABEL:fig:Meth_ge3, we observe the $ge$/3 Rabi
oscillation of the coupler with a rate of $\Omega_{\text{ge/3}}/2\pi=0.12$
MHz. The resonant frequency of the $ge$/3 oscillation is at $2.408$ GHz, which
corresponds to an effective Zeeman shift of $\Delta_{\text{Z},c}/2\pi=-21$MHz.
Using Eq. 22 and 23, we obtain $\lvert\frac{\phi_{c}}{\phi_{d}}\rvert\approx
2.5e-3$, corresponding to a drive asymmetry of $l<1\%$. It is possible that
the drive asymmetry is a little higher in the actual experiment, due to a
deviation from our assumption of perfectly symmetric junction energies, but we
expect this deviation to be small (Supplementary Note 7).

### Decoherence in the dual-rail subspace

We can derive the expected evolution and the beamsplitter fidelity from the
time evolution of the density matrix in the dual-rail subspace (Supplementary
Note 4):

$\displaystyle\rho_{\text{DR}}(t)$
$\displaystyle=\frac{1}{2}e^{-\kappa_{1}t}\;\times$
$\displaystyle\begin{bmatrix}1+e^{-\kappa_{\varphi}t}\cos\left(\Omega
t\right)&i\sin\left(\Omega t\right)e^{-\kappa_{\varphi}t}\\\
-i\sin\left(\Omega
t\right)e^{-\kappa_{\varphi}t}&1-e^{-\kappa_{\varphi}t}\cos\left(\Omega
t\right)\end{bmatrix},$ (24)

where we have used $\Omega=2g_{\text{BS}}$ for simplicity. The upper left
element measures the probability of finding a photon in Bob, which is Eq. 6 in
the main text.

To calculate the beamsplitter fidelity, we can set the evolution time to be
$t_{\text{BS}}=\frac{\pi}{4g_{\text{BS}}}$ and calculate the overlap with a
lossless evolution ($\rho_{|\psi^{\pm}\rangle}^{(0)}$) that has no leakage or
dephasing:

$\displaystyle\mathcal{F}(t_{\text{BS}})$
$\displaystyle=\text{Tr}\left[\rho_{\text{DR}}(t_{\text{BS}})\times\rho_{\text{DR}}^{(0)}(t_{\text{BS}})\right]$
$\displaystyle=\frac{1}{2}\left[e^{-\kappa_{1}t}+e^{-\left(\kappa_{1}+\kappa_{\varphi}\right)t}\right]$
$\displaystyle\approx 1-\kappa_{\text{BS}}t_{\text{BS}}.$ (25)

Here $\kappa_{\text{BS}}=\kappa_{1}+\frac{\kappa_{\varphi}}{2}$, which gives
us Eq. 7 in the main text.

Figure 7:  Measured driven decoherence. Decoherence rates are extracted from
short sections of the long-time evolution under both drive tones, as a
function of drive amplitude. The sideband collision of coupler and Bob (see
Supplementary Note 5) clearly limits the fidelity at high amplitudes, and we
operate on the boundary of this collision (grey dashed line), where
$\kappa_{\text{BS}}$ is limited by $\kappa_{1}$.

fig]fig:Meth_kappaBS

In order to efficiently extract $\kappa_{1}$ and $\kappa_{\phi}$ from the
experiment, we create a single photon in Bob, and measure Bob’s population
under the beamsplitter interaction within two sections of time:
$t\in[0,\>\delta t]$, and $t\in[T,T+\>\delta t]$. For small $\delta t$ and
large $T$, the measurement result within these time windows can be
approximated as

$\displaystyle P_{0}(t)$
$\displaystyle=A\left(1+\cos{(2g_{\text{BS}}t+\phi_{0})}\right)+B,$
$\displaystyle P_{T}(t)$
$\displaystyle=Ae^{-\kappa_{1}T}\left(1+e^{-\kappa_{\varphi}T}\cos{(2g_{\text{BS}}t+\phi_{0})}\right)+B.$
(26)

Here, $A$ and $A+B$ represent the amplitude and the offset of the oscillation,
which are close to but not strictly equal to 0.5 in the presence of state
preparation and measurement (SPAM) errors. From fitting to $P_{0}(t)$ we
extract $A$ and $B$, which are then used in the fitting of $P_{T}(t)$ to find
$\kappa_{1}$ and $\kappa_{\varphi}$. We choose $T=25\mu$s and $\delta t\approx
20t_{\text{BS}}$ to extract the decay rates at different drive amplitudes,
shown in LABEL:fig:Meth_kappaBS.

## Data availability

The data presented in this study is available at DOI:
https://doi.org/10.6084/m9.figshare.23589144.v1 and more detailed source data
is available from the corresponding authors upon request.

## Code availability

All computer codes used in this study are available from the corresponding
authors upon request.

## References

  *   * (1) Weedbrook, C., Pirandola, S., García-Patrón, R., Cerf, N.J., Ralph, T.C., Shapiro, J.H., Lloyd, S.: Gaussian quantum information. Rev. Mod. Phys. 84, 621–669 (2012). https://doi.org/10.1103/RevModPhys.84.621
  * (2) Ma, R., Saxberg, B., Owens, C., Leung, N., Lu, Y., Simon, J., Schuster, D.I.: A dissipatively stabilized mott insulator of photons. Nature 566(7742), 51–57 (2019). https://doi.org/10.1038/s41586-019-0897-9
  * (3) Hung, J.S.C., Busnaina, J.H., Chang, C.W.S., Vadiraj, A.M., Nsanzineza, I., Solano, E., Alaeian, H., Rico, E., Wilson, C.M.: Quantum simulation of the bosonic creutz ladder with a parametric cavity. Phys. Rev. Lett. 127, 100503 (2021). https://doi.org/10.1103/PhysRevLett.127.100503
  * (4) Schmidt, S., Koch, J.: Circuit qed lattices: Towards quantum simulation with superconducting circuits. Annalen der Physik 525(6), 395–412 (2013). https://doi.org/10.1002/andp.201200261
  * (5) Roushan, P., Neill, C., Megrant, A., Chen, Y., Babbush, R., Barends, R., Campbell, B., Chen, Z., Chiaro, B., Dunsworth, A., Fowler, A., Jeffrey, E., Kelly, J., Lucero, E., Mutus, J., O’Malley, P.J.J., Neeley, M., Quintana, C., Sank, D., Vainsencher, A., Wenner, J., White, T., Kapit, E., Neven, H., Martinis, J.: Chiral ground-state currents of interacting photons in a synthetic magnetic field. Nature Physics 13(2), 146–151 (2016). https://doi.org/10.1038/nphys3930
  * (6) Reagor, M., Paik, H., Catelani, G., Sun, L., Axline, C., Holland, E., Pop, I.M., Masluk, N.A., Brecht, T., Frunzio, L., Devoret, M.H., Glazman, L., Schoelkopf, R.J.: Reaching 10 ms single photon lifetimes for superconducting aluminum cavities. Applied Physics Letters 102(19), 192604 (2013). https://doi.org/10.1063/1.4807015
  * (7) Eickbusch, A., Sivak, V., Ding, A.Z., Elder, S.S., Jha, S.R., Venkatraman, J., Royer, B., Girvin, S.M., Schoelkopf, R.J., Devoret, M.H.: Fast universal control of an oscillator with weak dispersive coupling to a qubit. Nature Physics 18(12), 1464–1469 (2022). https://doi.org/10.1038/s41567-022-01776-9
  * (8) Milul, O., Guttel, B., Goldblatt, U., Hazanov, S., Joshi, L.M., Chausovsky, D., Kahn, N., Çiftyürek, E., Lafont, F., Rosenblum, S.: A superconducting quantum memory with tens of milliseconds coherence time (2023). https://arxiv.org/abs/2302.06442
  * (9) Chuang, I.L., Leung, D.W., Yamamoto, Y.: Bosonic quantum codes for amplitude damping. Phys. Rev. A 56, 1114–1125 (1997). https://doi.org/10.1103/PhysRevA.56.1114
  * (10) Cochrane, P.T., Milburn, G.J., Munro, W.J.: Macroscopically distinct quantum-superposition states as a bosonic code for amplitude damping. Phys. Rev. A 59, 2631–2634 (1999). https://doi.org/10.1103/PhysRevA.59.2631
  * (11) Leghtas, Z., Kirchmair, G., Vlastakis, B., Schoelkopf, R.J., Devoret, M.H., Mirrahimi, M.: Hardware-efficient autonomous quantum memory protection. Phys. Rev. Lett. 111, 120501 (2013). https://doi.org/10.1103/PhysRevLett.111.120501
  * (12) Teoh, J.D., Winkel, P., Babla, H.K., Chapman, B.J., Claes, J., de Graaf, S.J., Garmon, J.W.O., Kalfus, W.D., Lu, Y., Maiti, A., Sahay, K., Thakur, N., Tsunoda, T., Xue, S.H., Frunzio, L., Girvin, S.M., Puri, S., Schoelkopf, R.J.: Dual-rail encoding with superconducting cavities (2022). https://arxiv.org/abs/2212.12077
  * (13) Ofek, N., Petrenko, A., Heeres, R., Reinhold, P., Leghtas, Z., Vlastakis, B., Liu, Y., Frunzio, L., Girvin, S.M., Jiang, L., Mirrahimi, M., Devoret, M.H., Schoelkopf, R.J.: Extending the lifetime of a quantum bit with error correction in superconducting circuits. Nature 536(7617), 441–445 (2016). https://doi.org/10.1038/nature18949
  * (14) Sivak, V.V., Eickbusch, A., Royer, B., Singh, S., Tsioutsios, I., Ganjam, S., Miano, A., Brock, B.L., Ding, A.Z., Frunzio, L., Girvin, S.M., Schoelkopf, R.J., Devoret, M.H.: Real-time quantum error correction beyond break-even. Nature 616(7955), 50–55 (2023). https://doi.org/10.1038/s41586-023-05782-6
  * (15) Ni, Z., Li, S., Deng, X., Cai, Y., Zhang, L., Wang, W., Yang, Z.-B., Yu, H., Yan, F., Liu, S., Zou, C.-L., Sun, L., Zheng, S.-B., Xu, Y., Yu, D.: Beating the break-even point with a discrete-variable-encoded logical qubit. Nature 616(7955), 56–60 (2023). https://doi.org/10.1038/s41586-023-05784-4
  * (16) Gao, Y.Y., Lester, B.J., Chou, K.S., Frunzio, L., Devoret, M.H., Jiang, L., Girvin, S.M., Schoelkopf, R.J.: Entanglement of bosonic modes through an engineered exchange interaction. Nature 566(7745), 509–512 (2019). https://doi.org/10.1038/s41586-019-0970-4
  * (17) Tsunoda, T., Teoh, J.D., Kalfus, W.D., de Graaf, S.J., Chapman, B.J., Curtis, J.C., Thakur, N., Girvin, S.M., Schoelkopf, R.J.: Error-detectable bosonic entangling gates with a noisy ancilla. PRX Quantum 4, 020354 (2023). https://doi.org/10.1103/PRXQuantum.4.020354
  * (18) Joshi, A., Noh, K., Gao, Y.Y.: Quantum information processing with bosonic qubits in circuit QED. Quantum Science and Technology 6(3), 033001 (2021). https://doi.org/10.1088/2058-9565/abe989
  * (19) Campagne-Ibarcq, P., Zalys-Geller, E., Narla, A., Shankar, S., Reinhold, P., Burkhart, L., Axline, C., Pfaff, W., Frunzio, L., Schoelkopf, R.J., Devoret, M.H.: Deterministic remote entanglement of superconducting circuits through microwave two-photon transitions. Phys. Rev. Lett. 120, 200501 (2018). https://doi.org/10.1103/PhysRevLett.120.200501
  * (20) Axline, C.J., Burkhart, L.D., Pfaff, W., Zhang, M., Chou, K., Campagne-Ibarcq, P., Reinhold, P., Frunzio, L., Girvin, S.M., Jiang, L., Devoret, M.H., Schoelkopf, R.J.: On-demand quantum state transfer and entanglement between remote microwave cavity memories. Nature Physics 14(7), 705–710 (2018). https://doi.org/10.1038/s41567-018-0115-y
  * (21) Bergeal, N., Vijay, R., Manucharyan, V.E., Siddiqi, I., Schoelkopf, R.J., Girvin, S.M., Devoret, M.H.: Analog information processing at the quantum limit with a josephson ring modulator. Nature Physics 6(4), 296–302 (2010). https://doi.org/10.1038/nphys1516
  * (22) Frattini, N.E., Vool, U., Shankar, S., Narla, A., Sliwa, K.M., Devoret, M.H.: 3-wave mixing josephson dipole element. Applied Physics Letters 110(22), 222603 (2017). https://doi.org/10.1063/1.4984142
  * (23) Sivak, V.V., Shankar, S., Liu, G., Aumentado, J., Devoret, M.H.: Josephson array-mode parametric amplifier. Phys. Rev. Applied 13, 024014 (2020). https://doi.org/10.1103/PhysRevApplied.13.024014
  * (24) Zakka-Bajjani, E., Nguyen, F., Lee, M., Vale, L.R., Simmonds, R.W., Aumentado, J.: Quantum superposition of a single microwave photon in two different ’colour’states. Nature Physics 7(8), 599–603 (2011). https://doi.org/10.1038/nphys2035
  * (25) Sirois, A.J., Castellanos-Beltran, M.A., DeFeo, M.P., Ranzani, L., Lecocq, F., Simmonds, R.W., Teufel, J.D., Aumentado, J.: Coherent-state storage and retrieval between superconducting cavities using parametric frequency conversion. Applied Physics Letters 106(17), 172603 (2015). https://doi.org/10.1063/1.4919759
  * (26) Gao, Y.Y., Lester, B.J., Zhang, Y., Wang, C., Rosenblum, S., Frunzio, L., Jiang, L., Girvin, S.M., Schoelkopf, R.J.: Programmable interference between two microwave quantum memories. Phys. Rev. X 8, 021073 (2018). https://doi.org/10.1103/PhysRevX.8.021073
  * (27) Lecocq, F., Ranzani, L., Peterson, G.A., Cicak, K., Simmonds, R.W., Teufel, J.D., Aumentado, J.: Nonreciprocal microwave signal processing with a field-programmable josephson amplifier. Phys. Rev. Applied 7, 024028 (2017). https://doi.org/10.1103/PhysRevApplied.7.024028
  * (28) Noh, T., Xiao, Z., Jin, X.Y., Cicak, K., Doucet, E., Aumentado, J., Govia, L.C.G., Ranzani, L., Kamal, A., Simmonds, R.W.: Strong parametric dispersive shifts in a statically decoupled two-qubit cavity qed system. Nature Physics (2023). https://doi.org/%****␣SQUID_beamsplitter_v2.bbl␣Line␣625␣****10.1038/s41567-023-02107-2
  * (29) Zhou, C., Lu, P., Praquin, M., Chien, T.-C., Kaufman, R., Cao, X., Xia, M., Mong, R.S.K., Pfaff, W., Pekker, D., Hatridge, M.: Realizing all-to-all couplings among detachable quantum modules using a microwave quantum state router. npj Quantum Information 9(1), 54 (2023). https://doi.org/10.1038/s41534-023-00723-7
  * (30) Sank, D., Chen, Z., Khezri, M., Kelly, J., Barends, R., Campbell, B., Chen, Y., Chiaro, B., Dunsworth, A., Fowler, A., Jeffrey, E., Lucero, E., Megrant, A., Mutus, J., Neeley, M., Neill, C., O’Malley, P.J.J., Quintana, C., Roushan, P., Vainsencher, A., White, T., Wenner, J., Korotkov, A.N., Martinis, J.M.: Measurement-induced state transitions in a superconducting qubit: Beyond the rotating wave approximation. Phys. Rev. Lett. 117, 190503 (2016). https://doi.org/10.1103/PhysRevLett.117.190503
  * (31) Zhang, Y., Lester, B.J., Gao, Y.Y., Jiang, L., Schoelkopf, R.J., Girvin, S.M.: Engineering bilinear mode coupling in circuit qed: Theory and experiment. Phys. Rev. A 99, 012314 (2019). https://doi.org/10.1103/PhysRevA.99.012314
  * (32) Lescanne, R., Villiers, M., Peronnin, T., Sarlette, A., Delbecq, M., Huard, B., Kontos, T., Mirrahimi, M., Leghtas, Z.: Exponential suppression of bit-flips in a qubit encoded in an oscillator. Nature Physics 16(5), 509–513 (2020). https://doi.org/10.1038/s41567-020-0824-x
  * (33) Koch, J., Yu, T.M., Gambetta, J., Houck, A.A., Schuster, D.I., Majer, J., Blais, A., Devoret, M.H., Girvin, S.M., Schoelkopf, R.J.: Charge-insensitive qubit design derived from the cooper pair box. Phys. Rev. A 76, 042319 (2007). https://doi.org/10.1103/PhysRevA.76.042319
  * (34) Grimm, A., Frattini, N.E., Puri, S., Mundhada, S.O., Touzard, S., Mirrahimi, M., Girvin, S.M., Shankar, S., Devoret, M.H.: Stabilization and operation of a kerr-cat qubit. Nature 584(7820), 205–209 (2020). https://doi.org/10.1038/s41586-020-2587-z
  * (35) Puri, S., St-Jean, L., Gross, J.A., Grimm, A., Frattini, N.E., Iyer, P.S., Krishna, A., Touzard, S., Jiang, L., Blais, A., Flammia, S.T., Girvin, S.M.: Bias-preserving gates with stabilized cat qubits. Science Advances 6(34) (2020). https://doi.org/10.1126/sciadv.aay5901
  * (36) Chirolli, L., Burkard, G., Kumar, S., DiVincenzo, D.P.: Superconducting resonators as beam splitters for linear-optics quantum computation. Phys. Rev. Lett. 104, 230502 (2010). https://doi.org/10.1103/PhysRevLett.104.230502
  * (37) Lecocq, F., Claudon, J., Buisson, O., Milman, P.: Nonlinear coupling between the two oscillation modes of a dc squid. Phys. Rev. Lett. 107, 197002 (2011). https://doi.org/10.1103/PhysRevLett.107.197002
  * (38) Kamal, A., Clarke, J., Devoret, M.H.: Gain, directionality, and noise in microwave squid amplifiers: Input-output approach. Phys. Rev. B 86, 144510 (2012). https://doi.org/10.1103/PhysRevB.86.144510
  * (39) Yamamoto, T., Inomata, K., Watanabe, M., Matsuba, K., Miyazaki, T., Oliver, W.D., Nakamura, Y., Tsai, J.S.: Flux-driven josephson parametric amplifier. Applied Physics Letters 93(4), 042510 (2008). https://doi.org/10.1063/1.2964182
  * (40) Strand, J.D., Ware, M., Beaudoin, F., Ohki, T.A., Johnson, B.R., Blais, A., Plourde, B.L.T.: First-order sideband transitions with flux-driven asymmetric transmon qubits. Phys. Rev. B 87, 220505 (2013). https://doi.org/%****␣SQUID_beamsplitter_v2.bbl␣Line␣875␣****10.1103/PhysRevB.87.220505
  * (41) McKay, D.C., Naik, R., Reinhold, P., Bishop, L.S., Schuster, D.I.: High-contrast qubit interactions using multimode cavity qed. Phys. Rev. Lett. 114, 080501 (2015). https://doi.org/10.1103/PhysRevLett.114.080501
  * (42) Lu, Y., Chakram, S., Leung, N., Earnest, N., Naik, R.K., Huang, Z., Groszkowski, P., Kapit, E., Koch, J., Schuster, D.I.: Universal stabilization of a parametrically coupled qubit. Phys. Rev. Lett. 119, 150502 (2017). https://doi.org/10.1103/PhysRevLett.119.150502
  * (43) Dykman, M.I., Marthaler, M., Peano, V.: Quantum heating of a parametrically modulated oscillator: Spectral signatures. Phys. Rev. A 83, 052115 (2011). https://doi.org/10.1103/PhysRevA.83.052115
  * (44) Wang, C., Gao, Y.Y., Reinhold, P., Heeres, R.W., Ofek, N., Chou, K., Axline, C., Reagor, M., Blumoff, J., Sliwa, K.M., Frunzio, L., Girvin, S.M., Jiang, L., Mirrahimi, M., Devoret, M.H., Schoelkopf, R.J.: A schrödinger cat living in two boxes. Science 352(6289), 1087–1091 (2016). https://doi.org/10.1126/science.aaf2941
  * (45) Nigg, S.E., Paik, H., Vlastakis, B., Kirchmair, G., Shankar, S., Frunzio, L., Devoret, M.H., Schoelkopf, R.J., Girvin, S.M.: Black-box superconducting circuit quantization. Phys. Rev. Lett. 108, 240502 (2012). https://doi.org/10.1103/PhysRevLett.108.240502
  * (46) Minev, Z.K., Leghtas, Z., Mundhada, S.O., Christakis, L., Pop, I.M., Devoret, M.H.: Energy-participation quantization of josephson circuits. npj Quantum Information 7(1), 131 (2021). https://doi.org/10.1038/s41534-021-00461-8
  * (47) Krastanov, S., Albert, V.V., Shen, C., Zou, C.-L., Heeres, R.W., Vlastakis, B., Schoelkopf, R.J., Jiang, L.: Universal control of an oscillator with dispersive coupling to a qubit. Phys. Rev. A 92, 040303 (2015). https://doi.org/10.1103/PhysRevA.92.040303
  * (48) Geerlings, K., Leghtas, Z., Pop, I.M., Shankar, S., Frunzio, L., Schoelkopf, R.J., Mirrahimi, M., Devoret, M.H.: Demonstrating a driven reset protocol for a superconducting qubit. Phys. Rev. Lett. 110, 120501 (2013). https://doi.org/10.1103/PhysRevLett.110.120501
  * (49) Zhang, Y., Curtis, J.C., Wang, C.S., Schoelkopf, R.J., Girvin, S.M.: Drive-induced nonlinearities of cavity modes coupled to a transmon ancilla. Phys. Rev. A 105, 022423 (2022). https://doi.org/10.1103/PhysRevA.105.022423
  * (50) Chapman, B.J., de Graaf, S.J., Xue, S.H., Zhang, Y., Teoh, J., Curtis, J.C., Tsunoda, T., Eickbusch, A., Read, A.P., Koottandavida, A., Mundhada, S.O., Frunzio, L., Devoret, M.H., Girvin, S.M., Schoelkopf, R.J.: High-on-off-ratio beam-splitter interaction for gates on bosonically encoded qubits. PRX Quantum 4, 020355 (2023). https://doi.org/10.1103/PRXQuantum.4.020355
  * (51) Knill, E., Leibfried, D., Reichle, R., Britton, J., Blakestad, R.B., Jost, J.D., Langer, C., Ozeri, R., Seidelin, S., Wineland, D.J.: Randomized benchmarking of quantum gates. Phys. Rev. A 77, 012307 (2008). https://doi.org/10.1103/PhysRevA.77.012307
  * (52) Emerson, J., Alicki, R., Życzkowski, K.: Scalable noise estimation with random unitary operators. Journal of Optics B: Quantum and Semiclassical Optics 7(10), 347 (2005). https://doi.org/10.1088/1464-4266/7/10/021
  * (53) Boone, K., Carignan-Dugas, A., Wallman, J.J., Emerson, J.: Randomized benchmarking under different gate sets. Phys. Rev. A 99, 032329 (2019). https://doi.org/10.1103/PhysRevA.99.032329
  * (54) Proctor, T.J., Carignan-Dugas, A., Rudinger, K., Nielsen, E., Blume-Kohout, R., Young, K.: Direct randomized benchmarking for multiqubit devices. Phys. Rev. Lett. 123, 030503 (2019). https://doi.org/10.1103/PhysRevLett.123.030503
  * (55) Wu, Y., Kolkowitz, S., Puri, S., Thompson, J.D.: Erasure conversion for fault-tolerant quantum computing in alkaline earth rydberg atom arrays. Nature Communications 13(1), 4657 (2022). https://doi.org/10.1038/s41467-022-32094-6
  * (56) Krantz, P., Kjaergaard, M., Yan, F., Orlando, T.P., Gustavsson, S., Oliver, W.D.: A quantum engineer's guide to superconducting qubits. Applied Physics Reviews 6(2), 021318 (2019). https://doi.org/10.1063/1.5089550
  * (57) Blais, A., Grimsmo, A.L., Girvin, S.M., Wallraff, A.: Circuit quantum electrodynamics. Reviews of Modern Physics 93(2) (2021). https://doi.org/%****␣SQUID_beamsplitter_v2.bbl␣Line␣1200␣****10.1103/revmodphys.93.025005
  * (58) You, X., Sauls, J.A., Koch, J.: Circuit quantization in the presence of time-dependent external flux. Phys. Rev. B 99, 174512 (2019). https://doi.org/10.1103/PhysRevB.99.174512
  * (59) Kohler, S., Dittrich, T., Hänggi, P.: Floquet-markovian description of the parametrically driven, dissipative harmonic quantum oscillator. Phys. Rev. E 55, 300–313 (1997). https://doi.org/10.1103/PhysRevE.55.300
  * (60) Verney, L., Lescanne, R., Devoret, M.H., Leghtas, Z., Mirrahimi, M.: Structural instability of driven josephson circuits prevented by an inductive shunt. Phys. Rev. Appl. 11, 024003 (2019). https://doi.org/10.1103/PhysRevApplied.11.024003
  * (61) Shillito, R., Petrescu, A., Cohen, J., Beall, J., Hauru, M., Ganahl, M., Lewis, A.G.M., Vidal, G., Blais, A.: Dynamics of transmon ionization. Phys. Rev. Appl. 18, 034031 (2022). https://doi.org/10.1103/PhysRevApplied.18.034031
  * (62) Yan, F., Campbell, D., Krantz, P., Kjaergaard, M., Kim, D., Yoder, J.L., Hover, D., Sears, A., Kerman, A.J., Orlando, T.P., Gustavsson, S., Oliver, W.D.: Distinguishing coherent and thermal photon noise in a circuit quantum electrodynamical system. Phys. Rev. Lett. 120, 260504 (2018). https://doi.org/10.1103/PhysRevLett.120.260504
  * (63) Clerk, A.A., Devoret, M.H., Girvin, S.M., Marquardt, F., Schoelkopf, R.J.: Introduction to quantum noise, measurement, and amplification. Rev. Mod. Phys. 82, 1155–1208 (2010). https://doi.org/10.1103/RevModPhys.82.1155
  * (64) Heeres, R.W., Vlastakis, B., Holland, E., Krastanov, S., Albert, V.V., Frunzio, L., Jiang, L., Schoelkopf, R.J.: Cavity state manipulation using photon-number selective phase gates. Phys. Rev. Lett. 115, 137002 (2015). https://doi.org/10.1103/PhysRevLett.115.137002
  * (65) Zhang, Y., Curtis, J.C., Wang, C.S., Schoelkopf, R.J., Girvin, S.M.: Drive-induced nonlinearities of cavity modes coupled to a transmon ancilla. Phys. Rev. A 105, 022423 (2022). https://doi.org/10.1103/PhysRevA.105.022423
  * (66) Johansson, J.R., Nation, P.D., Nori, F.: QuTiP: An open-source python framework for the dynamics of open quantum systems. Computer Physics Communications 183(8), 1760–1772 (2012). https://doi.org/10.1016/j.cpc.2012.02.021
  * (67) Campagne-Ibarcq, P., Eickbusch, A., Touzard, S., Zalys-Geller, E., Frattini, N.E., Sivak, V.V., Reinhold, P., Puri, S., Shankar, S., Schoelkopf, R.J., Frunzio, L., Mirrahimi, M., Devoret, M.H.: Quantum error correction of a qubit encoded in grid states of an oscillator. Nature 584(7821), 368–372 (2020). https://doi.org/10.1038/s41586-020-2603-3
  * (68) Gertler, J.M., Baker, B., Li, J., Shirol, S., Koch, J., Wang, C.: Protecting a bosonic qubit with autonomous quantum error correction. Nature 590(7845), 243–248 (2021). https://doi.org/10.1038/s41586-021-03257-0
  * (69) Ambegaokar, V., Baratoff, A.: Tunneling between superconductors. Phys. Rev. Lett. 10, 486–489 (1963). https://doi.org/10.1103/PhysRevLett.10.486
  * (70) Wood, C.J., Gambetta, J.M.: Quantification and characterization of leakage errors. Phys. Rev. A 97, 032306 (2018). https://doi.org/10.1103/PhysRevA.97.032306

## Acknowledgments

We thank J. Venkatraman and and X. Xiao for helpful conversations on parasitic
nonlinear resonances. We are grateful to J. Curtis, A. Koottandavida and I.
Tsioutsios for technical assistance and providing useful code. We thank B.
Chapman for providing parametric amplifiers used in this experiment, and A.
Read for help with DC-flux line wiring. We are grateful to C. Wehr and H. M.
Moseley for help with designing substrate clamps and K. Chou for advice on
beamsplitter control. We thank A. Miano for help simulating the dc-flux
spectroscopy of the SQUID. We also thank P. Winkel, M. Devoret, B. Chapman, S.
de Graaf, S. Xue, J. D. Teoh, T. Tsunoda, S. Chakram and J. Huang for useful
discussions. This research was sponsored by the Army Research Office (ARO)
under grant nos. W911NF-16-1-0349, W911NF-18-1-0212 and W911NF-22-1-0053. The
views and conclusions contained in this document are those of the authors and
should not be interpreted as representing the official policies, either
expressed or implied, of the Army Research Office (ARO) or the US Government.
Fabrication facilities use was supported by the Yale Institute for Nanoscience
and Quantum Engineering (YINQE) and the Yale SEAS Cleanroom.

## Author contributions

Y.L., A.M. designed and conducted the experiment, under the supervision of
R.J.S.. Y.L., A.M. and S.G. designed the experimental device and developed the
differential drive scheme. Y.L. fabricated the DC-SQUID and ancilla transmon
chip devices with help from S.G. and L.F.. A.M., Y. L. and J.G. conducted
measurements and analyzed the data. J.G., A.M. and Y.L. implemented the dual-
rail qubit randomized benchmarking. Y.Z. and J.C. provided theoretical support
for analyzing parametric processes and randomized benchmarking, respectively.
S.M.G. provided important theoretical insights during the writing of the
manuscript. A.M., Y.L., J.G. and R.J.S. wrote the manuscript with feedback
from all co-authors.

## Competing Interests

L.F. and R.J.S. are founders and shareholders of Quantum Circuits, Inc (QCI).
S.M.G. receives consulting fees and is an equity holder in QCI. The remaining
authors declare no competing interests.

Supplementary Information

##  Supplementary Note 1 — Hamiltonian of the driven DC SQUID

The two-mode nature of a DC-SQUID has been well-studied in literature
SQUIDModesLecocq2011 ; KamalSQUIDModes . In this section, we derive the joint
Hamiltonian for these two modes in the presence of an AC-flux modulation. The
SQUID circuit is shown in Fig. S8, where we include a small but finite linear
inductance $L_{3}$ and capacitance $C_{3}$ that connects the two junctions.
The flux variables across the junctions and the linear inductance are labeled
as $\Phi_{1,2,3}$, where the first two directly correspond to the Josephson
phases $\hat{\theta}_{1,2}$ in the main text. From this general formalism, we
will derive the parity-protected Hamiltonian as a special case, with
$\hat{\theta}_{c}$ and $\phi_{d}$ representing the coupler and actuator
respectively. We will also derive the other limit of the common-driven SQUID
and explicitly establish the Hamiltonian under a continuous spectrum of drive
asymmetry.

Figure S8: The circuit diagram of a flux-driven SQUID.

We start with the general circuit Lagrangian for the driven SQUID in Fig. S8

$\displaystyle\mathscr{L}=\frac{1}{2}\sum_{i=1}^{3}C_{i}\dot{\Phi}_{i}^{2}+\sum_{i=1}^{2}E_{Ji}\cos{\frac{\Phi_{i}}{\phi_{0}}}-\frac{\Phi_{3}^{2}}{2L_{3}},$
(27)

with a constraint from the flux quantization relationship,

$\Phi_{1}-\Phi_{2}+\Phi_{3}=\Phi_{\text{ext}}+n\Phi_{0},$ (28)

where $\Phi_{0}\;(\phi_{0})$ is the (reduced) flux quantum. Therefore, we may
rewrite Eq. 27 as

$\displaystyle\mathscr{L}$
$\displaystyle=\frac{1}{2}C_{1}\left(\dot{\Phi}_{c}+\dot{\Phi}_{d}\right)^{2}+\frac{1}{2}C_{2}\left(\dot{\Phi}_{c}-\dot{\Phi}_{d}\right)^{2}$
$\displaystyle+E_{J_{1}}\cos{\frac{\Phi_{c}+\Phi_{d}}{\phi_{0}}}+E_{J_{2}}\cos{\frac{\Phi_{c}-\Phi_{d}}{\phi_{0}}}$
$\displaystyle+\frac{1}{2}C_{3}\left(2\dot{\Phi}_{d}-\dot{\Phi}_{\text{ext}}\right)^{2}-\frac{\left(2\Phi_{d}-\Phi_{\text{ext}}\right)^{2}}{2L_{3}}.$
(29)

Here, we have defined the “differential mode” $\Phi_{d}$, and the “common
mode” $\Phi_{c}$, as

$\Phi_{c,d}=\frac{1}{2}\left(\Phi_{1}\pm\Phi_{2}\right).$ (30)

Under the assumption that the linear inductance $L_{J3}$ and the capacitance
$C_{J3}$ are small, the differential mode frequency is much higher than that
of the common mode. This means that, even under circuit asymmetry such that
$C_{1}\neq C_{2}$ and $E_{J_{1}}\neq E_{J_{2}}$, $\Phi_{c}$ and $\Phi_{d}$
remain good approximations for the normal-mode flux variables of the circuit.
From this, the equation of motion for $\Phi_{d}$ can be approximated as

$\displaystyle\frac{d}{dt}\frac{\partial\mathscr{L}}{\partial\dot{\Phi}_{d}}-\frac{\partial\mathscr{L}}{\partial\Phi_{d}}\approx$
$\displaystyle
C_{d\Sigma}\ddot{\Phi}_{d}+L_{d\Sigma}^{-1}\Phi_{d}-2C_{3}\ddot{\Phi}_{ext}-2\frac{\Phi_{\text{ext}}}{L_{3}}=0,$
(31)

where $C_{d\Sigma}=C_{1}+C_{2}+4C_{3}$ and
$L_{d\Sigma}=\left(L_{J1}^{-1}+L_{J2}^{-1}+4L_{3}^{-1}\right)^{-1}$ are the
total capacitance and inductance of the differential mode. When the drive
frequency is much lower than the differential mode frequency, the steady state
solution of Eq. 31 is given by

$\left<\Phi_{d}\right>\approx\Phi_{\text{ext}}/2.$ (32)

This can also be seen from the fact that in Eq. 29, the last term is a high
energy term that requires its numerator to be a “frozen” degree of freedom,
$\left<2\dot{\Phi}_{d}-\dot{\Phi}_{ext}\right>$=0. Plugging this back in to
Eq. 29 yields

$\displaystyle\mathscr{L}$
$\displaystyle=\frac{1}{2}C_{1}\left(\dot{\Phi}_{c}+\frac{\dot{\Phi}_{ext}}{2}\right)^{2}+\frac{1}{2}C_{2}\left(\dot{\Phi}_{c}-\frac{\dot{\Phi}_{ext}}{2}\right)^{2}$
$\displaystyle+E_{J_{1}}\cos{\frac{\Phi_{c}+\frac{\Phi_{\text{ext}}}{2}}{\phi_{0}}}+E_{J_{2}}\cos{\frac{\Phi_{c}-\frac{\Phi_{\text{ext}}}{2}}{\phi_{0}}}.$
(33)

In reality, the differential flux modulation may also be contaminated by stray
voltage coupling of the common mode to the flux line, as well as by the non-
uniformity in the spatial distribution of the time-dependent flux that creates
uneven voltage drops (electromotive forces) across the shunting capacitors
(see Methods “Fine-tuning a differential drive” for a detailed analysis).
Without loss of generality, the Hamiltonian can be written as

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=4E_{C}\hat{n}_{c}^{2}+\epsilon\left(t\right)\hat{n}_{c}$
$\displaystyle-$ $\displaystyle
E_{J_{1}}\cos{\left(\hat{\theta}_{c}+\frac{\Phi_{\text{ext}}}{2\phi_{0}}\right)}-E_{J_{2}}\cos{\left(\hat{\theta}_{c}-\frac{\Phi_{\text{ext}}}{2\phi_{0}}\right)},$
(34)

where $\epsilon\left(t\right)$ represents the effective common mode drive
strength, and $\hat{n}_{c}$, $\hat{\theta}_{c}$ are the Cooper-pair number
operator and the superconducting phase operator of the common mode,
respectively. In the displaced frame, this Hamiltonian becomes (for symmetric
SQUID of $E_{J_{1}}=E_{J_{2}}=E_{J}/2$)

$\displaystyle\hat{\mathcal{H}}_{\text{disp}}$
$\displaystyle=4E_{C}\hat{n}_{c}^{2}+E_{J}\frac{\hat{\theta}_{c}^{2}}{2}$
$\displaystyle-
E_{J}\cos_{\text{NL}}{\left(\hat{\theta}_{c}+\phi_{1}\right)}-E_{J}\cos_{\text{NL}}{\left(\hat{\theta}_{c}+\phi_{2}\right)}$
$\displaystyle=4E_{C}\hat{n}_{c}^{2}-E_{J}\cos\phi_{c}\cos\phi_{d}\cos\hat{\theta}_{c}$
$\displaystyle+E_{J}\sin\phi_{c}\cos\phi_{d}\sin\hat{\theta}_{c}-E_{J}\phi_{c}\hat{\theta}_{c},$
(35)

where $\cos_{\textup{NL}}x=\cos{x}+\frac{x^{2}}{2}$, and
$\phi_{c,d}=\frac{1}{2}\left(\phi_{1}\pm\phi_{2}\right)$ are the displacements
in the phases of the common and the differential mode. From eq. 34 and eq. 35,
we obtain these displacements as $\phi_{d}=\Phi_{\text{ext}}/2\phi_{0}$, and
$\phi_{c}=\epsilon(t)\omega_{d}\hbar^{-1}/\left(\omega_{d}^{2}-\omega_{c}^{2}\right)$,
where $\omega_{c}$ and $\omega_{d}$ are the frequencies of the coupler and the
drive. A pure differential drive corresponds to $\phi_{c}=0$, where Eq. 35
reduces to

$\hat{\mathcal{H}}_{DDS}=4E_{C}\hat{n}_{c}^{2}-E_{J}\cos\phi_{d}\cos{\hat{\theta}_{c}},$
(36)

which is the parity-protected Hamiltonian (Eq. 1) in the main text. In
contrast, $\phi_{d}=0$ happens for a purely charge-driven transmon,

$\displaystyle\hat{\mathcal{H}}_{trans}$
$\displaystyle=4E_{C}\hat{n}_{c}^{2}-E_{J}\cos\phi_{c}\cos\hat{\theta_{c}}$
$\displaystyle+E_{J}\sin\phi_{c}\sin\hat{\theta}_{c}-E_{J}\phi_{c}\hat{\theta}_{c}.$
(37)

With Eq. 35, we can express the strengths of some of the parametric processes
in terms of the displacement amplitudes, $\phi_{c}$ and $\phi_{d}$. For
example, the beamsplitter interaction arises from the $\cos\hat{\theta}_{c}$
modulation in the third line of Eq. 35. For purely RF drives where $\phi_{d}$
contains no DC component, both the beamsplitter rate and the Zeeman shift are
proportional to $\phi_{c}^{2}+\phi_{d}^{2}$ under weak drive strengths. On the
other hand, the last line of Eq. 35 gives rise to the odd parity interactions,
with the leading order term being the ‘$ge/3$’ interaction,

$\displaystyle\hat{\mathcal{H}}_{ge/3}$ $\displaystyle\approx
E_{J}\left(\phi_{c}-\frac{\phi_{c}^{3}}{6}\right)\left(1-\frac{\phi_{d}^{2}}{2}\right)\hat{\theta}_{c}-E_{J}\phi_{c}\hat{\theta}_{c}$
$\displaystyle\approx-\frac{E_{J}}{2}\phi_{c}\left(\phi_{d}^{2}+\frac{\phi_{c}^{2}}{3}\right)\theta_{c,\text{zpf}}\left(\hat{c}^{\dagger}+\hat{c}\right),$
(38)

which we used to experimentally characterize the residual drive asymmetry
(Methods “Experimentally characterizing residual drive asymmetry”).

##  Supplementary Note 2 — Parity protection and drive-frequency engineering

In this section, we discuss the importance of parity protection and drive-
frequency engineering in the differentially-driven-SQUID (DDS), by contrasting
it to a charge-driven transmon. In the latter’s case (Eq. 37), only the ‘even
parity’ part of the Hamiltonian ($\hat{\theta}_{c}^{2n}$, coming from the
second term) that allows even numbers of mode quanta to interact, is useful
for the bilinear beamsplitting process ($\propto\hat{\theta}_{c}^{2}$), The
odd-parity terms ($\hat{\theta}_{c}^{2n+1}$, from the second line) introduce a
set of nonlinear processes that are not only unhelpful, but can be actively
harmful, leading to an incoherent excitation of the coupler that dephases the
beamsplitter process, or an unintended exchange of sensitive quantum
information with other modes in the system or the environment.

We are able to address both of these issues by utilizing the orthogonality of
the coupler and actuator in the differentially-driven SQUID (DDS, Eq. 36),
which offers a two-fold advantage through its parity protection and simpler
drive engineering.

First, the forbidden linear coupling between the drive and the coupler avoids
its displacement that converts dephasing into heating. This can be intuitively
understood as high-frequency coupler dephasing combining with the drive to
provide excitation at the coupler frequency. This process should be sensitive
only to the spectrum of coupler frequency noise at the detuning between the
drive and the coupler mode, and thus driving far-detuned from the coupler is
also helpful.

Second, incoherent excitations may also result from driven resonant processes
in combination with coupler decay. It is therefore imperative to avoid these
resonances, the strongest of which is a squeezing of the coupler that is
resonant near the coupler mode frequency (specifically, at half the g-f
transition frequency). Intuitively, this can be thought of as a coherent
excitation from the ground to the second excited state, combined with an
incoherent decay back to the first excited state, thus leading to an over-all
incoherent heating. This is a common issue in both the charge-driven transmon
as well as the DDS, since squeezing and higher-order even-parity interactions
are not forbidden by Eq. 36. However, in the DDS, the orthogonality between
the actuator and the coupler greatly facilitates the engineering of the drives
that are far-detuned from the coupler, suppressing both the squeezing and the
dressed dephasing.

Finally, the selection rule enforced by the parity-protected Hamiltonian
should forbid any odd-parity exchanges for modes participating in the couple,
which is especially useful since strong drives can cause frequency shifts of
the resonance condition. This includes resonances that could either cause
undesired exchanges between the cavity ($\hat{a}$) and the coupler
($\hat{c}$), like $\phi_{d}\;\hat{a}\hat{c}^{\dagger 2}+\text{h.c.}$ which was
a limiting factor in Yvonne2018_PRX_ProgrammableInterference , or resonances
that cause direct excitations of the coupler through terms like
$\phi_{d}^{3}\;\hat{c}^{\dagger}+\text{h.c.}$.

Which of these effects dominate in the experiment depends strongly on the
frequencies of interest and the spectrum of environmental noise that couples
to the system. Here, we take a broader approach by using numerical simulations
to illustrate the clear difference between a common and a differential drive
for a DC-SQUID, in the presence of one or both of the drive tones. We expect
the full dynamics in the presence of both cavities and a colored spectrum of
environmental noise to be significantly more complicated.

Figure S9:  Floquet simulations demonstrating advantages of the
differentially-driven SQUID. a, The heating rate of the coupler from the
Floquet ground state to the Floquet excited state, as a result of two drive
tones ($\omega_{d_{1}}=2\pi\times 2.772$GHz, $\omega_{d_{2}}=2\pi\times
3.003$GHz) with various amount of drive asymmetry. Plotting the heating rate
normalized by the spectral density of the dephasing
$S_{\omega_{c}\omega_{c}}[\Delta_{d_{2},c}]$ reveals the strong dependence of
the coupler heating on the drive asymmetry, with a suppression by several
orders of magnitude as we move from large asymmetries (orange) to small
asymmetries (teal) that are practically achievable (the black arrows indicate
the value of asymmetry in this experiment, and in the case of a single-
junction transmon, respectively). b, A two-tone Floquet-Markov simulation
illustrating the choice of drive frequencies for the differentially driven
SQUID, at a fixed drive amplitude. The driven coupler heating rate is plotted
(magenta) as a function of frequency of drive tone 1 (with the other drive
frequency fixed to satisfy the beamsplitter resonance condition), in the
presence of intrinsic coupler decay. We see that even for a purely
differential drive, the amplitude-dependent squeezing (here resonant at
$\omega_{gf/2}=6.94$ GHz) can convert decay into heating (magenta). We place
the drives far red-detuned from the coupler (grey dashed line), avoiding both
the coupler squeezing and weaker nearby resonances, like the
$\phi_{d}^{4}\;\hat{c}^{\dagger^{2}}+\text{h.c.}$ process. c, The impurity of
the driven steady-state of the coupler, under a single drive tone, for the
common and differentially driven cases. The blue and white regions imply a
cold and hot coupler respectively, while the red regions are a result of
chaotic behaviour or coupler ionization. The differentially-driven SQUID shows
a significantly larger available drive space, where one can drive at higher
amplitudes at almost every drive frequency without ionizing the coupler. Some
residual limiting features in this plot are also clearly recognizable, like
coupler squeezing (white dashed line is a guide to the eye), and the
$\phi_{d}^{2}\hat{c}^{\dagger^{4}}+\text{h.c.}$ process, both of which follow
a downward trend due to the Zeeman shift.

fig]floquet_sims

Our simulations use well-established Floquet-Markov methods HanggiFloquet ;
Verney2019 ; Yaxing_PRA_BilinearModeCoupling . The Floquet basis provides an
appropriate basis for analyzing the driven coupler, with the Floquet ground
state smoothly mapping to the undriven ground state. Intrinsic coupler decay
and dephasing can result in an effective excitation out of this ground state,
represented by a ‘hopping’ rate to other Floquet states
Yaxing_PRA_BilinearModeCoupling . This can cause incoherent shifts in the
beamsplitting strength and resonance condition, and thus dephase the
beamsplitting process. Quantifying the coupler heating rate provides a simple
way to gauge the unwanted effects of the drive, with the coupler excitation
providing an estimate for the expected coupler-induced infidelity of the
process.

We simulate the coupler in its undriven eigenbasis, with the first 50
eigenstates. We calculate the Floquet modes $|\alpha_{n}\rangle$, which for
small amplitudes adiabatically map to the $n$-th undriven eigenstate (Fock
state $|n\rangle$). We calculate the transition rates $W_{mn}$ from
$|\alpha_{m}\rangle$ to $|\alpha_{n}\rangle$ via

$\displaystyle W_{mn}$
$\displaystyle=\sum_{k=-\infty}^{+\infty}S_{FF}\left[\Delta_{nmk}\right]$
$\displaystyle\times|\int_{0}^{T}\frac{e^{-ik\omega_{d}t}}{T}\langle\alpha_{m}\left(t\right)|\hat{O}|\alpha_{n}\left(t\right)\rangle
dt|^{2}.$ (39)

Here, $\omega_{d}$ is the greatest common divisor of the frequencies of the
two drives, $\omega_{d1}$ and $\omega_{d2}$, and $T=2\pi/\omega_{d}$ is the
least common multiplier of their periods. $\hat{O}$ (in the simulation taken
to be the charge operator) is the system operator that is coupled to the bath,
through $\hat{\mathcal{H}}_{int}/\hbar=\delta F\left(t\right)\hat{O}$, where
the effect of bath is modelled by the noise parameter $\delta
F\left(t\right)$. The fluctuations in this parameter are captured by
$S_{FF}\left[\omega\right]$, the noise spectral density:

$S_{FF}\left[\omega\right]=\int_{-\infty}^{+\infty}d\tau
e^{i\omega\tau}\left<\delta F\left(\tau\right)\delta F\left(0\right)\right>.$
(40)

The argument of this spectral density in Eq. 39,
$\Delta_{nmk}=\omega_{m}-\omega_{n}-k\omega_{d}$, represents all the possible
frequencies at which the transition of interest is coupled to the bath.

We start with simulating the effects of dressed coupler dephasing, by looking
at the lab-frame Hamiltonian of the coupler in the presence of both drive
tones and incoherent coupler frequency fluctuations (LABEL:floquet_simsa). We
introduce a useful quantity, the drive asymmetry $l$,

$l=\frac{\lvert\phi_{c}\rvert}{\sqrt{\lvert\phi_{c}\rvert^{2}+\lvert\phi_{d}\rvert^{2}}},$
(41)

for which 0 corresponds to a fully differential drive, and 1 to a fully common
mode drive. As a function of $l$, we calculate $W_{01}\times t_{\text{SWAP}}$,
which measures the coupler excitation after a single Alice-Bob swap operation.
For a fair comparison, we also re-scale drive amplitudes to maintain a total
beamsplitter rate of $g_{BS}/2\pi\approx 2$MHz, corresponding to
$t_{\text{SWAP}}=125$ns.

We find, as predicted by theory, that this effect is primarily sensitive to
the noise spectral density of coupler frequency fluctuations
($S_{\omega_{c}\omega_{c}}[\Omega]$) at frequencies close to the detuning of
the drive and the coupler mode ($\Omega\sim\Delta_{d_{1,2},c}$). In flux
tunable devices, this noise spectrum is expected to decrease inversely with
increasing frequency, suggesting that with a reasonable large drive detuning
($\sim 4.25$ GHz for this experiment), the effect of flux-noise induced
dressed dephasing should be negligible. We show here however, that that
experimentally achievable drive symmetries suppress dressed dephasing by
several orders of magnitude regardless of the dephasing rate
(LABEL:floquet_simsa).

Through a similar Floquet-Markov calculation with both tones, we also
calculate an effective coupler heating rate arising from a dressing of
intrinsic coupler decay (LABEL:floquet_simsb), specifically for the
differentially driven Hamiltonian. Here, we find that despite the parity
protection, there exists a significant heating rate within $\sim 1$ GHz of the
resonant squeezing, which occurs near the coupler frequency. This heating can
intuitively be understood as a combination of a
$|g\rangle\rightarrow|f\rangle$ coherent transition due to squeezing, and a
$|f\rangle\rightarrow|e\rangle$ incoherent transition due to coupler decay.
Choosing a large drive-detuning strongly suppresses the squeezing strength and
avoids this issue. We also find weaker resonances at frequencies that are at a
half ($\phi_{d}^{4}{c^{\dagger}}^{2}$) and a third
($\phi_{d}^{6}{c^{\dagger}}^{2}$) of the resonant squeezing process, but we
are able to utilize the freedom in choosing our drive frequencies to evade
them.

Finally, to provide an intuition for the combined effect of the odd-order
processes forbidden by the parity protection, we show the marked difference in
the allowed ‘drive space’ for the common and the differentially driven SQUID
(LABEL:floquet_simsc). We define the drive space to be the allowed frequencies
and amplitudes of a single drive tone that can drive the converter while
keeping the coupler in its driven ground state. Since a large portion of
either plot shows chaotic behavior and coupler ionization TransmonIonization ,
we use the impurity of the steady-state density matrix as our performance
metric, which circumvents the need for a precise mapping between the Floquet
states and the undriven states. Since these are steady state calculations, the
impurity of the coupler in regions where it is not ionized should not directly
be interpreted as the heating rate for short time scales. We observe that the
differential drive’s allowed frequency and amplitude range is much larger than
those allowed for the common drive, clearly demonstrating the former’s
advantage in this simplistic scenario. This simulation only contains a single
drive tone and simulates just the coupler mode, but we expect this advantage
in drive space to hold even in the presence of the coupled storage modes and
multiple drives, since a large fraction of the additional parasitic processes
will also be forbidden for the parity-protected Hamiltonian.

##  Supplementary Note 3 — Suppressing the back-action of the buffer cavity

While we have discussed the elimination of a linear coupling between the
buffer mode and the coupler in our Methods, quantum fluctuations in the buffer
mode might still dephase the coupler even in the presence of a purely non-
linear coupling. This can be seen by writing the buffer-SQUID Hamiltonian in
the dressed ladder-operator basis,

$\displaystyle\hat{\mathcal{H}}/\hbar=\omega_{c}\hat{c}^{\dagger}\hat{c}+\omega_{B}\hat{B}^{\dagger}\hat{B}$
$\displaystyle-\frac{E_{J_{1}}}{\hbar}\cos_{\textup{NL}}\left(\hat{\theta}_{c}+k\hat{\theta}_{B}\right)-\frac{E_{J_{2}}}{\hbar}\cos_{\textup{NL}}\left(\hat{\theta}_{c}-k\hat{\theta}_{B}\right)$
$\displaystyle\approx\omega_{c}\hat{c}^{\dagger}\hat{c}+\omega_{B}\hat{B}^{\dagger}\hat{B}-\chi_{Bc}\hat{B}^{\dagger}\hat{B}\hat{c}^{\dagger}\hat{c}.$
(42)

Here, $\omega_{B}$, $\hat{B}$ and $\hat{\theta}_{B}$ are the frequency, the
ladder operator, and the phase operator of the buffer mode.
$k=\frac{1}{2}M_{Bc}L_{B}^{-1}$ measures the participation of the buffer-mode
phase in the superconducting phases across the SQUID junctions, with $L_{B}$
and $M_{Bc}$ being the self-inductance of the buffer mode and the mutual
inductance between the buffer mode and the SQUID loop, respectively. The
cross-Kerr term emerges from the differential coupling with a strength of
$\chi_{Bc}\approx
E_{J}\theta_{c,\text{zpf}}^{2}\theta_{B,\text{zpf}}^{2}k^{2}$, and translates
the photon-number fluctuations in the buffer mode to coupler dephasing noise.
For thermal-photon distribution and coherent-photon distribution, the
dephasing rates are given by Fei2018

$\displaystyle\Gamma_{\phi}^{\text{th}}$
$\displaystyle=\chi_{Bc}^{2}\frac{\overline{n}_{\text{th}}}{\kappa_{B}},$
$\displaystyle\Gamma_{\phi}^{\text{coh}}$
$\displaystyle=\chi_{Bc}^{2}\frac{\overline{n}_{\text{coh}}\frac{\kappa_{B}}{2}}{\left(\omega_{d}-\omega_{B}\right)^{2}+\left(\frac{\kappa_{B}}{2}\right)^{2}},$
(43)

where $\overline{n}_{\text{th}}$ and $\overline{n}_{\text{coh}}$ represent the
mean thermal-photon number and coherent-photon number of the buffer mode,
respectively. On the other hand, the beamsplitter rate is only linearly
dependent on the cross-Kerr strength,
$g_{BS}\propto\chi_{Bc}\overline{n}_{\text{coh}}$. This allows us to suppress
these dephasing rates while maintaining the beamsplitting rate, simply by
reducing the cross-Kerr strength, and increasing the buffer drive strength (to
increase $\overline{n}_{\text{coh}}$). We calculated using a HFSS-EPR
simulation EPRPaper  that our design has $\chi_{Bc}\approx 2\pi\times 0.1$kHz.
which corresponds to negligible dephasing rates of
$\Gamma_{\phi}^{\text{th}}\ll 1$Hz, $\Gamma_{\phi}^{\text{coh}}\ll 1$kHz. With
a drive power of $P\approx-50$dBm at the buffer cavity drive port, we can
realize a differential drive amplitude of $\phi_{d}\approx 0.2\Phi_{0}$, that
drives a sufficiently fast beamsplitter for this experiment.

##  Supplementary Note 4 — Decoherence in the dual-rail subspace

In this section, we derive the evolution of the dual-rail density matrix
(Methods Eq. 24), and discuss how the dual-rail decay and dephasing rates
($\kappa_{1}$ and $\kappa_{\varphi}$) are related to fluctuations in the
undriven system. Without loss of generality, the rotating frame Hamiltonian of
the driven dual-rail qubit is given by

$\displaystyle\hat{\mathcal{H}}(t)/\hbar$
$\displaystyle=\left(g_{\text{BS}}+\delta
g(t)\right)\left(\hat{a}^{\dagger}\hat{b}+\hat{a}\hat{b}^{\dagger}\right)$
$\displaystyle+\delta\omega_{a}(t)\hat{a}^{\dagger}\hat{a}+\delta\omega_{b}(t)\hat{b}^{\dagger}\hat{b}$
$\displaystyle+\delta f_{a}(t)\hat{a}e^{-i\omega_{a}t}+\delta
f_{b}(t)\hat{b}e^{-i\omega_{b}t}+h.c.$ (44)

Here, $\delta g$ represents the fluctuation of the beamsplitter rate, which
could arise from a dispersion of beamsplitting strength with respect to
beamsplitting strength. $\delta\omega_{a,b}$ are the frequency fluctuations of
the Alice and Bob cavities, which could arise from their cross-Kerr to the
coupler mode, or from sources like a hot ancilla or other intrinsic dephasing.
Imperfections in control electronics, including slower drifts over time, could
also contribute to these amplitude or frequency fluctuations. Finally, $\delta
f_{a,b}$ are the cavities’ couplings to environmental decay channels. These
noises ($\delta x=\\{\delta g,\delta\omega,\delta f\\}$) are characterized by
their spectral density (eq. 40), which are related to the decay and dephasing
rates of the individual cavities,

$\kappa^{a,b}_{\downarrow}=S_{f_{a,b}f_{a,b}}[\omega_{a,b}],\;\kappa^{a,b}_{\phi}=S_{\omega_{a,b}\omega_{a,b}}[0]/2.$
(45)

To analyze how eq. 44 gives rise to decoherence in the dual-rail subspace, we
work in the eigenbasis of the beamsplitter interaction,
$\lvert\text{Alice},\text{Bob}\rangle=\lvert 00\rangle$,
$|\psi^{\pm}\rangle=\left(|01\rangle\pm|10\rangle\right)/\sqrt{2}$, and
$\lvert 11\rangle$, and rewrite this Hamiltonian in the matrix form,

$\displaystyle\hat{\mathcal{H}}(t)/\hbar=$
$\displaystyle\begin{bmatrix}0&\delta f_{-}(t)&\delta f_{+}(t)&0\\\ \delta
f_{-}(t)^{*}&\delta\Omega_{-}(t)&\frac{\delta\omega_{-}(t)}{2}e^{-2ig_{\text{BS}}t}&-\delta
f_{-}(t)\\\ \delta
f_{+}(t)^{*}&\frac{\delta\omega_{-}(t)^{*}}{2}e^{2ig_{\text{BS}}t}&\delta\Omega_{+}(t)&\delta
f_{+}(t)\\\ 0&-\delta f_{-}(t)^{*}&\delta
f_{+}(t)^{*}&\delta\omega_{+}(t)\end{bmatrix},$ (46)

where

$\displaystyle\delta\Omega_{\pm}$
$\displaystyle=\delta\omega_{a}(t)/2+\delta\omega_{b}(t)/2\pm\delta g(t),$
$\displaystyle\delta f_{\pm}$ $\displaystyle=\left(\delta
f_{a}(t)e^{-i\omega_{a}t}\pm\delta
f_{b}(t)e^{-i\omega_{b}t}\right)\frac{e^{\mp ig_{\text{BS}}t}}{\sqrt{2}},$
$\displaystyle\delta\omega_{\pm}$
$\displaystyle=\delta\omega_{a}(t)\pm\delta\omega_{b}(t).$

This allows the calculation of the transition rates among the four eigenstates
using time-dependent perturbation theory Clerk2010 , from which we obtain the
decay rates from $|\psi^{\pm}\rangle$ to $\lvert 00\rangle$ as
$\kappa_{\psi^{\pm}}=\left\\{S_{f_{a}f_{a}}[\omega_{a}\pm
g_{\text{BS}}]+S_{f_{b}f_{b}}[\omega_{b}\pm g_{\text{BS}}]\right\\}/2$. The
relaxation between $\lvert\psi_{\pm}\rangle$ due to the cavity frequency
fluctuation $\delta\omega$, as well as the beamsplitter rate fluctuation
$\delta g$, lead to dephasing in the dual-rail subspace at a rate
$\kappa_{\varphi}=\left\\{\bar{S}_{\omega_{a}\omega_{a}}[2g_{\text{BS}}]+\bar{S}_{\omega_{b}\omega_{b}}[2g_{\text{BS}}]\right\\}/4+S_{gg}[0]$,
with
$\bar{S}\left[\omega\right]=\left\\{S\left[\omega\right]+S\left[-\omega\right]\right\\}/2$
being the symmetrized spectral density. At small $g_{\text{BS}}$ and weak
dissipation rates, and assuming the cavity frequency fluctuations are
uncorrelated, one can see how these rates are related back to cavity decay and
dephasing rates,
$\kappa_{\psi^{\pm}}=\kappa_{1}=(\kappa_{\downarrow}^{a}+\kappa_{\downarrow}^{b})/2$,
$\kappa_{\varphi}=(\kappa_{\phi}^{a}+\kappa_{\phi}^{b})/2+S_{gg}[0]$.

From here, it is straightforward to write down the evolution of the dual-rail
qubit’s density matrix, in the basis of $\lvert\psi^{\pm}\rangle$,

$\displaystyle\rho_{|\psi^{\pm}\rangle}(t)=\frac{1}{2}e^{-\kappa_{1}t}\begin{bmatrix}1&e^{i\Omega
t-\kappa_{\varphi}t}\\\ e^{-i\Omega t-\kappa_{\varphi}t}&1\end{bmatrix},$ (47)

where $\Omega=2g_{\text{BS}}$. This can be easily transformed to $\lvert
10\rangle$ and $\lvert 01\rangle$ basis (denoted as $\rho_{\text{DR}}$),

$\displaystyle\rho_{\text{DR}}(t)=\frac{1}{2}e^{-\kappa_{1}t}\times$
$\displaystyle\begin{bmatrix}1+e^{-\kappa_{\varphi}t}\cos\left(\Omega
t\right)&i\sin\left(\Omega t\right)e^{-\kappa_{\varphi}t}\\\
-i\sin\left(\Omega
t\right)e^{-\kappa_{\varphi}t}&1-e^{-\kappa_{\varphi}t}\cos\left(\Omega
t\right)\end{bmatrix},$ (48)

which becomes Eq. 24 in Methods.

Figure S10:  Coupler-cavity sideband collision on increasing drive amplitude
a, Frequency stack of the driven system. Due to the Zeeman shift, the coupler
frequency is pushed towards the Alice and Bob cavities. The parametric
modulation also creates sidebands of the coupler that are spaced by the two-
tone detuning. These sidebands are coupled to the cavities through detuned
iSWAP interactions, that get closer to resonance at stronger drive amplitudes.
b, Simulation of the frequencies of the Bob cavity and the second (red-
detuned) sideband of the coupler as a function of the drive amplitude. A clear
avoided crossing is seen around drive amplitude of 0.14$\pi$ when the coupler
sideband collides into Bob, corresponding to a resonant iSWAP interaction. c,
The Floquet-Markov simulations showing the reduction of the cavities’
lifetimes with respect to the drive amplitude. Alice’s and Bob’s $T_{1}$s are
obtained from two separate simulations, where only Alice or Bob mode is
included along with the coupler, respectively. In both simulations, decay
channels were set to produce similar undriven lifetimes of Alice, Bob and the
coupler to the experimental values. d, A comparison of the mean of the decay
rates (black circles) obtained from the simulation in c to the experimentally
obtained decay rate of the dual-rail qubit (red crosses, same as main text
Fig. 7). The curves show qualitative agreement, but the simulations only
capture processes that involve the individual cavities and the coupler, and
therefore a quantitative agreement is not expected.

fig]sideband_collision

##  Supplementary Note 5 — Limits on performance from coupler-storage sideband
collision

We now provide a hypothesis for the limitations on fidelity and achievable
beamsplitting strength that we find in our device. We see in experiment
(Methods Fig. 7) that the effective beamsplitter decoherence rate remains
stable as a function of drive amplitude, until it suddenly rises at
$|\phi_{d}|\sim 0.13\pi$. This is accompanied by a super-quadratic dependence
of $g_{\text{BS}}$ on the drive amplitude. A hypothesis for this effect is an
inflation of the interaction between the cavities and the coupler, mediated by
the beamsplitter drives. This could arise from a drive-induced coupler-cavity
sideband interaction, that is then magnified by the coupler’s Zeeman shift in
the regime of strong drives. Here, we investigate the strength of such a
’sideband collision’, first focusing on Bob, which is closer in frequency to
the coupler. By combining Eqs. 11,12, and 15 in the main text, one obtains a
photon exchange interaction between the coupler and Bob in the rotating frame:

$\displaystyle\frac{\hat{\mathcal{H}}_{bc}^{(\text{RWA})}}{\hbar}$
$\displaystyle\approx\sum^{\infty}_{n=1}\omega_{c}J_{n}\left(|\phi_{d_{1}}|\right)J_{n}\left(|\phi_{d_{2}}|\right)\cos\left(n\Delta_{d}t\right)$
$\displaystyle\times\beta_{b}\left(e^{i\Delta_{bc}t}\;\hat{b}^{\dagger}\hat{c}+\text{h.c.}\right).$
(49)

Here,
$\Delta_{bc}=\omega_{b}-\omega_{c}+\omega_{\text{Z},b}-\omega_{\text{Z},c}$
represents the detuning between Bob and coupler after accounting for Zeeman
shifts. When the modulation frequency satisfies $\Delta_{d}=\Delta_{bc}/n$, a
resonant ‘iSWAP’ interaction occurs between Bob and the coupler, which we
interpret as the n-th sideband collision. This can be intuitively understood
as a sideband interaction FirstOrderSideband , as illustrated in
LABEL:sideband_collisiona. The two-tone flux modulation creates coupler
sidebands spaced by the tone detuning $\Delta_{d}$ which, when pushed down by
the Zeeman shift, start interacting strongly with the storage cavities. We
identify from Floquet simulations that around a drive amplitude of $0.14\pi$,
the second sideband of the coupler crosses the Bob mode
(LABEL:sideband_collisionb), thereby limiting Bob’s lifetime due to an
inherited decay. We expect Alice to show similar behaviour, encountering the
coupler’s third sideband at a similar Zeeman shift, but we expect this process
to be higher order and therefore weaker.

To study the implications of this increased hybridization, we perform Floquet-
Markov simulations to find predicted $T_{1}$’s for both cavities as a function
of drive strength (LABEL:sideband_collisionc). These decay rates are extracted
from two separate simulations, with each simulation including the coupler mode
and only a single cavity mode, but should be sufficient to capture the iSWAP
interactions. We therefore do not necessarily expect the simulations to show
quantitative agreement with experiment, as the experiment could be limited by
processes that simultaneously involve both cavities, or the coupler’s readout
mode. However, for a qualitative test of this hypothesis, we compare the
predicted mean cavity decay rates $\kappa_{1}$ to those found in experiment
(main text Fig. 7), and find a reasonable agreement
(LABEL:sideband_collisiond). This issue of increased hybridization could
potentially be avoided by redesigning the coupler driven frequency at the
operating point, for example by placing it red-detuned to the cavities, such
that the cavities are limited by a higher-order sideband.

##  Supplementary Note 6 — Cavity nonlinearity induced by the driven SQUID

Figure S11:  Drive-induced cavity Kerr. a,The ratio between the beamsplitter
rate, $g_{BS}$, and the self-Kerr of Bob, $K_{B}$, plotted against the
beamsplitter rate. In the original circuit (red dots and blue crosses), this
ratio is limited by the Zeeman shift as well as the sideband collision issue.
By rearranging the frequency stack such that the coupler is further detuned
from the cavities ($f_{c}=$7.75GHz, adjusting $g_{ac}$ and $g_{bc}$
accordingly to maintain the coupler-cavity hybridization), this ratio
increases by almost a factor of 5 (orange dots), at the operating point where
$g_{BS}/2\pi\approx 2$MHz. Further optimization is possible by arraying the
device into 3 SQUIDs (green dots), which roughly increases this ratio by a
factor of $3^{2}=9$. b, Swap operation infidelity for various coherent states
as a function of drive-induced Kerr, at fixed $g_{BS}/2\pi=2.16$MHz. Static
Kerr is set to be the same as the current device (red), and the arrayed device
with 3 SQUID (blue). Cavity Kerr leads to infidelity through both an undesired
photon-number-dependent phase accumulation as well as a shift in the
beamsplitting resonance condition for each photon manifold. When
$K_{b}/g_{BS}\sim 1$, this amounts to photon exchange only occurring in a
single photon manifold. As $K_{b}/g_{BS}$ is reduced the swap fidelity of
coherent states reduces quadratically with this ratio (highlighted by the
black line corresponding to zero static Kerr), until it saturates to an
infidelity limit imposed by $K_{b}^{\text{static}}/g_{BS}$. The red and green
dashed lines indicate the drive-induced $K_{b}/g_{BS}$ of the current device,
and the arrayed SQUID at higher frequency, respectively.

fig]fig:drive_induced_Kerr

In addition to the increased cavity decay rates due to the sideband collision,
large drive amplitudes on the SQUID coupler may also induce more nonlinearity
to the cavities. We investigate the magnitude of this induced nonlinearity,
it’s effect on the beamsplitter, and strategies to mitigate the driven
nonlinearity in this section.

A natural contribution to the static Kerr of the cavities comes from the SQUID
coupler, which is in the transmon regime and has a self-Kerr of
$K_{c}/2\pi\approx 125$MHz. As a result of the linear dispersive coupling, the
storage cavities weakly hybridize with the coupler and inherit self-Kerrs,
$K_{a(b)}\approx\left(\frac{g_{a(b)c}}{\Delta_{a(b)c}}\right)^{4}K_{c}$. We
directly measure this self-Kerr by first displacing the cavity to
$|\alpha\rangle$, followed by a delay time of $t$ that incurs a phase space
rotation proportional to $\Delta+\alpha^{2}\frac{K_{b}}{2}$, where $\Delta$ is
the detuning between the cavity and drive frequencies. We then use a second
displacement pulse with opposite phase that attempts to bring the cavity back
to vacuum, with full revival determined by $\Delta$, $K_{b}$ and $t$. By
performing a measurement conditioned on the cavity being in the vacuum state
using the ancilla qubit and repeating this measurement as a function of
$\alpha$, we are able to quantify this revival and extract both the Kerr and
detuning SNAPPaperWithKerrMeasurement .

The cavity self-Kerrs are measured to be $K_{a}/2\pi\approx 5$kHz and
$K_{b}/2\pi\approx 15$kHz when the coupler is undriven. When the coupler is
driven, the Zeeman shift pushes the coupler frequency lower, bringing it
closer to the cavities and therefore imposing a higher amount of nonlinearity
to the cavities. This is accompanied by the sideband collision discussed in
the previous section, which also causes a stronger hybridization between the
coupler and the cavities. At the operating point where a beamsplitter rate of
$\sim 2.2$MHz is achieved, we measure the self-Kerr of Bob cavity to be around
$128$kHz. We expect the self-Kerr of Alice to remain lower than Bob’s, since
the sideband it encounters is higher-order. While this aggravated nonlinearity
is relatively harmless in the dual-rail subspace central to our investigation,
it could introduce a coherent error in beamsplitter operations that involve
large photon numbers, and optimizing the DDS’ performance by reducing the
drive-induced Kerr is an important future direction.

There exist multiple solutions that are fully compatible with our
architecture, as demonstrated in LABEL:fig:drive_induced_Kerr. For instance,
since the dominant source of cavity nonlinearity is from the Zeeman shift and
the sideband collision, we could rearrange our mode frequencies such that the
coupler is detuned further from the cavities. Thus, at the same drive
amplitude (that produces the same $g_{BS}$), both the coupler and its sideband
stay further away from the cavities, which is helpful for decreasing the
Kerr$/g_{BS}$ ratio. We could also redesign the coupler as an array of dc
SQUIDs, which would suppress the coupler Kerr by a factor $N^{2}$, where $N$
is the number of the SQUID loops. Finally, the cavity Kerrs can also be
dynamically canceled through additional microwave drives on the coupler, as
has been shown in Yaxing_Kerr , or corrected via additional unitaries like
SNAP SNAPPaperWithKerrMeasurement .

To place the quantity $g_{BS}/K_{b}$ into context, and to get a sense of how
significantly this nonlinearity affects higher-photon manifolds, we examine
the limits on swap fidelity imposed by this coherent error. Using time-domain
QuTiP qutip  simulations, we first calibrate a swap pulse at our operating
point of $g_{BS}/2\pi=2.16$ MHz optimized for the single-photon manifold, as
done in experiment. We then apply this pulse to a coherent state in the
presence of both static and drive-induced Kerr (scaling all driven Kerrs
proportionally), computing the final state fidelity
$|\langle\psi_{ideal}|\psi_{final}\rangle|^{2}$ to the ideal swapped coherent
state $|\psi_{ideal}\rangle$.

LABEL:fig:drive_induced_Kerrb shows this swap infidelity as a function of
drive-induced $K_{b}/g_{BS}$ for different coherent states and static Kerrs.
When the drive-induced Kerr is on the order of the beamsplitting rate, the
coherent state swap fidelity suffers drastically as the beamsplitter becomes
increasingly off-resonant for all but the single-photon manifold. As the
driven Kerr is reduced, the infidelity reduces quadratically for decreasing
$K_{b}/g_{BS}$ until eventually saturating to a background infidelity limit
imposed by the static Kerr.

These simulations demonstrate that by implementing modifications compatible
with our design, we can achieve a substantial improvement in performance.
Arraying the coupler ($N=3$) and shifting to a higher frequency would lead to
a Kerr-limited swap fidelity of more than $99.9\%$ for the $|\alpha=2\rangle$
coherent state, and above $99\%$ for $|\alpha=3\rangle$. Note that this effect
could be mitigated to some extent by calibrating the beamsplitter in photon
manifold corresponding to the mean coherent state photon number. While
applications with much more stringent linearity requirements exist
GKPExperiment ; VladGKPBeyondBreakeven , significant progress towards
realizing hardware-efficient error correction has been made with protocols
using states with mean photon numbers as low as the states simulated here
BeyondBreakeven2016 ; ChenWangAQEC ; BinomialBeyondBreakeven . However for
most applications to high fidelity control of logical qubits, including an
additional unitary resolving this coherent error will be useful.

##  Supplementary Note 7 — DC flux calibration and junction asymmetry

Figure S12:  DC flux and junction asymmetry. a, Spectroscopy of the coupler
readout resonator while tuning source voltage of the DC flux coil. This curve
provides an initial estimate of the zero flux point as well as the conversion
between coil voltage and Flux Quantum. b, Coupler frequency vs. DC flux. The
data is fit to a model that takes both linear inductance of the SQUID loop and
hybridization between coupler, Alice, and Bob into account. c,  Junction
asymmetry as backed out from room temperature resistance measurements from
SQUID devices on multiple wafers. These devices were fabricated using the same
methods as (and contemporaneous with) the device demonstrated in the main text
(section Supplementary Note 11).

fig]flux_and_LJ

Operating at the DC flux sweet spot of the SQUID is not only desirable for
combating flux noise, but is also critical to keeping our aforementioned
parity-protection. To ensure that we operate at $\Phi_{\text{DC}}=0$, we have
a coil built in to our package to apply a static magnetic field through the
SQUID loop. This coil is made out of thin NbTi wire wrapped around a copper
spool, similar to the magnet in many resonant JPA designs, and is placed
sufficiently far from the converter and cavities, so as to not spoil any
high-Q modes. The coil is connected to a DC power supply, where we can tune
the source Voltage as desired.

The spectroscopy of the coupler readout serves as an initial estimate of the
zero flux point, where we observe a periodic tuning of the coupler readout
frequency as the DC source is swept, due to the repulsion of the coupler mode.
(LABEL:flux_and_LJa).

To obtain a more precise and more direct measurement of the zero-flux point,
we then perform coupler spectroscopy while sweeping the now-calibrated DC
source (LABEL:flux_and_LJb). The coupler frequency is measured near zero flux
as well as further on the slope. As the coupler frequency sweeps into
resonance with Alice and Bob, the modes strongly hybridize, obscuring the
coupler resonance. The measured coupler frequency is well fit to a model that
takes both the linear inductance in the SQUID loop and the hybridization with
Alice and Bob into account, yielding the simulated points in
LABEL:flux_and_LJb.

Another requirement of the differentially-driven Hamiltonian (Eq. 1 in the
main text) is that the Josephson energies of the two junctions in the SQUID
are the same, which in the practical sense means the junctions should have
identical size. While this is the case in design, imperfections in the
fabrication processes, especially in the e-beam writing step, results in an
inevitable deviation of the $E_{J}$ of the two junctions. To quantify the
typical value of the $E_{J}$ asymmetry, we measure the room-temperature
resistance of the Josephson junctions, which is related to their Josephson
energies through the Ambegaokar-Baratoff relation AB_relation . Through this
way, we measured the $E_{J}$ asymmetry, which we define as
$\left(E_{J_{1}}-E_{J_{2}}\right)/\left(E_{J_{1}}+E_{J_{2}}\right)$, for the
SQUID devices on four different wafers that were fabricated around the same
time (LABEL:flux_and_LJc). From these measurements, We observe a negligibly
small $E_{J}$ asymmetry of 0.3%$\pm$1.3%.

##  Supplementary Note 8 — Calibrating a beamsplitter pulse

In order to calibrate the beamsplitter that is used in the RB experiment, we
must accurately tune up the pulses on the RF controls. While fitting a chevron
(Main text Fig. 4b) can approximately find the resonance condition,
calibrating the system to the degree required for fidelities above 99.9%
requires more precise protocols. Using a four-wave mixing process to engineer
the beamsplitter interaction also incurs drive-induced Zeeman shifts of the
cavity frequencies that are on the same order as the beamsplitter rate,
further complicating the calibration. If these shifts were identical for Alice
and Bob, then the beamsplitting resonance condition would not change, but any
relative difference between these shifts will manifest as a shift in the
beamsplitting resonance condition. Most notably, this means that the resonance
condition changes as a function of drive amplitude, complicating the process
of calibrating a pulse with finite ramps.

Figure S13:  Beamsplitter Pulse Calibration. a, Sweeping drive amplitude and
detuning, we perform $2N_{\text{swaps}}$ consecutive beamsplitter pulses with
identical phase, leading to bands that correspond to plausibly correct
calibration values. b, By performing the same number of beamsplitter pulses as
(a), but changing the phase of every two beamsplitters relative to the prior
two, we can effectively engineer a sequence equivalent to
$U_{\text{SWAP}}Z(\varphi)U_{\text{SWAP}}\dots Z(\varphi)U_{\text{SWAP}}$. For
a value of $\theta_{R}=1.3$ radians (empirically found), we recover bands that
now trace paths roughly perpendicular to those found from the (a) protocol. c,
The correct calibration setting would pass both protocols, so we take the
overlap, or product, of (a) and (b) to recover a region that the correct
calibration is in. d-f, By increasing $N_{\text{swaps}}$, the bands get
tighter individually, but also more densely spaced. This leads to multiple
regions in the overlap plot (f), but only one region has the same center as
regardless of $N_{\text{swaps}}$, in this case the center region. By
iteratively narrowing the scope and increasing $N_{\text{swaps}}$, we can
refine our calibration ever finer. g, Trajectories of states initialized on
the poles of the dual-rail Bloch Sphere for the correct calibration. The pulse
we calibrate for is detuned from the coherent evolution resonance condition,
but the detuning combines with Zeeman shifts during the ramps (which appear as
Z rotations on the Bloch sphere) to perfectly compensate and give the desired
unitary. h, Calibration vs $N_{\text{swaps}}$. At each value of
$N_{\text{swaps}}$, we fit the overlap region to a 2D Gaussian and extract
center frequency (brown) and amplitude (green) as well as the fit error. By
dividing these values by the final drive tone detuning
($\Delta_{d}=237.26367\text{ MHz}\pm 50\text{ Hz}$) and final amplitude
($\sqrt{\phi_{d_{1}}\phi_{d_{2}}}=0.1214367\pm 4\times 10^{-7}$) respectively,
we can plot the calibration normalized to the final parameter set as a
function of $N_{\text{swaps}}$ i, To track the significance of calibration
drifts, we run the protocol for $N_{\text{swaps}}=501$ every 15 minutes over a
24 hour period, finding no more than a 0.03% fluctuation over this timescale.
The frequency calibration has no observable drifts.

fig]Pulse_Cal_Fig

Within the dual-rail qubit, Eq. 18 in the main text reduces to

$\displaystyle\hat{\mathcal{H}}_{\text{BS}}/\hbar$
$\displaystyle=((\Delta_{ab}-\Delta_{d})+\Delta_{\text{Z},ab}\varepsilon(t))\frac{\hat{\sigma}_{Z}}{2}$
$\displaystyle+g_{\text{BS}}\;\varepsilon(t)\;(e^{i\varphi_{\text{BS}}}\hat{\sigma}_{+}+e^{-i\varphi_{\text{BS}}}\hat{\sigma}_{-}),$

where $\varepsilon(t)\in[0,1]$ is a normalized envelope function corresponding
to the square of the pulse shape, as both pulses are ramped simultaneously.

In this Bloch sphere picture, the polar angle $\theta(t)$ of a state starting
at the North Pole is given by

$\displaystyle\theta(t)$
$\displaystyle=2\int_{0}^{t}g_{\text{BS}}\;\varepsilon(t^{\prime})dt^{\prime}.$

The azimuthal phase drift during the operation is then given by

$\displaystyle\delta\phi(t)$
$\displaystyle=\int_{0}^{t}\sin(\theta(t^{\prime}))\frac{((\Delta_{ab}-\Delta_{d})+\Delta_{\text{Z},ab}\;\varepsilon(t^{\prime}))}{2}\;dt^{\prime}$
(50)

There exists a static detuning $\Delta_{d}$ such that, at a time
$\tau_{\text{BS}}$ given by $\theta(\tau_{\text{BS}})=\frac{\pi}{2}$ we can
have complete cancellation of the azimuthal phase drift
$\delta\phi(\tau_{\text{BS}})=0$. Such a pulse takes the trajectory of a
detuned process during the steady state before arriving back to the correct
final position on the Bloch sphere during the ramp down
(LABEL:Pulse_Cal_Figg).

To find the correct parameter set that achieves this, we can prepare a single
photon in Bob, apply a candidate beamsplitter pulse a total of $4n+2$ times,
and measure P(Bob = 0) to see whether the photon is successfully swapped out
of Bob. Specifically, we fix the total pulse time and envelope shape while
simultaneously sweeping tone 1 drive frequency (and thus $\Delta_{d}$) and
drive strength (and thus $g_{\text{BS}}$ and $\Delta_{\text{Z},ab}$) of the
candidate pulse to scan the parameter space. If the candidate pulse does
indeed implement a beamsplitter, then we will have implemented a total of
$N_{\text{swaps}}=2n+1$ swap operations on the photon, leaving Alice with the
photon and Bob in vacuum.

To gain some insight as to what we should expect when sweeping these
parameters, we consider an adapted Rabi model with an amplitude-dependent
detuning. Under this model, we find the probability $P_{\text{
B$\rightarrow$A}}$ to oscillate as:

$\displaystyle P_{\text{ B$\rightarrow$A}}$
$\displaystyle\approx\sin^{2}\left(2\sqrt{\Omega^{2}+(\Delta_{0}+\alpha\Omega)^{2}}t\right),$

where $\Omega$ and $\Delta_{0}$ are drive strength and detuning respectively,
and $\alpha$ is some coefficient relating drive strength to frequency shift.
For a fixed time $T$, we would expect contours of constant $P_{\text{
B$\rightarrow$A}}$ to be of the form:

$\displaystyle\Omega^{2}+(\Delta_{0}+\alpha\Omega)^{2}$
$\displaystyle=\left(\frac{(2n+1)\pi}{4T}\right)^{2}.$

This equation, when varying $\Omega$ and $\Delta_{0}$, results in ellipses
with different bands corresponding to different values of $n$. As seen in
LABEL:Pulse_Cal_Figa,d, a sweep of the parameter space does indeed recover
these bands, with different bands corresponding to different odd integer
fractions of the correct band. For example, a sequence that used $(4(5)+2)$
pulses but only implemented 3 total swaps would produce a band at
$\frac{3}{5}$ of the radius of the correct band.

To narrow down precisely where along these bands the correct detuning to
properly null eq. 50 is, we run a complementary pulse sequence. This
complementary sequence consists of pairs of candidate beamsplitter pulses,
with the phase $\varphi_{\text{BS}}$ of each pair shifted by an angle
$\theta_{R}$ relative to the previous pair. This relative angle between pulse
pairs results in a rotation of the bands as shown in LABEL:Pulse_Cal_Figb,e,
and a value of $\theta_{R}=1.3$ radians makes the new bands approximately
orthogonal to the unrotated version. The correct set of calibration parameters
successfully swaps the photon under both of these pulse sequences, so by
taking the overlap (product) of these 2D sweeps we can get patches of possible
parameter sets as seen in LABEL:Pulse_Cal_Figc,f.

Since the density of the contour lines is set by the number of pulses in the
sequence, increasing $N_{\text{swaps}}$ not only makes the bands thinner
individually, but also makes them appear closer together. Putting all of this
together, we can start at a low value of $N_{\text{swaps}}$ where the bands
are spaced far apart and the overlap region is wide. Then, by refining the