Random access quantum information processors using multimode circuit quantum electrodynamics

Qubit connectivity is an important property of a quantum processor, with an ideal processor having random access—the ability of arbitrary qubit pairs to interact directly. This a challenge with superconducting circuits, as state-of-the-art architectures rely on only nearest-neighbor coupling. Here, we implement a random access superconducting quantum information processor, demonstrating universal operations on a nine-qubit memory, with a Josephson junction transmon circuit serving as the central processor. The quantum memory uses the eigenmodes of a linear array of coupled superconducting resonators. We selectively stimulate vacuum Rabi oscillations between the transmon and individual eigenmodes through parametric flux modulation of the transmon frequency. Utilizing these oscillations, we perform a universal set of quantum gates on 38 arbitrary pairs of modes and prepare multimode entangled states, all using only two control lines. We thus achieve hardware-efficient random access multi-qubit control in an architecture compatible with long-lived microwave cavity-based quantum memories.

The device is heat sunk via an OFHC copper post to the base stage of a Bluefors dilution refrigerator (10-30 mK). The sample is surrounded by a can containing two layers of µ-metal shielding, thermally anchored using an inner close fit copper shim sheet, attached to the copper can lid. The schematic of the cryogenic setup, control instrumentation, and the wiring of the device is shown if Supplementary Figure 1. The device is connected to the rest of the setup through three ports: a charge port that applies qubit and readout drive tones, a flux port for shifting the qubit frequency using a DC-flux bias current and for applying RF sideband flux pulses, and an output port for measuring the transmission from the readout resonator. The charge pulses are generated by mixing a local oscillator tone (generated from an Agilent 8257D RF signal generator), with pulses generated by a Tektronix AWG5014C arbitrary waveform generator (TEK) with a sampling rate of 1.2 GSa/s, using an IQ-Mixer (MARQI MLIQ0218). The charge drive pulses are amplified at room temperature (30 dB), and subsequently combined with the readout drive pulse, generated from a second Agilent 8257D RF signal generator, which is also controlled by digital trigger pulses from the TEK. The combined signals are sent to the device, after being attenuated a total of 60 dB in the dilution fridge, using attenuators thermalized to the 4K (20 dB) and base stages (10 + 30 dB). The charge drive line also includes a lossy ECCOSORB CR-117 filter to block IR radiation, and a low-pass filter with a sharp roll-off at 6 GHz, both thermalized to the base stage. The flux-modulation pulses are directly synthesized by a Tektronix AWG70001A arbitrary waveform generator (50 GSa/s) and attenuated by 20 dB at the 4 K stage, and bandpass filtered to within a band of 400 MHz -3 GHz at the base stage, using the filters indicated in the schematic. The DC flux bias current is generated by a YOKOGAWA GS200 low-noise current source, attenuated by 20 dB at the 4 K stage, and low-pass filtered down to a bandwidth of 3 kHz using a home built C-L-C π filter with a coil inductor (superconducting wire wound on a spool machined out of Vim Var core iron) and NPO capacitors. The DC flux bias current is combined with the flux-modulation pulses at a bias tee thermalized at the base stage. The state of the transmon is measured using the transmission of the readout resonator, through the dispersive circuit QED readout scheme [1]. The transmitted signal from the readout resonator is passed through a set of cryogenic circulators (thermalized at the base stage) and amplified using a HEMT amplifier (thermalized at the 4 K stage). Once out of the fridge, the signal is filtered (tunable narrow band YIG filter with a bandwidth of 80 MHz) and further amplified. The amplitude and phase of the resonator transmission signal are obtained through a homodyne measurement, with the transmitted signal demodulated using an IQ mixer and a local oscillator at the readout resonator frequency. The homodyne signal is amplified (SRS preamplifier) and recorded using a fast ADC card (ALAZARtech).

Supplementary Note 2 -Device design and fabrication
The CPW resonators in the array have a center pin width of 12 µm and a gap width of 6 µm. They are coupled to each other via interdigitated capacitors, where each side of the capacitor has 6 digits that are 107 µm long, 6 µm wide, and spaced by 6 µm. The capacitor coupling the array to the qubit is identical to the intra-array capacitors to minimize disorder of the resonators in the array. The transmon is capacitively coupled to ground via CPW capacitors on either side of the SQuID, with a center pin width of 20 µm and gap of 10 µm. The SQuID is a 20 µm by 10 µm loop, with two square junctions that are 170 nm and 125 nm wide. The flux bias line 6 µm from the SQuID is dipolar, with 25 µm long and 2 µm wide wires on each end. The ground plane of the chip has an array of 5 µm wide square holes spaced by 50 µm for flux vortex pinning [2]. The linear elements of the full circuit design have been simulated with a commercial 3D finite element analysis software (ANSYS HFSS).
The device (shown in Figure 1 in the main text) was fabricated on a 430 µm thick C-plane sapphire substrate. The base layer of the device, which includes the majority of the circuit (excluding the Josephson junctions of the transmon), consists of 100 nm of aluminum deposited via electron-beam evaporation at 1 Å/s, with features defined via optical lithography and reactive ion etch (RIE) at wafer-scale. The wafer was then diced into 7x7 mm chips. The junction mask was defined via electron-beam lithography with a bi-layer resist (MMA-PMMA) in the Manhattan pattern, with overlap pads for direct galvanic contact to the optically defined capacitors. Before deposition, the overlap regions on the pre-deposited capacitors were milled in-situ with an argon ion mill to remove the native oxide. The junctions were then deposited with a three step electron-beam evaporation and oxidation process. First, an initial 35 nm layer of aluminum was deposited at 1 Å/s at an angle of 29 • relative to the normal of the substrate, parallel azimuthally to one of the fingers in the Manhattan pattern [3] for each of the junctions. Next, the junctions were exposed to 20 mBar of high-purity O 2 for 12 minutes for the first layer to grow a native oxide. Finally, a second 120 nm layer of aluminum was deposited at 1 Å/s at the same angle relative to the normal of the substrate, but orthogonal azimuthally to the first layer of aluminum. After evaporation, the remaining resist was removed via liftoff in N-Methyl-2-pyrrolidone (NMP) at 80 • C for 3 hours, leaving only the junctions directly connected to the base layer, as seen in the inset of Figure 1 of the main text. After both the evaporation and liftoff, the device was exposed to an ion-producing fan for 15 minutes, in order to avoid electrostatic discharge of the junctions.
The device is mounted and wirebonded to a multilayer copper PCB microwave-launcher board. Additional wirebonds connect separated portions of the ground plane to eliminate spurious differential modes. The device chip rests in a pocketed OFHC copper fixture that presses the chip against the launcher board. Notably, the fixture contains an additional air pocket below the chip to alter 3D cavity modes resulting from the chip and enclosure, shifting their resonance frequencies well above the relevant band by reducing the effective dielectric constant of the cavity volume.

