Circuit quantum electrodynamics enables the combined use of qubits and oscillator modes. Despite a variety of available gate sets, many hybrid qubit-boson (i.e. qubit-oscillator) operations are realizable only through optimal control theory, which is oftentimes intractable and uninterpretable. We introduce an analytic approach with rigorously proven error bounds for realizing specific classes of operations via two matrix product formulas commonly used in Hamiltonian simulation, the Lie–Trotter–Suzuki and Baker–Campbell–Hausdorff product formulas. We show how this technique can be used to realize a number of operations of interest, including polynomials of annihilation and creation operators, namely
Today, many quantum computing architectures are homogeneous, with the same type of qubit used throughout the device. From devices made of superconducting qubits [1][2][3] to ion trap qubits [4], prior work largely focuses on linking qubits of the same type together in faulttolerant ways. However, there is emerging work [5][6][7][8][9] studying the use of heterogeneous quantum computers that leverage two or more types of quantum architectures (e.g. qubits and oscillator modes). Heterogeneous devices hold promise because they can be tailored for speci昀椀c physical simulation problems, which would be especially useful in applications like material discovery [10], molecular simulation [11,12], topological models [13] or lattice gauge theories [14].
In particular, hybrid qumode-qubit models [15] hold some advantages: for example, microwave qumodes have long lifetimes and large accessible Hilbert spaces, making them attractive targets for quantum error correction [15]. Introducing qumodes also enables new physical gates, such as the Mølmer-Sørensen gate [16] while at the same time enabling new forms of transduction between qubits and qumodes [17]. Oscillator interactions have unique features, like nonlinearities, which are challenging to simulate even with homogeneous quantum architectures [8].
Ef昀椀ciently compiling logical operations to physical pulses is a critical, but computationally expensive task. In theory, pulse design techniques like optimal control theory (OCT) [18][19][20] can produce pulses that implement arbitrary quantum transformations on a hybrid qumodequbit system. These techniques have been applied to a variety of physical systems, including NMR [19], superconducting transmon qubits [20], and qumodes [21][22][23]. In practice, OCT is Table 1. Overview of techniques for synthesizing particular unitary transformations and the quantum gates needed. Our formulas allow the manipulations of broad classes of Hamiltonian block encodings, denoted B A , B B . Each row contains the formula used, the target to approximate, the preconditions, and a reference to the location of the precise statement of the performance of the method. Most formulas use 'conjugated' hamiltonians (UHU † ) which can be achieved via conjugation of the exponential. The formula provided denotes hybrid gates with σ i terms, and single-qubit gates are capitalized (e.g. S, X, H). The bounds on the number of gates depend on the accuracy required of the approximation and are given in the corresponding theorems.
Target Preconditions Reference BCH(itσ i B, itσ i A) exp(t 2 [A, B])1 A, B Hermitian theorem 2.4 BCH(itσ j A, itσ k B) exp(-i t 2 σ i {A, B}) A, B Hermitian equation ( 18)
[A, B] = 0, AB = (AB) † equation ( 32)
computationally intensive and in昀氀exible, requiring pulses to be recompiled on a case-by-case basis. Furthermore, OCT is almost always uninterpretable, yielding only a pulse which performs the desired operation without providing any physical intuition. Experimentally, these limitations prevent the high-昀椀delity realization of quantum algorithms; theoretically, the complexity of our quantum circuits often inappropriately ignores the classical cost of required compilation.
Inspired by recent experimental progress [24,25], we introduce an extensible control scheme for a universal, hybrid qumode-qubit quantum computer (table 1). Speci昀椀cally, we show how block-encoded operators can be manipulated using the Lie-Trotter-Suzuki and Baker-Campbell-Hausdorff matrix product formulas. We thus enable the creation of instruction set architectures (ISAs) that can be analytically compiled to experimentally available gate sets. Prior art, namely Jacobs [26], has studied similar techniques to compile operations; we generalize these techniques to work in a variety of domains, including in settings with multiple qumodes and more exotic operators, and prove concrete error bounds on these techniques.
We develop two parallel approaches, one which primarily uses the creation and annihilation operators which we refer to as based on 'Fock methods,' and the other primarily relying on position and momentum operators which we refer to as based on 'phase-space methods. ' We demonstrate that both methods can be used to generate an ISA for qumode-qubit devices. Because our methods operate in both Fock and phase-space, we can achieve transformations that are natively described in either picture; for example, exponentials of polynomials of annihilation and creation operators for Hamiltonian simulation and state preparation in Fock space and controlled parity and beam splitter operators in phase space. Our methods obtain almostlinear asymptotic scaling.
Furthermore, we use the previously mentioned formulas to realize a number of operations of interest, including polynomials of annihilation and creation operators, namely a p a † q for integer p, q. These block-encoded operations are crucial for quantum signal processing (QSP) Figure 1. A qubit-qumode Hilbert space. On the left, the qubit is a typical twodimensional Hilbert space. On the right, the qumode consists of Fock states of value at most Λ = 4. and certain problems in quantum simulation. Finally, we give examples for the Hamiltonian of a nonlinear material and applications to key unsolved problems in quantum simulation such as the Fermi-Hubbard model. While these approaches are expensive in terms of raw gate counts, because they are analytic they provide intuition into synthesizing the gate robustly. Furthermore, these gate sequences can be used as a starting point for optimal control methods, helping to avoid cold start issues.
In this section, we introduce the hybrid qumode-qubit architecture we operate on and the matrix product formulas we will use.
We aim to provide a generic toolbox to build unitary transformations on hybrid qumodequbit devices. Such devices are common in quantum systems, spanning circuit quantum electrodynamics (superconducting qubits coupled to microwave photons) to ion-trap quantum computing (for which the mechanical modes of oscillation are coupled to atomic qubits). The challenge is that fundamentally different insights are needed to compile unitaries in the hybrid setting over typical binary-based approaches.
We review the mathematical properties of qumode quantum mechanics that are needed to understand the basic operations considered for the ISA architecture that we consider. Speci昀椀cally, we present an analytic ISA based on the Lie-Trotter-Suzuki (Trotter) and Baker-Campbell-Hausdorff (BCH) decompositions for decomposition of gates of the form U = e i Ĥσ j , where Ĥ is a Hermitian operator composed of phase-space operators and Pauli gates. Before jumping into the speci昀椀c details of our gate operations, we need to review the basics of qumode quantum mechanics as well as the mathematical results needed to use these qumode operations to compile a given unitary transformation.
We 昀椀rst produce a mathematical description of a qumode, then describe operations which can be performed on the qumode and qubit (昀椀gure 1). Qumodes store bosonic states. Bosons are are already commonplace in quantum computing experiments: photons (energy quanta of the electromagnetic 昀椀eld) are used in photonic chips, cavity QED, and hybrid circuit QED, while phonons (quanta of mechanical vibrations) are used to couple ion-spin qubits in ion traps. However, we are interested in using the qumode as an explicit computational resource, rather than as a conduit for entangling operations or source of noise.
There are two different bases that are commonly used to describe the state of the qumode:
1. Phase-space representation, where operators are written in terms of position (x) and momentum (p) operators 2. Fock-space representation, where operators are written in terms of qumode creation (a † ) and annihilation (a) operators.
In the phase-space representation, the computational basis corresponds to the strength of the electric 昀椀eld in the case of photons (or equivalently the position of a mechanical oscillator for vibrational systems). We refer to this with the operator x and we have that for any x ∈ R, x|x⟩ = x|x⟩. We also use the corresponding operator for momentum p. This describes the magnetic 昀椀eld for a photonic system. In practice, cutoffs are imposed on the values of the 昀椀eld and further discretization error on the gates and the outputs prevents arbitrary precision readout; however, for simplicity we ignore the latter issue in order to provide a simpler if less realistic computational model and ignore the issue that even when cutoffs are imposed the vector space does not strictly form a Hilbert space without also including spatial discretization.
In the Fock-space representation, we track the number of bosons (number of photons or the energy level of the harmonic oscillator for the vibrational case) in the computational basis. In this representation the computational basis is de昀椀ned to be an eigenvector of the boson number operator n|n⟩ = n|n⟩, where n = a † a is the number operator and a and a † add and remove a boson from the system, respectively. Formally this spectrum is countably in昀椀nite, but after truncation it forms a 昀椀nite-dimensional Hilbert space and thus can be thought of as a qudit. For example, assuming a cutoff Λ = 3 on the boson number n
Here P Λ is the projector onto the subspace of the cavity containing at most Λ photons
Observe that this Hilbert space has dimension Γ + 1 corresponding to the Fock states from |0⟩ to |Γ⟩. Truncation of the Hilbert space is required to provide error bounds, otherwise the remainder terms become unde昀椀ned. Provided that an appropriate cutoff is picked for the system, the discrepancies between the truncated and untruncated systems will often be negligible. For notational clarity, we assume a cutoff of Λ for all further equations and assume the annihilation and creation operators implicitly have the projectors P Λ .
To incorporate the qubit's state, we take the tensor product of the qubit and qumode Hilbert spaces, so that the state space is H 2 ⊗ H Λ+1 . This means that, for example, a computational basis state will be of the form |q⟩ ⊗ |m⟩ where the state |q⟩ here can be thought of as the union of the qubits in the system, and |m⟩ represents a qumode state where the system is either in position x = m for the phase-space encoding or has m photons if the Fock-space encoding is used.
2.1.2. Representing operations and measurements via block encodings.
Depending on whether the qumode is described in Fock or phase space, there exist two pairs of complementary qumode operators. In Fock space, the operators are creation and annihilation operators a † , a; in phase-space, the operators are position and momentum x, p.
Our techniques can be applied in both spaces because Fock and phase-space operators have the equivalencies
and commutation relations
To illustrate the qumode-qubit ISA, we de昀椀ne three types of operations with examples:
1. Qubit-exclusive: these include typical qubit gates like the Pauli (X, Y, Z) gates,
and Hadamard (H) and phase (S) gates
2. Qumode-exclusive: we assume that linear optical operations (which are at most quadratic in the 昀椀eld operators) can be performed on the qumode. This includes the displacement operations e αa † +α * a for ³ ∈ C, phase delays (or phase-space rotations) e -iαa † a , and squeezing operations e α(a † ) 2 -α * a 2 . Qumodes can also be entangled with other qumodes via beamsplitter operations e (αa † b-α * ab † ) where b is the creation operator acting on a different qumode. In this work, we primarily focus on the single-qumode case. The multi-mode case with beamsplitters is discussed in [5,14] but without the rigorous convergence bounds provided in the present work. 3. Qumode-qubit entangling: there are several entangling operations between individual qumodes and qubits that widely appear in the circuit QED literature. Two common operations we consider are the conditional displacement operation [25] e -iσ z ⊗(αa † +αa) and the Selective Number-dependent Arbitrary Phase (SNAP) gate [24] e -iσ z ⊗ ∑ n αn Pn , where Pn = |n⟩⟨n| is the projector onto the nth Fock state.
Note that, for clarity, we use uppercase letters for qubit-exclusive gates and lowercase letters for qumode-exclusive or qumode-qubit gates. For example, we use S, X, H for qubit-exclusive gates and σ i notation for Paulis in hybrid gates.
Our compilation strategy describes operations in terms of qubit-exclusive gates and the S 1 gate, a primitive qumode-qubit entangling gate (de昀椀nition 2.1). S 1 is a useful 'block-encoding' primitive to compose complex gates because it embeds a 昀椀rst-order Fock operator in the offdiagonal blocks. 'Block-encoded' matrices refer to how 'blocks'/submatrices of a larger matrix can be expressed as an existing matrix. Block encodings are frequently used in quantum algorithm design [27,28]. Definition 2.1 (S 1 primitive gate). For any t > 0 and any positive integer cutoff Λ, we de昀椀ne S 1 to be the unitary acting on the Hilbert space H 2 ⊗ H Λ+1 that has the following representation as a block matrix
Note that S 1 can itself be decomposed into conditional displacements (see appendix A) or implemented directly via OCT [29][30][31][32]. Also note that the block encoding can also be expressed as the sum of qumode-qubit tensor products
We also consider the broader class of block-encoded Hamiltonians:
Definition 2.2 (Block encodings). For a qumode operator A acting upon H Λ+1 , we denote the joint block-encoded Hamiltonian to be
so that the subscript is the upper right block and the lower left block is the transpose and complex conjugate to preserve Hermiticity. I.e., B A is Hermitian for any A, and thus is a suitable Hamiltonian. Note that B A can describe a qumode-qubit Hamiltonian acting on
In this block-encoding notation, observe that S 1 = exp itB a † . Throughout this work, we consider increasingly exotic A matrices that can be created via polynomials of qumode operators.
Finally, to simplify notation, we oftentimes write tensor products implicitly as follows
This notation is used extensively to abbreviate the action of σ i on H 2 and M (which is comprised of qumode operators) on H Λ+1 . We further assume that the qubit can be measured directly, but the qumode can only be measured by entangling it with a qubit and reading out the state of the qubit to obtain a single classical bit of information about the qumode state [5].
Matrix product formulas describe the behavior of products of matrix exponentials (namely e A e B ). These formulas are well-known in Hamiltonian simulation [33][34][35] and are used to approximate a discretized version of the time evolution operator e -iHt using the Trotter formula. Thus, they often have rigorous error bounds that describe how e -iHt can be approximated given constituent e -iHjt . In our setting, e -iHt will be the operation we seek to implement using the S 1 gate.
As stated above, we use two product formulas: the BCH formula and the Trotter formula (昀椀gure 2). The BCH formula is used to create a commutator (or anticommutator) of operators. The Trotter formula is used to add these commutators and anticommutators together. We introduce the informal theorems below: Theorem 2.3 (Informal Trotter theorem from [34]). Suppose we may implement e Mλ and e Nλ for arbitrary λ ∈ R and anti-Hermitian M, N. Then, a pth order Trotter formula has the error scaling
requiring no more than 4 • 5 p-1 exponentials.
Theorem 2.4 (Informal BCH theorem from [36]). Suppose we can implement the operators e Mλ , e Nλ for λ ∈ R and anti-Hermitian M, N. Then, a BCH formula of order p has the error scaling
requiring no more than 8 • 6 p-1 exponentials.
These product formulas are de昀椀ned recursively and are comprised of sequences of e Mλ , e Nλ gates evaluated at varying values of λ. While the sequences may be long, as there is an exponential dependence on p, it is important to note that these error bounds are known to be loose, so experimental implementations may bene昀椀t from improved accuracy.
Finally, note in our notation these formulas take the logarithm matrix and output sequences of matrix exponentials which approximate a desired exponential. Additionally, note that Hamiltonians can be 'conjugated', i.e. we can transform a Hamiltonian from H to UHU † for any unitary U via:
Conjugation, Trotter, and BCH give us immense 昀氀exibility for the Hamiltonian manipulations which can be achieved (table 1).
In this section, we formalize our technique to build polynomials of block encodings in both Fock and phase space. This is notable because it shows how Hermiticity is not a requirement for our compilation scheme: whereas phase-space operators x and p are Hermitian, Fock operators a, a † are decidedly non-Hermitian. We 昀椀rst show how our hybrid qumode-qubit architecture allows for the synthesis of anticommutators of Hermitian operators and, by proxy, matrix products of phase-space operators. We then use similar techniques with non-Hermitian operators such as Fock-space operators to manipulate block encodings of matrices. Finally, we contextualize these methods with asymptotic error bounds, providing theoretical analyses of our proposed techniques.
Our qumode-qubit architecture uses BCH to natively implement anticommutator-like exponentials by exploiting the qubit's Hilbert space. In conjunction with Trotter, these anticommutators are used to build larger block encodings; we generate nontrivial transformations like beamsplitters and the Hong-Ou-Mandel effect in section 4.
To begin, we recall the qumode-qubit commutators
where ϵ ijk is the Levi-Civita symbol, and i, j, k ∈ {x, y, z}. We use these relations, as well as the Pauli product identity
to decompose the anticommutator σ i {A, B} in terms of a qumode-qubit commutator
Note this assumes A, B ∈ H Λ+1 are Hermitian so that exp[i σ j A] is unitary and commutes with σ i , as is the case when A, B are mode-only operators. Thus, by using a hybrid qubit-cavity operation of the form exp i Aσ j , exp i Bσ k , the BCH formula can convert commutators into anticommutators.
Finally, because 1 2 [A, B] + 1 2 {A, B} = AB, we use the Trotter formula to produce
assuming we may implement exp 1 2 σ i [A, B] via a traditional BCH formula.
In this section, we show how to achieve A q for an arbitrary qumode operator A. We also show that our techniques work in a multi-qumode setting. This extends the prior techniques, which require A to be Hermitian. When A is a Fock space operator, this corresponds to realizing arbitrary powers of a and a † , which is known to generate a universal set of operations on the qumode. This is useful for the simulation of nonlinear materials, which naturally lead to terms that are polynomial in a, a † , as well as QSP [27].
Our method again uses the qubit coupling to induce a phase in the compound qumodequbit system, similar to the previous section. We begin with block encodings as described in equation (11). To manipulate the block encodings, we begin by recognizing that qubit-only operations can modify the exponential via 'conjugation,' i.e.
Thus, given any block encoding, we can also create the auxiliaries
recalling that S is a qubit phase gate. Applying BCH yields the commutators
These commutators themselves can be conjugated. Recall that HZH = X and
Thus, using only B A (t) gates, we can approximate B A 2 . We now lift this procedure to produce exp itB AB given exp itB A , exp itB B for commuting A, B ∈ H Λ+1 . Observe the following commutators (whose exponentials we can implement via BCH)
Provided that [A, B] = 0 and via conjugation, the equality simpli昀椀es to
Via Trotter, we can directly implement the sum:
We select τ = t 2 to obtain the desired time and conjugate by σ x to produce the desired matrix. This procedure is described in algorithm 1 and thus allows us to approximate exp i tB AB . This process can be repeated iteratively, assuming AB commutes with B; for example, if A = B = a, then this process can be used to produce higher powers a k , (a † ) k .
Our formulas require [A, B] = 0 to build higher order polynomials. This requirement is tolerable, as we still may achieve a broad class of transformations, including homogeneous polynomials of a or a † . Our formulas can also be extended to more general cases where the operators commute, e.g. when the synthesized unitary operates on two different qumodes, as in the conditional beamsplitter (a gate that acts as a beamsplitter controlled on an ancillary qubit).
Algorithm 1 is an extension of the prior commutator approaches in phase space because the σ i = -i 2 [σ j , σ k ] relation is natively expressed in the algorithm; i.e. if we have B A = B B = B x = exp itxσ x , the 'Left' term vanishes and the'Right' term is the commutator we would apply.
Finally, in algorithm 2, we demonstrate how to implement AB if AB = (AB) † . This process places AB in the upper left block, which is useful to exactly execute exp itAB by preparing the qubit state to |0⟩, but it prevents the process from being executed recursively. Simplifying the commutator from equation ( 29) 昀椀nds
The prior description of our algorithm assumes errorless product formulas. However, the BCH and Trotter formulas introduce errors which must be accounted for, especially when applying our algorithm recursively. In this section, we cite the error scaling of the general addition algorithm described in algorithm 1 and the multiplication algorithm described in algorithm 2.
The proofs and full results are included in appendix B.
Input:
for some constant c and order p A , p B ⩾ 1 where [A, B] = 0. Then, the application of algorithm 1 will yield the scaling
with C TOTAL = max(∥AB∥ , ∥BA∥ , C 2 BCH ) and C BCH = max(∥A∥ , ∥B∥ , c), using no more than 1.07 • 30 q exponentials, where q = max(⌈ min(p1,p2)-1 2 ⌉, 1).
for constant c and p A , p B ⩾ 1 where AB = (AB) † . Then, algorithm 2 has the error
with C = max(∥A∥ , ∥B∥ , c), using no more than 8 • 6 q-1 exponentials where q = max(⌈ min(pA,pB)-1 2 ⌉, 1).
While the asymptotic error analysis suggests that the cost of this method is onerous, we note that the product formulas often have overly pessimistic error scaling and operation counts [37]. In the applications below, we provide numerical simulations which suggest our technique is more readily implementable than theory suggests.
In this section, we show how our technique is a powerful tool for analytically realizing desired operations. This technique succeeds both for Hamiltonian simulation problems and general control problems. In particular, we show how the aforementioned physical intuition for a desired transformation is often suf昀椀cient to produce an approach to create desired operations.
We include applications in both phase and Fock space:
1. Phase-space techniques are demonstrated to be useful in the case where displacements (e (αa † +α * a) = e i αx for ³ real or e αp for ³ imaginary) are the only experimentally available gates. We produce the controlled parity operator e itσ z a † a (section 4.2); the beamsplitter e -itσ z (a † b+ab † ) (section 4.4); gates for two encodings of universal control of the restricted span{|0⟩ , |1⟩} Hilbert space (appendix E); and gates for simulation of Fermi-Hubbard lattice dynamics using the two lowest Fock states of the cavity (appendix F), including samesite, hopping, controlled-beamsplitter (appendix C.2), and FSWAP gates. 2. Fock space techniques are shown to be useful assuming compilation to S 1 (de昀椀nition 2.1) and single-qubit operations to produce polynomials of annihilation and creation operators, (namely a p a † q for integer p, q). We demonstrate how polynomials of these operators can be used in Hamiltonian simulation (e.g. with χ (3) nonlinear materials, appendix D.2) and state preparation (appendix D.3).
As a simple application, let us consider the case of simulating a χ (3) nonlinear material. These interactions commonly occur in nonlinear optics and appear when the index of refraction for a material varies linearly with the intensity of the electromagnetic 昀椀eld. Such interactions can be modeled for a single qumode using the expression
Our goal here is to examine the cost of a simulation of such a Hamiltonian in our model for time t and error tolerance ϵ and to determine the parameter regimes within which a hybrid simulation using our techniques could provide an advantage with respect to a conventional qubit-based simulation of the Hamiltonian.
Each of the terms can be approximated using formulas from table 1. The ωa † a term requires an embedding of Hermitian a † a, so that
The second order term is treated in the same way, noting that B (a † ) 2 can be produced via table 1, so that
Thus, via the BCH formula, we can block-encode the two Hamiltonian terms. Trotterizing allows us to block-encode the entire Hamiltonian into the upper-left quadrant. Thus, by setting the qubit to |0⟩, we can approximate the Hamiltonian. The error scaling is as follows and is proven in appendix D.2:
Hamiltonian:
(i.e. a Hamiltonian with a Kerr non-linearity). Let t be the evolution time and ϵ be the target error tolerance. For any positive integer q we can approximate an exponential of the blockencoded Hamiltonian with error at most ϵ in the operator norm using re O(q) S 1 operations where r ∈ Ω (Λ 4 t) 1+1/(q-3 4 )
.
This shows that we can perform a simulation of the dynamics within error ϵ using a number of operations within our instruction set that scales near-linearly with the evolution time and subpolynomially with ϵ. Further, this approach requires no ancillary memory and can be done with a single qumode and a qubit. In contrast, a qubit-only device would require a polylogarithmic number of qubits in Λ.
It is worth noting that in this case the ancillary qubit is not being used directly in the model. Instead it is being used to control the dynamics and generate the appropriate nonlinear interaction between the photons present in the model.
We now demonstrate how the approach can extend beyond problems in Hamiltonian simulation. We begin with an example of the technique for control: In particular, we seek to perform a nondestructive measurement of the qumode in which we project the information into the qubit [11,38].
To construct such a nondestructive measurement, we seek to implement e itnσ z where n = a † a is the number operator. If we could implement this gate for arbitrary t, we could perform phase estimation on the qubit to nondestructively project the qumode into a 昀椀xed number of bosons. This could be done by setting t suf昀椀ciently small so that tΛ ⩽ 2π is calculable with phase estimation. Alternatively, for t = π, this operation checks the parity of the qumode and applies an RZ gate for odd parities. We employ the instruction set in the phase-space representation to synthesize the in昀椀nitesimal conditional rotation gate
for k = x, y, z. We rewrite n in terms of the phase-space operators by recognizing:
Applying equations ( 3), (4), and (45) yields
such that the gate is expressed via the Trotter decomposition as the product of exp
Given the Pauli commutator relation, the 昀椀rst commutator is
and the second commutator is
such that both terms are amenable to BCH decomposition, and the in昀椀nitesimal conditional rotation is composed with a gate-depth lower bound of nine. To perform an error analysis, we may directly apply the error scaling of BCH and Trotter to 昀椀nd:
Theorem C.1. Suppose we can implement e itxσ i , e itpσ i without error. Then, we may approximate exp itB x2 +p 2 with arbitrary error scaling p as
where C = max( x2 + p2 , ∥x∥ 2 , ∥p∥ 2 ) and using no more than 4 • 5 p 2 -1 4 exponentials. We then provide numerics in 昀椀gure 3. As expected, the wavefunction initialized in the second excited state of the cavity and the ground state of the associated qubit has an autocorrelation function that oscillates with phase exp(2it). Dynamics are well reproduced with 2000 time steps for a 昀椀nal time of 20 with cutoff Λ = 14. Note the units are arbitrary in the absence of de昀椀nition of the cavity frequency ω, with the only units de昀椀ned by setting the reduced Planck constant to unity ℏ = 1arbitrary units. The close agreement between the BCH-synthesized and exact gates is supported by the error scaling after a single gate application computed for time step t, which features a power law scaling in agreement with the predicted error scaling for both BCH and Trotter decompositions. We can obtain a similar decomposition with Fock-space operators. Observe that the MULT subroutine applied to the exp itB a , exp itB a † using theorem 3.2 can also yield
when B a 's implementation is error-free. Note that aa † = a † a + 1, so the block encoding that is actually applied is
Thus, our Fock-space methods would also achieve the same transformation, albeit requiring a phase and RZ correction.
Consider the case where we seek to prepare Fock state |k⟩ b on the qumode. On qubit devices, preparing integer states is trivial, assuming the qubits represent logic in a binary fashion. However, hybrid boson-qubit devices natively implement exponentials of the phase-space or Fock-space operators. Thus, preparing |k⟩ b directly can often be challenging. Existing work [39] enables the preparation of states, but are less granular in their control of the qumodequbit state. We present a method that only swaps the |1⟩|0⟩ ↔ |0⟩|k⟩ states for integer k, otherwise leaving the initial state intact.
To begin, we aim to implement (a † ) k on the vacuum. It is suf昀椀cient to approximate
which we call the 'unprotected' state preparation operator. Selecting appropriate t yields precisely the desired behavior, which gives the following result:
Theorem D.5 (Unprotected state preparation). For k ⩽ Λ, we can take t = (2n + 1)
While T k (t) performs the desired transformation, it may incur unwanted side effects if the starting state is of the form |1⟩ ⊗ |b⟩ for b > 0. We can use our same approach to produce the following operation: Theorem D.6. Consider the Fock preparation unitary P k with the form
When t = (2n + 1) π
√ k! , we have that P k performs our desired state preparation
We claim that we can approximate this unitary with P k,p where
The proofs of theorems D.5 and D.6 are provided in appendix D.3. Though this subroutine appears expensive, numerical results suggest it is far more implementable than theory would suggest. In the following simulations, we apply the above technique but always use a second-order symmetrized BCH formula and second-order (symmetrized) Trotter formula. This amounts to 480 exponentials for the unprotected case and 960 exponentials for the protected case. The synthesized gates are provided in 昀椀gures 4 and 5.
As can be seen in 昀椀gure 6, the protected gate has higher compilation errors than the unprotected gate. This is likely due to the fact that the protected formula requires an additional application of the Trotter product formula. This will further incur commutator error over the unprotected gate.
We also analyze the error scaling as the order of the BCH formulas used increases. Figure 7 describes the error-resource tradeoff as the Trotter step within each BCH formula increases. Observe that the error decays as higher order formulas are used or Trotter step size is reduced, as expected. While modest depths have relatively high in昀椀delity, in theory our compilations can achieve arbitrary accuracy at the level of both entire gates and matrix individual elements. The operations we seek to realize need not act on a single qumode; in fact, our techniques are extensible to hybrid setups with multiple qumodes or qubits. Consider the conditional (controlled-phase) beam splitter
This gate naturally pertains to certain lattice gauge theories [14] and gives rise to exponential SWAP (eSWAP) [40] and controlled-SWAP (cSWAP) gates for state puri昀椀cation and SWAP tests, when paired with an ordinary (uncontrolled) beam splitter [41]. The argument in phasespace representation is
Figure 6. Absolute element-wise error when compiling T 2 gate via scheme: left is when the compilation target is the unprotected gate, while right is when the target is the protected gate. The protected gate has higher synthesis errors. such that the gate is decomposed in terms of a Trotter expansion as the product of two exponential terms exp
According to the Pauli commutation relation, the 昀椀rst commutator is
and the second is
with the following error scaling:
Theorem C.2. Assume we may implement e i txmσ j , e i tpmσ j for m ∈ {1, 2}; i.e. we may implement the qubit-conditional position shifts and momentum boosts on either qumode without error. Then, we may approximate B x1 x2+ p1 p2 with arbitrary error scaling p as
2 ) and using no more than 4 • 5
p 2 -1 4 exponentials.
The two exponential terms are decomposed via the BCH formula for a lower-bound gate depth of eight. Results are shown in 昀椀gure 8 for 15 states per cavity with a shared qubit over a 昀椀nal time of π/2 with 200 equal time steps, where the system is initially in the 昀椀rst excited state of each cavity and the ground state of the shared qubit |11g⟩. As expected for the conditional beam splitter, the gate exhibits the Hong-Ou-Mandel effect, in which the occupation of cavity 1 oscillates between the 昀椀rst excited mode and a superposition of the ground and the second excited states of the cavity. The BCH-synthesized results closely agree with that of the original gate, with no visible leakage beyond the physical states (the lowest three states of the cavity) into the working space under the time duration studied. As for the conditional rotation gate, the relative error of the BCH-synthesized gate computed for a single time step of length t was found to scale according to a power law with the time step, in accordance with the analytic result for Trotterization and BCH decomposition.
Our main contribution in this paper is a systematic approach to synthesizing unitary dynamics on a hybrid quantum computer that has access to both qubit and qumode operations. Such gate sets naturally model systems such as cavity quantum electrodynamics systems and ion trap-based quantum computers. Our main innovation here is the development of high-order analytic formulas that can be used to place bounds on the complexity of implementing arbitrary unitary operations on such a hybrid device. Speci昀椀cally, we see that these methods are capable of achieving subpolynomial scaling with the inverse error tolerance (1/ϵ) and allow us to implement arbitrary nonlinearities in the 昀椀eld operators in the generator of the unitary that we wish to implement at low cost asymptotically. In particular, we focus on using a construct known as a block-encoded creation operator as our fundamental construct and show numerically highly accurate approximations to the exponential of a block encoding of the square of the creation operator. Further, we study the Hong-Ou-Mandel effect and observe that the synthesized operations used in our construction can have negligible error with respect to our target precision. Beyond compilation, our work thereby provides intuition when compiling a variety of gates. For example, our methods give mathematical understanding for the recent successes of echoed conditional displacement (ECD) gate sets for quantum computing [25].
While this work enables better analytic control of qubit-qumode systems, there remain many open questions. In particular, lowering the resource cost of these compilations, especially in noisy environments, is key to practical implementations. We note that for both the Hong-Ou-Mandel effect and the block encoding of (a † ) 2 that thousands of gate operations are needed to achieve in昀椀delities of 10 -3 or smaller. This makes such sequences impractical for near-term applications where the gate in昀椀delities are on the order of 1%. Ideally, these analytic sequences can be used as an initial sequence which can be further optimized via existing optimal control techniques. This would reduce warmstart challenges and potentially improve convergence.
Another concern arises from the nature of the qubit system, namely that creating qubitqumode interactions requires truncating the qubit system to a two-level system. For example, transmon qubits are typically truncated to two levels [15,42,43]; in reality, these transmons can experience leakage to higher states, meaning that the qubit may itself exhibit qumodelike properties when poorly calibrated. While leakage is a practical concern, it can be somewhat avoided via well-de昀椀ned pulse shaping methods (e.g. DRAG [44] pulses). These methods allow one to approach the control speed limit set by the anharmonicity of the level structure.
Several avenues of approach exist that could be used to improve upon these results: 昀椀rst, note that our technique is connected to ideas from QSP [27,45]. Thus, QSP could be applied to approximate polynomials of qumode operators; this could potentially improve the scaling with respect to the error tolerance from the subpolynomial scaling currently demonstrated to polylogarithmic scaling. Alternatively, our approach could seed gradient-descent optimization procedures for control such as GRAPE [19] at the pulse level or numerical optimization of parameterized gates at the SNAP [46] or ECD [25] instruction level. These locally optimized sequences may then prove to be either better, or more understandable, than existing gradient-optimized pulse sequences for control of such systems. Designing pulses is, in practice, sensitive to device properties and noise environments, and presents a rich area for further research.
We demonstrate how to obtain the S 1 operator with the native gates present in the dispersive coupling regime of the Jaynes-Cummings model. Astute readers may recognize the S 1 operator as being a Jaynes-Cummings Hamiltonian when rewritten in the form
It is thus unsurprising that we may implement these types of operations on a circuit QED device that exhibits a similar native Hamiltonian.
The conditional displacement operator with magnitude ³ is written:
For a Fock state |n⟩, a † a|n⟩ = n|n⟩; therefore e ia † aθ a † e -ia † aθ = e iθ a † . Taking ³ = ³ * , we have
This operation can be built using single-qubit operations on a controlled displacement gate with additional phase delays on the oscillator. Both are linear optical operations or singlequbit operations, which we expect to be inexpensive in our computational model. Next, note that e i(α(a † +a))
Thus to O(³ 2 ) the block-encoded creation operation can be constructed using single-qubit, controlled-displacement, and linear-optical operations through S 1 ≈ e i(π/2)a † a e i(α(a † +a))⊗σ y e -i(π/2)a † a e i(α(a † +a))⊗σ x . (66)
Figure 9 describes the circuit.
To assess the error scaling of our approach, we must consider three sources of error: the underlying implementation error from using approximations of B A (t), the error from BCH, and the error from Trotter. In appendix B.1, we show that the BCH and Trotter formulas can still be applied on exponentials that have error. Then, we use these formulas in appendix B.2 to produce the error bounds for addition. Finally, in appendix B.3, we produce the error bounds for multiplication.
We begin by formally stating the Trotter and BCH formulas when there is no implementation error:
Theorem B.1 (BCH Product Formula (Theorem 2 from [36])). Let A and B be bounded complex-valued matrices and assume without loss of generality t ∈ R + is assumed for the purposes of asymptotic analysis to be in o (1). We then de昀椀ne
We further de昀椀ne BCH p,k recursively for p ⩾ 2 and odd k ⩾ 1:
with the constants
This recursive formula has the error scaling, where µ = max(∥A∥ , ∥B∥
and uses 8 • 6 p-1 exponentials when k = 1 and 4 • 6 p-1 exponentials otherwise.
. Let {H j : j = 1 . . . m} be a set of M bounded Hermitian operators acting on a Hilbert space of dimension 2 n and assume without loss of generality that t ≥ 0. For H = H j , the error in the Trotter-Suzuki formulas of order k and timestep r obeys the error bound
where τ = ∥H∥ t and
using no more than 2m5 k-1 r exponentials. We de昀椀ne the kth order Trotter formula as
where p k = (4 -4 1/(2k-1) ) -1 and k > 1. The relation has the base case
e Hjλ/2
This implies the following corollary:
Corollary B.1. If r = 1, i.e. there is no time stepping, the Trotter formula exhibits the error scaling
using no more than 2m5 k-1 exponentials.
Both of these formulas, however, assume that our implementation of exponentials occurs without error. However, if we would like to apply our technique recursively, our matrix product formulas must account for implementation error, i.e. be able to use primitives that themselves may have error. Thus, we restate both Trotter and BCH when the operators have asymptotic error:
Lemma B.1 (BCH under implementation error). Suppose there are some ideal operators U 1 (t), U 2 (t) that are exponentials of some anti-Hermitian matrix, i.e.:
We seek to build exp t 2 [A 1 , A 2 ], the exponential of the commutator of the matrices. Also suppose that may approximate U 1 (t), U 2 (t) with U 1,p1 (t), U 2,p2 (t) with the error scaling
for some p 1 , p 2 ⩾ 1. Call the logarithms of Ũ1,p1 , Ũ2,p2 to be Ã1,p1 (t), Ã2,p2 (t). Then, applying a qth-ordered BCH formula, where q = max(⌈ min(p l ,p2)-1 2 ⌉, 1) on the implementable U 1,p1 (t), U 2,p2 (t) can still approximate the commutator exponential by applying theorem B.1:
where C = max(∥A 1 ∥ , ∥A 2 ∥ , c). This procedure uses 8 • 6 q-1 total exponentials.
Proof. Recognize that we may decompose the error involved in implementing the commutator exponential into two sources: the error incurred from the BCH formula intrinsically and the implementation error from the realizable terms. Thus, by the triangle inequality
We begin with the LHS term. By theorem B.1,
where C BCH = max(∥A 1 ∥ , ∥A 2 ∥). For the RHS, recall by Box 4.1 of Nielsen and Chuang [47] that implementation errors accumulate at most linearly; thus, we can sum over the 8 • 6 q-1 operations used by BCH. By symmetry of the BCH formula, we apply the U 1,p1 (t), U 2,p2 (t) exponentials precisely 4 • 6 q-1 times. Thus,
Setting q = max{⌈ min(p l ,pr)-1 2 ⌉, 1} and recalling µ ⩾ 1, we observe
as desired.
Given lemma B.1's assumptions, theorem B.2's assumptions, and assuming ∥A 1 + A 2 ∥ ⩾ 1, an operator can be constructed by applying a qth-ordered Trotter formula Trotter 2q by setting q = max(⌈ min(p l ,p2)-1 2 ⌉, 1) so that
where C = max(∥A 1 + A 2 ∥ , c) and using no more than 4 • 5 q-1 operator exponentials.
Proof. Similarly, we may apply the triangle inequality in order to determine a bound by separating the error accrued into the intrinsic Trotter error and the implementation error
To analyze the LHS, which represents the Trotter error, corollary B.1 provides a bound
To analyze the RHS, which represents the implementation error, recall by box 4.1 of Nielsen and Chuang [47] that the error accrues linearly. Furthermore, the number of operations in Trotter is no more than 4 • 5 q-1 operations in total. Therefore, we only apply each constituent operation at most 2 • 5 q-1 times. Thus, a loose upper bound can be written as
It is suf昀椀cient to set q = max(⌈ min(p1,p2)-1 2 ⌉, 1) so that
as desired.
We apply the above results to produce the error analysis of algorithm 1:
Theorem 3.1. Suppose we have approximations B A (t), B B (t) with the error scaling
for some constant c and order p A , p B ⩾ 1 where [A, B] = 0. Then, the application of algorithm 1 will yield the scaling
with C TOTAL = max(∥AB∥ , ∥BA∥ , C 2 BCH ) and C BCH = max(∥A∥ , ∥B∥ , c), using no more than 1.07 • 30 q exponentials, where q = max(⌈ min(p1,p2)-1 2 ⌉, 1).
Proof. Our proof proceeds by applying the above theorems upon our operations. By setting q = max(⌈ min(p1,p2)-1 2 ⌉, 1), lemma B.1 implies that:
where C BCH = max(∥A∥ , ∥B∥ , c) and p BCH = min(p A , p B ). We'll need to set τ = t 2 to achieve the desired time evolution. Additionally, Pauli conjugation has no impact on the error scaling. We call these formulas 'Left' and 'Right'
where
Finally, we use lemma B.2, implying that a Trotter formula with order q has the error scaling
where
To bound the number of operations used, we recognize that the Trotter formula requires at most 4 • 5 q-1 commutators, each of which requires 8 • 6 q-1 constituent operators. Thus, the total number of operators required is bounded as
This implies the following corollary for the annihilation/creation operators:
Assume we can implement the following k l , k r th order approximations of S with error scaling p l , p r :
S kr,pr (t) -S kr (t
with c ⩾ Λ max(k l ,kr)/2 . Then, we can implement higher-order operators with comparable t scaling
using no more than 1.07 • 30 q S k l , S kr operators.
Proof. To synthesize the block encoding of higher-order annihilation/creation operators, we directly apply theorem 3.1. We quantify the error by bounding the block-encoding norm. Note that:
A kth order block-encoded operator has a bounded norm
Thus, the constant C BCH is bounded, as C BCH ⩽ max(Λ k l /2 , Λ kr/2 , c) ⩽ max(Λ max(k l ,kr)/2 , c) = c by hypothesis. Next, observe that ∥AB∥ , (AB) † ⩽ Λ (k l +kr)/2 . Thus, C TOTAL ⩽ max(Λ (k l +kr)/2 , c 2 ) = c 2 . Therefore, the 昀椀nal error scaling will be upper bounded by O((c 2 t) min(p l ,pr)/2 ).
Algorithm 2's error scaling follows directly from the BCH formula: Theorem 3.2. Suppose we have approximate block encodings B A , B B with the error
for constant c and p A , p B ⩾ 1 where AB = (AB) † . Then, algorithm 2 has the error
with C = max(∥A∥ , ∥B∥ , c), using no more than 8 • 6 q-1 exponentials where q = max(⌈ min(pA,pB)-1 2 ⌉, 1).
Proof. We can directly apply lemma B.1 using the B A , B B operators. When applied, we 昀椀nd that
where C = max(∥A∥ , ∥B∥ , c) and q = max(⌈ min(pA,pB)-1 2 ⌉, 1). Then, by taking τ = t 2 , we yield
By counting the number of exponentials in the result via lemma B.1, we 昀椀nally 昀椀nd that the number of exponentials needed is at most 8 • 6 q-1 .
In the following section, we derive error bounds for the two phase-space applications described: the conditional rotation gate and the controlled-phase beam splitter. Note that the cutoff approach we employ to bound ∥a∥ , a † is more complex for position and momentum operators.
To obtain error bounds, we leave all expressions in terms of ∥x∥ , ∥p∥. We leave a more concrete bound that can be found by applying a cutoff upon both x, p simultaneously to future work.
Theorem C.1. Suppose we can implement e itxσ i , e itpσ i without error. Then, we may approximate exp itB x2 +p 2 with arbitrary error scaling p as
where C = max( x2 + p2 , ∥x∥ 2 , ∥p∥ 2 ) and using no more than 4 • 5 p 2 -1 4 exponentials. Proof. Begin by directly applying theorem B.1 to identify the error scaling. This implies that
Without loss of generality, select σ i = σ y and σ j = σ z so that [σ i , σ j ] = 2i σ x . Then, by selecting τ = t 2 , the BCH formula is an approximation of exp itB x2 . Thus,
Similarly, for p2 ,
We apply lemma B.2 and observe
where C = max( x2 + p2 , ∥x∥ 2 , ∥p∥ 2 ) and q = ⌈ p 2 -1 4 ⌉. This requires no more than 4 • 5 q-1 operator exponentials, thus implying
C.2. Controlled-phase beam splitter gate Theorem C.2. Assume we may implement e i txmσ j , e i tpmσ j for m ∈ {1, 2}; i.e. we may implement the qubit-conditional position shifts and momentum boosts on either qumode without error. Then, we may approximate B x1 x2+ p1 p2 with arbitrary error scaling p as
2 ) and using no more than 4 • 5
p 2 -1 4 exponentials.
Proof. We may take a similar approach as above. Applying theorem B.1 yields
where we set ∥x∥ = max(∥x 1 ∥ , ∥x 2 ∥) and ∥p∥ = max(∥p 1 ∥ , ∥p 2 ∥). We again take τ = t 2 so that
Applying lemma B.2 gives
where
) and q = ⌈ p 2 -1 4 ⌉. This requires no more than 4 • 5 q-1 operator exponentials, thus implying
We now introduce a series of techniques that allow us to realize polynomials of Fock-space operators. We 昀椀rst begin in appendix D.1 by identifying the error scaling of an arbitrary order Fock-space block encoding, i.e. B a k . In appendix D.2, we show how the techniques can be used to simulate the Jaynes-Cummings Hamiltonian, which itself is a polynomial of Fockspace operators. In appendix D.3, we demonstrate how this technique can be extended beyond simulation into realizing more general operators, such as a unitary for state preparation.
We seek to demonstrate the following result:
Theorem D.1. For positive integer k ⩾ 1 and timestep t ∈ R, we seek to implement the target block encoding T k (t) de昀椀ned as
For any ϵ > 0 and p > 1 there exists an implementable unitary operation T k,p of order p such that
Algorithm 3. POWER(k, t, p) .
Input: k = 2 ℓ for nonnegative integer ℓ, timestep t > 0, order p > 1
and the number of applications of S 1 (t) needed to implement the operation scales in r • n 1.6 30 np 420 n 2 p/2 6 log 2 n+1 , (124
where r ∈ Θ (Λ k/2 t) 1+1/(p-1)
.
Because we can add two lower-order block encodings via algorithm 1, we can exploit a binary expansion to achieve arbitrary orders (e.g. (a † ) 9 = (a † ) 2 3 a † ). Thus, our 昀椀rst task is to demonstrate the implementation of these block encodings with orders that are a power of two. This is achievable through the recursive POWER algorithm:
This builds to the following result:
Theorem D.2. For any t ⩾ 0, p ⩾ 1, and 昀椀xed k = 2 ℓ for some ℓ ⩾ 1 we have that the unitary implemented by algorithm 3, POWER acting on
using no more than 6 log 2 k • 420 kp/2 unitary S 1 operators.
To bound the error of this algorithm, we begin by identifying the implementation error of the second-order formula, i.e. S 2 , the 昀椀rst operator with implementation error: Lemma D.1 (Implementing second-order block encodings). Suppose we can implement the following operation without error (as de昀椀ned in de昀椀nition 2.1 and subject to a bosonic cutoff)
Then, we can approximate S 2 (t) to the pth order, i.e. implement S 2 such that
using no more than 6 • 14 p S 1 (t) operations.
Proof. Note that, if S 1 (t) is errorless, then we only need to account for error incurred by the BCH and Trotter formulas. By employing a pth order BCH formula we can produce commutator exponentials with error O((Λ 1/2 τ ) 2p+1 ) by theorem B.1 and fact B.1.
We again set τ = t 2 so that the error scales in at worst O((Λt) p+ 1 2 ). Then, we apply a Trotter formula lemma B.2 of order ⌈ p 2 ⌉ so that
with C ⩽ max(Λ, Λ 2/2 ) = Λ. Our worst case error scaling is then O((Λt) p+ 1 2 ). This requires no more than 2 • 2 • 5 ⌈ p 2 ⌉-1 of the commutators, each of which required 8 • 6 p-1 昀椀rst order operations. Thus, the cost scales in no more than
total number of S 1 operations.
This base case allows us to analyze the performance of algorithm 3:
Proof of theorem D.2. We demonstrate the bounds inductively. The base case (ℓ = 1) holds via lemma D.1. For the inductive hypothesis, we assume that, for any p ′ ⩾ 1 and k = 2 ℓ , we may implement POWER(k, t, p ′ ) as
To demonstrate the inductive step, we seek to apply corollary B.2 directly to the implementable operators from the inductive hypothesis. Thus, we set p ′ = 2p so that POWER (2k, t, p)
our desired error scaling. By corollary B.2, we require an adder of order max(⌈ p ′ -1 2 ⌉, 1) ⩽ p + 1 2 . Thus, the adder requires 1.07
Since COST(2, kp) ⩽ 6 • 14 kp by lemma D.1, the number of S 1 operations is upper bounded by COST (2k, p) ⩽ 6 log 2 k+1 • 420 kp =⇒ COST (k, p) ⩽ 6 log 2 k • 420 kp/2 . (136) underlying operators. By assuming symmetry of the Trotter formula for ADD, each addition requires 1 2 1.07 • 30 q of the underlying operator. Thus, we can obtain a bound on the number of applications required of each fundamental POWER operator as log 2 n s=1 1 2 1.07 • 30 2 s-1 p+1/2 ⩽ 3 log 2 n 30 np ⩽ n 1.6 30 np (139)
because the sth layer of ADD requires constituent operators of order 2 s p, so q ⩽ 2 s-1 p + 1 2 . Finally, considering the total cost by adding up the cost of the individual POWER operators multiplied by the number of applications required yields
as desired. Now, we seek to 昀椀nalize the number of operations required in terms of ϵ. Recognize that we may use timeslicing to reduce the error arbitrarily. Note that: Lemma D.2. Suppose we may implement T k,p (t), an approximation of T k (t) with p > 1 such that
Then, by timeslicing the approximation, we can produce T r k,p (t) where
where T r k (t) requires r ∈ Θ (Λ k/2 t) 1+1/(p-1)
applications of the T k (t) operator.
Proof. We de昀椀ne the timeslicing of T k (t) as applying T k (t/r) operator r times
To 昀椀nd a Taylor expansion for T k,p (t/r) r , note the explicit form for T k,p (t)
where ∥∆(t)∥ ∈ O(1). This allows us to express the re昀椀ned operator as
Thus, when we analyze the error,
To bound the implementation error by ϵ, i.e. ϵ ∈ O (Λ k/2 t) p r p-1
, we should select r as
as desired.
Finally, we demonstrate our original theorem statement, which allows us to create a bound on the number of S 1 operations necessary to achieve an arbitrarily ordered operator:
Proof of theorem D.1. By theorem D.3, we can perform a single Trotter step of timestep t r using n 1.6 30 np 420 n 2 p/2 6 log 2 n+1 S 1 operations. Thus, the total number of S 1 operations required scales in r • n 1.6 30 np 420 n 2 p/2 6 log 2 n+1 ,
where, by lemma D.2, it is suf昀椀cient to set r ∈ Θ (Λ k/2 t) 1+1/(p-1)
.
Hamiltonian:
(i.e. a Hamiltonian with a Kerr non-linearity). Let t be the evolution time and ϵ be the target error tolerance. For any positive integer q we can approximate an exponential of the blockencoded Hamiltonian with error at most ϵ in the operator norm using re O(q) S 1 operations where r ∈ Ω (Λ 4 t) 1+1/(q-3 4 )
ϵ 1/(q-3 4 )
.
Proof. We 昀椀rst show that the two Hamiltonian terms are implementable separately. Then, via Trotter, we combine them and perform an error analysis. In particular, we hope to embed the Hamiltonian such that we approximate the operator
i.e. where the Hamiltonian is embedded in the upper left hand block. Thus, when applied to a system with the |0⟩ qubit, this amounts to implementing exp i tH on the mode.
We begin by embedding the a † a term. Notice that a † a Hermitian; thus, we can apply algorithm 2 on the a, a † block encodings. By the error analysis in theorem B.1 and the bound on the norm from fact B.1,
By setting τ = ωt 2 , we can achieve the a † a term:
requiring 8 • 6 q-1 total exp i tB a † operations (which can be reduced to S 1 operations via conjugation).
We now tackle the second order term. Recall that we can block-encode (a † ) 2 and (a) 2 via algorithm 3. We then apply algorithm 2 to (a † ) 2 , a 2 to yield the desired upper-left block encoding. Namely, we can implement exp i tB (a † ) 2 with the following error scaling:
using no more than 6 • 14 p S 1 operations, we can apply algorithm 2 to 昀椀nd
by setting ℓ = ⌈ p-1 2 2 ⌉ ⩽ p 2 + 1 and thus using 8 • 6 ℓ-1 exponentials. When τ = κt 4 , MULT i κt 4 B (a † ) 2 , i κt 4 B (a) 2 -exp i t κ 2 a † 2 (a) 2 0 0 (a) 2 a † 2 ∈ O Λ 4 κt p 2 + 1 4 , (157) using no more than 8 • 6 ℓ-1 • 6 • 14 p ⩽ 48 • 6 p/2 • 14 p ⩽ 48 • 35 p total S 1 operators. Then, we may set p = 2q so that, given no more than 48 • 35 2q total S 1 operations, we can implement the BCH formula with error scaling O((Λ 4 κt) q+ 1 4
). Finally, we apply the Trotter formula to the two subterms via lemma B.2. De昀椀ne the approximate matrix exponentials as follows:
so that,
where the constant factor can be obtained by observing the Hamiltonian norm is bounded via the triangle inequality. Now, by setting s = ⌈ 1 2 (q -3 4 )⌉ ⩽ q 2 + 5 8 , we can obtain the desired error scaling. This formula requires no more than 4 • 5 s-1 total operations; by symmetry, we can assume each of the BCH formulas only must be applied 2 • 5 s-1 times. Thus, the total number of S 1 operations is no more than
To produce an ϵ scaling, we apply lemma D.2 to the Trotterized operator, implying that we require the following r scaling for 昀椀xed q:
where the total number of S 1 operations is no more than
We 昀椀rst need to demonstrate the connection between block-encoded powers of annihilation/creation operators and state preparation. First, observe that the ideal block encoding would allow for initialization from the vacuum:
where the matrix products have well-de昀椀ned forms
When t √ k! = (2n + 1) π 2 for n ∈ N, the |1⟩ ⊗ |0⟩ term vanishes and we are left with the |0⟩ ⊗ |k⟩ Fock state, as desired.
While this result allows us to prepare the |k⟩ Fock state, it also will incur unwanted transformations on starting states other than the vacuum (|1⟩ ⊗ |n⟩, n ̸ = 1). By applying the BCH formula, we can isolate this operation so that it only operates on the |1⟩ ⊗ |0⟩ term. In particular, we argue: Theorem D.6. Consider the Fock preparation unitary P k with the form
When t = (2n + 1) π
√ k! , we have that P k performs our desired state preparation
We claim that we can approximate this unitary with P k,p where
using no more than 4 • 5 q-1 T k,p subroutines.
Proof. The general construction of the operator emerges from the use of a Trotter formula in conjunction with a phase rotation gate. We begin by de昀椀ning the rotation operator:
Definition D.7. Call R Z0 the phase-昀氀ip operator acting on some set of modes B to be
i.e. only 昀氀ip the phase for the vacuum. This operator is implementable using a 0-controlled cavity-conditioned qubit rotation gate.
Then, because R Z0 is self-adjoint, we can conjugate T k (-t) without error as
where speci昀椀c left-and right-hand (2|0⟩⟨0| -1) terms vanish given the annihilation/creation operators. We then apply the Trotter formula upon
as desired.
To compute the error scaling, recall our result from theorem D.3, which states the error scaling of T k,p (and, respectively,
so that, by lemma B.2,
when q = max(⌈ p-1 2 ⌉, 1) and using no more than 4 • 5 q-1 operator exponentials. Thus, we can set q = p+1 2 ⩾ max(⌈ p-1 2 ⌉, 1) so that we use no more than 4
Thus, our approximate operators can be applied to yield the same result with high probability:
Theorem D.8. We can prepare the |0⟩ ⊗ |k⟩ with probability at least 1 - ¶ using no more than r ∈ Θ (Λ k/2 t) 1+1/(p-1)
(179) S 1 operators.
Proof. Begin by identifying the ϵ precision necessary to yield a failure probability less than ¶. A suf昀椀cient condition would be that
Observe that our idealized operator has a success probability; thus, we seek to demonstrate that
Because the probability of measuring |0⟩ ⊗ |k⟩ lies in [0, 1], the above inequality holds when
We recognize that we can lower bound the norm
by requiring P k -P k,p ∞ ⩽ 1. This allows us to produce a lower bound on the original LHS
Thus, it is suf昀椀cient to hold
We apply lemma D.2 to theorem D.6 so that the time-sliced P r k,p has
applications of P k,p (t/r). The S 1 bound follows from a similar analysis to theorem D.1 applied to the result from theorem D.6.
To further demonstrate the ef昀椀cacy of the instruction set, we employ the approach to encode a qubit in a cavity either via generation of effective Pauli gates or imposition of an effective Hubbard interaction in the Jaynes-Cumming Hamiltonian. In this sense, the techniques presented here are analogous to those in [48], in that we use our results to effectively truncate the quantum information to a two-dimensional subspace despite the fact that the natural dynamics of the systems causes the quantum information to leak from this space into the larger Hilbert space of the cavity. Error bounds for these complex operations are considered outside of the scope of the present work and are left for future research.
For universal control in the restricted span {|0⟩ , |1⟩} Hilbert space, we generate three effective Pauli operators σ x eff , σ y eff , and σ z eff that produce Pauli rotations in the lowest two modes of the cavity, with minimal leakage to higher energy states. The form of the effective Pauli operators is determined by expressing the standard Pauli operators
in terms of creation and annihilation operators truncated to the 昀椀rst two Fock states
which yields
To reduce leakage into higher energy states, we ensure the creation operator â † eff only acts on the ground state |0⟩ and the annihilation operator âeff only acts on the 昀椀rst excited state |1⟩ with the projector
where n is the number of photons in the cavity and where only the span {|0⟩ , |1⟩} states are populated. Since the operator is a projector, it obeys the relation
such that the effective Pauli gates are
Consider the in昀椀nitesimal σ x -rotation gate in the span {|0⟩ , |1⟩} Fock space
Expression of the exponent in terms of phase-space operators equations (3), ( 4) and (45) gives
The gate is therefore given by Trotter for each of three terms: exp
where Â1 corresponds to a position displacement and B1 corresponds to the y-conditional rotation gate. The argument of the second term is already in the form of a commutator, such that
where Â2 corresponds to an unconditional momentum boost, and B2 corresponds to the zconditional rotation gate. Lastly, the third term already belongs to the ISA and needs no further decomposition.
The in昀椀nitesimal σ x -rotation gate in the span {|0⟩ , |1⟩} Fock space is therefore composed of a product of nine rotation and displacement gates or 21 displacement gates.
The in昀椀nitesimal σ y -rotation gate in the span {|0⟩ , |1⟩} Fock space is determined analogously
where the exponent of Â2 corresponds to a conditional momentum shift and the exponent of B2 is a y-conditional rotation gate. The in昀椀nitesimal σ y -rotation gate in the span {|0⟩ , |1⟩} Fock space therefore has a lower bound gate depth of nine displacement and rotation gates or 21 in displacement gates.
The in昀椀nitesimal σ z -rotation gate in the span {|0⟩ , |1⟩} Fock space is U span{0,1},z = e iλ 2 σ z eff σ z (224) = e -λ 2 (I-2â † (I-n)â)σ z , (
whose argument in terms of ladder operators is
Given the ladder operator commutator equation (6),
the relationship between the fourth-order ladder operator term and the number operator is
The argument of the exponential in terms of number operators is then
To further demonstrate the power of our strategy, we employ the approach to simulate fermionic dynamics on qubit-qumode systems. We detail the required operations, with error estimates an area for future work. We consider the Fermi-Hubbard lattice Hamiltonian
The kinetic energy term TFH describes the nearest-neighbor interaction for hopping of a single spin between two sites with hopping parameter J and spin σ given annihilation operators {ĉ j,σ } and creation operators ĉ † j,σ for sites {j}. The potential energy term VFH describes the samesite interaction, which gives the energetic unfavorability of a spin up ↑ and spin down ↓ coexisting on the same site i, where nj,σ gives the number of spin σ particles on site j. According to fermion statistics, no more than a single particle of a given spin can exist on a single site.
Each individual set of qumodes represents either a spin up or spin down particle on a single lattice site, for direct comparison to the qubit-based schemes of [49][50][51]. Each qmode set is connected to a qumode set that represents the same site of opposite spin to facilitate computation of the potential energy VFH , as well as to qumode sets of the same spin on neighboring sites to facilitate computation of the kinetic energy TFH . Qumode sets are also connected along Jordan-Wigner strings to take into account fermionic statistics.
The |0⟩ qumode represents absence of a spin and the |1⟩ qumode represents presence of a spin. Within each cavity, only the qumodes in span {|0⟩, |1⟩} are considered, as in section (E), which prevents leakage into unphysical high-energy qumodes. At the end of each operation, the qumode must be in either the |0⟩ or |1⟩ mode and the transmon state must also be in the ground state |g⟩, which provides an error syndrome and therefore a degree of error detection not employed in qubit-based representations of the Fermi-Hubbard lattice.
Propagation of any combination of up and down spins is simulated with three gates that operate on two sets of qumodes. The 昀椀rst two gates-the same-site and hopping gates-are de昀椀ned as the propagator of the same-site and hopping Hamiltonians, respectively. The samesite term of the Hamiltonian for site i is (244)
The diagonal propagator is recognized as both the conditional cross-Kerr interaction and the controlled-phase (CPHASE) gate in the reduced subspace span {|0⟩ , |1⟩}. The hopping term of the Hamiltonian for each σ spin in sites i, (i + 1) is
where the latter expression employs the commutator relationship of the annihilation and creation operators. The hopping Hamiltonian for the speci昀椀ed mapping is then the off-diagonal matrix
which gives the hopping propagator U hop = e -iH hop τ U hop = 1 0 0 0 0 cos (tτ ) i sin (tτ ) 0 0 i sin (tτ ) cos (tτ ) 0 0 0 0 1
which is recognized as a conditional controlled-phase beamsplitter restricted to span {|0⟩ , |1⟩}, or equivalently a Givens/iSWAP-like gate in the reduced span {|0⟩ , |1⟩} subspace [50,51]. The 昀椀nal gate of the three-gate set incorporates the fermionic statistics of the spins via the fermionic SWAP (FSWAP) gate [49,51]. The state of each set of qumodes is swapped with one of its neighbors with inclusion of a phase where both spins are present in neighboring sets as
1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 -1
which is recognized as the product of a conditional rotation gate and a beam-splitter on 3D cQED systems. Finally, initial states are prepared by the universal set of gates in span {|0⟩ , |1⟩} detailed in section E.
We consider the in昀椀nitesimal conditional cross-Kerr gate
which is also employed in Gottesman-Kitaev-Preskill codes encoded in qubit-qumode systems [52].
The argument is expressed in terms of a commutator according to the Pauli commutation relation equation (16) [A, B]λ 2 = iλ 2 n1 n2 σ z (251)
⌉, 1) and p l , p r are the orders of the