Hardware-efficient quantum random access memory with hybrid quantum acoustic systems
Connor T. Hann Chang-Ling Zou Yaxing Zhang Yiwen Chu Robert J. Schoelkopf Steven M. Girvin Liang Jiang





Departments of Applied Physics and Physics, Yale University, New Haven, Connecticut 06511, USA


Key Laboratory of Quantum Information, CAS, University of Science and Technology of China, Hefei, China.


Departments of Applied Physics and Physics, Yale University, New Haven, Connecticut 06511, USA


Department of Physics, ETH Zürich, 8093 Zürich, Switzerland


Departments of Applied Physics and Physics, Yale University, New Haven, Connecticut 06511, USA


Departments of Applied Physics and Physics, Yale University, New Haven, Connecticut 06511, USA


Departments of Applied Physics and Physics, Yale University, New Haven, Connecticut 06511, USA


Hybrid quantum systems in which acoustic resonators couple to superconducting qubits are promising quantum information platforms.
High quality factors and small mode volumes make acoustic modes ideal quantum memories, while the qubit-phonon coupling enables the initialization and manipulation of quantum states. 
We present a scheme for quantum computing with multimode quantum acoustic systems, and based on this scheme, propose a hardware-efficient implementation of a quantum random access memory (QRAM). 
Quantum information is stored in high-Q phonon modes, and couplings between modes are engineered by applying off-resonant drives to a transmon qubit.
In comparison to existing proposals that involve directly exciting the qubit, 
this scheme can offer a substantial improvement in gate fidelity for long-lived acoustic modes. 
We show how these engineered phonon-phonon couplings can be used to access data in superposition according to the state of designated address modes—implementing a QRAM on a single chip. 

Hardware-efficient quantum random access memory with hybrid quantum acoustic systems
    Liang Jiang
    December 30, 2023
====================================================================================








Introduction.—The coupling of superconducting qubits to microwave resonators, termed circuit quantum electrodynamics (cQED) <cit.>, constitutes one of today's most promising quantum computing architectures. 
Microwave modes provide good quantum memories <cit.>, while superconducting nonlinearities enable the initialization <cit.>, manipulation <cit.>, readout <cit.>, and protection <cit.> of quantum states encoded in microwave photons. However, long microwave wavelengths could potentially limit the scalability of cQED systems. On-chip resonators face trade-offs between compactness and quality factor <cit.>, and microwave modes with millisecond coherence have thus far only been demonstrated in large 3D cavities <cit.>.

Recently, coherent couplings between superconducting qubits and acoustic resonators have been demonstrated in a remarkable series of  experiments <cit.>.
These so-called circuit quantum acoustodynamic (cQAD) systems (Fig. <ref>) possess many of the  advantageous properties of cQED systems, e.g., superconducting qubits can be used to generate arbitrary superpositions of acoustic Fock states <cit.> and perform phonon-number resolving measurements <cit.>. Yet relative to electromagnetic modes, acoustic modes can provide dramatic benefits in terms of size and coherence times. The velocities of light and sound differ by five orders of magnitude, and short acoustic wavelengths enable the fabrication of ultra-compact phononic resonators <cit.>. Furthermore, acoustic modes can be exceptionally well-isolated from their environments—quality factors in excess of 10^10 were recently demonstrated in GHz frequency phononic crystal resonators <cit.>.
A variety of applications for such platforms have been proposed, including quantum transduction <cit.>, entanglement generation <cit.>, and quantum signal processing <cit.>, but surprisingly the direct use of cQAD systems for quantum computing has received relatively little attention, with the notable exception of Ref. <cit.>.



In this work, we propose a hardware-efficient and scalable quantum computing architecture for multimode cQAD systems. 
Quantum information is stored in high-quality acoustic modes, 
and interactions between modes are engineered by applying off-resonant drives to an ancillary superconducting transmon qubit.
During these operations, the transmon is only virtually excited, so the effects of transmon decoherence are mitigated. This is a crucial property, since the transmon's decoherence rate can exceed that of the phonons by orders of magnitude.
In comparison to existing proposals that involve directly exciting the transmon <cit.>, this virtual approach can offer substantial improvement in gate fidelity for long-lived phonons. This scheme is also directly applicable to multimode cQED <cit.>. 

Furthermore, to demonstrate the benefits that the proposed cQAD architecture affords in hardware efficiency, we propose an implementation of a quantum random access memory (QRAM) <cit.>. A classical RAM is a device that can query a database. Given an address j as input, the RAM outputs the element D_j stored at position j in the database. 
Analogously, a QRAM is a device that, when provided with a superposition of addresses, returns a correlated superposition of data,

    ∑_j=1^Nα_j|j⟩_a|0⟩_b∑_j=1^Nα_j|j⟩_a|D_j⟩_b,

where the subscripts a and b denote the address and output qubit registers, respectively. 
The ability to perform such queries efficiently, i.e. in log N time, is a prerequisite for a variety of quantum algorithms that provide speedups over their classical counterparts <cit.>. 
However, building a QRAM is highly non-trivial; even a small-scale QRAM of the sort described in Ref. <cit.> has yet to be experimentally demonstrated. 
One major challenge is that, to query a database of size N, a QRAM requires O(N) quantum resources <cit.>. 
Hardware-efficiency is thus crucial for QRAM queries of large datasets, and in our implementation this efficiency is enabled by the on-chip integration of superconducting circuits with compact acoustic resonators.
Out proposal both reduces physical resource requirements and provides a roadmap for a near-term demonstration of QRAM.

Note that “random access quantum memories” <cit.>—memories that can take a single classical address j as input and return a corresponding qubit |ψ_j⟩—cannot generally perform operation (1) and should not be confused with QRAM, which can use superpositions of addresses to query multiple memory elements.






Quantum computing in cQAD.—In multimode cQAD, a transmon qubit is piezoelectrically coupled to a collection of acoustic modes. These modes can be supported in 
bulk acoustic wave (BAW) <cit.> or surface acoustic wave (SAW) <cit.> resonators, or in an array of phononic crystal (PC) resonators <cit.> (Fig. <ref>). 
Quality factors of ≈ 10^5, 10^8, and 10^10 have been measured at GHz frequencies in SAW <cit.>, BAW <cit.>, and PC resonators <cit.>, respectively, and 
the transmon can be simultaneously coupled to many high-Q modes on a single chip <cit.>.
These systems can be described by the Hamiltonian 

    H    =ω_qq^†q-α/2q^†q^†qq
       +∑_k(ω_km_k^†m_k+g_k q^†m_k+ g_k^* q m_k^†)+H_d

Here, q and m_k denote the annihilation operators for the transmon and phonon modes, respectively. The transmon is modeled as an anharmonic oscillator with Kerr nonlinearity α and is coupled to the k^th phonon mode with strength g_k (typically a few MHz <cit.>). In combination with external drives on the transmon H_d = ∑_j Ω_j q^† e^-iω_j t+ H.c., this coupling provides the basic tool to initialize, manipulate, and measure phononic qubits <cit.>. For example, itinerant photon-encoded qubits sent to the system can be routed into a particular phonon mode via pitch-and-catch schemes <cit.>.

Interactions between phonon modes can be engineered by applying off-resonant drives to the transmon, and we use these interactions to implement a universal gate set for phononic qubits. 
The main idea is that the transmon's Kerr nonlinearity enables it to act as a four-wave mixer <cit.>, so phonons can be converted from one frequency to another by driving the transmon.
For example, phonons can be converted from frequency ω_A to ω_B by applying two drive tones whose frequencies ω_1,2 satisfy the resonance condition ω_2-ω_1=ω_B-ω_A, see Fig. <ref>. This driving gives rise to an effective Hamiltonian  H=g_v^(1)m_Am_B^†+H.c., where g_v^(1)=-2αg_A/δ_Ag_B^*/δ_BΩ_1^*/δ_1Ω_2/δ_2(1-β^(1)). 
Here, δ_j ≡ω_j - ω_q, and β^(1) is a drive-dependent correction (see supplementary material <cit.> for derivations). 
Evolution under this coupling for a time π/2g_v^(1) implements a  gate, which exchanges the states of modes m_A and m_B, while evolution for a time π/4g_v^(1) implements a 50:50 beamsplitter operation <cit.>. 




Three-mode interactions can be similarly engineered (Fig. <ref>). 
Applying a single drive tone with frequency ω_1=ω_A+ω_B-ω_C gives rise to the effective Hamiltonian H=g_v^(2)m_Am_Bm_C^†+H.c., where g_v^(2)=-2αg_A/δ_Ag_B/δ_Bg_C^*/δ_CΩ_1^*/δ_1(1-β^(2)) <cit.>. 
This three-mode interaction can be used to implement a controlled phase () gate for qubits encoded in the |0,1⟩ phonon Fock states <cit.>. 
To perform a  gate between qubits in modes A and B, mode C is used as an ancilla and initialized in |0⟩. 
Evolving for a time π/g_v^(2) then enacts the mapping |110⟩_ABC→|001⟩→ -|110⟩, while leaving all other initial states unaffected. 
The state |11⟩_AB acquires a relative geometric phase, thereby implementing the  gate. 