Supplementary Note 3 -Multimode quantum memory Hamiltonian
Our multimode quantum memory implementation uses the eigenmodes of a linear array of n = 11 identical, strongly coupled superconducting resonators [4], as shown in Figure 1 of the main text. The array is described by the where ν r is the resonance frequency of the identical resonators, g r is the tunnel-coupling between neighboring resonators, andĉ † j (ĉ j ) is the operator that creates (annihilates) photons in the resonator at spatial index j. The coupling (g r ∼ 250 MHz) between neighboring resonators is larger than the disorder in the resonator frequencies (ν r = 6.45±0.1 GHz). Thus, the single-photon eigenmodes of this circuit are 11 distributed "momentum" states of the array, with eigenfrequencies in a band from 6 to 7 GHz and mode spacing varying from 50 to 150 MHz. The eigenmodes and eigenfrequencies of the tight-binding Hamiltonian of Supplementary Equation (1) are [4]: where |ψ k and ν k are the kth eigenstate and eigenfrequency, respectively, and |1 j is the state with a single photon in the jth resonator of the array and with all other resonators in the ground state. The coupling between the transmon and a given eigenmode is given by:

Supplementary Note 4 -Parametric control Hamiltonian
The frequency of the transmon is tunable using the magnetic flux threading the SQuID loop of the transmon, controlled by passing a current through a nearby flux line. For a sinusoidally modulated flux, the transmon |g − |e transition frequency is: is the mean qubit frequency during the flux modulation. The relation between the frequency ( ) and flux ( Φ ) modulation amplitudes and the DC-shift of the transmon frequency during flux-modulation are: respectively. The frequency is shifted from its bare value due to the non-linear flux dependence of the transmon frequencies, and is quadratic in the modulation amplitude. We obtain a simple description for parametric control of this system by considering the Hamiltonian of Equation (1) in the main text and restricting to the lowest two transmon levels. The Hamiltonian then reduces to: The lowest two levels form a qubit, whose frequency is modulated over time by using the flux bias. In this work, we focus on iSWAP interactions between the trasmon and the resonator mode. As a result, we modulate the transmon frequency near the the difference frequencies of the memory modes and transmon. The (b kσ− + c.c.) terms in Supplementary Equation (7) can be therefore be dropped in the rotating-wave approximation. When the modulation frequency is resonant with the detuning between the eigenmode and the transmon, this realizes effectively resonant interactions and stimulated vacuum Rabi oscillations of a single photon between the transmon and the mode. These sidebands manifest in a rotating frame defined by the transformation [5,6]: In this rotating frame, the Hamiltonian is transformed to: where ∆ k = ν k −ν ge , is the detuning between the qubit and the k th eigenmode. When ν sb = ∆ k , we obtain resonant first-order sideband transitions between the transmon and mode k, described by: We perform universal operations on the multimode memory using iSWAP operations between a given mode and both the |g − |e and |e − |f transmon transitions, with the latter allowing the realization of entangling gates between arbitrary eigenmodes. The minimal description of our gate operations on the multimode memory therefore involves three transmon levels, with the parametric control of the eigenmodes described by an extension of the Hamiltonian of Supplementary Equation (9) to a single qutrit coupled to the harmonic memory modes. In addition to sideband transitions, the Hamiltonian also includes dispersive shifts arising from photons in the memory modes, due to the bare coupling between the eigenmodes and the transmon. We ignore the correction to the dispersive shift due to the modulation amplitude dependence of the bare term (∝ J 0 ( m /2ν sb )), whose lowest order contribution is quadratic in /2ω m . These effects are described by the following simplified Hamiltonian: In the above, 2ν sb . α k (t) for α ∈ {ge, ef } are the strengths of time-dependent parametric frequency modulation tones addressing mode k and ∆ α k = ν k −ν α is the detuning between mode k and the frequency of the corresponding transmon transition frequency α ∈ {ge, ef }. Ω α (t) are the strengths of the transmon charge drives and χ e,f k are the dispersive shifts of the |e and |f levels resulting from the addition of a photon in mode k. In addition to the dispersive shift, there are second-order tunneling terms of the form b † l b k for l = k arising from the virtual hopping of photons between different eigenmodes via the transmon. These terms are of the same order as the dispersive shift, but their effect can be ignored since they correspond to off-resonant tunneling (∼ 1 MHz) between non-degenerate levels (spaced by ∼ 100 MHz). We note that there is also a shift (DC-offset) of the qubit frequency during the flux modulation, arising from the non-linear flux-frequency relation of the transmon. Given that the flux pulses used in the experimental sequences in this work are sequential, we include their effect as: with additional cross terms being present if different flux tones were turned on simultaneously. We can further simplify the Hamiltonian above by ignoring off-resonant terms. If the transmon charge and flux-modulation tones are of the form q,sb cos (ω q,sb t + φ q,sb ), and we consider near resonant operations with a single eigenmode k, the drive phases (φ q,sb ) enter the effective Hamiltonian as: Reducing to 2 × 2 subspaces over which each of these terms act, and taking the top row to be the state with the higher transmon level, and with θ(t) = 2Ω α t and 2g sb,α t for the sideband and qubit drives, we obtain: MHz, while the SQuID loop junction asymmetry, (E J1 − E J2 )/(E J1 + E J2 ) = 0.1. These parameters correspond to maximum and minimum qubit frequencies of 5.84 GHz and ∼ 2 GHz, respectively. The experiments in this work were typically performed with the transmon biased between 3.9 − 4.7 GHz (see Supplementary Figure 2a). This frequency band is ∼ 2 GHz away from the eigenmodes of the resonator array. As a result, photons in the multimode manifold cause relatively small dispersive shifts of the transmon frequency. Additionally, the slope of the flux-frequency curve in this regime allows for sufficiently large frequency modulation amplitudes, while limiting sensitivity to flux noise to maintain transmon coherence. The transmon qubit state is probed using a capacitively coupled CPW readout resonator. The frequency and the quality factor of the readout resonator are ν read = 5.255 GHz and Q = 15000, and the coupling to the qubit is g read = 47 MHz. For the typical transmon frequency range, we obtain single-shot readout fidelities between 0.3 − 0.85 using dispersive [1] and high-power [7] circuit QED readout schemes. The readout signal is calibrated by appending a sequence with no pulse, and one with a transmon |g − |e π pulse at the of each set of experimental sequences. Upon averaging over 1000-2000 experiments, the readout signal results in a visibility of ∼ 99%, limited by the fidelity of the single qubit gates, and consistent with the RB data.
The coherence of the transmon is characterized by standard lifetime (T 1 ) and Ramsey (T * 2 ) experiments. The measured T 1 of the |e state and T * 2 of the |g − |e transition are show as a function of the |g − |e transition frequency in Supplementary Figure 2b and 2c. The T 1 is found to show a slight decrease with increasing frequency in this regime, explained partially by increased Purcell loss from coupling to the readout resonator (see Supplementary Figure 2b). The T 1 at a given flux bias slowly varies with time (over the course of weeks) by ∼ 25 − 30%. The T * 2 is found to increase with frequency, consistent with reduced sensitivity to flux noise due to the decreasing slope of the frequency-flux curve transform. We note that T * 2 showed no improvement from reducing the cutoff frequency of the external cryogenic low-pass filter on the DC flux bias from 100 to 3 kHz. At ν q = 4.36 GHz, the T 2 obtained from a spin-echo experiment with a single π pulse is 3.7 µs. This time could be increased to ∼ 14 µs using a dynamical decoupling sequence (Carr-Purcell-Meiboom-Gill with 61 pulses) [8]. We note that the measured T * 2 jumped from 400 ns to 1.2 µs, 2-3 weeks following cooling the fridge to the base temperature (20 mK), coincident with a shift and stabilization of the applied current corresponding to a flux quantum. The coherence of the |f level is characterized by analogous lifetime and Ramsey experiments. The lifetime of the |f level at ν ge q = 4.3325 GHz is T 1,ef = 3.7 µs while the phase coherence time is T * 2,ef = 1.2 µs.

Supplementary Note 6 -Stimulated vacuum Rabi oscillations with multiple transmon levels
The eigenmodes of the resonator array are probed by sideband spectroscopy using the RF flux bias (see Supplementary Figure 1). The flux-modulation pulses are directly synthesized using a Tektronix AWG70000A arbitrary waveform generator (AWG) with a sampling rate of 50 GSa/s. The target pulse waveform envelopes are typically square, with Gaussian edges (σ = 10 − 20 ns, 2σ cutoff) added to reduce the pulse bandwidth and thus minimize crosstalk between modes. The protocol used to the measure the |g − |e sideband spectrum is shown in Figure 2 of the main text and reiterated in Supplementary Figure 3a, with the corresponding measured spectrum shown in with the corresponding |g − |e sideband spectrum is shown in Supplementary Figure 4a. The stimulated vacuum Rabi chevrons for some modes are found to be distorted from the expected shape [6] (see mode 0 and 10). These distortions are due to resonances in the transmission profile of the flux bias line, as seen in the DC-offset spectroscopy. The memory modes appear as avoided crossings in a Ramsey experiment on |g − |e transition, due to interference of Ramsey fringes arising from the DC-offset and resonant stimulated vacuum Rabi oscillations. At a given frequency, the DC-offset shows a quadratic dependence on the modulation amplitude seen in Supplementary Figure 4(c), as expected from equation (6). The amplitude profile of the transfer function T (ν) of the flux bias line is obtained from the DC-offset at fixed AWG drive voltage, with |T (ν)| ∝ |δν DC |.
For short pulse durations and large stimulated vacuum Rabi rates, the bandwidth of the pulse becomes commensurate with the frequency scale over which the transfer function of the flux bias varies significantly, causing distortion of the flux pulses. However, this effect is corrected using the knowledge of the transfer function of the flux bias. The complete complex transfer function (characterizing the amplitude and the phase of the flux bias distortion) is extracted only from the amplitude of the transfer function, by assuming that the response of the line is causal and enforcing the Kramers-Kronig relations [9]. We enforce these relations here by fitting the amplitude of the measured transfer function to the functional form: We account for flux-pulse distortion by modifying the pulses generated by the AWG to account for the transfer function of the flux bias line. The AWG waveform used to generate a given pulse f (t) = Re[f c (t)] at the location of the qubit is:

Supplementary Note 8 -Coupling between the transmon and the multimode memory
The rate of stimulated vacuum Rabi oscillations are related to the modulation strength and the bare coupling according to Supplementary Equation (10) of the main text. We can therefore extract the bare couplings from the measured sideband Rabi oscillation rates, particularly since the strength of the modulation can be independently calibrated from spectroscopy of the DC-offset of the transmon frequency (see Supplementary Note 7).
Instead, we measure the eigenmode-state dependent dispersive shift of the transmon frequency. The shift for each mode k is measured with a transmon Ramsey interference experiment conducted after loading a photon into mode k, according to the protocol shown in Supplementary Figure 7 The qubit frequency is stabilized to 50 kHz prior to each measurement and corresponds to the indicated error bar. b, The transmon-mode couplings g k , extracted from the measured dispersive shift using equation (21), with the errors inferred from those in χ.
oscillation frequency ν osc and the Ramsey frequency ν Ram according to δ i = ν osc − ν Ram , and is plotted as a function of mode number in Supplementary Figure 7(b). We extract the coupling rate g k from the measured dispersive shift χ k , which are related by [10]: where α is the transmon anharmonicity and ∆ k = ν q − ν k is the detuning between the transmon and mode k. The g k 's extracted from this expression are shown in Supplementary Figure 7(c). The bare coupling rates extracted from the stimulated vacuum Rabi rate and from the dispersive shift are found to be consistent.

Supplementary Note 9 -Hamiltonian tomography
We use Hamiltonian tomography [11] to extract the 2N − 1 parameters of a chain of N nearest-neighbour coupled resonators. We assume Hamiltonian for this chain is given by equation (1)  the individual resonator frequencies (ν r,i ) and tunnel couplings (g r,i ): We extract these parameters using the frequencies (ν k ) and couplings to the transmon (g k ) of the eigenmodes of the array (2N numbers). The coupling to the transmon is proportional to the amplitude of the memory-mode wavefunction at the edge resonator (g eff,k ∝ |φ k 1 |), where the creation operator for eigenmode k isb † k = i φ k iâ † i . The bare frequencies and tunnel-couplings of the resonator are then extracted by iteratively solving the Schrödinger equation starting from the transmon end of the chain, while imposing the constraints from wavefunction normalization , as shown below: The bare frequencies and tunnel couplings thus extracted are shown in Supplementary Figure 8a and 8b. We infer that two of the normal modes (see Supplementary Figure 3) being extremely weakly coupled is likely a result of coupler 7 and 9 being defective.

Supplementary Note 10 -Coherence of the multimode memory
The coherence times of the memory modes are characterized through protocols analogous to those for the transmon, with the qubit pulses sandwiched with a pair of transmon-mode iSWAP pulses to transfer the quantum state between the transmon and the mode. The results of T 1 and T * 2 measurements for mode 1 of the memory are shown in Supplementary Figure 9a   The frequency of the iSWAP pulse acting on a particular mode is obtained by choosing the frequency corresponding to maximum contrast of the stimulated vacuum Rabi chevrons, such as those in Supplementary Figure 3. The amplitude and pulse bandwidth of the flux-modulation pulses are optimized to maximize the oscillation rate, while minimizing cross talk with neighboring modes and sideband transitions across other transmon levels. The length of an iSWAP pulse is obtained using fits of stimulated vacuum Rabi oscillations such as in Figure 2c of the main text. To achieve high fidelity gate operations, we also calibrate and correct phase errors arising during the sideband pulses.
The main phase error in the flux-modulation pulse is due to the DC-shift of the transmon frequency during fluxmodulation (δν DC in equation (6)). The frequency of the center of the stimulated vacuum Rabi chevron is detuned from the difference frequency between the mode and the relevant transmon transition by −δν DC . Since the flux-pulse frequency is set to the center of the chevron, the clock (rotating frame) of the drive generator is shifted from the the frame of the Hamiltonian of equation 7. If the drive generator idles on resonance, there is an additional phase that accrues during that time. In the Ramsey experiment measuring the coherence time (T * 2 ) of the modes (see Supplementary Figure 9a (top)), the accrued phase shifts the frequency of the Ramsey fringes by δν DC . We can then account for the misalignment of clocks by advancing the phase of the subsequent pulse by 2πδν DC τ . This correction can be easily implemented by keeping the drive clock aligned with the bare qubit-resonator system when the flux pulse is off, and incrementing the drive frequency by δν DC during the iSWAP pulse to bring it into resonance with the DC-shifted frame.
Fixing the drive clock to be in sync with the Hamiltonian of Supplementary Equation (7) results in the absence of idle-time dependent phase errors. We additionally need to calibrate an additional dynamical phase (σ z error) that occurs due to the change in the qubit frequency during the ramp up of the flux pulse. This phase is manifest in a rotating frame corresponding to the instantaneous qubit frequencyν ge (t) in equation (8). Repeating the transformation of equation (9) with a time-dependent qubit frequency results in an additional term: If we consider a square flux pulse with modulation amplitude corresponding to a DC-offset of ν DC and pulse duration of t π , the additional term in the Hamiltonian results in a dynamical phase of πν DC t π . This error is calibrated using the sequence shown in Supplementary Figure 10a  Supplementary Figure 10 | a, Circuit for calibrating phase of the iSWAP gate, where we sweep the phase δφ 2 added and subtracted from the iSWAP pulses used to load and unload the state to the memory mode. We introduce an additional offset phase (φ π 2 ) which corrects σz errors occurring in the qubit pulse, ensuring that the mode state at the end of the first sideband pulse is free from all σz errors. The optimal phase for the iSWAP pulses is obtained by minimizing P (δφ). b, iSWAP phase calibration for mode 9 of the multimode memory, optimal phase (92.5 ± 0.1 • ) indicated by the black dashed line.
iSWAP phase (φ π−π ), we add (subtract) φ π−π 2 to every iSWAP pulse for loading (unloading) an excitation into each of the memory modes. Subsequent |g − |e iSWAP pulses in all circuit diagrams include this phase correction, and are represented byπ when represented in an equation.