A variety of other operations can be similarly implemented. Single- and two-mode squeezing can be implemented by driving the transmon at appropriate frequencies, and phase shifts can be imparted by tuning the drive phases during  operations. Together, these two- and three-mode interactions are universal <cit.>. In the remainder of this work, however, we focus on the beamsplitter, , and  operations, as these are the only operations we require to implement a QRAM.

Note that in BAW and SAW resonators, phonon mode frequencies are approximately uniformly spaced, i.e. ω_j+1-ω_j=ν, where ν is the free spectral range. 
This uniform spacing can lead to problematic degeneracies in the resonance conditions above. Nonuniform mode spacing is thus necessary to enable selective coupling, and in <cit.> we describe several ways to engineer nonuniformity in BAW and SAW systems.
As shown in Fig. <ref>, one approach is to couple the phonons to an external mode, such as a microwave resonator, so that the resulting hybridization shifts mode frequencies <cit.>. 
In <cit.>, we also introduce a metric, Δν, to quantify the nonuniformity. Roughly speaking, Δν is the scale at which the mode spacing varies.






Gate fidelities.—During the gates described above, the transmon is never directly excited; instead, it is only virtually excited, so infidelity attributable to transmon decoherence is suppressed.
These virtual gates can thus provide great advantage in cQAD systems, where transmon decoherence is likely to be the limiting factor. 
This is in contrast to existing proposals <cit.>, in which gates between resonator mode qubits are implemented by swapping information directly into the transmon using resonant interactions of the form g_d(q^†m+qm^†), which can be engineered, e.g., by modulating the transmon's frequency. 
In the following, we compare the predicted fidelities of the virtual gates proposed here and the direct gates considered in Refs. <cit.>.

In a multimode architecture, there exists a fundamental tradeoff between decoherence and spectral crowding. Slower gates are more prone to decoherence, while faster gates have reduced frequency resolution and can disrupt other modes. 
The infidelities of the direct and virtual gates, respectively 1-ℱ_d and 1-ℱ_v, can be approximated as a sum of contributions from these two effects <cit.>,

    1-ℱ_d   ≈ (κ+γ)c_dπ/2g_d+(g_d/ν)^2,
    
    1-ℱ_v   ≈κ̅_γc_vπ/2g_v+(g_v/Δν)^2,

where κ and γ are the bare phonon and transmon decoherence rates, and c_d,v are constants accounting for the durations of each gate (c_v=1 for , and c_v=2 for , as these gates have durations π/2g_v and π/g_v respectively. As discussed in Ref. <cit.>, c_d=5 for  and c_d=4 for .) 

The first terms in Eq. (<ref>) and Eq. (<ref>) account for decoherence. During direct gates, information spends roughly equal time in the phonon and transmon modes, so the total decoherence rate is κ+γ. 
During virtual gates, the total decoherence rate, κ̅_γ, is κ̅_γ=(κ_γ^A+κ_γ^B) for , and κ̅_γ=(κ_γ^A+κ_γ^B+κ_γ^C)/2 for .
Here, κ_γ^j=κ+γ(g_j/δ_j)^2(1+β^(γ)) denotes the dressed decay rate of mode j, which includes a contribution from the inverse Purcell effect <cit.> and a drive-dependent correction β^(γ) <cit.>. The second term in each expression accounts for spectral crowding. 
The probability of accidentally exciting another mode scales as (g_d/ν)^2 in the direct case, and as (g_v/Δν)^2 in the virtual case. We note that the infidelities (<ref>,<ref>) are defined with respect to the ideal two-qubit unitaries; see <cit.> for derivations. 

The competition between decoherence and spectral crowding results in an optimal coupling rate <cit.>. 
By adjusting the drive strengths, g_d,v can be tuned to their respective optima. The optimal infidelities are 

    1-ℱ_d   ≈3/2[ c_d π(κ+γ)/√(2)ν]^2/3,
    
    1-ℱ_v   ≈3/2[c_v πκ̅_γ/√(2)Δν]^2/3.

While transmon and phonon decoherence contribute equally to 1-ℱ_d, transmon decoherence only makes a small contribution to 1-ℱ_v via the inverse Purcell effect, wherein γ is suppressed by a factor of (g/δ)^2≪1. The virtual gates can thus be expected to attain higher fidelities when there is a large disparity between γ and κ, i.e. for sufficiently long-lived phonon modes. Indeed, ℱ_v>ℱ_d whenever κ_γ≲(κ+γ)Δν/ν, provided the optimal coupling rates can be reached.



In Fig. <ref>, we plot the optimal infidelities of direct and virtual gates as a function of κ and γ for currently feasible experimental parameters. 
The comparison reveals that virtual gates can be performed with high fidelity (>99%) given long-lived phonons, and that virtual gates attain higher fidelities than direct gates in the same regime. 
Indeed, realistic improvements in phonon coherence are likely to bring near-term devices into this ℱ_v ≫ℱ_d regime (Fig. <ref>c,f). 

We briefly note other factors relevant to the comparison of direct and virtual gates.
Multi-phonon encodings. Direct gates require that qubits be encoded in the |0,1⟩ phonon Fock states, while virtual operations are compatible with multi-phonon encodings, including some bosonic quantum error-correcting codes <cit.>. 
Parallelism. Direct gates must be executed serially, while virtual gates can be executed in parallel by simultaneously applying the requisite drives (though care should be taken to ensure that the additional drives do not bring spurious couplings on resonance).
Speed. Virtual gates are inherently slower than direct gates, with realistically attainable coupling rates g_v/2π∼ 10-100 <cit.>. 





QRAM Implementation.—To illustrate the advantages of cQAD systems, we propose an implementation of a QRAM <cit.>. As defined by Eq.(<ref>), a QRAM is a device which can query a database with an address in superposition. 
The ability to perform such queries efficiently is a prerequisite for a variety of quantum algorithms, including Grover's search <cit.>, matrix inversion <cit.>, and various proposals in the field of quantum machine learning <cit.>. 
While demanding hardware and connectivity requirements have thus far precluded an experimental demonstration of a QRAM, our proposed cQAD implementation is naturally hardware-efficient. Indeed, a small-scale cQAD QRAM can be implemented with a single multimode resonator.



The elementary building block of our QRAM implementation is a quantum router, shown in Fig. <ref>(a). 
The router directs an incoming qubit into different output modes conditioned on the state of a routing qubit. When the routing qubit is in state |0⟩(|1⟩), an incoming qubit |ψ⟩ in the top mode is swapped to the left(right) mode. The routers are implemented using the operations described above: the routing circuit contains a  and a controlled- gate, the latter of which is implemented using  and beamsplitter operations (Fig. <ref>b).  

To implement a QRAM, a collection of routers is arranged in a binary tree, with the outputs of routers at one level acting as inputs to routers at the next (Fig. <ref>c). To query the database at the bottom of the tree, qubits from the address register are routed sequentially into the tree, with earlier address qubits controlling the routing of later ones in a “bucket-brigade” scheme <cit.>. 
A “bus” qubit then follows the path paved by the address qubits and extracts the data, after which it is routed back out of the tree and into the output register.
Finally, to disentangle the address and routers, the address qubits are routed back out of the tree. 
Since all routing operations are quantum-controlled, preparing the address register in superposition allows access to the data in superposition, thereby implementing operation (<ref>). Further details are provided in <cit.>.

We highlight three appealing properties of this cQAD-based QRAM.

    
  * Hardware-efficiency. Hundreds of phonon modes can simultaneously couple to a transmon on a single chip <cit.>. Thus, the hardware and fabrication cost of a cQAD-based QRAM can be drastically reduced in comparison to cavity- <cit.> or circuit-QED <cit.> implementations. 
    
  * Scalability.  
    It is not necessary to control all routing through a single transmon; since only adjacent routers are coupled, different regions of the tree can be controlled and implemented independently. 
    For example, the QRAM can be built out of several modules, where each module comprises a group of routers controlled by one transmon. 
    The phononic modes in each module could be supported in physically separate resonators, or multiple transmons could be simultaneously coupled to the same multimode resonator to give access to a large bandwidth of modes, potentially spanning several GHz <cit.>. 
    
  * Error resilience. Because our implementation follows the bucket-brigade model, it inherits a favorable log N error scaling <cit.>. In particular, the scaling argument of Ref. <cit.> directly applies to phonon-loss errors: the query infidelity scales as 1-ℱ∼εlog N, where ε is the phonon loss probability. Remarkably, one can show that the infidelity scales logarithmically for arbitrary independent, incoherent errors, such as phonon loss, dephasing, and heating <cit.>. 


Discussion.—We have proposed a quantum computing architecture for multimode cQAD and an implementation of a QRAM based on it. 
The implementation is hardware-efficient, owing to the compactness of multimode cQAD systems that is enabled by small acoustic wavelengths. We emphasize that hardware efficiency is not only crucial for scaling to large system sizes, but that it is also particularly advantageous for near-term experiments. Indeed, a small-scale QRAM can be implemented even with just a single multimode resonator. In the long term, the use of bosonic quantum error correcting codes <cit.> and compatible logical gates <cit.> to implement a QRAM that is both fault-tolerant and hardware-efficient is an intriguing direction for future research. These ideas can also be directly applied to multimode cQED.

To be viable, our scheme requires long phonon coherence times (1/κ≫ 1/g_v). 
Though loss due to intrinsic material processes like phonon-phonon scattering or two level systems limit phonon coherence, such mechanisms should not prevent access to this regime. 
Indeed, both BAW and PC quality factors can approach 10^10 before encountering such limits <cit.>, corresponding to κ/2π∼ 1Hz at GHz frequencies.
Additionally, intrinsic sources of phonon dephasing are not expected <cit.>. 
Phonon decoherence in current cQAD experiments is thus likely dominated by extrinsic mechanisms that can be mitigated with improved fabrication techniques, though the extent to which the coupling to superconducting circuits may limit phonon coherence in cQAD is an important open question. 

Acknowledgements.—We thank Vijay Jain, Prashanta Kharel, Patricio Arrangoiz-Arriola, Alex Wollack, Hong Tang, Peter Rakich, and Oskar Painter for helpful discussions. C.T.H. acknowledges support from the NSF Graduate Research Fellowship Program (DGE1752134). We acknowledge support from the ARL-CDQI (W911NF-15-2-0067, W911NF-18-2-0237), ARO (W911NF-18-1-0020, W911NF-18-1-0212), ARO MURI (W911NF-16-1-0349), AFOSR MURI (FA9550-15-1-0015), DOE (DE-SC0019406), NSF (EFMA-1640959, DMR-1609326), and the Packard Foundation (2013-39273).

apsrev4-1




















Supplemental Material: Hardware-efficient quantum random access memory with hybrid quantum acoustic systems








§ VIRTUAL COUPLING RATES

In this section, we study the virtual coupling rates

    g_v^(1)    =-2αξ_1^*ξ_2λ_Aλ_B^*(1-β^(1)),
    
    g_v^(2)    =-2αξ_1^*λ_Aλ_B^*λ_C(1-β^(2)).
 
Below, we define the notation, derive these expressions, and discuss the importance of the corrections β^(1,2) for cQAD systems. Then, in order to verify the accuracy of these expressions, we compare them to numerical results obtained using the Floquet theory methods of Ref. <cit.>.  



 §.§ Derivation of the virtual coupling rates

To derive the expressions (<ref>) and (<ref>), we begin with the multimode cQAD Hamiltonian (Eq. (2) of the main text) and perform a unitary transformation defined by U_1 = exp i H_0 t, where H_0 = ω_q q^† q+∑_k ω_k m_k^† m_k. Thus,

    H=∑_j(Ω_jq^†e^-iδ_jt+H.c.)+∑_k(g_km_kq^†e^-iδ_kt+H.c.)-α/2q^†q^†qq,

where δ_k = ω_k-ω_q is the detuning of the k^th phonon mode, while δ_j = ω_j-ω_q and Ω_j are the detuning and the strength of the j^th drive tone, respectively.
In the spirit of Ref. <cit.>,
we first perform unitary transformations to eliminate the qubit-phonon
couplings and drive terms then consider the effects of the anharmonicity. For notational convenience, we introduce the dimensionless
parameters λ_k≡ g_k/δ_k and ξ_j≡Ω_j/δ_j.
To leading order in λ_k≪1, the unitary that eliminates
the couplings is U_2=exp∑_k(λ_k^*m_k^†qe^iδ_kt-H.c),
and that which eliminates the drives is U_3=exp∑_j(ξ_j^*qe^iδ_jt-H.c).
The combined effect of these two transformations is to enact the mapping
q→ q+∑_jξ_je^-iδ_jt+∑_kλ_km_ke^-iδ_kt≡ Q,
so that the Hamiltonian becomes

    H=-α/2Q^†Q^†QQ.

Note that we have neglected linear terms of the form
(Ω_j^*λ_km_ke^i(δ_j-δ_k)t+H.c.). This omission is justified in the RWA provided that |δ_j-δ_k|≫λ_kΩ_j,
i.e. that the drives are sufficiently far detuned from any modes in
which we are interested.
For the moment, we also neglect frequency (Stark) shifts of the phononic eigenmodes—we consider their effects in Sec. (<ref>). 

When two drive tones are applied whose frequencies satisfy the resonance condition ω_2-ω_1=ω_B-ω_A,
the Hamiltonian (<ref>) contains a resonant beamsplitter-type coupling, g_v^(1)m_Am_B^†+H.c.,
where

    g_v^(1)=-2αξ_1^*ξ_2λ_Aλ_B^*.

Similarly, when a single drive tone is applied with frequency[In this section we consider the case ω_A < ω_B < ω_1 < ω_C, which nicely highlights the similarities between g_v^(1) and g_v^(2). The derivations proceed analogously for other cases, such as the case of ω_C<ω_A < ω_B < ω_1 shown in Fig. 2 of the main text.] ω_1=ω_A+ω_C-ω_B,
the Hamiltonian contains a resonant three-mode coupling
g_v^(2)m_Am_B^†m_C+H.c.,
where

    g_v^(2)=-2αξ_1^*λ_Aλ_B^*λ_C.








 §.§ Corrections to the virtual coupling rates


The Hamiltonian (<ref>) contains many terms beyond
just the resonant terms discussed above (see Table <ref>). 
Most of these terms are rapidly-rotating and can be neglected in the RWA assuming dispersive coupling
(λ≪1) and weak drives (ξ≪1).
However, other terms can produce corrections β^(1,2)
to the coupling rates.
In this section, we first calculate these corrections to leading order
in λ and ξ. Then, we derive nonperturbative contributions
associated with the AC Stark shift. 

1.4




The leading order contribution to β^(1,2) is zeroth order
in both λ and ξ. The only terms in the Hamiltonian (<ref>)
which contribute to β^(1) and β^(2) at this order
are, respectively, 

    [-α(
    q^†2ξ_2λ_Am_A
    +q^†2ξ_1λ_Bm_B)
    e^-i(δ_B+δ_1)t+H.c.]
     and [-α(
    q^†2λ_Am_Aλ_Cm_C
    +q^†2ξ_1λ_Bm_B)e^-i(δ_B +δ_1)t+H.c.]

The corrections from these terms can be calculated via standard perturbation
theory, 

    β^(1,2)=α/δ_B+δ_1+α.

For the  operation, where the drives are far-detuned, this correction
is typically negligible. However, for the  operation, this correction
can significantly reduce the coupling rate since δ_1,δ_B
can be comparable to α. We note that the expression for β^(1)
matches the leading order expression derived in Ref. <cit.>. 

Contributions to β^(1,2) at higher orders in λ can
be neglected since we have assumed the dispersive regime, λ≪1.
Contributions at higher orders in ξ can be systematically calculated
with perturbation theory in principle, but such calculations quickly
become tedious. Here, we employ an alternative approach. We consider
the AC Stark shift type terms, -2α∑_j|ξ_j|^2Q^†Q,
and compute their contributions to β^(1,2) nonperturbatively
by working in a rotating frame. 

Let S denote the qubit's AC Stark shift. In the frame where the
qubit mode rotates at its Stark-shifted frequency, ω̃_̃q̃=ω_q+S,
the system Hamiltonian is

    H=-Sq^†q+∑_j[Ω_jq^† e^-iδ̃_jt+H.c.]+∑_k[g_k m_kq^†e^-iδ̃_kt+H.c.]-α/2q^†q^†qq

where δ̃=ω-ω̃_̃q̃. Performing unitary
transformations analogous to those above eliminates the coupling and
drive terms so that H=-SQ̃^†Q̃-α/2Q̃^†Q̃^†Q̃Q̃,
where Q̃=q+∑_jξ̃_je^-iδ̃_̃j̃t+∑_kλ̃_km_ke^-iδ̃_kt.
Here, ξ̃_j=Ω_j/δ̃_j and λ̃_k=g_k/δ̃_k.
The Stark shift terms can then be cancelled by setting[
This equation determines S implicitly; to leading order in the
drives, S=-2α∑_j|Ω_j|^2/δ_j^2. However,
the Hamiltonian (<ref>) contains the terms (αξ_1,2q^†2qe^-iδ_1,2t+H.c.),
which also contribute to S at this order. Employing perturbation
theory, one finds S=-2α∑_j|Ω_j|^2/δ_j(δ_j+α),
which matches the leading order calculation in Ref. <cit.>. This latter expression is used in the numerics throughout this work.] S=-2α∑_j|ξ̃_j|^2. 
In the frame where the Stark shift terms are eliminated, one finds
modified expressions for the corrections,

    β^(1)    =1-δ_1δ_2δ_Aδ_B/δ̃_1δ̃_2δ̃_Aδ̃_Bδ̃_B+δ̃_1/δ̃_B+δ̃_1+α
    β^(2)    =1-δ_1δ_Aδ_Bδ_C/δ̃_1δ̃_Aδ̃_Bδ̃_Cδ̃_B+δ̃_1/δ̃_B+δ̃_1+α.

Hence, the coupling rates are


    g_v^(1)    =-2αξ̃_1^*ξ̃_2λ̃_Aλ̃_B^*δ̃_̃B̃+δ̃_̃1̃/δ̃_̃B̃+δ̃_̃1̃+α
    
    g_v^(2)    =-2αξ̃_1^*λ̃_Aλ̃_B^*λ̃_̃C̃δ̃_̃B̃+δ̃_̃1̃/δ̃_̃B̃+δ̃_̃1̃+α.