Supplementary Note 12 -Randomized benchmarking
We characterize the fidelity of single-mode operations using randomized benchmarking (RB) [12]. As described in the main text (see Figure 3), single-mode Clifford gates are realized by sandwiching single-qubit Clifford gates (C i ) with transmon-mode iSWAP pulses.C i = Uπ sb C i Uπ sb (φ=π) = C i To load the excitation to and from the transmon, we use σ z error corrected sideband iSWAP pulses that are 180 • out of phase with each other, so that the mode Cliffords are mapped directly from their transmon qubit counterparts. The Cliffords are generated by concatenating an operator each from {0, π 2 y , π y , − π 2 y } and {0, π 2 x , π x , − π 2 x , π 2 z , − π 2 z }, to generate all 24 elements of the single qubit Clifford group. The circuit showing the sequence used for RB of the modes is shown in Supplementary Figure 11.

C1 C2 Cinv
Supplementary Figure 11 | Randomized benchmarking (RB) characterizes the average random gate fidelity by acting random sequences of single-mode Clifford gates of increasing length, inverting the sequence, and then measuring the qubit state as a function of length. The sequence is acted on the multimode system initialized in the ground state, with the mode occupation error ( ) measured at the end of each sequence. The RB fidelity is extracted from the decay of the occupation fidelity (1 − ) as a function of the length of the benchmarking sequence. For the data shown in Figure 3 of the main text, we use sequence lengths corresponding to those in [12,13] and average over 32 random sequences.
The RB fidelity (p) is extracted by fitting the decay curves to the form Ap m + B, where m is the sequence length. We estimate the coherence limit to the RB fidelity to be: where t π sb,k and T 1,k are the iSWAP durations and lifetimes, respectively, of mode k, and p q is the RB fidelity of the transmon. The experimentally measured fidelities approach these coherence limits ( Figure 3 of the main text).

Supplementary Note 13 -Two mode gates using sideband transitions
The level diagram describing the relevant multimode states and transitions involved in the CZ gate are shown in Supplementary Figure 12a. We break the 2π |e − |f sideband pulse used to provide a conditional phase between the transmon and the second mode (|e1 → − |e1 , see Figure 4a of main text) into two π pulses. Control of the relative phases between these pulses allows for the correction of additional phase errors arising from the dispersive shift and the realization of an arbitrary controlled phase gate.
We obtain a CNOT gate by inserting an |e − |f transmon charge π pulse (π ef q ) between the two |e − |f sideband iSWAP pulses. This allows for mapping |e0 to |e1 (and vice versa) via the |f 0 (qubit-mode CNOT), which again becomes a mode-mode CNOT gate when sandwiched with two |g − |e sideband iSWAP pulses. The pulse sequence, energy level diagram, and relevant transitions for the CNOT gate are shown in Supplementary Figure 12b.
Slight modifications of these pulse sequences allow the realization of other two-mode gates such as the mode-mode CY and SWAP gates. The pulse sequences (without corrections from the dispersive shift) for realizing these two-mode gates are summarized in Supplementary Table 1.

Supplementary Note 14 -Phase error corrections to two-mode gates
The previous discussion of two-mode gates only involved resonant first-order sideband transitions and ideal transmon charge drive pulses. This idealized description is corrected by additional terms in the Hamiltonian of equation (1) (1) (2) (4) Supplementary Figure 12 | a and b, Energy level diagrams showing the multimode states and pulses involved in CZ and CNOT gates, respectively. The e iπ phase factor arising from sideband SWAP operations between |e01 ←→ |f 00 is modified by the dispersive shift (χ) of |e01 level (red). The additional phase arising from the dispersive shift is corrected by adjusting the phase of an ef sideband pulse on mode 1. CY gates are realized by adjusting the phase of the π ef q pulse.