These expressions have the same form as above, but with the replacements
δ→δ̃, i.e. detunings are now defined
relative to the qubit's Stark-shifted frequency. It follows that there also exists a Stark shift correction β^(γ) to the inverse-Purcell enhancement 

    κ_γ
    = κ +γ (g/δ)^2(1+β^(γ)),

where β^(γ)=(δ/δ̃)^2-1, i.e. κ_γ=κ+γ(g/δ̃)^2.
 These corrections are important whenever the drives are strong enough that the qubit's Stark shift becomes comparable to the drive
or mode detunings. Expressions (<ref>), (<ref>),
and (<ref>) are used to produce the plots in
this work.







 §.§ Comparison with numerical Floquet calculation


To assess the accuracy of the expressions (<ref>) and
(<ref>), we compare with numerical calculations of the
coupling rates using the methods developed in Ref. <cit.>. First, we briefly summarize the main results
of that work. The authors consider the process of engineering a bilinear
interaction between two microwave cavity modes that are mutually coupled
to a transmon qubit. Treating the couplings as a perturbation, they
calculate the linear response of the driven transmon. This perturbative
treatment is justified in the dispersive regime. They show that g_v^(1)
can be calculated in terms of a susceptibility matrix χ^(1)(ω_A,ω_B;ω_1,ω_2),
that describes the response of the driven transmon at frequency ω_A
to a weak probe field at ω_B, when subject to drives at
ω_1 and ω_2. The susceptibility can then be computed
numerically to all orders in the drive amplitudes using Floquet theory.
The authors find good quantitative agreement between their theoretical
predictions and experimental results, even for strong drives (ξ>1).

This approach can be directly applied to calculate g_v^(1).
To calculate g_v^(2), we analogously define a higher-order
susceptibility matrix χ^(2)(ω_A,ω_B,ω_C;ω_1)
that captures the response of the transmon at frequency ω_A
to weak probes at ω_B and ω_C, when subject to
a drive at ω_1. Rather than computing χ^(2) directly,
which can be numerically tedious, we note that χ^(2) can be
computed in terms of χ^(1). In the calculation of χ^(1,2),
the drives and probes are treated identically at the Hamiltonian level;
both the drive and probe terms are of the form H=f_j q^† e^-iω_jt+H.c..
For the drives, f_j=Ω_j, while for the probes, f_j=gm_j.
Since the susceptibility is calculated to all orders in the drive fields
but only to leading order in the weak probe fields, going beyond
leading order does not change the result in the limit where the field
f_j is weak. Weak probes and weak drives are thus interchangable,
the only difference being a matter of interpretation. It follows that

    χ^(2)(ω_A,ω_B,ω_C;ω_1)=χ^(1)(ω_A,ω_B;ω_1,ω_C).

This equivalence holds for g_C≪δ_C, which is the same limit
that was already assumed to justify the perturbative treatment. Thus,
the numerical procedure for calculation of g_v^(1) can also
be straightforwardly applied to calculate g_v^(2).

In Figs. <ref>(a) and (b), we calculate g_v^(1,2) numerically as described above, and we compare the results with the analytical expressions (<ref>)
and (<ref>). Good agreement is observed for weak drives
(ξ≲0.4 for the parameters used in the plots). Discrepancies
emerge at stronger drives, but this is expected because the corrections
are obtained perturbatively. In Fig. <ref> (c),
(d), the coupling rates are plotted as a function of δ_A
to make apparent the importance of the AC Stark shift corrections.
Due to the Stark shift, the corrected expressions and numerics are
both red-shifted relative to the uncorrected expressions. Were the
corrections not included, this relative shift would result in a systematic
overestimation of the coupling rates for blue-detuned phonon modes.



The AC Stark shift is responsible for the interesting non-monotonic
behavior of expressions (<ref>) and (<ref>)
with ξ. Intuitively, this behavior is explained by the fact that
the Stark shift causes the qubit to move away from the phonon
modes in frequency space. This reduces the participation of the phonons
in the qubit mode, therby reducing the coupling rate. When optimizing
g_v so as to minimize the  or  infidelity,
the non-monotonicity effectively restricts the drive amplitudes to
the range ξ≤ξ_crit., where ξ_crit.
is the value of ξ for which g_v is maximal.
For example, in Fig. 3 of the main text, the virtual couplings rates are restricted to |g_v^(1)|/2π < 100 and |g_v^(2)|/2π < 25, in accordance with the maximal values of g_v^(1,2) obtained in Fig. <ref>.
Good agreement
between the analytics and numerics is observed for ξ≲ξ_crit.,
validating the use of expressions (<ref>) and (<ref>). 

The numerical calculations in Fig. <ref> (a)
suggest that g_v^(1) could be further increased by allowing
ξ>ξ_crit.. To be conservative, however, we choose
not to exploit this possibility; strong drives can induce deliterious processes
such as the multiphoton resonances described in Ref. <cit.>. 

This comparison illustrates the importance of the corrections derived
above and confirms that the virtual coupling rates are well-described
by expressions (<ref>) and (<ref>) for
the drive strengths considered in this work. 







§ ENGINEERING NONUNIFORM MODE SPACING


As discussed in the main text, nonuniform mode spacing is necessary in order to ensure that the resonance conditions are nondegenerate, i.e. to ensure that a given pair or triple of modes can be selectively coupled. In this section, we first formalize the meaning of nonuniform, then present several schemes for engineering nonuniformity in BAW and SAW systems—those in which modes are otherwise approximately uniformly spaced.
For concreteness, we also provide example schematics for BAW and SAW devices with engineered nonuniformity. Note that phononic crystal resonators are not generally plagued by such degeneracies, since mode frequencies can be controlled by engineering the geometry of each individual phononic resonator. 



As shown in Fig. <ref>, a set of modes is nonuniformly spaced if there exist mode pairs {i,j} and {k,ℓ} for which ν_ij≠ν_kℓ, where ν_ij=|ω_i-ω_j| is the frequency spacing between modes i and j. 
In the context of multimode coupling, it is useful to quantify this nonuniformity as follows. 
Let 𝒮 denote the set of all modes that are used to store quantum information, and let 𝒫 denote the set of all mode pairs that one chooses to couple (we provide examples below). 
The connectivity of the system is then described by a graph with vertices 𝒮 and edges 𝒫. 
As a practically relevant measure of the nonuniformity, we define the quantity 

    Δν = min_{i,j}∈𝒫[min_{k∈𝒮,ℓ}≠{i,j}|ν_ij-ν_kℓ|],

which lowerbounds the frequency selectivity of two-mode couplings. 
Explicitly, the beamsplitter resonance condition, ω_2-ω_1 = ω_B-ω_A for a pair of modes {A,B}∈𝒫 is detuned from all other beamsplitter resonance conditions involving any mode in 𝒮 by at least Δν. 
Highly selective virtual couplings thus require g_v/Δν≪ 1. Note that since Δν depends on the choices of 𝒮 and 𝒫, there can exist a tradeoff between selectivity and the effective size and connectivity of the system.
The definition of Δν can be straightforwardly generalized to the case of three-mode couplings.

Whether a given pair or triplet of modes can be selectively coupled depends on the structure of the nonuniformity, and 
in this regard it is convenient to classify different sorts of nonuniformity according to properties of ν_j,j+1. We study two such classes in the examples below: point defect nonuniformities, for which ν_j,j+1 is constant except in the vicinity of a single defect, and periodic nonuniformities, for which ν_j,j+1 is periodic. Of course, other classes exist, but we focus on these two classes since instances can readily be engineered in cQAD systems.



 §.§ External mode hybridization

A point defect nonuniformity can be created by coupling the phonons to some external mode, such as a microwave resonator. 
As demonstrated in Ref. <cit.>, and sketched in Fig. <ref>(a,b), the resulting mode hybridization can significantly shift phonon mode frequencies within some bandwidth 𝒟 of the external mode. 
The nonuniformity Δν is dictated by the magnitude of these frequency shifts. For example, frequency shifts of order 1 were demonstrated in Ref. <cit.>.

This class of nonuniformity can enable selective coupling: selective two-mode coupling is possible if one or both involved modes lie in 𝒟, and selective three-mode coupling is possible if two of the three involved modes lie in 𝒟. Hence, the set 𝒮 can include arbitrarily many modes, but the set 𝒫 can only include mode pairs with at least one mode in 𝒟. While modes outside of 𝒟 cannot be directly coupled to one another, information from these modes can instead be swapped into modes in 𝒟, manipulated, and swapped back. 
Note that the coherence of the external mode should be comparable to that of the phonons, lest the hybridization result in a significant increase in effective decay rates, and in general there may exist a tradeoff between increased nonuniformity and enhanced decay.






 §.§ Two phonon mode families


Another approach is to create a periodic nonuniformity by simultaneously coupling the transmon to two families of phonon modes <cit.> with different free spectral ranges (FSRs). 
While modes within each family are uniformly spaced, the FSR difference causes the spacing between modes from different families to vary, as shown in Fig. <ref>(c,d). 
This nonuniformity enables two modes from different families to be selectively coupled, but because of the periodicity, selectivity is only guaranteed over a finite bandwidth smaller than one period. 
With two mode families, a set 𝒮 containing ≈ν/Δν modes can be found wherein any two modes from different families can be selectively coupled with Δν =  |ν_1-ν_2|, where ν_1,2 are the FSRs of the two families.

By itself, the use of two mode families does not enable selective three-mode coupling[At least two modes from any three come from the same family, and since the modes in each family are uniformly spaced, there necessarily exists another set with the same resonance condition.], but this limitation can be circumvented by coupling the transmon to one or more external modes. 
For example, the BAW devices of Refs <cit.> are housed in microwave cavities, and coupling the transmon to a high-Q cavity mode can enable selective three-mode coupling between the cavity and any pair of modes in 𝒮. 
In a SAW device, the additional mode could come from another SAW resonator or a microwave resonator. The transmon itself could even serve as the external mode, but gate fidelities would then be directly limited by transmon coherence.
Ideally, the coherence of the external mode should be comparable to that of the phonons, lest it limit gate fidelity.



 §.§ Composite resonators

Yet another approach is to employ a composite acoustic resonator, in which phonons propagate in media with different indices of refraction [Fig. <ref>(e,f)]. Reflections at the interfaces can give rise to a periodic modulation of the FSR <cit.>. 
As in the case of two mode families, this periodic nonuniformity can enable selective two-mode coupling within a finite bandwidth 𝒮, though the magnitudes of both 𝒮 and Δν depend on the nature of the modulation. 

Whether selective three-mode coupling within 𝒮 is feasible depends on the of the specific nature of the FSR modulation. In cases where it is not already possible, selective three-mode coupling can be enabled by either coupling the transmon to some external mode, as previously described, or alternatively by restricting the bandwidth over which the transmon-phonon coupling is appreciable. 
For example, if the transmon-phonon coupling is only appreciable within 𝒮, as in Fig. <ref>(e), then selective three-mode coupling is possible since the system contains an effectively finite number of nonuniformly spaced modes. 
In SAW systems, the coupling bandwidth can be tuned by changing the number of fingers in the interdigitated transducer[Because SAW resonators have finite bandwidth, care should be taken to avoid coupling to unconfined modes. This problem can be solved in general by engineering the transmon-phonon coupling bandwidth to lie within the SAW resonator bandwidth. The size of both bands can be tuned by varying the number of fingers in the respective interdigitated transducers <cit.>. ] <cit.>. 
In BAW systems, the coupling bandwidth can be similarly tuned by changing the electromechanical transducer's geometry. 
For instance, in a transducer comprised of alternating layers of piezoelectric and non-piezoelectric materials, the spacing, thickness, and number of such layers could be chosen so that the coupling has a narrow response centered at a particular frequency, as in a Bragg reflector.



 §.§ Example schematics



For concreteness, in Fig. <ref> we provide example schematics for SAW and BAW devices in which nonuniformity is engineered according to the strategies described above. Fig. <ref>(a) shows a SAW device that exploits the external mode hybridization strategy. A SAW resonator is fabricated on a piezoelectric substrate, and coupling between the transmon and the phononic modes is enabled by an interdigitated capacitor. A superconducting coplanar waveguide resonator is also coupled to the phononic modes, and the hybridization of the phononic modes with the resonator mode creates the necessary nonuniformity.

Fig. <ref>(c) shows a BAW device that exploits the two mode families strategy. The device is based on those demonstrated in Refs. <cit.>; a three-dimensional (3D) transmon is housed inside a microwave cavity, and thin disks of piezoelectric material (transducers) fabricated in the transmon's pads enable the transmon to couple to BAW modes in the substrate. Two modifications have been made relative to the devices in Refs. <cit.>. First, an additional transducer has been added so that the transmon simultaneously couples to two families of modes. Second, the thickness of the substrate beneath one of the transducers has been reduced so that the two families have different FSRs. The microwave cavity mode, which dispersively couples to the transmon, provides the external mode necessary to enable selective three-mode couplings. We note that other elements, e.g. a separate readout resonator for the transmon, can be integrated into 3D architectures in such a way that the transmon can be driven and measured without involving the cavity mode <cit.>. 







§ DETAILED DESCRIPTION OF THE CQAD QRAM

In this section we provide a more complete description of the QRAM proposed in the main text. The operation of the quantum quantum routers is first discussed in detail, and then schemes for extracting data from either a classical or a quantum database are presented.



 §.§ Operation of a quantum router


During a QRAM query, the operation of each quantum router can be divided ino four stages: initialization, downstream routing, upstream routing, and extraction. In the initialization stage [Fig. <ref>(a)], an incoming address qubit in the top mode is stored in the routing mode at the vertex so that it can control the routing of subsequent qubits. Initialization is performed simply by swapping the states of the top and routing modes. In the downstream routing stage [Fig. <ref>(b)], the router directs incoming qubits to one of two output modes conditioned on the state of the address qubit. As shown in the corresponding circuit, this routing operation can be implemented using a  gate and a controlled- gate. Note that the action of the controlled- gate is trivial when the routing qubit is in |0⟩. Similarly,  the action of the  gate is trivial when the routing qubit is in |1⟩, since the affected modes are both in |0⟩ after the controlled- has been performed.

Once all downstream routing is complete and the data has been transferred to the bus, the inverse operations (upstream routing and extraction) must be performed in order to disentangle the address and bus qubits from the routers.  Both the  and controlled- unitaries are their own inverses, so the upstream routing [Fig. <ref>(c)] and extraction [Fig. <ref>(d)] operations can be implemented by time-reversing the downstream routing and initialization operations, respectively. 





 §.§ Database access schemes


In a QRAM, the address and bus registers must necessarily be composed of qubits. However, the database itself can contain either classical or quantum data, depending on the intended application. How the database is accessed depends on the type of data. In the classical case, data can be directly copied into the bus, but in the quantum case the no-cloning theorem prevents one from copying quantum data. Instead, the bus can be entangled or swapped with the data, which generally leaves the address-bus system entangled with the database after a query. Additionally, while a QRAM should in general be able to both read from and write to the database, in many applications a quantum read-only memory (qROM) is sufficient <cit.>. This reduced functionality can be used to simplify the database access scheme. Below we present schemes for reading and writing both quantum and classical data, as well as a simplified scheme for read-only classical data.





Figure <ref> illustrates procedures for accessing classical and quantum databases. 
To read quantum data from the database [Fig. <ref>(a)], we introduce a so-called pointer qubit, initially prepared in the state |1⟩.
The purpose of this qubit is to indicate which database entries should be extracted. After all address qubits have been routed into place, the pointer qubit follows their path to the bottom of the tree. 
Modes at the bottom of the tree then serve as controls in controlled- operations that extract the desired data qubit. 
In each branch of the superposition, only the database entry adjacent to the pointer qubit is swapped out. 
The data, pointer, and address qubits are then routed out of the tree in sequence to complete the query. To write quantum data to the database[Writing either classical or quantum data to the database generally leaves the database in a superposition, unless entries are written one by one.], the entire procedure can be run in reverse. 
Note that, if a query is error-free, the pointer qubit remains in the state |1⟩ throughout and is disentangled from the rest of the system at the end of the query. 
Measuring this qubit after the query could thus provide a useful means of detecting errors.  

In the classical case, a data bit D_j can be encoded in a mode as the Fock state |D_j⟩. The operations described above for the quantum case then suffice to read or write classical data. 
However, since the read operation involves swapping the bus and data qubits, data is removed from the database during each query. 
In the classical case, data can instead be copied to the bus [Fig. <ref>(b)], which ensures that queries leave the database undisturbed. 
To copy data, the bus qubit is initially prepared in the state (|0⟩+|1⟩)/√(2)≡|+⟩, and routed to the bottom of the tree.
The desired classical bit is then copied to the bus by performing a  operation between the data and bus qubits; if |D_j⟩=|0⟩, the bus is unaffected, but if |D_j⟩=|1⟩, the  operation flips |+⟩ to |-⟩, hence encoding the bit in the |±⟩ basis. Importantly, this operation does not disturb other modes, since these phase shifts leave the state |0⟩ unaffected. 
The bus and address qubits are then routed out of the tree to complete the query.

In the classical read-only case [Fig. <ref>(c)], it is not necessary to encode data in the Fock states of phonon modes.
Instead, the entire database can be stored in a purely classical memory array so that it is not susceptible to decoherence. 
To query such a database, the bus is again prepared in |+⟩ and routed to the bottom of the tree. The desired classical bit is copied to the bus by applying phase shifts to each mode at the bottom of the tree in accordance with the database entries.
If D_j=1, a phase shift of -1 is applied to the j^th mode, and if D_j=0, no phase shift is applied. 
As discussed in the main text, phase shifts can be imparted to a mode by tuning the relative driving phases during a consecutive sequence of two  operations with an ancillary mode. 







§ DERIVATION OF GATE FIDELITY EXPRESSIONS



Given a quantum process ℰ constituting a noisy implementation of a gate U, the gate fidelity is typically defined as 

    ℱ = ∫ dψ⟨ψ|U^†ℰ(ψ) U |ψ|,⟩