Two-mode gate Pulse Sequence
CZ j,k π ge sb,j + π ef sb,k + π ef sb,k + π ge sb,j (φ = π) CX j,k π ge sb,j + π ef sb,k + π ef q,y + π ef sb,k + π ge sb,j (φ = π) CY j,k π ge sb,j + π ef sb,k + π ef q,x + π ef sb,k + π ge sb,j (φ = π) SWAP j,k π ge sb,j + π ef sb,k + π ge sb,k + π ef sb,k + π ge sb,j Supplementary  For the case of the dispersive shift, the corrections to the target unitaries depend on the quantum state of the multimode memory and result in a transmon-mode ZZ error. If we ignore photons in the rest of the memory, the effect of the dispersive shift on ther modes involved in a two-mode entangling gate can be inferred from the energy level diagram of Supplementary Figure 12. The dispersive shift causes the |e10 and |e01 levels to be shifted (red) from their bare values (black). As a result, sidebands resonant with |e00 ←→ |g10 are off-resonant from |e01 ←→ |g11 .
The dispersive shift results in a phase and population error in the |e01 state. For the dispersive shifts χ and iSWAP times t π used in this work, this population error is ∝ (χt π ) 2 and, at worst, ∼ 5% over the course of a CZ gate. This error is uncorrected and factors into the total gate error. The phase error on other hand is ∝ (χt π ) and results in a more significant effect (see Supplementary Figure 13a). Given that the |e01 state affected by the dispersive shift is selectively addressed by the |e − |f sideband pulses used in the gate (see Supplementary Figure 12 is calibrated and corrected by adjusting the relative phase between these pulses. The experimental protocols used for the gate calibration and phase error correction for the CZ gate are described in Supplementary Note 15. The state dependent phases arising in the gate can be calculated by considering the effective Hamiltonian of equation (14) in the 8 × 8 subspace of levels relevant for the gates and shown in Supplementary Figure 12: Here, the multimode state is labeled |n j , n k and the phases of equation (16) have been absorbed into the g's and Ω's (which are time dependent), i.e., The |e −|f sideband pulses act only on one transition and are unaffected by the state dependent shift when considering only two modes. We chose the |e − |f sideband freqeuncy to be resonant with the |f 00 and the dispersively shifted |e01 level. In the rotating frame of equation (26), this corresponds to |e − |f first-order sidebands acquiring the following time-dependence: g i,ef is proportional to the envelope of the |e − |f sideband pulse, and δ k and δ j are the dispersive shifts of |e0 j 1 k and |e1 j 0 k respectively. We compute the action of the CZ and CNOT gate sequences by evolving the Hamiltonian above, with time dependent coefficients and phases as per Supplementary Table 2. In these pulse sequences, only one of the drive terms is on at any given time and the corresponding unitaries obtained upon integration of the Schrödinger equation are generalizations of those in Supplementary Equation (18), with corrections arising from the dispersive shift. The effective unitary thus

Supplementary Note 15 -CZ gate phase calibration sequences
Here, we describe protocols used to calibrate and correct each of the additional phases arising in the CZ gate. The phases of the iSWAP pulses used in the CZ gate are defined below, where φ 1,2 are the controlled phases: j is the control mode and k is the target mode of the CZ gate, with the states labeled |n j , n k , andπ indicating iSWAP pulses for which the DC-offset σ z error is corrected within the pulse. From equation (29), we see that only two relative phase adjustments (φ a − φ d and φ b − φ c ) are required to correct the dispersive shift error. Here, we adjust these relative phases by controlling φ 1 = φ b and φ 1 = φ d , while leaving φ a and φ c fixed. We measure each phase error through Ramsey experiments with initial states that are appropriate superpositions of the basis states, as indicated in Supplementary Table 3 The phase error from the dispersive shift is obtained by preparing the system in the state |ψ p = |10 + |11 . We measure the relative phase between the basis states after applying the CZ gate, using the sequence in Figure 14a. The dispersive shift results in the |11 acquiring an additional phase in the preparation (φ p ), gate (φ g ) and measurement segments (φ m ). Similar additional phases also accrue during the gate (φ g ) and the measurement (φ m ) segments. We sweep the phase (δφ) of the first |e − |f sideband pulse of the CZ gate (see Supplementary Equation (31)). The phase that maximizes the final measured transmon population provides the total added phase, φ ds t = φ p + φ g + φ m . CZ j,k (φ ds t , 0) is a combination of an ideal CZ gate and the Cφ gate that cancels additional phases arising from the dispersive shift over the entire sequence. We isolate the state-dependent phase error arising only during the CZ gate by adding a second Cφ gate to the previous sequence, as shown in Supplementary Figure 14b (δφ) of the (first) |e − |f sideband pulse of the second Cφ gate, with φ 1 = φ t ds for the first Cφ gate. Given the same preparation and measurement sequences, the SPAM phases are corrected by construction by the first Cφ gate. We find the phase δφ that minimizes the population of the transmon, thus realizing a CZ gate that flips the sign of the |11 state. The CZ k,j (φ g , 0) gate therefore is corrected for phases from dispersive shifts occurring during the gate sequence.
We obtain a fully corrected CZ gate by correcting the relative phase between the {|00 , |01 } and {|10 , |11 } manifolds. These state manifolds have a relative phase resulting from the transmon frequency DC-offset occurring during the |e − |f sidebands of the CZ gate. We correct this additional phase by adjusting the phase of the final |g − |e sideband pulse of the CZ gate (φ 2 in equation (31)), using the experimental sequence of Supplementary  Figure 14c. The resulting CZ j,k (φ ds g , φ DC ) gate is therefore corrected of errors from dispersive shifts and qubit fluxmodulation DC-offsets.
The phase error resulting from dispersive shifts due to off-resonant first-order sidebands are measured by acting the CZ gate on the |00 + |01 state. The CZ gate nominally does not change this state, and we correct this phase error using subsequent qubit pulses as shown in Supplementary Figure 14d. This phase is significant only for gates between modes with spectral spacing near the anharmonicity of the transmon.

Supplementary Note 16 -Two-mode quantum state tomography
Reconstructing the density matrix of an arbitrary two-qubit state requires the measurement of all possible two-qubit correlations { XI , XX . . . ZZ }, i.e.; These correlators can be measured through Ramsey interferometry, as described in the main text [14]. We equivalently measure all the necessary correlations with the aid of the single and two-mode gate operations prior to measuring the state of the transmon. A sideband iSWAP pulse (π sb ) on the |g − |e transition, along with single qubit rotations alone can be used to measure all single-mode correlators The entanglement information is present in two-mode correlators, We measure these correlators by acting two-mode gates before measuring a single-mode correlators. For instance, the XX correlator of a given state (|ψ i ) is measured by acting CX gate prior to the measurement of XI . In the Heisenberg picture [15], the transformation is shown below: Here |ψ i is the two-mode state to be measured and |ψ f = U CX |ψ f is the state obtained following action of the CX gate. A summary of pulse sequences used for the measurement of each of the correlations required for two-mode tomography are shown in Table 4. We extract the correlators and construct the density matrix of the two-mode state from the measured transmon  population P ij at the end of the sequence for each correlator C ij using; In general, fast measurement and reset [16] of the transmon would allow us to perform sequential measurements of two-mode correlations using the transmon without requiring mode-entangling gate operations. For each mode, we would map the mode state to the transmon with an iSWAP, measure the transmon, and reset it to the ground state. The transmon state could be reset with an iSWAP back to the measured mode or to an auxiliary mode. The transmon can subsequently be used to measure the next mode. Additionally, we can perform Wigner tomography [17] of the multimode chain through direct measurements of the multimode fields and parametric amplification. These techniques pose more stringent conditions on the measurement fidelity and speed and are beyond the scope of this work.
We account for spurious phases arising in the Bell state tomography sequence by varying the phase of the final sideband pulse used in each correlator measurement. The results of such phase sweeps for the |Φ + and |Ψ + Bell states are shown in Supplementary Figure 20a and b. We note that for these states, the only correlators that are functions of the final sideband phase are XX, XY, YX, YY and the |Φ + and |Ψ + Bell states give opposite answers for measurements of the two-mode parity ZZ. The optimal sideband phase that accounts for the additional phases shifts are indicated by the dashed lines. We extract the density matrices from the results of these measurements using equation (34). The state fidelities for the two states are calculated from the overlap of the ideal ρ id and measured ρ m density matrices: F p = Tr(ρ id ρ m ).

Supplementary Note 17 -Process tomography of two-mode gates
Process tomography of a two-qubit gate consists of quantum state tomography after acting the gate on a set of 16 linearly independent input states that form a basis for representing an arbitrary two-qubit density matrix [18]. Process tomography of two-mode gates therefore consists of a set of 240 measurement sequences of the form shown in Supplementary Figure 16.
Supplementary Figure 16 | Process tomography sequence for two-mode gates, broken down into preparation (red), gate (green), and tomography (blue) segments. For the preparation sequence, we use qubit rotations R1,2 = I, Ry π 2 , Rx π 2 , Rx (π) and DC-offset corrected sideband iSWAP pulses with an additional − π 2 phase, such that the target multimode state at the end of the preparation sequence is |ψi = R1 ⊗ R2 |0j0 k . We measure the density matrix of the gate outputs for given input density matrices ρi = |ψi ψi| using the state tomography protocols of Supplementary Note 16 and the sequences in Supplementary Table 4, corresponding to σ1,2 = {I, X, Y, Z} in the sequence shown above. We note that the iSWAP gate acts on mode k for some of the correlators, and that the tomography sequences that measure single-mode correlators have no additional two-mode gate, corresponding to σ1 = 1. The qubit and iSWAP operations that are indicated by the dashed lines have errors arising from the dispersive shift.
The gate calibration protocols described in Supplementary Note 15 for the CZ gate, and analogous protocols for the CX gate, correct phase errors due to dispersive shifts and the transmon DC-offset from flux modulation during the gate. We additionally correct errors arising from the dispersive shift during the state preparation and tomography (SPAM) segments of the various process tomography sequences. These errors occur in the qubit and iSWAP operations indicated by dashed lines in Supplementary Figure 16. The dispersive shift causes amplitude and phase errors in the transmon and sideband pulses. We again correct only phase errors to first-order in χ/Ω sb . These controlled-phase errors can be formally incorporated as CΦ gates at the end of the preparation sequence (CΦ p ) and prior the to measurement sequence (CΦ m ). These additional gates are concatenated into the gate and tomography sequences as shown in Supplementary Figure 17. CZ and Cσ 1 are chosen to give phase-corrected CZ and Cσ 1 gates when concatenated with CΦ p and CΦ m , respectively. The preparation error is corrected through an added phase (φ p ) in the first |e − |f sideband of the first gate, while the tomography error is corrected through an added phase (φ m ) in the second |e − |f sideband of the last gate of the sequence. We thereby correct the sequence to first-order in the dispersive-shift error. The sideband phases are chosen in this manner in order to correct errors in both the |10 and |11 states (see equation (29), (30)). The phase errors depend on the duration of the qubit and sideband pulses used in the sequence. In the absence of loss, they can be calculated based on the dispersive shift and pulse shapes.
We calibrate the additional phase errors through process tomography of the Identity (1) gate (idling for 10 ns). We find the optimal phases by sweeping the added controlled-phases of the Cφ and Cσ gates, and comparing results for corresponding correlators with and without CX/CY gates (such as XX and XI, respectively) as shown in Supplementary Figure 18.
This scheme allows us to isolate state preparation and measurement errors (φ p and φ m ). In the Heisenberg picture, working backward from the transmon measurement, we consider how correlators are modified by the dispersive shift and the correcting two-mode gates. As an example, for the prepared state |xx (R 1,2 = Y π 2 ), the expected values of the correlators XI and XX are: where φ p and φ m are the phase errors of the |11 state (relative to the other computational basis states) in state preparation and measurement, respectively. Finding and correcting φ p and φ m amounts to choosing φ c,1 and φ c,2 such that XI = XX = 1.
The additional phases only depend on the shape of the qubit and sideband pulse waveforms. As a result, we can calibrate φ p and φ m for all 240 sequences using a total of 13 unique experiments. We can then extract the full process matrix by measuring at the optimal angles obtained from the calibration experiments. We check that the validity of the calibrations by also additionally sweeping the phase of the final sideband pulse. In order to reduce SPAM error from decoherence, we combine the state preparation and measurement correction gates (as indicated by the arrows in Supplementary Figure 18) during process tomography of the 1 and CZ gates, noting that CΦ φc,1 commutes with both of them.
We perform process tomography of the CZ gate by inserting it in place of the 1 in Supplementary Figure 18, after calibrating the tomography axes. A two-mode gate is fully characterized by the completely positive map E; A m =B i ⊗B j , withB i ∈ {I, X, Y, Z}, forms a basis of operators acting on a two-mode state ρ. χ mn is the process matrix characterizing the two-mode gate, and is extracted from the measured output density matrices ρ out j for 16 linearly independent input density matrices (ρ j ) as shown below: ⇒ λ jk = Tr[ρ k ρ out j ] = mn β mn jk α mn .
Supplementary Equation (39) is directly inverted to obtain the process matrix α mn . We do not impose the completeness condition, mn χ mnÂmÂ † n = 1 as a constraint. This constraint arises from the probabilities of states in the relevant two-mode space summing to 1. This is not necessarily the case when there are several memory modes. The process fidelities are extracted from the measured (χ m ) and ideal process matrices (χ t ) using F p = Tr [χ m χ t ].

a b
Supplementary Figure 19 | a, Experimentally measured correlators after correcting for phase errors arising during state preparation and measurement for process tomography of the CZ gate between mode j = 2 and k = 6. b, Process matrix extracted from the resulting measurements by inverting equation (39).

Supplementary Note 18 -Preparation of entangled states
We use a slight variant of the scheme described in the main text to prepare multimode entangled states. We prepare Bell states [19] between two modes with the following protocol: starting with the transmon in its excited state, we swap half of the excitation via a sideband pulse ( √ iSWAP) to the first mode. This creates a |Ψ + Bell state between the first mode and the transmon: We can rotate this state into the |Φ + Bell state by flipping the transmon state via its charge bias control: |Φ + = 1 √ 2 (|g000...0 + |e100...0 ) .
For either of the states in equations (40) and (41), we can create the corresponding Bell states between two arbitrary modes by simply swapping the transmon state to another mode. To extend this protocol and create a n-mode Greenberger-Horne-Zeilinger (GHZ) state [20], we again utilize the anharmonicity of the transmon and map the population of |e to |f . This allows us to transfer this excitation to the second mode via a sideband of the |e − |f transition without disturbing the population in the ground state: These last two pulses can be repeated for each of the remaining eigenmodes before finally swapping the transmon back to the nth mode to complete the GHZ state: |ψ GHZ = 1 √ 2 (|g000...0 + |g111...1 ) .
We measure 3-qubit correlators by including an additional two-mode gate to the pulse sequences in Table 4 for measuring 2-qubit correlations. The XXX correlator is measured through the experimental sequence shown in Supplementary Figure 21a. We reduce tomography error by compiling consecutive |g − |e iSWAP pulses between successive CNOT gates, as well as the CNOT gate and the final iSWAP pulse. Additionally, we mitigate loss arising from transmon T * 2 by inserting a π pulse prior to the final π/2 qubit pulse, analogous to a spin-echo. The remaining correlators comprising W are measured by changing the phases of the CX gates and the final π/2 pulse. The results of the measurement of W, as a function of the time τ between the π and π/2 pulses are shown in Supplementary Figure 21b. The echo results in a revival of W at a time τ ∼ 530 ns. At the point indicated by the red star, W = 2.35 ± 0.04, with the indicated error that arising from readout visibility (4×visibility error). While the echo mitigates errors from transmon dephasing, the measured value of W is a lower bound that includes loss from tomography. The tomography sequence prior to the echo includes 5 sideband and 2 qubit π pulses. Accounting for the loss in these pulses based on the RB data for the |g − |e iSWAPs and qubit pulses gives W > 2.65. This is still a conservative estimate for the witness, since the |e − |f iSWAPs have a lower fidelity than the corresponding |g − |e iSWAPs. a π ge sb,j + π ef sb,k + π ef q + π ef sb,k + π ef sb,l + π ef q + π ef  Figure 21 | (a) Gate sequence for measurement of the XXX correlator with a transmon echo. As indicated by the pulse sequence, consecutive |g − |e iSWAP pulses on the same mode in the red shaded region have been compiled and removed. (b) GHZ entanglement witness measured as a functon of time τ following the echo π pulse. The blue shaded region indicates the EPR bounds for W, corresponding to the states being biseperable. The memory modes used in the measurement are j = 6, k = 8 and l = 2.