where the integral is over the uniform measure dψ on state space, with normalization ∫ dψ = 1.  For an n-qubit gate, the states |ψ⟩ are elements of the Hilbert space of the n qubits involved in the gate. With this definition, crosstalk with modes not involved in the gate is quantified implicitly via its deleterious effects on qubits involved in the gate, and the fidelity is independent of whether other modes are used to store quantum information.

Other definitions are possible in the context of multimode systems. For example, given a multimode system in which M modes are used to store qubits, one can consider an alternate definition of the fidelity,

    ℱ̃ = ∫ dψ⟨ψ|Ũ^†ℰ̃(ψ) Ũ |ψ|,⟩

where Ũ = U ⊗1^⊗(M-n) is the tensor product of the ideal n-qubit gate and the identity operation on the remaining M-n qubits, ℰ̃ is the corresponding noisy implementation of the gate, and the states |ψ⟩ are elements of the Hilbert space of all M qubits in the system. With this definition, crosstalk is accounted for explicitly via deviations of ℰ̃ from the ideal operation Ũ, and the fidelity generally depends on M.  


In this section, we compute the direct and virtual gate fidelities according to both of the metrics, ℱ and ℱ̃, which we refer to as the local and global fidelities, respectively.
The local fidelity ℱ is commonly used in the literature in the context of multi qubit architectures <cit.> and has the advantage that it is largely independent of implementation differences between BAW, SAW, and phononic crystal devices. 
In contrast, the global fidelity ℱ̃ is a more comprehensive metric, but it depends sensitively on implementation details, making it challenging to analyze in general. 
In the main text we thus compute the local fidelity so as to facilitate a straightforward and implementation-independent comparison with other architectures.
For concreteness, though, we compute the global infidelity for a particular implementation based on a BAW device below, and we show that the main claim—that virtual gates can offer higher fidelities for long-lived phonon modes—holds when using either metric. A closely related discussion can be found in Ref. <cit.>.




 §.§ Local infidelity

We now derive the expressions for the local infidelity 1-ℱ given in the main text (Eqs. 3 and 4).
There are two dominant contributions to the infidelity: decoherence and spectral crowding (crosstalk). The contributions to the direct and virtual gate infidelities from decoherence are, respectively,

    (κ+γ)t_d   and  κ̅_γ t_v,

where t_d = (c_d π/2 g_d) and t_v = (c_v π/2 g_v) are the total gate times, and κ+γ and κ̅_γ are the total decoherence rates for the direct and virtual gates respectively. For two-mode couplings, the total decoherence rate is κ̅_γ = κ_γ^A + κ_γ^B, while for three-mode couplings the rate is κ̅_γ =( κ_γ^A + κ_γ^B + κ_γ^C)/2. The factor of 1/2 in the latter expression results from averaging over the gate duration; when there is a phonon each in modes A and B, the rate is κ_γ^A + κ_γ^B, but once these two phonons have been converted into a single phonon in mode C the rate is κ_γ^C.  
Note that we have neglected the inverse Purcell enhancement in the direct case under the assumption that this enhancement is negligible relative to γ. 


The presence of other modes also contributes to the infidelity, regardless of whether these other modes are used to store quantum information. When performing a gate, transitions between the modes involved in the gate and the other modes are driven off-resonantly, and we approximate the spectral crowding infidelity as the probability that one of these unwanted transitions occurs. This probability is computed in Ref. <cit.> for direct gates: assuming a set of uniformly spaced modes with free spectral range ν, the infidelity is approximately[Following Ref. <cit.>, we neglect a constant prefactor of order 1 on the right hand side, with the justification that (g_d/ν)^2 is actually a pessimistic upper bound; the spectral crowding infidelity can be reduced by smoothly ramping up the drives.]

    ∑_n(g_d/δ_n)^2≈(g_d/ν)^2,

where the sum on the left runs over all unwanted transmon-mode transitions, each detuned by successive multiples of ν, i.e. δ_n = {±ν,± 2ν ,…}. Importantly, this infidelity is independent of the total number of modes in the system. Similarly, for virtual gates the spectral crowding infidelity is approximately

    ∑_n(g_v/δ_n')^2≈(g_v/Δν)^2,

where the sum runs over all unwanted virtual couplings that affect modes involved in the gate. More precisely, n indexes all unwanted two-mode (three-mode) couplings which involve at least one mode from the set {A,B} ({A,B,C}), and δ_n' is the resonance condition detuning between the n^th unwanted coupling and the desired coupling. For example, in the two-mode case, these detunings are of the form δ'_n = (ω_B-ω_A) - (ω_i - ω_j), with either i or j ∈{A,B}. In performing the sum we have assumed that unwanted couplings are detuned by successive multiples of Δν, i.e. δ_n' = {±Δν, ± 2Δν ,…}. In general the δ_n' depend on the specific structure of the nonuniformity, but the scaling (g_v/Δν)^2 holds regardless.

Summing the contributions from decoherence and spectral crowding gives the infidelity approximations used in the main text,

    1-ℱ_d   ≈ (γ+κ_γ)[c_dπ/2g_d]+(g_d/ν)^2,
    
    1-ℱ_v   ≈κ̅_γ[c_vπ/2g_v]+(g_v/Δν)^2.

Note that these first-order expressions are valid only in the limit 1-ℱ≪ 1. When plotting the infidelity, we include higher-order terms to ensure the output is restricted to [0,1]; the above expressions are of the form k t + ε, but we plot corresponding higher-order expressions 1- e^- k t(1-ε).



 §.§ Global infidelity

We now compute the global infidelities 1-ℱ̃ for direct and virtual gates. As before, decoherence and spectral crowding are the two dominant contributions to the infidelity. For a multimode system in which M≥2 modes are used to store qubits, the total decoherence contribution to the direct gate infidelity 1-ℱ̃_d is

    [ (κ+γ) +(M-2)κ]t_d,

since the two qubits involved in the gate experience a total decoherence rate of approximately κ + γ, while the remaining M-2 qubits are left idling in their respective resonator modes during the gate. 
The total decoherence contribution to the virtual gate infidelity 1-ℱ̃_v is κ̅_γ t_v, where κ̅_γ is the total decoherence rate, which now includes contributions from all M modes,

    κ̅_γ = 
         
           Mκ_γ   for two-mode coupling,
    (M-1/2)κ_γ   for three-mode coupling.

The three-mode coupling expression contains (M-1/2) as opposed to M because the conversion of two phonons into a single phonon and back reduces the average decoherence rate.  


Next, consider the contribution from spectral crowding. 
Both direct and virtual gates require frequency selectivity, and spectral crowding infidelity results from the limited ability to resolve a desired transition from others which are nearby in frequency space. In order to compute the global spectral crowding infidelity, it is thus necessary to specify the full distribution of transition frequencies. Because these distributions can vary from platform to platform, however, it is challenging to make general statements about the spectral crowding infidelity. Computing this infidelity for the multitude of possible implementations is beyond the scope of this work, so in the following we consider a particular implementation based on a BAW or SAW device as an illustrative example. We stress that the following discussion applies only to this particular implementation. 

Consider a multimode cQAD device in which the acoustic modes are supported in a BAW or SAW resonator. In the context of direct gates, we consider a single multimode resonator, in which modes are approximately uniformly spaced with FSR ν. In the context of virtual gates, nonuniform mode spacing is required, and we consider the case where this nonuniformity is engineered through the use of two mode families with different FSRs, |ν_1-ν_2|=Δν≠ 0, as discussed in Sec. <ref>.  

The spectral crowding contributions to the direct and virtual infidelities, 1-ℱ̃_d and 1-ℱ̃_v, can be respectively written as

    S_d(g_d/ν)^2   and   S_v(g_v/Δν)^2,

where the coefficients S_d,v remain to be determined. The spectral crowding contribution to 1-ℱ̃_d is given by

    ∑_n(g_d/δ_n)^2 = π^2/3(g_d/ν)^2,

where δ_n = {±ν,± 2ν ,…} as before. Thus, S_d = π^2/3.
Notice that unwanted transmon-phonon transitions contribute to the direct gate infidelity regardless of whether qubits are encoded in the other modes, cf. the spectral crowding term in 1-ℱ_d.

We compute the spectral crowding contribution to 1-ℱ̃_v numerically by generating a large array of modes and directly computing the sum ∑_k'(g_v/δ_n')^2. Here, n indexes all unwanted virtual couplings that affect at least one of the M qubits (cf. the local infidelity, where only unwanted couplings that affected at least one of the qubits involved in the gate were considered).
We take the average over all possible desired transitions among the M modes to obtain S_v.  As shown in Fig. <ref>(a), we find that that S_v increases with M for fixed ν/Δν. Note that for the particular implementation discussed in Sec. <ref>, S_v is the same for both engineered two- and three-mode couplings.

Summing the contributions from decoherence and spectral crowding gives

    1-ℱ̃_d   ≈[ (κ+γ) +(M-2)κ]t_d
        + S_d(g_d/ν)^2
        ,
    
    1-ℱ̃_v   ≈ 
         Mκ̅_γ t_v
        +S_v(g_v/Δν)^2
    .

As before, g_d,v may be tuned to minimize the infidelities, and we plot the minimal infidelities for the cases of M=2 and M=10 in Fig <ref>(b,c). 
For the M=2 case, the plots look very similar to the local infidelity plots in Fig. 3 of the main text; the only difference results from the fact that we keep track of the prefactors S_d and S_v in the global infidelity. For the M=10 case, reductions in both direct and virtual gate fidelities are observed relative to the M=2 case, as expected. The reductions are more pronounced in the virtual case, due in part to the fact that S_v scales with M while S_d does not.   
This scaling difference is a property of the specific BAW/SAW implementation that we consider here, and in fact we have chosen this example to illustrate that such differences can exist. Because of this scaling difference, this example effectively constitutes a worst-case scenario for the comparison of direct and virtual gate fidelities.  
Despite this, there are still large regions of parameter space where virtual gates offer higher fidelities.




Two additional remarks are necessary. First, we stress that the spectral crowding infidelities can vary from platform to platform, so the above results are specific to the BAW/SAW implementation that we have considered. In phononic crystal devices, for example, S_d,v will have different a dependence on M. In these devices the mode frequencies must be distributed within a fixed bandgap, so the effective FSR scales with 1/M, and hence S_d scales with M^2 <cit.>. It is possible to arrange the mode frequencies so that S_v also scales with M^2. As a result, the comparison between direct and virtual gates would be less affected by the choice of M in such devices. 

Second, the above analysis illustrates that for large M, it will be impractical to control all M qubits using only a single transmon. A more scalable approach is to construct larger devices out of smaller modules, where each module contains its own acoustic modes and transmon. In this context, the optimal number of modes per module will typically to be closer to 10, rather than to 100 or 1000.  
For example, the performance of a module can be quantified using the quantum volume <cit.>,

    V = max_M [min(M, d(M)) ]^2,

where d(M) = 1/M(1-ℱ̃) is the maximal circuit depth[This assumes no parallelization. Unlike direct gates, virtual gates can be performed in parallel, and parallel operation can increase the quantum volume.]. For the BAW/SAW device discussed above, we find that the M which optimizes the quantum volume satisfies M≤12 across the entire parameter space shown in the plots of Fig. <ref>(b,c).
Similarly, in Ref. <cit.>, M≈15 was found to optimize the quantum volume of multimode cQAD devices based on phononic crystal resonators, assuming realistic future device parameters. 










§ EFFECTS OF PHONON MODE FREQUENCY SHIFTS


Due to hybridization with the transmon, the acoustic modes in multimode cQAD systems can experience frequency shifts caused by transmon drives (AC Stark shifts), or the presence of phonons (Kerr interactions). In this section, we discuss the effects of such shifts on the two- and three-mode resonance conditions given in the main text. 



 §.§ AC Stark shifts

The AC Stark shift terms (Table I, Row 2),

    -2α∑_k|ξ_k|^2|λ_j|^2m_j^†m_j≡  S_jm_j^†m_j,

describe a drive-dependent shift of the phonon mode frequencies. In order to engineer resonant interactions in the presence of these shifts, the drive strengths and frequencies must be carefully chosen so that interactions between the Stark-shifted modes are brought on resonance <cit.>. For example, a resonant beamsplitter-type coupling between modes A and B can be engineered by applying two drive tones tones whose frequencies ω_1,2 satisfy,

    ω_2-ω_1=(ω_B+S_B)-(ω_A+S_A).

Since the shifts S_j are calculable, it is straightforward to ensure that this modified resonance condition is satisfied. For example, one can fix the drive amplitudes (which, when the drives are far-detuned, approximately fixes the S_j), then tune the frequencies to satisfy the modified condition.

As discussed in the main text, nonuniform mode spacing is necessary to ensure that couplings between modes can be selectively engineered. It is crucial that this nonuniformity persists in the presence of the phonon modes' Stark shifts, lest the resonance conditions become degenerate. In general, resonance conditions will remain nondegenerate provided that the Stark shifts are small relative to the nonuniformity, S≪Δν. Since the Stark shifts are comparable in magnitude to the engineered two-mode couplings, S∼ g_v^(1), the same degree of nonuniformity needed to guarantee highly selective couplings, Δν≫ g_v^(1), also suffices to ensure that the Stark shifts do not create additional degeneracies in the resonance conditions. 



 §.§ Self- and Cross-Kerr interactions

We next consider the effects of and self- and cross- Kerr interactions (Table I, row 5),

    -α/2|λ_j|^4m_j^†m_j^†m_jm_j   ≡ -χ_jj/2m_j^†m_j^† m_jm_j, 
    and     
        -2α|λ_j|^2|λ_k|^2m_j^†m_k^†m_jm_k   ≡ -χ_jkm_j^†m_jm_k^†m_k,    (j≠ k).

Below, we discuss the effects of these interactions on engineered two-mode and three-mode couplings.
Unlike the AC Stark shifts, the effects of these terms on the resonance conditions cannot be compensated in general. 
However, we show that the infidelity associated with these terms can be made negligible relative to the infidelity associated with decoherence and spectral crowding. 



  §.§.§ Effects of Kerr on two-mode couplings

For an engineered two-mode interaction H^(1)= (g_v^(1) m_A^† m_B + H.c.), we first consider the effects of cross- and self-Kerr terms involving only modes A and B. Recall that we consider encodings with ≤1 phonon per mode, so only the initial state |1,1⟩_A,B is affected by these interactions.  From this initial state, the coupling H^(1) induces evolution within the subspace {|11⟩,|20⟩,|02⟩}. In this subspace, the full Hamiltonian H = (H^(1) + Kerr terms) is

    H =
        [         -χ_AB  √(2) g_v^(1)  √(2) g_v^(1); √(2) g_v^*(1)         -χ_AA             0; √(2) g_v^*(1)             0         -χ_BB ].

We define the infidelity in terms of the overlap between the final state obtained by evolution under this Hamiltonian, exp(-i H t)|11⟩, and the ideal final state in the absence of these Kerr terms, exp(-i H^(1) t)|11⟩,

    1-ℱ = 1-|⟨11|e^i H^(1) te^-i H t|11||⟩^2.

The infidelity can be approximated by making two assumptions: (i) that χ≪ g_v^(1), and (ii) that χ_AA≈χ_BB≈χ_AB/2. The first assumption is justified by the fact that the transmon-induced self and cross-Kerr interactions, which are fourth order in λ = g_j/δ_j ≪1, are typically much weaker than the engineered coupling g_v^(1), which is only second order in λ. For example, in Ref. <cit.>, χ_AB/2π≈ 1 kHz, while g_v^(1) / 2π≈ 30 kHz. The second assumption is justified in the regime where λ_A ≈λ_B, i.e. when the transmon-phonon hybridization is comparable for modes A and B. This will be the case, for example, when the bare couplings are comparable, g_A ≈ g_B, and the relative difference in the detunings is small, |δ_B-δ_A|/δ_A,B≪ 1.
Under these assumptions, one finds that the infidelity scales quadratically with the small parameter (χ_AB/g_v^(1)),

    1-ℱ≈1/64sin^4(2 g_v^(1) t)(χ_AB/g_v^(1))^2 + O[(χ_AB/g_v^(1))^4].

We provide numerical estimates of this infidelity below. 

In multimode systems, cross-Kerr interactions with other modes j∉{A,B} can also contribute to the infidelity. The effect of such interactions is to modify the two-mode resonance condition,

    ω_2-ω_1=ω_B-ω_A+ D^(1),

where D^(1)=∑_j(χ_Aj-χ_Bj)n_j, and n_j∈{0,1} denotes the number of phonons in mode j. Since the n_j are not known a priori, such shifts cannot generally be compensated and contribute to the gate infidelity. When the detuning of the drives from the proper resonance condition is small, D^(1)≪ g_v^(1), the associated contribution to the infidelity is approximately

    1-ℱ≈ (D^(1)/g_v^(1))^2.


To numerically estimate the infidelities typically associated with these cross- and self-Kerr interactions, we consider the experiment performed in Ref. <cit.>, in which a beamsplitter-type coupling was engineered between two microwave cavity modes via their mutual coupling to a transmon[Though the experiment used electromagnetic modes, rather than acoustic modes, the engineered coupling and transmon-induced Kerr interactions are expected to have similar magnitudes since they arise via the same mechanisms.]. In that work, engineered couplings as large as g_v^(1)/2π = 30kHz were demonstrated, and Kerr interactions were measured to be χ/2π≲ 1kHz.  With these parameters, the infidelity contribution (<ref>) is < 10^-4. 


To estimate the contribution (<ref>), we note that the cross-Kerr differences (χ_Aj-χ_Bj) will typically be much smaller than the cross-Kerr interactions themselves. 
For example, when the frequency separation between the modes satisfies ν≪δ_A,B, 
and assuming g_A≈ g_B, 

    (χ_Aj-χ_Bj)≈ 2 ν/δ_Aχ_Aj≪χ_Aj.

For the parameters in Fig. 3 of the main text, typical cross-Kerr differences are ∼100 Hz .
The contribution (<ref>) also depends on the total number of phonons, which will usually be order 10 or fewer (see Sec. <ref>). Thus, typically D^(1)/2π≈ 100 Hz - 1 kHz, and the infidelity contribution (<ref>) is 10^-5 - 10^-3. 
We conclude that, for two-mode couplings, the infidelity associated with self- and cross-Kerr interactions is typically negligible relative to the infidelity associated with decoherence and spectral crowding. 




  §.§.§ Effects of Kerr on three-mode couplings

Next we study the effects of Kerr interactions on engineered three-mode couplings. Given an engineered three-mode interaction H^(2)= (g_v^(2) m_A m_B m_C^† + H.c.), we first consider the effects of cross- and self-Kerr terms involving only modes A, B, and C.  As discussed in the main text, the initial states we consider have the form |n_A≤1,n_B≤1,n_C=0⟩, and of these states only |1,1,0⟩ is affected nontrivially by such Kerr interactions. This is also the only initial state on which the engineered coupling H^(2) acts nontrivially. Thus, the effects of cross-Kerr and self-Kerr interactions between the modes A, B, and C can be fully compensated simply by replacing (ω_A +ω_B)→ (ω_A +ω_B-χ_AB) in the resonance condition. 

Cross-Kerr interactions with other modes j∉{A,B,C} yield the modified resonance condition,

    ω_1=ω_A+ω_B-ω_C+ D^(2)

where D^(2)=∑_j(χ_Cj-χ_Aj-χ_Bj)n_j. Since n_j is not known a priori these shifts cannot generally be compensated, and the associated contribution to the infidelity would be approximately (g_v^(2)/D^(2))^2. This contribution could be significant because D^(2) depends on the magnitudes of the cross-Kerr couplings themselves, rather than their differences, cf. D^(1). 

However, since all operations described in this work preserve the total phonon number[While couplings of the form (m_Am_Bm_C^†+H.c.) do not preserve phonon number, the gates we describe only require cyclic evolutions under such couplings, which do preserve total phonon number.], the detuning D^(2) can potentially be reduced by subtracting off the average cross-Kerr shift. Let χ̅ denote the average cross-Kerr coupling across all populated modes, and let n_tot denote the total number of phonons in the system. If the drive frequency is chosen to be ω_A + ω_B - ω_C - χ̅(n_tot-2), the detuning becomes[Because the coupling acts nontrivially only on initial state |11⟩_AB, the total number of phonons in modes j∉{A,B,C} can be assumed to be n_tot-2.]

    D'^(2)= D^(2) +  χ̅(n_tot-2) = ∑_j[
        (χ_Cj-χ̅)
        -(χ_Aj-χ̅)
        -(χ_Bj-χ̅)
        ]n_j.

D'^(2) depends only on deviations in cross-Kerr couplings from the average; when these deviations are small, the infidelity is suppressed. We estimate that the detuning D'^(2)/2π lies in the range ≈ 100Hz - 1kHz, so the associated infidelity is (g_v^(2)/D'^(2))^2 ≈ 10^-4 - 10^-2. This infidelity can be even further suppressed using  mode-dependent compensations (e.g. by leveraging the knowledge that, say, λ_A>λ_B>λ_C). We thus conclude that, with appropriate compensations to the drive frequency, the infidelity associated with self- and cross-Kerr interactions can be made negligible relative to the infidelity associated with decoherence and spectral crowding.


In summary, when engineering two-mode couplings, the drive frequencies should be chosen to satisfy the resonance condition

    ω_2 - ω_1 = ω_B - ω_A + S_B - S_A.

This choice fully compensates the effects of AC Stark Shifts. While the effects of self- and cross-Kerr interactions cannot generally be compensated, the resulting detuning is only a function of small cross-Kerr differences.
Similarly, when engineering three-mode couplings, the drive frequency should be chosen to satisfy the resonance condition

    ω_1 = ω_A + ω_B - ω_C + S_A + S_B - S_C - χ_AB - χ̅ (n_tot-2).

This choice fully compensates the effects of AC Stark shifts, as well as cross- and self-Kerr terms between the modes A, B, and C. Cross-Kerr interactions with other modes generally contribute to the infidelity. However, the associated detuning from the resonance condition (<ref>) is only a function of small deviations in cross-Kerr, and the infidelity can be further suppressed by using appropriate mode-dependent compensations.








§ UNIVERSALITY 

In the main text, we focus on the beamsplitter, , and  operations, as these are the only operations we require to implement a QRAM. For completeness, we now show that engineered two- and three-mode interactions can be used to implement a universal gate set. We refer the reader to Refs. <cit.> for details.

Consider the encoding of a qubit into two phononic modes, A and B, as

    |0̅⟩ = |1⟩_A|0⟩_B,  |1̅⟩ = |0⟩_A|1⟩_B,

where |0̅, 1̅⟩ denote the qubit states, and |0,1⟩_A,B denote phonon Fock states. Note that this encoding differs from the single-mode encoding discussed in the main text; this dual-rail encoding—commonly considered in optical quantum computing <cit.>—is the simplest that allows for universal computation in our scheme. With this encoding, engineered two-mode interactions of the form

    H^(1) = g_v^(1) m_A^† m_B + H.c.

can be used to implement arbitrary single-qubit rotations: for g_v^(1) real (imaginary), H^(1) generates rotations around the X (Y) axis of the Bloch sphere. The phase of g_v^(1) is controllable via the phases of the drives used to engineer the interaction. 

Given two qubits, respectively encoded in modes {A_1,B_1} and {A_2,B_2} with the encoding (<ref>), and an ancillary mode C initialized in |0⟩_C, we can implement a  gate using the engineered three-mode interaction,

    H^(2)=g_v^(2)m_B_1m_B_2m_C^†+H.c.

Evolution under H^(2) for a time π/g_v^(2) imparts a -1 phase to the |1̅1̅⟩ state, while leaving all other initial states unaffected, thereby realizing the CZ gate. Arbitrary single-qubit rotations and CZ are universal. Note that these engineered two- and three-mode interactions can similarly be used to implement a universal set of logical gates for certain bosonic quantum error-correcting codes <cit.>.







§ PARALLELISM AND LIMITATIONS 

One appealing property of virtual gates is that they can be performed in parallel. Parallelism, together with the potential for a high degree of connectivity, makes the proposed architecture highly efficient. In this section, we discuss how virtual gates can be performed in parallel and consider practical limitations on the degree of parallelization possible. 

We first provide a simple example to illustrate how gates may be performed in parallel. Say one wants to implement a  gate between modes A and B in parallel with a  gate between modes C and D. A collection of drives must be applied so that the requisite couplings, (m_A^†m_B+H.c.) and (m_C^†m_D+H.c.) are brought on resonance simultaneously. Assuming |ω_B-ω_A|≠|ω_D-ω_C|, this requires a minimum of three distinct drive tones. For example, one can apply three drives whose frequencies satisfy the resonance conditions

    ω_2-ω_1   =ω_B-ω_A, 
    ω_3-ω_2   =ω_D-ω_C,

where for simplicity we have the neglected compensations for Stark shifts and cross-Kerr interactions (see Sec. <ref>).
Importantly, in order to ensure that only the intended interactions are brought on resonance, we must also have |ω_3-ω_1|≉|ω_j-ω_k| for any modes j,k which are used to store quantum information. Provided that this condition is satisfied, the  gates are implemented in parallel without disturbing other modes, and the operation in parallel is not fundamentally different from serial operation. 

The restriction |ω_3-ω_1|≉|ω_j-ω_k| illustrates an important general point: when applying multiple drives to perform gates in parallel, all combinations of the drive frequencies must be carefully checked in order to verify that only the desired couplings are brought on resonance. More specifically, for d drives applied simultaneously, d(d-1)/2 pairwise combinations must be checked to ensure that only the desired two-mode couplings are resonant. Furthermore, each of the d frequencies must also be individually checked to ensure that only the desired three-mode couplings are resonant. As the number of drives grows, it is likely to become increasingly challenging to avoid accidentally bringing some unwanted coupling on resonance. The need to avoid unwanted resonant couplings can thus limit on the degree of parallelization possible. This limit generally depends on the specific structure of the nonuniformity, the desired connectivity, and the total number of modes. 

An additional practical limit may come from the need to avoid drive-induced heating of the transmon. As discussed in Ref. <cit.>, 
the applied drives can facilitate an effective “heating” of the transmon via its natural dissipation and dephasing processes, and this heating can limit the operation fidelity.  With two drive tones, the rate of this drive-induced heating is <cit.>[There is a factor of 2 missing in Eq. (48) in Ref. <cit.>; expression (<ref>) is the correct expression. ]

    γ[ 
        |αξ_1^2/2δ_1+α|^2 +
        |αξ_2^2/2δ_2+α|^2 +
        |2αξ_1ξ_2/δ_1+δ_2+α|^2
        ] + 
        2γ_ϕ^(hf)[ |ξ_1|^2 + |ξ_2|^2 ]

to leading order in ξ_1,2. Here, γ is the transmon loss rate, and γ_ϕ^(hf) is the pure dephasing rate of the transmon due to noise at high frequencies[While a typical Ramsey experiment probes the noise spectrum near zero frequency, it is the spectral content at high frequencies ∼δ_1,2 that is relevant to drive-induced heating via dephasing. ], which was measured to satisfy γ_ϕ^(hf)≪γ in Ref. <cit.>.
For the weak drives considered in this work, the infidelity associated with this additional heating is negligible relative to that from the inverse Purcell effect. As more drives are added, however, this effective heating rate will increase and could limit operation fidelity.


