basis to the basis $\\{|\phi_{\mathbf{s}}\rangle\\}$ required for use of the
non-Abelian state transfer protocol. Specifically, we seek to map from the
computational basis state $|r\rangle$ to a corresponding state
$|\phi_{\mathbf{s}}\rangle$. To specify an appropriate correspondence, recall
that in this context the string $\mathbf{s}\in\\{+1,-1\\}^{n+a}$ is associated
with a position
$q_{\mathbf{s}}=\lambda(\sum_{j=1}^{n+a-1}s_{j}2^{j-1}-s_{n+a}2^{n+a-1})$ (see
Eq. (550)), that takes discrete values in the range
$[-\lambda(2^{n+a}-1),\lambda(2^{n+a}-1)]$ with spacing $2\lambda$. It is
therefore suitable to index each such possible value by an integer $s$,
defined as

$\displaystyle s:$
$\displaystyle=\tfrac{1}{2\lambda}\big{(}q_{\mathbf{s}}+\lambda(2^{n+a}-1)\big{)}\in\\{0,1,...,2^{n+a}-1\\}.$
(578)

Thus, we assume access to a unitary change of basis operation $\mathcal{T}$
that acts as $\mathcal{T}|s\rangle=|\phi_{\mathbf{s}}\rangle$, such that its
application to Eq. (577) generates the state

$\displaystyle\sum_{s=0}^{2^{n+a}-1}C_{s}|\phi_{{s}}\rangle,$ (579)

where we are now indexing by the associated integer $s$ to simplify analysis.
As shown in Ref. [466], $\mathcal{T}$ can be easily implemented as a product
of local operations:

$\mathcal{T}=(HX)^{\otimes(n+a-1)}\cdot H_{n+a}\cdot Z_{1}Z_{n+a}.$ (580)

Step 3: The next step uses the non-Abelian state transfer protocol to transfer
the state to an oscillator as per Eq. (552). To achieve this, we prepare the
oscillator in the sinc state
$|0,2\lambda\rangle_{B}=\tfrac{1}{\sqrt{2\lambda}}\int dq\
\text{sinc}\big{(}\tfrac{\pi q}{2\lambda}\big{)}|q\rangle_{B}$, and apply
$U_{\text{st}}^{(n+a)}(2\lambda)$ to the joint system, producing the state

$\displaystyle|\mathbf{0}\rangle_{Q}\otimes\sum_{s=0}^{2^{n+a}-1}C_{s}\frac{1}{\sqrt{2\lambda}}\int
dq\
\text{sinc}\big{(}\tfrac{\pi(q-q_{\mathbf{s}})}{2\lambda}\big{)}|q\rangle_{B}.$
(581)

The infidelity suffered in this step is quantified by the expression in Eq.
(554).

Step 4: Next, we symmetrize this state about $q=0$ by applying a displacement
operation $D(\lambda)$. Using the identities $s=2^{n}\cdot k+x$ and
$q_{s}=2\lambda(s-\tfrac{1}{2}2^{n+a}-\tfrac{1}{2})=2\lambda(2^{n}\cdot[k-\tfrac{1}{2}2^{a}]+x-\tfrac{1}{2})$,
this yields the wave function

$\displaystyle\sum_{x=0}^{2^{n}-1}\sum_{k=-2^{a}/2}^{2^{a}/2-1}\frac{1}{\sqrt{2^{a}}}c_{\mathbf{x}}\frac{1}{\sqrt{2\lambda}}\text{sinc}\big{(}\tfrac{\pi(q-2\lambda(2^{n}\cdot
k+x))}{2\lambda}\big{)}.$ (582)

Step 5: Subsequently, we apply the Fourier gate to the oscillator, which
enacts a continuous Fourier transform on the wave function. Noting the Fourier
transform relation $\frac{1}{\sqrt{2\pi}}\int
dq^{\prime}e^{iqq^{\prime}}\frac{1}{\sqrt{2\lambda}}\text{sinc}\big{(}\tfrac{\pi(q^{\prime}-b)}{2\lambda}\big{)}=\sqrt{\frac{\lambda}{\pi}}1_{|q|\leq\pi/2\lambda}e^{iqb}$,
where $1_{|q|\leq\pi/2\lambda}$ is an indicator function (i.e. a box function
in this context), this transforms the above wave function to

$\displaystyle\sqrt{\frac{\lambda}{\pi}}1_{|q|\leq\pi/2\lambda}\sum_{x=0}^{2^{n}-1}c_{j}e^{iq2\lambda
x}\left[\frac{1}{\sqrt{2^{a}}}\sum_{k=-2^{a}/2}^{2^{a}/2-1}e^{iq2\lambda
2^{n}k}\right]$ (583) $\displaystyle\qquad=:\tilde{\psi}(q).$

Step 6: After free evolution, we perform the inverse of Steps 4 and 5, but now
with reciprocal spacing
$2\lambda^{\prime}=\tfrac{2\pi}{2^{n+a}}\tfrac{1}{2\lambda}$ chosen to produce
the QFT. We first enact the inverse of step 5, which consists of applying the
displacement $D(-\lambda^{\prime})$ to map the wave function to
$\tilde{\psi}(q+\lambda^{\prime})$.

Step 7: Next, we perform the inverse of Step 4, by transferring the oscillator
state to $n+a$ qubits by using
$U_{\text{st}}^{(n+a)}(2\lambda^{\prime})^{\dagger}$. As per Eq. (VI.1.2),
this maps the above state to

$\displaystyle\sum_{s=0}^{2^{n+a}-1}\int dq\
\tilde{\psi}(q+q_{s}^{\prime}+\lambda)\text{sinc}(\tfrac{\pi
q}{2\lambda})|\phi_{\mathbf{s}}\rangle_{Q}|q\rangle_{B}=$ (584)
$\displaystyle\sum_{s=0}^{2^{n+a}-1}\sum_{x=0}^{2^{n}-1}\sqrt{\frac{\lambda}{\pi}}\int_{-\pi/2\lambda}^{\pi/2\lambda}dq\
c_{\mathbf{x}}e^{i(q+q_{s}^{\prime}+\lambda)2\lambda x}$
$\displaystyle\qquad\quad\times\Big{[}\frac{1}{\sqrt{2^{a}}}\sum_{k}e^{i(q+q_{s}^{\prime}+\lambda)2\lambda
2^{n}k}\Big{]}\text{sinc}(\tfrac{\pi
q}{2\lambda^{\prime}})|\phi_{\mathbf{s}}\rangle_{Q}|q\rangle_{B},$

where $q_{s}^{\prime}=2\lambda^{\prime}s-\lambda^{\prime}(2^{n+a}-1)$. Note
that we have not used any approximations to simplify this expression, as was
done in Sec. VI.1.2.

As in the previous QFT protocol, the factor
$\sum_{k=-2^{a}/2}^{2^{a}/2-1}e^{i(q+q_{s}^{\prime}+\lambda)\Delta 2^{n}k}$ in
brackets is the Dirichlet kernel, which in the limit of large $a$ behaves as
the sum of delta functions
$\sum_{l\in\mathbb{Z}}\delta(\frac{(q+q^{\prime}_{s}+\lambda^{\prime})\Delta
2^{n}}{2\pi}-l)$, with an overall error $\mathcal{O}(1/2^{a})$ at finite $a$.
One can take this limit to drastically simplify this state; the full technical
analysis of this step is presented in Ref. [466]. The final result is that the
state of the system can be simplified to the product state

$\displaystyle\Bigg{[}\sum_{l=0}^{2^{n}-1}\Big{(}\frac{1}{\sqrt{2^{n}}}\sum_{x=0}^{2^{n}-1}c_{j}e^{i2\pi
lx/2^{n}}e^{-i\pi
x/2}\Big{)}|\phi_{s=2^{a}l}\rangle_{Q}\Bigg{]}|\tilde{q}=0\rangle_{B}.$ (585)

where $|\tilde{q}=0\rangle_{B}$ is a state concentrated around $q=0$, and the
infidelity suffered in reaching this expression is $\mathcal{O}(1/2^{a})$. The
oscillator and qubits are now decoupled, and the qubits state is nearly the
quantum Fourier transform evaluated in the basis $|\phi_{s=2^{a}l}\rangle$.
Again, the presence of the phase $e^{i2\pi lx/2^{n}}$ resembling the QFT
arises from the continuous-discrete correspondence of the Fourier transform
discussed earlier.

Step 8: To extract the quantum Fourier transform from this state, the last
step is to transform from the $\\{|\phi_{s}\rangle\\}$ basis back to the
computational basis using $\mathcal{T}^{\dagger}$. This acts as
$\mathcal{T}^{\dagger}|\phi_{s=2^{a}l}\rangle=|2^{a}l\rangle=|l\rangle|0\rangle^{\otimes
a}$, such that its application to the previous state produces

$\displaystyle\left[\sum_{l}\Big{(}\frac{1}{\sqrt{2^{n}}}\sum_{x=0}^{2^{n}-1}c_{\mathbf{x}}e^{i2\pi
lj/2^{n}}e^{-i\pi x/2}\Big{)}|l\rangle|0\rangle^{\otimes
a}\right]_{Q}\otimes|\tilde{q}=0\rangle_{B}$ (586)
$\displaystyle~{}~{}~{}=\left[\sum_{x=0}^{2^{n}-1}e^{-i\pi
x/2}c_{\mathbf{x}}U_{\text{QFT}}\big{(}|x\rangle\big{)}|0\rangle^{\otimes
a}\right]_{Q}\otimes|\tilde{q}=0\rangle_{B}$

The remaining state on the first $n$ qubits is nearly the QFT of the original
$n$-qubit state $|\psi\rangle_{Q}$, up to a phase factor $e^{-i\pi x/2}$, just
as in the previous QFT protocol.

Again, this phase factor arises because this protocol really implements a
symmetric QFT, which can be mapped to the usual QFT at the expense of
incurring a factor $e^{-i\pi x/2}$. As in the QFT previous protocol, this
phase factor is inconsequential and can be removed by appending $|00\rangle$
to the initial state, such that $e^{-i\pi x/2}=1$ for all $x$ with nonzero
coefficients $c_{\mathbf{x}}$. This protocol then outputs the state
$U_{\text{QFT}}(|\psi\rangle|00\rangle)=|+\rangle|+\rangle
U_{\text{QFT}}(|\psi\rangle)$, from which the QFT can be extracted.

Performance: Consider first the gate complexity of this QFT protocol, the bulk
of which is due to the use of the non-Abelian state transfer protocol.
Therefore, as per the performance bounds in Sec. VI.1.2, the gate count scales
as $\mathcal{O}(n+a)$. However, these gates require a total displacement
$\mathcal{O}(\lambda 2^{n+a})$, which as before, implies a time complexity of
$\mathcal{O}(\lambda 2^{n+a})$ when the displacements are implemented with a
fixed coupling between the qubits and oscillator, but this can be reduced with
a sufficiently tunable and strong coupling.

Next, consider the fidelity of this protocol. The first source of infidelity
occurs in Step 3 due to the usage of the non-Abelian QSP state transfer
protocol. The infidelity of this operation is quantified in Eq. (554); a
careful analysis of this scenario, presented in Ref. [466], indicates that the
resulting infidelity scales as $\mathcal{O}(1/2^{a})$. Moreover, the other
source of infidelity arises in Step 6 when taking the large $a$ limit, which
also contributes an infidelity $\mathcal{O}(1/2^{a})$. Lastly, we note that no
infidelity is suffered in the inverse state transfer step, as we do not use
any approximation to simplify this step. Therefore, summing these
contributions to the infidelity, we find the overall fidelity of this QFT
protocol is

$F=1-\mathcal{O}(1/2^{a}).$ (587)

###  VI.3 Quantum Random Walks

Quantum interference can cause non-trivial effects which are hard to observe
using deterministic processes. Quantum random walks (QRWs) use multiple rounds
of measurements (whose results are non-deterministic) to observe such
phenomena and bring out some of the novel features of quantum mechanics.
Introduced in Ref. [471], QRWs describe the trajectory of a particle with
internal degrees of freedom moving on a discrete lattice. One possible
realization is to use an oscillator in a coherent state as the particle and
represent its internal degrees of freedom using an auxiliary qubit. In this
setting, QRWs involve repeated conditional displacements (CD$(\epsilon)$, see
Box III.7) of the oscillator, each followed by a multi-dimensional coin toss
which decides the basis in which the auxiliary qubit is measured. In other
words, the multi-dimensional coin toss decides the auxiliary qubit rotations
about the vectors on the $XY$ plane of the Bloch sphere ($R_{\phi}(\theta)$,
see Table A.1) before measuring in the $Z$ basis. Thus, the phase-space ISA
can be exploited to straightforwardly implement QRWs in hybrid CV-DV quantum
computing architectures.

Ref. [471] showed that QRWs could be used to study non-trivial effects in
quantum mechanics that arise due to quantum interference. For example, let us
conditionally displace an oscillator by a very small amount $|\epsilon|\ll 1$.
Using the definitions from Eq. (483) (but now including the normalization
factors for the cat states), we can write,

$\begin{split}\mathrm{CD}(\epsilon)\ket{0}\ket{+}&=\left(\rule{0.0pt}{10.33327pt}\ket{\epsilon}\ket{0}\
+\ \ket{-\epsilon}\ket{1}\right)/\sqrt{2}\\\
&=\frac{\mathcal{N}_{\text{even}}}{2}\ket{\mathcal{C}_{0,\epsilon,+}}\ket{+}\
+\
\frac{\mathcal{N}_{\text{odd}}}{2}\ket{\mathcal{C}_{0,\epsilon,-}}\ket{-}\\\
&\approx\frac{\mathcal{N}_{\text{even}}}{2}\ket{0}\ket{+}\ +\
\frac{\mathcal{N}_{\text{odd}}}{2}\ket{1}\ket{-}~{}.\end{split}$ (588)

The final equation uses an approximation in which we have equated a small even
(odd) cat state with Fock state $\ket{0}$ ($\ket{1}$). It is important to note
that the definitions used here for cat states
$\ket{\mathcal{C}_{0,\epsilon,\pm}}$ include their normalization prefactors,
giving rise to the unequal amplitudes for each term in the last two equations.
More specifically, the probability of each Fock state depends on the
normalization constants of even ($\mathcal{N}_{\text{even}}$) and odd
($\mathcal{N}_{\text{even}}$) cat states of size $\epsilon$,

$\displaystyle\mathcal{N}_{\text{odd}}$
$\displaystyle=\sqrt{2(1-e^{-2\epsilon^{2}})}$ (589)
$\displaystyle\mathcal{N}_{\text{even}}$
$\displaystyle=\sqrt{2(1+e^{-2\epsilon^{2}})}.$ (590)

The fidelity of even and odd cat states with the Fock states $\ket{0}$ and
$\ket{1}$ states improves with decreasing value of $\epsilon$ as long as
$\epsilon\neq 0$, as shown by Eq. (591) and (592).

$\displaystyle F_{0}$ $\displaystyle=|\langle
0|\mathcal{C}_{0,\epsilon,+}\rangle|^{2}=4\frac{e^{-\epsilon^{2}}}{\mathcal{N}^{2}_{\text{even}}},\quad$
$\displaystyle\alpha\in\mathbb{R}~{},$ (591) $\displaystyle F_{1}$
$\displaystyle=|\langle
1|\mathcal{C}_{0,\epsilon,-}\rangle|^{2}=4\epsilon^{2}\frac{e^{-\epsilon^{2}}}{\mathcal{N}_{\text{odd}}^{2}},\quad$
$\displaystyle\alpha\in\mathbb{R}~{}.$ (592)

Note that, for $\epsilon\rightarrow 0$, the small odd cat has a high overlap
with the Fock $\ket{1}$ state (see Eq. (592)) in contrast to the small even
cat which has a large overlap with $\ket{0}$ (see Eq. (591)) largely
resembling the initial vacuum state. Thus, Eq. (588) shows that there can be a
non-trivial unit increase in the photon number of the oscillator upon qubit
measurement after an extremely small conditional displacement. We show the
marginal distribution in position $|\psi(x)|^{2}$ and the phase-space quasi-
probability distribution $\mathcal{W}(x,p)$ in Fig. 42 for pictorial
intuition. A qubit measurement in the $X$-basis collapses the oscillator onto
one of the two states shown in this figure.

Let us discuss the consequences of this on a position basis. The even cat
(left) resembles a vacuum and has its peak centered at $x=0$ while the odd cat
(right) has two peaks centered at $x=\pm 1$. If the qubit measurement
collapses the oscillator onto the odd cat state, the position measurement
would yield $|x|=1$ with maximum probability, hence yielding a net
displacement of $|\delta x_{-}|\sim 1$. Here the subscript
${}^{\prime}-^{\prime}$ denotes displacement upon position measurement of the
oscillator state conditioned on the qubit measurement yielding the $\ket{-}$
state.

Thus, it is evident that upon measuring the position of the oscillator there
is a non-zero probability that this measurement yields an effective
displacement of $\delta x_{-}\sim\mathcal{O}(1)$ in position, much larger than
the original displacement of $\epsilon$. This is a quantum mechanical effect
that is only possible due to a local phase of $\pi$ incurring a minus sign
between the two slightly displaced coherent states $\ket{\pm\epsilon}$. The
effect of this minus sign is to subtract the wave functions corresponding to
the two displaced coherent states. While highly non-trivial, the collapse of
the oscillator to this ${}^{\prime}-^{\prime}$ superposition is also an
extremely rare event because the ratio,

$\displaystyle\mathcal{N}_{\text{odd}}/\mathcal{N}_{\text{even}}$
$\displaystyle=\sqrt{\frac{1-e^{-2\epsilon^{2}}}{1+e^{-2\epsilon^{2}}}}=\epsilon+O(\epsilon^{5}).$
(593)

vanishes as $\epsilon\to 0$.

Figure 42:  Marginal probabilities ($|\psi(x)|^{2}$) and Wigner distribution
($\mathcal{W}(x,p)$) of the oscillator states entangled with orthogonal qubit
states $\ket{\pm}$ after small conditional displacements
$\mathrm{CD}(\epsilon=0.071)$ as shown in Eq. (588). (Left) Oscillator state
entangled with $\ket{+}$ qubit state. This state closely resembles the vacuum
state $\ket{0}$ and occurs for qubit measurement $M_{X}=+1$ with a high
probability $p=0.995$. (Right) Oscillator state entangled with $\ket{-}$ qubit
state. This state is well-approximated by Fock state $\ket{1}$, and occurs for
qubit measurement $M_{X}=-1$ with a low probability $p=0.005$. The large
change in the position of the peaks, shown in the marginal probability
distribution, between the two states is achieved merely by a small
displacement of $\epsilon-0.071$. This phenomenon is an artifact of quantum
interference and can be observed using quantum random walks as described in
Ref. [471].

The probability of observing a large state change (displacement in phase
space) induced by the measurement back-action can be improved by changing the
qubit measurement axis, a task addressed by a quantum random walk. To see this
let us rewrite the state in Eq. (588) in terms of an arbitrary auxiliary qubit
basis $\vec{\sigma}_{\theta}=\cos{\theta}\sigma_{z}+\sin{\theta}\sigma_{x}$,
the $\pm 1$ eigenstates of which are represented by $\ket{\pm\theta}$.
Alternatively, we can rotate the qubit by
$R_{Y}\big{(}{\frac{\theta}{2}}\big{)}$ and then measure it in $Z$-basis. In
this case, the composite oscillator-qubit state just before measurement is
equal to (with implicit normalization),

$\displaystyle\ket{\psi}$
$\displaystyle=(\ket{\epsilon}+\tan{(\theta/2)}\ket{-\epsilon})\ket{+\theta}$
$\displaystyle~{}~{}~{}~{}+(\tan{(\theta/2)}\ket{\epsilon}-\ket{-\epsilon})\ket{-\theta}.$
(594)

Ref. [471] shows that under the assumption $\epsilon\rightarrow 0$, we can use
$e^{\pm i\epsilon\hat{p}}\approx I\pm i\epsilon\hat{p}$ and approximate the
final state in Eq. (594) effectively as the state after an asymmetric
conditional displacement $D_{\theta}(\delta x_{+},\delta x_{-})$ of vacuum.

$\displaystyle\ket{\psi}$ $\displaystyle\approx D_{\theta}(\delta x_{+},\delta
x_{-})\ket{0}_{\text{vac}}$ (595)
$\displaystyle=\sqrt{\mathcal{P_{+}}}\ket{\delta
x_{+}}\ket{+\theta}+\sqrt{\mathcal{P_{-}}}\ket{\delta x_{-}}\ket{-\theta},$
(596) where $\displaystyle\delta
x_{\pm}=\frac{1\mp\tan{(\theta/2)}}{1\pm\tan{(\theta/2)}}\epsilon,$ (597)
$\displaystyle D_{\theta}(\delta x_{+},\delta x_{-})$
$\displaystyle=e^{-i\delta x_{+}\hat{p}}\ket{+\theta}\bra{+\theta}+e^{-i\delta
x_{-}\hat{p}}\ket{-\theta}\bra{-\theta},$ (598)

and $\mathcal{P}_{\pm}$ denotes the probabilities of collapsing the oscillator
onto coherent states $\ket{\delta x_{\pm}}$. This approximation holds as long
as $|\delta x_{\pm}|\ll\Delta_{x}$, where $\Delta_{x}$ is the position
uncertainty of the initial state (here assumed to be $\Delta_{x}=\frac{1}{2}$,
in correspondence with an initial vacuum state).

Despite this restriction, it is possible to observe either $|\delta
x_{+}|\gg\epsilon$ or $|\delta x_{-}|\gg\epsilon$ with significant probability
by applying the modified qubit rotation $R_{Y}(\theta)$ in place of
$R_{Y}(\pi/2)$ prior to measurement. For example, choosing
$\tan{(\theta/2)}=1+\gamma$ such that $\epsilon/\Delta_{x}\ll|\gamma|\ll 1$
yields $\delta x_{-}=-2\epsilon/\gamma$ which satisfies the requirement
$\epsilon\ll|\delta x_{-}|\ll\Delta_{x}$. That is, the effective displacement
of the oscillator $|\delta x_{-}|$, upon collapse to $\ket{-\theta}$ qubit
state, is much larger than $\epsilon$ while being much smaller than the
initial position uncertainty $\Delta_{x}$. The probability of this event is
$P_{-}=\gamma^{2}/4$. Note that for a pre-measurement rotation of
$R_{Y}(\pi/2)$ such that $\gamma=0$, the probability of collapsing to
$\ket{\delta x_{-}}$ under these assumptions vanishes. In contrast, by
choosing a non-zero $\gamma$, and thus rotating qubit by
$R_{Y}(2\tan^{-1}{(1+\gamma)})$, the probability of observing a non-trivial
displacement of the oscillator will be enhanced. Ref. [471] shows that for
$\epsilon=0.01$ and $\gamma=-0.1$, this effect is apparent after repeating the
protocol for 10 rounds with post-selection upon $-1$ measurement outcomes, and
the qubit reinitialized to $\ket{+}$ after each measurement. The authors also
show another example of an interesting quantum mechanical observation using
photon numbers using the side-band ISA. The proposed random walk shows that
the photon number of the oscillator can be changed by a large amount even if
the interaction only allows for the addition or removal of a single photon.

Another interesting scenario is to carry out the above procedure but skip the
‘measurement + post-selection’ step. In this case, the iterative procedure
corresponds to a random walk for the composite system with $N$ pairs of
conditional displacements and random qubit rotations in the equatorial plane.
Such a sequence can be described by a unitary channel that produces a
superposition of $2^{N}$ different state trajectories. To this end, we note a
close similarity between QRWs and the phase-space trajectory illustrated in
Fig. 17(b). Ref. [472] experimentally demonstrated such quantum walks via
unitary channels composed entirely out of the phase-space ISA and used them to
explore a novel topological effect which they argued represents a potential
technique for the study of topological quantum materials.

This application of the phase-space ISA (see Sec. III.5.3) paves the way for
observations of novel quantum mechanical phenomena using hybrid CV-DV
architectures.

###  VI.4 Quantum Simulation

This section offers a brief overview of several hybrid CV-DV quantum
simulation methods which are introduced in a companion paper in Ref. [37], a
thorough manual on how to simulate fermions, bosons, and gauge-fields in CV-DV
hardware.

Quantum simulation aims to solve classically intractable physical problems
with controllable quantum devices and has been identified as a likely arena
for the first demonstration of practical advantage over existing classical
methods [473, 474, 475, 476, 477]. Exact classical simulation of quantum
systems is often intractable due to the exponential growth of Hilbert space
with system size. As a result, classical methods often resort to
approximations that involve either truncation or finite sampling methods to
remove the exponential computational overhead at the price of introducing
approximation errors that may or may not be acceptable for a given
application. Such approximations include the use of mean-field theory [478,
479] and semiclassical approximations of quantum phase space [480] which do
not capture particle-particle correlations. Other representative
approximations that truncate the Hilbert space such as many-body perturbation
theory, coupled cluster theory [481, 482], and tensor network-based methods
[483] face difficulties with strongly correlated systems as they are capable
of only limited treatment of many-body correlations. Separately, sampling-
based methods such as quantum Monte Carlo suffer from sign or phase problems
for frustrated spin systems and interacting fermionic matter [484, 485, 486].

As originally noted by Richard Feynman [487], quantum hardware, unlike
classical hardware, is naturally suited to the simulation of quantum
physics151515It is important to distinguish between two cases: (a) Simulation
of the time evolution of a quantum system with a specified Hamiltonian and a
specified starting state; and (b) Finding the ground state (or other
eigenstate) of a specified Hamiltonian. For most physical systems, the former
is in the complexity class and the latter is in which is a quantum analog of
and is believed to be exponentially more challenging than [1]. To that end,
there are two popular approaches: analog quantum simulation and digital
quantum simulation. The former corresponds to the _native_ realization of a
Hamiltonian in a controlled quantum system, while the latter corresponds to
the _synthesis_ of a Hamiltonian using gate-based instructions on a universal
quantum computer (or if not universal, at least a programmable quantum
simulator).

While both approaches have proven promising in recent years [488, 489], most
experiments have thus far focused on “single-species” quantum systems. For
example, demonstrations of digital qubit-based quantum simulation have
typically concerned either spin [490, 491] or fermionic [492] systems (though
we note the recent proposal for hybrid fermion-qudit processors in Ref.
[493]). Similarly, analog quantum simulators have successfully implemented
interacting spin [494, 495], all-bosonic [496, 57], and all-fermionic [497]
systems [477], yet only a few experiments have demonstrated the capability to
simulate both bosons and fermions within a single system [498, 499]. However,
multi-species quantum Hamiltonians are among the most important, and are
crucial to the development of fundamental science across a range of energy
scales, from the interactions between force-carrying bosons and fermionic
elementary particles underlying the standard model to the nuclear vibrational-
electron interactions underpinning quantum chemical phenomena such as
photosynthesis [500]. Consequently, this motivates the need for flexible
approaches to quantum simulation that enable hardware-efficient simulation of
different species within a single platform.

We focus on the problem of digital quantum simulation here as it is most
pertinent to the scope of this work, but will make a comparison to analog
quantum simulation where applicable. Turning first to conventional qubit-based
hardware, we note that it is possible to encode fermionic [78], bosonic [84,
163], and spin-$S$ [381] degrees of freedom using many qubits, in principle
enabling multi-species digital quantum simulation. In practice, however,
bosonic degrees of freedom and the corresponding bosonic operations are costly
to implement in all-qubit hardware, as discussed previously in Sec. III.6. For
example, the implementation of a single displacement gate on a qubit-encoded-
oscillator containing 8 bosons requires on the order of $10^{4}$ two-qubit
gates using the Fock-binary encoding [84] (see App. F). This inefficiency is
particularly problematic in cases where it becomes important to represent high
mean boson numbers with large density fluctuations, regardless of the choice
of encoding.

Hybrid oscillator-qubit platforms can simulate bosonic degrees of freedom
natively while retaining the capabilities of qubit-based hardware. Such hybrid
platforms show immense promise for the study of both single- and multi-species
models containing spin, fermionic, and bosonic degrees of freedom.
Furthermore, the universal ISAs presented in this work enable the
implementation of arbitrary Hamiltonian terms needed to simulate strongly
correlated matter and facilitate measurement of both local and nonlocal
observables, the latter playing a crucial role in the identification of exotic
phases of matter [501, 502].

In the remainder of this section, we extrapolate on the opportunities provided
by digital hybrid oscillator-qubit hardware toward problems in quantum
simulation. We begin by first describing high-level techniques that enable the
study of dynamics and ground states and allow for the measurement of
observables. Following this, we provide several explicit examples of
physically relevant models that include bosonic modes, spins, or some
combination of the two, and demonstrate how one can synthesize the Hamiltonian
terms of interest for each case using the techniques outlined in Section V.
Examples include the Rabi, Bose-Hubbard, and Heisenberg models, as well as a
$\mathcal{Z}_{2}$ lattice gauge theory model. Many of the specific compilation
ingredients discussed throughout are summarized in Table VI.1, which lists a
selection of common bosonic Hamiltonian terms and strategies for realizing the
unitaries they generate with native instructions. We then conclude with some
comments on the potential for simulation of fermion-boson problems in hybrid
oscillator-qubit platforms. For an in-depth exploration of synthesising
fermionic, bosonic, and gauge-field Hamiltonians, see Ref. [37], which in part
is based on the synthesisation methods in Ref. [243]; and Section V.

Term $H_{i}$ | Description | Synthesis of $e^{-iH_{i}t}$ | Example Models
---|---|---|---
$a_{i}^{\dagger}a_{i}$ |  Mode energy, Chemical potential | $R(\theta)$ |  Jaynes-Cummings [503], Bose-Hubbard [504]
$\hat{n}_{i}\left(\hat{n}_{i}-1\right)$ |  On-site interactions | $\textrm{SNAP}(\vec{\varphi})$ | Bose-Hubbard [504]
$e^{i\phi}\hat{a}_{i}^{\dagger}\hat{a}_{j}+\mathrm{h.c.}$ | Hopping | $\textrm{BS}(\theta,\varphi)$ |  Bose-Hubbard [504], Su–Schrieffer–Heeger [505]
$Z_{i}(a+a^{\dagger})$ | Spin-dependent displacement | $\textrm{CD}(\alpha)$ |  Rabi [506], Dicke [507], Hubbard-Holstein [508]
$Z_{i}Z_{j}(a+a^{\dagger})$ | Spin-spin-dependent displacement |  $\textrm{CD}(\alpha)$, SWAP, CP (Sec. V.3.1) | Photosynthetic complex [509]
$Z_{i}(e^{i\varphi}a^{\dagger}_{i}a_{j}+\mathrm{h.c.})$ | Spin-dependent hopping | $\textrm{CBS}(\theta,\varphi)$ | Bosonic $\mathbb{Z}_{2}$ LGT [510]
$a_{i}^{\dagger}a_{j}a_{k}$, $a^{\dagger}_{i}a^{\dagger}_{j}a_{k}a_{l}$, etc. | Higher-order terms |  $R(\theta)$, $D(\alpha)$, $\textrm{TMS}(r,\varphi)$, $\textrm{BS}(\theta,\varphi)$ \+ Approx. |  Bosonic SYK [511, 512]
$S^{z}_{i}S^{z}_{j}$, ($\vec{S}_{i}\cdot\vec{S}_{j})^{k}$, etc. | (Large)-spin interactions |  $R(\theta)$, $\textrm{BS}(\theta,\varphi)$ \+ Approx. |  AKLT [513]
$f_{i}^{\dagger}f_{j}$, $f_{i}^{\dagger}f_{j}^{\dagger}f_{k}f_{l}$, etc. | Fermionic matter |  Multi-qubit gates (Sec. V.3.2 ) |  Electronic structure [514], Hubbard-Holstein [508] Fermi-Hubbard [515]
Table VI.1: Common native and synthesized Hamiltonian terms in condensed
matter physics and quantum chemistry. ‘Higher-order’ terms refer to those
composed of $>2$ creation or annihilation operators multiplied together (i.e.,
polynomials of order 3 or higher in terms of $\hat{x}$ and $\hat{p}$).
$f_{i}^{\dagger},f_{i}$ are fermionic creation and annihilation operators,
realized via fermion-to-qubit mappings such as Jordan-Wigner or Bravyi-Kitaev
[516]. We note that there will, in general, be multiple avenues to synthesize
$e^{-iH_{i}t}$, and the aim here is to give a rough, non-exhaustive guide to
the native gates and techniques needed. Here, ‘Approx.’ refers to the use of
approximate compilation techniques, such as BCH and Trotter-Suzuki formulas
(see Sec. V.4). See Table LABEL:tab:gates-qubit-osc for definitions of common
hybrid gates listed in the Hamiltonian synthesis column. See Ref. [37] for a
thorough manual on how to simulate fermions, bosons, and gauge fields in CV-DV
hardware.

####  VI.4.1 Dynamics

Time dynamics involves the study of a system evolving under the time-evolution
operator $U(t)=e^{-iHt}$, where $H$ is some model Hamiltonian of interest. For
digital quantum simulation, the goal is therefore to synthesize and enact
$U(t)$ through decomposition into a discrete gate set. Many strategies exist
to achieve this, including Trotter-Suzuki and product formulas [517, 424,
242], linear-combination of unitaries [423, 443] and multi-product formulas
[445], quantum walks [518, 519, 520], qubitization [521], and various
randomized approaches [522, 446, 244]. For the sake of simplicity, here we
will focus on the first of these, but note that hybrid oscillator-qubit
hardware facilitates any of these methods. For further discussion of
Trotterization in the context of compilation, see Section V.4.1.

In the Trotter-Suzuki approach, the unitary time evolution operator $U(t)$
corresponding to the Hamiltonian $H=\sum_{j}H_{j}$ is written as

$\begin{split}U(t)&=e^{-i\sum_{j=1}^{M}H_{j}t}\\\
&=\left(\prod_{j=1}^{M}e^{-iH_{j}t/r}\right)^{r}+O\left(\sum^{M}_{j,k}\|[H_{j},H_{k}]\|t^{2}/r\right)\end{split}$
(599)

where $r$ is a positive integer introduced to decrease the error in the
approximation by increasing the number of timesteps, and $M$ is the total
number of Hamiltonian terms. This decomposition reduces the problem to the
implementation of $e^{-iH_{j}t/r}$ for each $j$, via either a single native
gate or through compilation using a particular instruction set (see Section
V). For example, a chemical potential term $H_{j}=e^{-i\mu\hat{n}_{k}(t/r)}$
can be implemented via a phase-space rotation $R(\mu t/r)$ (see Box III.2) ,
while high-order interactions such
$H_{j}=Z_{i}a_{i}^{\dagger}a_{i}a_{k}^{\dagger}a_{k}$ can be synthesized using
BCH methods (see the examples portion of V.4). We note that higher-order
variants of Eq. (599) exist, such as the second order (or symmetric) Trotter
formula that achieves an improved error scaling at the cost of higher gate
count [242].

####  VI.4.2 State preparation

Initial state preparation is a key component for both static and dynamic
simulations. For static problems, state preparation is used to construct a
prior quantum distribution over eigenstates that can be refined through
variational techniques [427] or quantum phase estimation [523]. In general,
the question of whether or not quantum computers can efficiently prepare high-
fidelity approximations to ground states of physically relevant systems
remains an open question; however, if a generic quantum algorithm existed for
preparing such states then a consequence would be that $\BQP=\QMA$. $\BQP$ and
$\QMA$ are the quantum analogs of $\P$ and $\NP$ for classical computation
[1].

For calculating dynamics, quantum states need to be provided that correspond
to the initial quantum state $\ket{\psi(0)}$ that we wish to evolve into
$\ket{\psi(t)}=U(t)\ket{\psi(0)}$. It is likely however that calculating
observables of thermal states relevant to material simulations is easier than
for ground states [524, 525, 526].

Digital methods capable of determining low energy states using only shallow-
depth quantum circuits are of particular interest for simulation of statics in
the NISQ-era, with the variational quantum eigensolver (VQE) [427] and the
quantum approximate optimization algorithm (QAOA) [428] being two examples
that are particularly well suited to near-term digital oscillator-qubit
hardware [453, 527]. Annealing, or adiabatic state preparation approaches can
also be used to solve such problems or as intuition for an ansatz for QAOA or
VQE [428, 528]. However, the cost of these methods scales with the energy gap
between the ground state and the first excited state. This means that while
adiabatic algorithms can be used to construct arbitrarily accurate
approximations to the quantum ground state for gapped systems, the fact that
the excitation gap typically shrinks exponentially with the system size causes
such approaches to be inefficient in general.

In the fault-tolerant regime, we do not usually assume that fundamental
limitations exist on gate depth and several qubits/qumodes. This allows us to
use alternative methods that have exact bounds on success probability and
fidelity (or provable speedups) that can be implemented on hybrid oscillator-
qubit hardware, with the instruction sets presented in this work enabling the
implementation of controlled-unitaries implementing $e^{-i\theta H}$ or
various energy filter functions. These include quantum signal processing (QSP)
[148], explored in depth in section 23 and discussed in [37] in the context of
a simple lattice gauge theory, as well as quantum phase estimation (QPE)
[523].

With some modifications, the above preparation methods can be used to prepare
either ground or excited states. For example, in the context of variational
methods, instead of minimizing $\langle\hat{H}\rangle$ to target the ground
state, we can instead minimize $\langle(E-\hat{H})^{2}\rangle$ with $E$ the
estimated energy of a target excited state. In the context of fault-tolerant
algorithms, QPE for example can be used to prepare excited states as long as
the initial trial state has non-negligible overlap with the target excited
state of interest[529]. In the case of QSP or qubitization, a filter function
centering at the excited state energy can be implemented to facilitate excited
state preparation [148].

####  VI.4.3 Observables

Extracting useful information from a quantum state requires the ability to
make measurements and estimate observables. These range from local
observables, such as the energy density of a particular state or a time-
ordered correlation function of a spin, to nonlocal observables crucial to the
identification of exotic phases of matter and topological properties [530,
501, 502]. For qubit systems, measurements only in the computational basis are
inadequate to determine the complex phase factors in superposition states that
determine the expectation values of observables that have off-diagonal
components. Thus if measurements in the computational basis are the only ones
possible, it is essential to be able to rotate the state before measurement to
map off-diagonal observables into the computational basis, thereby effectively
creating the ability to measure in other bases. Similarly, in CV systems one
needs to be able to make measurements in various bases (e.g., the Fock (boson
number) basis and in the position and momentum quadrature bases).

Using the gate sets and techniques presented in this work, a wide range of
observables can be measured in hybrid oscillator-qubit hardware. We note that,
particularly in the NISQ-era, error mitigation strategies [373, 374] such as
zero-noise extrapolation (ZNE) and probabilistic error cancellation (PEC) play
an important role in the determination of observables from noisy CV [375] and
DV [491] hardware. Separately, classical shadow methods [531] – a class of
techniques to estimate observables from a small number of measurements – have
garnered much recent interest in the literature, and have recently been
extended to CV systems [396, 395]. We leave a more complete description of
these techniques to the above-cited works, and here simply note their
importance for the estimation of observables.

In the following paragraphs, we provide a few illustrative (but non-
exhaustive) examples of how one can measure Hamiltonian terms of interest in
hybrid oscillator-qubit hardware.

Number operator. The occupation of a mode can be measured in multiple ways.
First, one can probe whether or not a mode is in a particular Fock state $n$
by, for example, using an SQR gate (Sec. III.4.3) and encoding the answer in
the state of a qubit. However, asking this question for all $n$ (up to some
cutoff $n_{\textrm{max}}$ is costly, several attempts that is linear in
$n_{\textrm{max}}$. Alternatively, it has been demonstrated experimentally
[90] that one can make a binary search using successive Hadamard tests with
the parity operator, where the $k$th round reveals the $k$th bit of $n$ (i.e.,
$\lfloor n/2^{k}\rfloor\textrm{ mod }2)$ via phase-kickback on the qubit [90].
This latter approach scales logarithmically in $n_{\textrm{max}}$. For
details, see App. G. The method whose details are discussed in App. G.3 is
recompiled into the phase-space ISA in Sec. E.1.2.

Either technique enables the measurement of operators that are functions of
$n$, such as the chemical potential $\langle n\rangle$ and onsite interaction
$\langle n^{2}\rangle$, for example, or the measurement of other operator
strings which may include $\hat{n}$ using (shot-by-shot) correlated
measurements, such as the Hamiltonian terms discussed in Sec. VI.4.6.

Field operators. We can directly synthesize measurements of the resonator
field operators $a=\hat{x}+i\hat{p}$ (here in Wigner units) via measurements
of $\hat{x}$ and $\hat{p}$ [296]. To achieve this, recall from Eqs. (448-450)
that a conditional momentum boost $\mathrm{CD}(ik)$ (Box III.7) can be viewed
either as a momentum boost $\pm k$ conditioned on $\sigma_{z}$ of the qubit or
as a qubit rotation about the $z$ axis by an angle $\theta=-4k\hat{x}$ that is
linearly dependent on the position of the oscillator. Taking the latter
viewpoint, by starting the ancilla in the $|+\rangle$ state and then measuring
$\sigma_{y}$ after application of the conditional momentum boost, one obtains
a one-bit estimate of $\sin(4k\hat{x})$. For sufficiently small $k$ such that
the sine can be linearized, this is essentially a homodyne measurement of
$\langle\hat{x}\rangle$. Using this technique (and the corresponding one for
measuring the other homodyne quadrature, $\langle\hat{p}\rangle$, it is
possible to construct from repeated shots the moments of the homodyne
quadrature distributions and thus expectation value of operators such as
squeezing $\langle a^{\dagger}a^{\dagger}+aa\rangle$ in much the same way that
has been achieved for flying microwave photons using traditional microwave
homodyne techniques (and huge amounts of careful averaging to eliminate
amplifier/detector noise) [532].

Hopping terms. Many bosonic models include terms of the form $T_{ij}\propto
a^{\dagger}_{i}a_{j}+a_{i}a^{\dagger}_{j}$ describing hopping of a single
particle from site $i$ to site $j$. Notably, such a term exactly generates a
beam-splitter gate (see Box III.4). To measure this term, we can note that
conjugation by an appropriate beam-splitter gate maps $T_{ij}$ onto the photon
number difference between sites $i$ and $j$ via what is effectively
interferometry,

$\displaystyle\textrm{BS}^{\dagger}\left(\frac{\pi}{2},-\frac{\pi}{2}\right)T_{ij}\textrm{BS}\left(\frac{\pi}{2},-\frac{\pi}{2}\right)$
$\displaystyle=$ $\displaystyle
a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}-a^{\dagger}_{j}a^{\phantom{\dagger}}_{j}.$
(600)

Consequently, measurement of $T_{ij}$ is achieved by first applying the beam-
splitter $\textrm{BS}^{\dagger}\left(\pi/2,-\pi/2\right)$ followed by
measurement of the number operator at sites $i$ and $j$, in close analogy to
the combination of a Hadamard gate and $Z$-basis measurement facilitating an
$X$-basis measurement for qubits.

Bosonic Interaction Terms. Another useful capability is the measurement of
interaction terms involving multiple modes. For example, consider an
interaction term of the form $V_{ij}\propto
a_{i}^{\dagger}a_{i}a_{j}^{\dagger}a_{j}$. To measure $V_{ij}$, we can proceed
in multiple ways. First, we can use the previously described strategy for
measuring the number operator to independently determine $n_{i}$ and $n_{j}$
and, over repeated shots, correlate the results to determine $\langle
n_{i}n_{j}\rangle$ (which, importantly, differs from $\langle
n_{i}\rangle\langle n_{j}\rangle$ when $n_{i}$ and $n_{j}$ are correlated).
While sufficient as a terminal measurement (i.e., at the very end of a
circuit), this strategy extracts more information than needed, collapsing each
mode onto a Fock state in the process. In situations where we wish to measure
$V_{ij}$ without disentangling modes $i$ and $j$ (for example, if we expect
$V_{ij}$ to be a conserved quantity and want to check this mid-experiment),
the conditional cross-Kerr gate synthesized in Eq. (440) provides an
alternative approach as it rotates an ancilla by an amount proportional to
$V_{ij}$. In analogy with the digital homodyne measurement discussed above,
one could measure the sine of the rotation angle to estimate the mean value of
$V_{ij}$. In principle, one could also use the phase estimation algorithm to
do an efficient binary search for the value of the integer $n_{i}n_{j}$,
however, in practice, this will require repeated fast and accurate synthesis
of the cross-Kerr unitary in Eq. (440) for a variety of different rotation
angles.

####  VI.4.4 One Qubit, One Mode: The Rabi Model

One of the simplest yet most profound models in quantum physics is the Rabi
model [506, 533] which, at the simplest level, describes the interactions
between a harmonic oscillator and a single spin-1/2 system. Its applications
are wide-ranging, from the fundamental description of quantum optics and
molecular physics to real-world technologies such as nuclear magnetic imaging.
Along with the Jaynes-Cummings model [503] (which, as discussed below, the
Rabi model reduces to in a certain parameter regime), it is widely used to
describe light-matter interactions under the dipole approximation and
consequently serves as the foundation for the fields of cavity and circuit
quantum electrodynamics – the latter of which forms the very basis for the
superconducting hardware considered in this paper [54, 27].

The Rabi model is described by the Hamiltonian

$H_{\textrm{Rabi}}=\hbar\omega_{m}a^{\dagger}a+\hbar\omega_{q}\sigma_{z}+g\left(a^{\dagger}+a\right)\sigma_{x},$
(601)

where $\omega_{m}$ and $\omega_{q}$ are the frequencies of the oscillator and
the two-level system, respectively, while $g$ is the coupling strength between
the two, here taken (via a gauge choice) to be real. The first and second
terms of $H_{\textrm{Rabi}}$ correspond to the free evolution of the
oscillator and qubit, while the final term describes interactions between the
two. This is to a good approximation of the natural physical Hamiltonian
natively available (and always ‘on’) in circuit QED settings. However, we may
wish to simulate a Rabi Hamiltonian with with different (programmable and
possibly time-varying) parameters than are natively available in circuit QED
setups, trapped ions, or neutral atoms. Doing so within the instruction sets
presented in this work is straightforward: the first two terms can be
implemented via oscillator and single-qubit rotations (in software via simple
rotating frame changes), respectively, while the last is exactly the generator
of a conditional displacement, $\textrm{CD}(\alpha)$, up to single qubit
rotations.

As discussed in App. A, in the parameter regime where
$g\ll\\{\omega_{m},\omega_{q}\\}$, one can make the replacement
$\sigma_{x}\to\sigma_{+}+\sigma_{-}$ and discard rapidly-rotating, non-energy-
conserving terms $a^{\dagger}\sigma_{+}$ and $a\sigma_{-}$ under what is known
as the rotating-wave approximation. This results in the comparatively simpler
Jaynes-Cummings model,

$H_{\textrm{JC}}=\hbar\omega_{m}a^{\dagger}a+\hbar\omega_{b}\sigma_{z}+g\left(a\sigma_{+}+a^{\dagger}\sigma_{-}\right),$
(602)

which reduces light-matter interactions to the number-conserving exchange of
mode and qubit excitations, e.g., absorption ($a\sigma^{+}$) and emission
($a^{\dagger}\sigma^{-}$) of a single photon. However, as one increases the
coupling such that $g\ll\\{\omega_{m},\omega_{q}\\}$ no longer holds, the
system enters the so-called ultra-strong coupling regime ($g/\gtrsim
0.1\\{\omega_{m},\omega_{q}\\}$) and, for even larger $g$, the deep strong
coupling regime ($g\gtrsim\\{\omega_{m},\omega_{q}\\}$). It is predicted that
these parameter regimes display exotic behavior such as ground state virtual
excitations, and as such have been the topic of both theoretical and
experimental investigations over the past decade [534, 535].

The Rabi model constitutes perhaps the simplest nontrivial physical model
well-suited for simulation in hybrid oscillator-qubit platforms and can be
easily simulated within the instruction sets presented here.

####  VI.4.5 Many Modes: The Bose-Hubbard Model and the Fractional Quantum
Hall Effect

Several models in condensed matter, optical, and nuclear physics are purely
bosonic and therefore well-suited for simulation using hardware with native
bosonic modes. For example, boson sampling, in which photons pass through a
linear-optical network before measurement in the number basis, is a candidate
for quantum advantage [536, 301, 304, 305, 464]. Furthermore, the Bose-Hubbard
model has become a testbed for phenomena such as the well-known superfluid-to-
Mott-insulator transition, whose first observation spurred forward the field
of quantum simulation with cold atoms in optical lattices [496]. Moreover, it
is conjectured that the Bose-Hubbard model hosts topological phases with novel
dynamics akin to the ones encountered in the fractional quantum Hall effect
[537] (FQHE) for electrons and which could be observed in circuit QED systems
[538, 539, 540, 541, 542].

While these models can often be implemented in analog quantum simulators, such
approaches are limited in terms of Hamiltonian terms that can be realized and
observables that can be measured. For example, the complex phases in the
hopping term that incorporate effects of an external magnetic field (required,
for example, in the description of the FQHE) are challenging to implement in
cold atoms due to their neutrality. Synthetic magnetic fields have been
employed in past work [543], but an FQHE-like state has only so far been
observed in a few small-scale experiments in cold atom [544] and
cavity/circuit QED platforms [545, 540, 546]. Furthermore, entanglement-based
measures, important for verifying the topological properties of candidate FQHE
states [502, 501], are challenging to obtain; cold atom experiments, for
example, have typically relied on measurement of occupations which supplies no
information regarding coherences [544] (though we note a recent proposal for
measuring arbitrary properties in analog simulators [547]).

Digital simulation with oscillator-qubit hardware presents advantages on all
of these fronts. Not only does the hardware natively support bosonic modes and
operations, but arbitrary Hamiltonian terms, such as magnetic field terms
crucial to the observation of the FQHE can be easily synthesized by adjusting
the phases of the beam-splitter couplings on each link of the lattice, using
the instruction set architecture presented in this work. The programmability
of such simulators allows one to easily tune Hamiltonian parameters and
synthesize additional terms. Finally, as previously discussed, digital
hardware offers the ability to measure a wide range of observables that
otherwise present a challenge in analog simulations.

With the above motivations in mind, we now turn explicitly to the Bose-Hubbard
model as an illustrative example and discuss its digital simulation in hybrid
oscillator-qubit hardware. Its Hamiltonian (for zero magnetic field) is given
by

$\begin{split}H_{\textrm{BH}}&=-\mu\sum_{i}\hat{n}_{i}+\frac{U}{2}\sum_{i}\hat{n}_{i}\left(\hat{n}_{i}-1\right)+H_{\textrm{hop}}\\\
H_{\textrm{hop}}&=-J\sum_{i}\left(a_{i}^{\dagger}a_{i+1}+\mathrm{h.c.}\right),\end{split}$
(603)

where $\mu$ is the chemical potential, $J$ the hopping strength, and $U$ the
interaction strength (with $U>1$ repulsive and $U<1$ attractive). Time
evolution $e^{-iH_{i}\Delta t}$ under each term $H_{i}$ for a duration $\Delta
t$ can be directly realized in the Fock-space ISA: $H_{\textrm{hop}}$
generates a beam-splitter $\textrm{BS}(-2J\Delta t,0)$, the chemical potential
term is realized using a phase-space rotation gate $R(-\mu\Delta t)$, and the
onsite interactions can be implemented using a SNAP gate (Box III.10) by
choosing $\varphi_{n}=Un(n-1)\Delta t/2$, such that

$\begin{split}\textrm{SNAP}(\varphi_{0},\varphi_{1},\varphi_{2},\ldots)&=e^{-i\sum_{n}Un(n-1)\Delta
t\ket{n}\bra{n}/2}\\\ &=e^{-i\sum_{n}U\hat{n}(\hat{n}-1)\Delta
t/2}.\end{split}$ (604)

As a simple illustration, Fig. 43b displays the Trotterized evolution of the
mode occupation implemented in Bosonic Qiskit [102].

Figure 43: Simulation of a four-site Bose-Hubbard model using Bosonic Qiskit.
(a) A schematic of the Bose-Hubbard model, illustrating nearest-neighbor
hopping with strength $J$ and on-site repulsion of strength $U>0$. This
competition between these Hamiltonian terms leads to a superfluid-to-Mott-
insulator phase transition [496]. (b) Occupancy as a function of time for
various choices of $J$ and $U$. In all cases, a time step of $\Delta t=0.2$
was used and the system was initialized to have two bosons at each edge site,
with intermediate sites vacant.

While here we have simplified the case where $U$, $\mu$, and $J$ are identical
across all sites, we note that additional possibilities arise by varying these
parameters across the lattice. For example, it has been shown in cold atomic
platforms that the 1D Bose-Hubbard model with the added disorder in the
chemical potential $\mu_{i}$ leads to many-body localization [548].
Separately, alternation of the hopping strength $J_{i}$ leads to an
(interacting) extension of the Su–Schrieffer–Heeger model, known to display
topological features such as edge states [505].

Moreover, magnetic fields and strong interactions in the Bose-Hubbard model
give rise to fractionally charged quasi-particles, long-ranged entanglement,
and anyonic exchange statistics [537, 549]. In this direction, the selection
of the beam-splitter phase $\phi$ allows one to implement hopping in a
magnetic field, thus facilitating potential investigation of the FQHE. In all,
this multitude of possibilities constitutes a rich and promising playground
for quantum simulation of many-body bosonic models using hybrid oscillator-
qubit hardware.

####  VI.4.6 Large Spin Models: The Transverse-field Ising Model

Several large-spin models, i.e., those with spin-$S$ degrees of freedom with
$S>1/2$, are of interest in condensed matter physics [513, 381] and chemistry
[550]. While individual qubits are naturally suited to represent spin-1/2
degrees of freedom, encoding a generic spin-$S$ system in qubit-based hardware
requires $\log(2S+1)$ qubits. Though not prohibitively expensive in terms of
qubit resources, this approach comes at an additional cost in the two-qubit
gate count required to implement both single- and multi-spin-$S$ operations,
prohibiting straightforward simulation, particularly on near-term devices.

In contrast, in hybrid oscillator-qubit hardware, one can make use of the
larger Hilbert space of the bosonic modes to represent $d$-level systems in a
hardware-efficient manner. For example, Ref. [47] demonstrated universal
control on the lowest-$d$ Fock levels of a single bosonic mode for arbitrary
$d$ with simple Jaynes-Cummings interactions, which is in principle flexible
enough to simulate large-spin models. Other promising avenues include qudit
error correction codes, such as the $d$-level generalization of GKP [187] and
cat codes [113].

While these strategies can be used to represent an arbitrary $d$-level system
(and, by extension, a spin-$S$ system with $d=2S+1$), here we discuss an
alternative approach which _directly_ implements a spin-$S$ and its defining
commutation relations in terms of simple combinations of bosonic operators.
This possibility was first noted by Julian Schwinger [551], and we therefore
refer to it as the Schwinger boson representation (SBR) henceforth.

The core idea is to first recognize that, for two oscillators with
annihilation operators $a$ and $b$, if one restricts to the subspace of total
photon number $N$ (consisting of $N+1$ total states), then the operators
$a^{\dagger}b$ and $ab^{\dagger}$ act exactly as spin-$S$ raising and lowering
operators $S^{+}$ and $S^{-}$ within this subspace, with the total spin being
given by $S=N/2$. For example, if we take $N=1$, this subsystem consists of
the Fock states $\ket{1,0}$ and $\ket{0,1}$, with $S^{+}=a^{\dagger}b$ and
$S^{-}=ab^{\dagger}$ acting exactly as raising and lowering operators for an
effective spin-1/2 system (with $d=2S+1$ possible states). Notably, this case
is precisely the idea behind the dual-rail qubit in a pair of oscillators (see
Section IV.4.3).

More generally, one can represent spin-$S$ systems of arbitrary $S$ by
restricting to the subspace of $N=2S$ photons distributed between a pair of
oscillators (representing a spin-$S$ in a single oscillator is also possible
[552, 553]). To see this, we can make the assignment

$\begin{split}S^{x}&=(a^{\dagger}b+ab^{\dagger})/2\\\
S^{y}&=(a^{\dagger}b-ab^{\dagger})/2i\\\
S^{z}&=(a^{\dagger}a-b^{\dagger}b)/2,\\\ \end{split}$ (605)

which exactly obeys the commutation relations

$[S^{x},S^{y}]=iS^{z}\,\,\,\mathrm{et\ cyc.}$ (606)

We also see that the beam-splitter (Box III.4) can be written as

$\mathrm{BS}(\theta,\phi)=e^{-i\theta[\cos(\phi)S_{x}-\sin(\phi)S_{y}]}.$
(607)

Beam-splitters alone therefore give us access to two generators of rotations
and allow for arbitrary SU(2) rotations of the two-mode qudit. We note that
the third SU(2) generator, $S^{z}$, can be implemented with phase-space
rotations (see Box III.2). For the case of a spin-$1/2$ ($N=1$), the two modes
form a two-level qubit and SU(2) rotations consequently provide universal
single-spin operations.

For $S>1/2$ ($N>1$), however, the most general unitary single-spin operation
lives in SU($d$) with $d=N+1>2$, and beam-splitters alone therefore do not
provide universal control. Instead, we need additional hybrid oscillator-qubit
operators to synthesize generators nonlinear in the spin operators. To that
end, one can recognize that any polynomial of $S^{z}$ is diagonal in the Fock
basis. Consequently, it is possible to realize nonlinear unitaries of the form
$U=\textrm{exp}\,(-i\theta f(S^{z}))$ via a single SNAP gate on mode $a$.
Here, $f(S^{z})$ is an arbitrary polynomial in $S^{z}$, and we have
additionally used the fact that
$S^{z}=(a^{\dagger}a-b^{\dagger}b)/2=a^{\dagger}a+N/2$. Together with beam-
splitter enabled $SU(2)$ rotations, this nonlinear gate provides universal
control for $SU(d)$. This is analogous to the universal control of high-spin
atoms in atomic physics utilizing the non-linear Zeeman effect [554, 354,
555].

Measurement of the spin-$S$ system in the $S^{z}$ basis has a natural
implementation in the Fock-space ISA – one must simply resolve the number
operator $a^{\dagger}a$ (or $b^{\dagger}b$), which can then be used to
determine $S^{z}$ (see Appendix G for details on photon number measurement).
By extension, this ability enables the measurement of the spin along any
direction, $\hat{n}\cdot\vec{S}$, by first rotating to the $\hat{z}$-basis
with an appropriate beam-splitter (much like a Hadamard gate prepending a
$Z$-basis measurement enables an $X$-basis measurement for a qubit). With this
primitive, one can more generally carry out quantum state tomography on a
spin-$S$ system by measuring $\hat{n}\cdot\vec{S}$ along a minimum of $4S+1$
distinct directions of $\hat{n}$ [556]. To reduce statistical uncertainty in
the state construction, however, it is highly beneficial to increase the
number of distinct measurement axes beyond the bare minimum [557]. Notably,
the tomographic protocols described here require only beam-splitters and
photon-number measurements within the SBR framework, independent of the value
of $S$.

As a simple example of how the SBR can be leveraged for quantum simulation,
let us consider the 1D transverse-field Ising model of spin-$S$. Its
Hamiltonian is

$\displaystyle
H_{\textrm{TFIM}}=-\sum_{i}S_{i}^{z}S_{i+1}^{z}-\lambda\sum_{i}S^{x}_{i},$
(608)

where the first and second terms correspond to spin-spin interactions and
coupling to a transverse field, respectively. To implement this model using
the SBR, we envision that the spin at site $i$ is composed of two oscillators
with corresponding annihilation operators $a_{i}$ and $b_{i}$. To simulate
dynamics according to this Hamiltonian, we can rely on Trotter-Suzuki methods
discussed in Sections V.4.1 (about compilation methods) and VI.4.1 (for a
dynamics-focused discussion), reducing overall compilation to that of the
spin-spin terms and transverse field terms alone.

The transverse field terms can be implemented via beam-splitter gates of the
form $\textrm{BS}(\theta,0)$, following the discussion surrounding Eq. (607).
The spin-spin term, on the other hand, is more challenging given that it is
quadratic in spin operators. To see how to proceed, we can first recast in
terms of bosonic operators,

$S_{i}^{z}S_{i+1}^{z}=\frac{1}{4}(a_{i}^{\dagger}a_{i}-b_{i}^{\dagger}b_{i})(a_{i+1}^{\dagger}a_{i+1}-b_{i+1}^{\dagger}b_{i+1})$
(609)

and can additionally note that the four above terms commute, such that

$\begin{split}e^{-i\theta
S_{i}^{z}S_{i+1}^{z}}&=e^{-i\frac{\theta}{4}\hat{n}_{a_{i}}\hat{n}_{a_{i+1}}}e^{-i\frac{\theta}{4}\hat{n}_{b_{i}}\hat{n}_{b_{i+1}}}\\\
&\times
e^{i\frac{\theta}{4}\hat{n}_{a_{i}}\hat{n}_{b_{i+1}}}e^{-i\frac{\theta}{4}\hat{n}_{b_{i}}\hat{n}_{a_{i+1}}}.\end{split}$
(610)

Taking advantage of the fact that $S=(\hat{n}_{a_{j}}+\hat{n}_{b_{j}})/2$ is
conserved, we can simplify this to

$\displaystyle e^{-i\theta S_{i}^{z}S_{i+1}^{z}}$
$\displaystyle=e^{-i\theta(\hat{n}_{a_{i}}-S)(\hat{n}_{a_{i+1}}-S)}.$ (611)

Therefore, synthesis of $\textrm{exp}(-i\theta S_{i}^{z}S^{z}_{i+1})$ reduces
to that of generators that are products of number operators of distinct
oscillators. While challenging to implement exactly, we can note that each
unitary is equivalent to the conditional cross-Kerr gate approximately
synthesized using BCH techniques in Eq. (440) if we make use of an ancilla
qubit in the state $\ket{0}$. Thus, the digital simulation of the transverse-
field Ising model within the gate sets presented here is straightforward.

The strategy of using the SBR to simulate spin systems is of course much
broader than this simple example, with possible applications of classical and
quantum spin models spanning the fields of magnetism [558], materials [559],
artificial intelligence [560], and even neural processing [561]. Furthermore,
because each spin is entirely encoded in a photon-number-preserving subspace,
the SBR can be combined with error detection strategies; similar to the dual-
rail qubit, one can detect photon loss errors by making measurements of the
joint-parity in the corresponding pair of oscillators [122, 124]. Thus,
studying spin models using the SBR in hybrid oscillator-qubit hardware not
only provides an avenue for direct study of a wide range of interesting models
but does so in a particularly promising NISQ-friendly manner by offering a
degree of noise resilience through post-selection.

####  VI.4.7 Many Qubits, Many Modes: $\mathbb{Z}_{2}$ Lattice Gauge Theory

A variety of interesting Hamiltonians include both qubits and bosonic modes
and are thus a natural fit for the hybrid oscillator-qubit hardware. Examples
include multi-spin-boson [562] models, bosonic lattice gauge theories [510],
the Jaynes-Cummings-Hubbard model describing many-body states of (dressed)
light [563, 564], the exciton-vibrational Hamiltonians governing energy
transfer and relaxation pathway in molecules and materials [565, 566].
However, such Hamiltonians are challenging to implement both on qubit-based
hardware, due to the inefficiency of implementing bosonic gates, and analog
quantum simulators, due to the multi-species nature of the models. Hybrid
oscillator-qubit platforms provide a promising hardware-efficient approach for
the quantum simulation of these Hamiltonians due to the ease with which multi-
oscillator, multi-qubit, and hybrid oscillator-qubit gates can be implemented.

Figure 44:  Implementation of the $\mathbb{Z}_{2}$ lattice gauge theory in
Bosonic Qiskit. Both panels display the average occupation at each site as a
function of time. The initial state is chosen such that there is one boson at
the central site, and all others vacant. For Gauss’s law, we fix $G_{j}=+1$.
In the left panel, we show results for $g/J=0.1$ (hopping-dominated) and
observe a light-cone spread of the boson. In the right panel, we show results
for $g/J=5$ (field-dominated) and observe that the boson is (mostly) trapped
at the central site. In both cases, a step size of 0.1$J$ was used.

Here, we take as an example a simple 1D bosonic $\mathbb{Z}_{2}$ lattice gauge
theory [510]. As in the previous sections, we aim to provide a high-level
discussion of the strategy to simulate this model. Readers interested in an
in-depth discussion and extension to $2$D should consult Ref. [37].

The $\mathbb{Z}_{2}$ lattice gauge model we consider describes bosonic matter
that hops along a 1D chain of sites alternated by qubits that represent a
$\mathbb{Z}_{2}$ gauge field living on the links of the lattice. The
Hamiltonian of this model is given by

$\begin{split}H_{\textrm{LGT}}&=H_{\textrm{hop}}+H_{\textrm{field}}\\\
H_{\textrm{hop}}&=-J\sum_{i}a^{\dagger}_{i}Z_{i}a_{i+1}+a_{i}Z_{i}a_{i+1}^{\dagger}\\\
H_{\textrm{field}}&=-g\sum_{i}X_{i},\end{split}$ (612)

where $a_{i}$ ($a_{i}^{\dagger}$) removes (adds) a particle at site $i$, and
$X_{i}$ and $Z_{i}$ are Pauli operators for the qubit at link $i$ (which
“connects” site $i$ and site $i+1$). To better understand this model (and the
competition between $H_{\textrm{hop}}$ and $H_{\textrm{field}}$ it describes),
it is helpful to first imagine all qubits are initialized in an $X$
eigenstate. $H_{\textrm{field}}$ describes the energy of the field and is
minimized by aligning all qubits in the $+X$ direction. In contrast,
$H_{\textrm{hop}}$ describes bosons hopping between sites, causing a phase-
flip (sending $+X$ to $-X$ and vice-versa) at each intermediary link, altering
the field energy. Consequently, the model displays drastically different
behavior in opposite parameter regimes, with $|J/g|\gg 1$ leading to a
“deconfined” phase where matter freely hops, and $|J/g|\ll 1$ leading to a
“confined” phase where the large field inhibits hopping [510].

Much like in electrodynamics, we can define an analog of Gauss’s law that
relates the number of bosons (or _charges_) at site $i$ with the product of
“left” and “right” field values $X_{i}X_{i+1}$ (akin to a field integral in
electrodynamics). In particular, we can identify this gauge symmetry with the
generators

$\hat{G}_{j}=X_{j}e^{i\pi a^{\dagger}_{j}a^{\phantom{\dagger}}_{j}}X_{j+1}.$
(613)

To see that the generators $\hat{G}_{j}$ correspond to a conserved quantity,
note that when a boson hops from site $j$ to site $j+1$, the factors of
$\hat{G}_{j}$ change as follows: (1) $a^{\dagger}_{j}a_{j}$ is decremented by
one, flipping the sign of the photon number parity $e^{i\pi
a^{\dagger}_{j}a_{j}}$, and (2) the presence of $Z_{j+1}$ in the hopping
amplitude flips the sign of $X_{i+1}$, leaving $\hat{G}_{j}$ unchanged
overall. In other words, this means that $[H_{\textrm{hop}},\hat{G}_{j}]=0$.
Noting that $\hat{G}_{j}$ trivially commutes with $H_{\textrm{field}}$,
$\hat{G}_{j}$ is therefore a conserved quantity. Such a feature is
advantageous for simulations on noisy hardware, as we can think of the
generators $\hat{G}_{j}$ as stabilizers in an error-detecting code in the
sense that we can perform error mitigation by post-selection of results that
obey the Gauss law constraints. For more discussion on this approach, see Ref.
[37].

Finally, we comment on the implementation of time evolution under this
Hamiltonian. Both $H_{\textrm{field}}$ and $H_{\textrm{hop}}$ can be
implemented directly with gates previously discussed: $H_{\textrm{field}}$ is
simply a single-qubit rotation, while $H_{\textrm{hop}}$ corresponds to a
controlled beam-splitter – see Eq. (395) and the surrounding text. As a
demonstration, we implement Trotterized time dynamics of $H_{\textrm{LGT}}$ in
Bosonic Qiskit for both hopping-dominated and field-dominated regimes in a
five-site system. We initialize a boson at the central site, with all others
vacant.

Fig. 44 displays the resulting time-evolved occupancy at each site. The
competing order between $J$ and $g$ is evident, as the hopping-dominated
regime gives rise to a light-cone spread of a boson initialized to the central
site, whereas the field-dominated regime prohibits hopping due to the
prohibitive cost to flip link qubits.

####  VI.4.8 Outlook: Fermion-Boson Problems

Classically intractable fermion-boson problems include non-equilibrium
dynamics in QCD relevant to heavy-ion collisions [567], calculation of
vibrational spectra of large anharmonic molecules coupled to an environment
for identification of unknown samples [476], and the study of phonon dynamics
in strongly correlated materials such as high-temperature superconductors
[568]. Hyrbid oscillator-qubit hardware is well-suited to such problems, as
the native oscillators naturally encode bosonic excitations, while fermionic
degrees of freedom can either be mapped to onto physical or logical qubits or
to high-$Q$ oscillators [569, 570], depending on the specifics of the problem
at hand. Furthermore, the instruction set architectures we have presented
offer universal control, such that arbitrary fermionic and boson gates can be
synthesized. For an in-depth manual for using CV-DV hardware for simulating
fermions and bosons, including gauge theories in 2D, see Ref. [37].

A wide range of encodings can facilitate the mapping of fermions onto qubits.
Well-known options include Jordan-Wigner [78], Bravyi-Kitaev [34], Verstraete-
Cirac [571], and Derby-Klassen [572] encodings, all of which enable the
encoding of fermionic degrees of freedom using either logical qubits derived
from CV-to-DV, DV-to-DV, or concatenated codes (see Section IV.1.4). For near-
term applications, one could also use the physical qubits themselves. Notably,
fermion-to-qubit mappings require non-local qubit-qubit gates which, in this
case, can be implemented at the physical layer using the oscillator-assisted
multi-qubit entangling gate techniques described in Section V.3.2. However, in
the case of superconducting hardware, it is important to note that transmon
qubit lifetimes are typically on the order of $\sim 100$ $\mu$s, whereas 3D
microwave cavity modes have far longer lifetimes on the order of $\sim 1$ ms
[23]. Therefore, even without proper error correction, it may be advantageous
to consider encoding the fermions into oscillators using, for example, the
lowest-lying Fock states of cavity modes or dual-rail/Schwinger-boson qubits.

As a simple demonstration, we now illustrate the compilation of the so-called
fermionic simulation gate [573, 574],

$\textrm{fSim}(\theta,\phi)=\left(\begin{matrix}1&0&0&0\\\
0&\cos(\theta)&-i\sin(\theta)&0\\\ 0&-i\sin(\theta)&\cos(\theta)&0\\\
0&0&0&e^{-i\phi}\end{matrix}\right),$ (614)

a particular application of which is the fermionic-SWAP or fSWAP gate,
equivalent to a SWAP gate up to an additional minus sign applied to the
$\ket{11}$ state to account for antisymmetric fermionic exchange statistics
(see Table A.1). Specializing to the case where we encode qubits in the
$\ket{0}$ and $\ket{1}$ Fock states of a pair of oscillators indexed $j$ and
$k$, we can compile the fSim gate as

$\begin{split}\textrm{fSim}(\theta,\phi)&=\textrm{exp}(-i\phi\sigma_{z}a_{j}^{\dagger}a_{j}a_{k}^{\dagger}a_{k})\times\textrm{BS}\left(2\theta,0\right),\end{split}$
(615)

where the first term is a cross-Kerr gate, synthesized using BCH methods in
Eq. (440) and controlled on an ancillary qubit prepared in the state
$\ket{0}$. Alternatively, one could replace this approximate cross-Kerr gate
with a $\textrm{CPHASE}(\phi)$ gate realized using the methods of Ref. [124].

As previously mentioned, the fermionic SWAP operator can be implemented using
an fSim gate,

$\begin{split}\textrm{fSWAP}&=\left[R\left(-\frac{\pi}{2}\right)\otimes
R\left(-\frac{\pi}{2}\right)\right]\times\textrm{fSim}\left(\frac{\pi}{2},0\right)\\\
&=\begin{pmatrix}1&0&0&0\\\\[-10.0pt] 0&0&1&0\\\\[-10.0pt]
0&1&0&0\\\\[-10.0pt] 0&0&0&-1\end{pmatrix},\end{split}$ (616)

where $R(\theta)$ corresponds to a phase-space rotation of the oscillator, and
we have furthermore displayed only the portion of the unitary acting within
the encoded two-qubit subspace.

Here we have focused on qubits encoded in Fock states $\ket{0}$ and $\ket{1}$
for a pair of cavities. We present fermionic swap networks in the transmon
qubits as well as their advantages in the dual rail encoding in Ref. [37].

##  VII Conclusions and Outlook

The primary contribution of our work is the introduction of explicit
instruction set architectures (ISAs) and abstract machine models (AMMs) for
hybrid CV-DV quantum processors composed of qubits and quantized bosonic
modes. We start from the physical behavior of bosonic states and
systematically develop bosonic gate models, and hybrid architectures including
error correction strategies, compilation theory and techniques, as well as
high-level proxy applications and algorithms. We have presented a few examples
illustrating how such hybrid quantum processors can be utilized to solve
practical problems, including DV-CV state transfer, the quantum Fourier
transform, bosonic quantum error correction, and quantum simulation of matter
with bosonic degrees of freedom. We show for the latter that hybrid digital
architectures may provide substantial computational advantages over existing
approaches that use only qubits because of the ability to avoid costly
arithmetic operations in favor of equivalent native bosonic operations.

The physical-level interactions that enable the proposed ISAs and AMMs can be
realized on quantum hardware containing both qubits and bosonic degrees of
freedom, and we have discussed here superconducting, trapped ion, and neutral
atom platforms. Dissipation and quantum noise on such hardware can be
corrected or suppressed with bosonic quantum error-correcting codes
(potentially further concatenated with qubit codes). These physical
interactions allow us to define several instruction sets summarized in Table
III.4, including among others, the Fock-space, phase-space, and sideband ISAs.
We have also shown explicitly how instructions from these ISAs can be composed
to achieve universal hybrid CV-DV quantum computation and further demonstrate
how to characterize the quantum device with a continuous variable analog of
gate set tomography.

Systematic compilation techniques can compose these universal gates into
useful primitives and algorithms that solve practical problems, as we
discussed in Sec. V. These compilation techniques include strategies for
compiling hybrid single-qubit oscillator unitaries, which are obtained by
generalizing the well-known Trotter-Suzuki and product formulas, quantum
signal processing, and linear combination of unitaries techniques from DV
systems to hybrid CV-DV systems. We presented exact analytical gate synthesis
methods for multi-qubit entangling gates and multi-oscillator entangling gates
in superconducting systems that are inspired by the Mølmer-Sørensen gate for
ion traps.

At the system level, we defined three key quantum AMMs: 1) qubit-centric, 2)
oscillator-centric, and 3) hybrid oscillator-qubit, which all have universal
computation capabilities. These AMMs can represent abstractions of the machine
for a high level of the stack that is user-facing, defining the mathematical
abstract Hilbert space where the computation happens. Alternatively, they can
be lower-level abstractions that are transpiler-visible, presenting to the
transpiler what physical degrees of freedom are available as computation
components and what their connectivity and allowed gates are. The formal
establishment of these AMMs should prove particularly useful for resource
estimation on hybrid CV-DV quantum processors.

We further examined applications of quantum computing, such as quantum
simulation, illustrating how quantum electrodynamics and related lattice gauge
Hamiltonians can be compiled into hybrid qubit/bosonic operations. We find
evidence of advantage for near-term applications of hybrid systems over
standard methods for purely qubit-based simulations. In addition, we presented
extensions of quantum signal processing to the non-abelian case for which the
qubit rotation angles are functions of the non-commuting position and momentum
variables of bosonic modes. We examined applications of quantum signal
processing within these architectures and discussed explicit approaches to
implement the quantum Fourier transform.

In conclusion, the ISAs and AMMs proposed here significantly expand our
understanding of how distinct quantum resources can be leveraged in a unified
computational framework and demonstrate the potential of DV-CV hybrid quantum
processors for the upcoming post-NISQ era of early-stage fault tolerance.

###  VII.1 Open Questions and Directions for Future Work

Despite these achievements, there are still many open questions ranging from
the fundamental side to more practical aspects of hybrid DV-CV quantum systems
that remain to be addressed in the future.

First and foremost, as we have seen in Sec. VI, hybrid quantum hardware opens
many possibilities for interesting applications. Identifying more use cases
where CV-DV hardware can deliver an advantage over qubit-only or classical
computers will be important. Examples of such applications likely include
quantum simulation of coupled nuclei-electron processes [575, 576] or light-
matter interactions [577, 578] relevant to chemical reactions and materials
science applications. As another example, the hybrid CV-DV quantum signal
processing framework established in the current work serves as a potentially
ideal tool for processing quantum signals. Classical signal processing
techniques from electrical engineering for filter design and classical CV-DV
conversion [579] can provide much intuition and guidance for further
developments in the quantum case.

And while the major quantum algorithms become united in the framework of
quantum signal processing and quantum singular value transforms, the newly
emerging hybrid CV-DV quantum signal processing techniques extend
possibilities to the non-abelian case, for which a complete theory remains to
be discovered. From a NISQ viewpoint, further development is required for
variational CV-DV algorithms beyond Sec. V.6 that can be run on such hybrid
processors[453]. Such a set of tools would dramatically increase the scope of
the signal-processing transformations that can be implemented on input
parameters and allow more powerful analytic compilation methods for general
multi-variable transformations based on unitary transformations of bosonic
modes. While we describe some QSP methods in this work, the presentation is
mathematically abstract, and higher-level CV-DV quantum programming languages
such as Bosonic Qiskit [102] are needed to simplify the programming process
for a broader community of users.

We have presented various CV-DV ISAs based on discrete sets of control
Hamiltonians but continuously parameterized gates. An open question in this
direction is to address the complexity of generating arbitrary unitaries in CV
systems. While the Trotter-Suzuki-based approaches considered here and in [69]
can provide sub-polynomial scaling with the error tolerance for the
synthesized gates, existing approaches fall short of the poly-logarithmic
scaling provided by circuit synthesis methods such as the Solovay-Kitaev
theorem [580, 581] which allows poly-logarithmic scaling for the error
tolerance when compiling over a discrete gate set. The development of better
CV gate synthesis methods may not only open up new possibilities in compiling
unitaries but also open up the door to the development of bosonic error
correcting codes where the physical and logical operations are natural bosonic
operations such as discrete displacement or squeezing operations.

The non-Abelian QSP techniques used in this work involve composite pulses in
phase space designed to cancel first-order errors. Generalizations of these
composite pulses similar to Ref. [410, 177, 444] are required to develop
optimal sequences for arbitrary CV Hamiltonian construction. Formally, this
task is equivalent to designing a constructive algorithm that can find optimal
non-abelian QSP sequences for arbitrary polynomials $f(\hat{x},\hat{p})$
[181]. In Refs. [582] the authors find a constructive algorithm to generate
arbitrary polynomials for the case of single-variable QSP. The generalization
of this scheme to QSP using two commuting variables was shown in Ref. [429].
Such a construction for non-abelian QSP extends to the case of multi-variable
QSP with non-commuting variables, and thus, remains an open question. The
multiple non-commuting variables for the case of the phase-space ISA will be
the numerous quantum phase space vectors $\hat{v}=\alpha\hat{x}+\beta\hat{p}$.
The examples from Ref. [435] discussed in Secs. V.4.2 and V.5 belong to a
subclass of this framework, QSP with just two non-commuting variables. A
hierarchy of various QSP frameworks developed with an eye towards the hybrid
phase-space instruction set can be found in Ref. [435].

A related aspect of algorithm design is to identify new useful universal
instruction sets for such hybrid quantum processors. In particular, we have
seen in Sec. II.2.4 that the stellar representation provides an intuitive
picture of how CV states and operations can be characterized in terms of the
roots of complex analytic functions and their associated time dynamics.
Designing an ISA that works directly within the stellar representation picture
would be very interesting and useful.

From an architectural perspective, we need more complete ISAs that include a
memory hierarchy and sophisticated quantum algorithmic and arithmetic logic
units (ALUs). In terms of memory, a mixed system of qubits and bosonic modes
provides opportunities by balancing the fast gates and high connectivity of
qubits with the long lifetimes and novel gate sets of bosonic modes. The
needed high connectivity can be obtained through quantum communication using
bosonic modes, for example by the cavity SWAP gates discussed in Sec. V.3
and/or through use of the LAQCC model discussed in Secs. IV.2 and IV.5.1 to
efficiently produce qubit GHZ states that can be used for quantum
teleportation.

At the simplest level, recent architectures [263] have proposed using bosonic
modes as small register memories, embedding them in architectures specifically
designed to match the connectivity requirements of algorithms. Recent work
with multi-mode bosonic systems has shown their utility as memories for
factoring [583] and unique error correction architectures [584]. This is
particularly important if proposed memory hierarchies are to be realized [261]
as in classical systems. Building upon these memory hierarchies, more
sophisticated quantum ALUs that take advantage of the hybrid CV-DV gates
(Table LABEL:tab:gates-osc and LABEL:tab:gates-qubit-osc) and the continuous-
variable nature of oscillator dynamics need to be developed. Sec. VI
demonstrated key quantum primitives including the Fourier transform and
Hamiltonian simulation, but creating more advanced primitives and proxy
applications, such as modular exponentiation and matrix arithmetic, remains to
be done. Instruction- and processor-level parallelism often used in modern
classical computers [38] could also be deployed based on a more sophisticated
ISA and memory hierarchy.

Qubit gate fidelities can be efficiently measured using randomized
benchmarking and related techniques. To date, however, there are no
corresponding systematic methods for CV operations, due to the absence of
unitary $2$-designs for CV systems. While we provide in Sec. IV.5.2 an
approach to gate set tomography using Wigner function tomography (see App.
G.6), this approach cannot be generalized to arbitrarily high dimensions
without making strong assumptions about the quantum channels being learned.
Such channel tomography estimates will be necessary for us to assess the
degree to which these systems will need quantum error correction.

Regarding the topic of noise and reliability, we have discussed examples of
bosonic error correction codes that encode discrete logical qubits into
oscillators. Creation and quantum control of such bosonic logical qubits
requires non-Gaussian operations that in turn require the presence of a
transmon or other auxiliary qubit or anharmonic element which typically has
poorer coherence properties than the oscillator they are controlling. This
raises fundamental fault-tolerant control and measurement challenges which are
only now beginning to be addressed [585, 129, 122, 169, 314, 315] and many
opportunities remain for further progress. An interesting future direction
would be to build fault tolerance into the cost function when designing
numerically optimized quantum control sequences.

To bring out the role of bosonic error correction in terms of fault tolerance,
we have discussed this direction in Sec. IV.1.4. We discuss how oscillators
can be useful in improving the cost of magic state distillation, a resource-
intensive component of fault-tolerant quantum computing (FTQC). So far, all
versions of quantum error correction have deemed magic state distillation to
be a necessary component to implement non-Clifford gates, an essential
requirement for quantum advantage. We pointed out in Sec. IV.1.4 that
oscillator-based codes can quadratically reduce the overhead cost associated
with this FTQC component by changing the noise profile of the underlying qubit
[363]. In addition, we highlight that there are efficient ways to realize an
error-corrected qudit in an oscillator. Qudits can be helpful in various ways
including reduction of the computation cost and one may expect that low error
qudits will be high in demand in the future. The ability to engineer high-
fidelity qudits via error correction in the oscillator also takes improvement
in magic state distillation several steps further. It has been shown that
qudits can directly improve the resource overhead of this preparation by
several orders of magnitude by directly enabling more efficient distillation
schemes [586, 587].

We have also briefly discussed the more challenging question of error
correction of oscillators (qua oscillators, and not as logical qubits). So far
there has been only one attempt at such a construction, given by Ref. [357],
which corrects errors under a classical displacement channel. Soon after, in
Ref. [360] it was shown that code families designed under this noise model
assumption do not have a threshold. Given that the true noise model of an
oscillator is dominantly photon loss, it remains an open question if there is
an oscillator code family that can protect oscillators against this more
realistic noise. Formally if we introduce a cutoff on the boson number the
problem reduces to qudit error correction for large qudit dimension $d$, but
it seems extremely challenging to efficiently preserve the oscillator states
and operator algebra for large $d$. If such a construction exists, it would be
a revolutionary step in reducing the complexity of many quantum simulations
and CV algorithms potentially yielding a quantum advantage. However, as
analogous ideas for classical computing have failed to deliver advantage thus
far [588, 589, 590] a complete investigation of all resources needed to
achieve a computational advantage in near-term quantum devices will likely be
needed to more fully understand the benefits and drawbacks of this approach.

Finally, while this effort seeks to build bridges translating between the
languages and worlds of physics and computer science, the core ideas of
compilation techniques need to evolve further beyond those presented here. We
have demonstrated how techniques such as the Trotter-Suzuki formula and a
linear combination of unitaries can be used to achieve effective compilation
of single qubit-oscillator operations and multi-qubit-oscillator entangling
gates (Sec. V). However, more work is required to develop additional
compilation strategies. Recently, ZX-calculus – a useful graphical tool for
reasoning about DV circuits [591, 592] – was formulated for CV quantum
circuits [593]. An extension to hybrid CV-DV circuits constitutes an
interesting future direction. Separately, as already noted, a complete theory
of non-abelian hybrid QSP does not yet exist, but its formulation would be
very helpful for analytical compilation. More generally, a complete formal
complexity-theoretic analysis of hybrid CV-DV circuit synthesis for tasks such
as state transfer and general unitary synthesis does not yet exist. In App. E
we present cross-compilation techniques among different pairs of ISAs, but
much work remains to be done to improve the efficiency and practicality of
cross-compilation.

At a fine-grain level, compilation techniques capable of identifying long idle
times and inserting SWAPS to memory are crucial. This may be performed in
conjunction with the use of flag qubits [594] or single-shot error mitigation
[595]. Furthermore, hybrid architectures will require transpilation techniques
that map programmatic elements with long idle times to qubits near memories.
Compilers should also be adapted to take advantage of the new sets of native,
low-cost gates afforded by bosonic hardware. More practically, it would be
very useful to create a software library or toolkit built on top of these
analytical and numerical compilation techniques that can perform application-
tailored resource estimation in terms of appropriate metrics such as the gate
count in hybrid ISAs.

Looking further into the future, optimization of quantum compilation using
quantum computers [456] is an exciting possibility. Another potentially
important direction that would apply to algorithm development for both
traditional and hybrid hardware would be the extension of the computation
model for quantum computers beyond algorithms based on purely unitary
operations and final state measurements to more general quantum channel
operations which could be created by including engineered dissipation and/or
mid-circuit measurements and feed-forward.

###### Acknowledgements.

This research was funded by the U.S. Department of Energy, Office of Science,
National Quantum Information Science Research Centers, Co-design Center for
Quantum Advantage (C2QA) under contract number DE-SC0012704 (Basic Energy
Sciences, PNNL FWP 76274). C2QA led this research. Shraddha Singh and Steven
Girvin acknowledge support by the Army Research Office (ARO) under Grant
Number W911NF-23-1-0051. Jasmine Sinanan-Singh was supported in part by the
Army Research Office under the CVQC project W911NF-17-1-0481. Jasmine Sinanan-
Singh and Yuan Liu were supported in part by NTT Research. Yuan Liu
acknowledges startup funding support from North Carolina State University.
Alexander Schuckert acknowledges funding from the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers, and
Quantum Systems Accelerator. John M. Martyn acknowledges support from the
National Science Foundation Graduate Research Fellowship under Grant No.
2141064. Micheline B. Soley was supported by the Yale Quantum Institute
Postdoctoral Fellowship and the Office of the Vice Chancellor for Research and
Graduate Education at the University of Wisconsin-Madison with funding from
the Wisconsin Alumni Research Foundation. The views and conclusions contained
in this document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the Army
Research Office (ARO), or the U.S. Government. The U.S. Government is
authorized to reproduce and distribute reprints for Government purposes
notwithstanding any copyright notation herein. The authors acknowledge helpful
discussions with Jacob Curtis, Michael DeMarco, Shruti Puri, Yongshan Ding,
Robert Schoelkopf, Cindy Regal, Baptiste Royer, Timothy Stavenger, and Ang Li,
and thank Luke Bell for helpful discussions and verification of some of the
key equations presented in the compilation section. External interest
disclosure: SMG receives consulting fees and is an equity holder in Quantum
Circuits, Inc.

##  Appendix A Physical Implementations

Here we present a physical-level discussion of quantum control of hybrid
systems consisting of bosonic modes and qubits, which serves as background
information for the hybrid gates presented in Sec. III.4. Table A.1 summarizes
all qubit-only gates available in superconducting circuits and ion trap
systems and Table A.2 presents the microscopic Hamiltonian terms freely
available in these systems.

|  | Gate Name | Parameters | Definition
---|---|---|---|---
Qubit-only |  Single- | Qubit $x,y$ rotation1 | $\theta\in[0,4\pi),\varphi\in[0,2\pi)$ | $R_{\varphi}\left(\theta\right)=\exp\left(-i\frac{\theta}{2}\sigma_{\varphi}\right)$
|  | $\sigma_{\varphi}=\sigma_{x}\cos\varphi+\sigma_{y}\sin\varphi$
Qubit $z$ rotation1 | $\theta\in[0,4\pi)$ | $R_{z}(\theta)=\exp\left(-i\frac{\theta}{2}\sigma_{z}\right)$
Multi- |  $ZZ$ interaction | $\theta\in[0,2\pi)$ | $ZZ(\theta)=\exp\left(-i\frac{\theta}{2}\sigma_{z}\otimes\sigma_{z}\right)$
CNOT (aka C$X$) |  | $\frac{I+\sigma_{z}}{2}\otimes I+\frac{I-\sigma_{z}}{2}\otimes\sigma_{x}$
CPHASE | $\varphi\in[0,2\pi)$ | $\frac{I+\sigma_{z}}{2}\otimes I+\frac{I-\sigma_{z}}{2}\otimes e^{i\varphi\frac{I-\sigma_{z}}{2}}$
C$Z$ | $\varphi=\pi$ | $\frac{I+\sigma_{z}}{2}\otimes I+\frac{I-\sigma_{z}}{2}\otimes\sigma_{z}$
$i$SWAP | $\theta\in[0,2\pi)$ | $\exp(i\frac{\theta}{2}(\sigma_{x}\sigma_{x}+\sigma_{y}\sigma_{y}))$
| $\theta=\pi/2$ | $\exp(i\frac{\pi}{4}(\sigma_{x}\sigma_{x}+\sigma_{y}\sigma_{y}))$
fSim2 | $\theta\in[0,2\pi),\varphi\in[0,2\pi)$ |  $i$SWAP($-\theta$) $\cdot$ CPHASE($\varphi$)
fSWAP |  | $i\left[R_{z}(\frac{\pi}{2})\otimes R_{z}(\frac{\pi}{2})\right]\textrm{fSim}(\frac{\pi}{2},0)$
Table A.1: Common qubit-only gates. The single-qubit gates are typically natively available. If not natively available, the multi-qubit gates can be readily synthesized through their coupling to the bosonic modes as discussed in Sec. V.3.2. 1In general, qubit rotations about any axis of the Bloch-sphere are given by $R_{\hat{n}}\left(\theta\right)=\exp\left(-i\frac{\theta}{2}\hat{n}\cdot\vec{\sigma}\right)$ where $\vec{\sigma}=\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)$ and $\hat{n}$ is a unit-vector, and this gate can be implemented directly on some physical qubits with x-y-z control. However, for many qubits (such as fixed-frequency transmons), only x-y control is native, and the z rotations are implemented in software [228], hence our choice to separate qubit rotations into either x-y or z rotations. 2fSim (Sec. VI.4.8) is the fermion simulation gate [573], while fSWAP (Sec. VI.4.8) is the fermionic SWAP gate, which swaps the state of two qubits and adds a minus sign to the state $\ket{11}$ to account for antisymmetric fermionic statistics.  Interaction | Hamiltonian | Platform
---|---|---
Dispersive Coupling | $a^{\dagger}a\sigma_{z}$ | Superconducting Circuits
Sideband Interaction | $a\sigma_{+}+a^{\dagger}\sigma_{-}$ | Trapped Ions, Superconducting Circuits, Neutral Atoms
Linear Drive | $\alpha a^{\dagger}-\alpha^{*}a$ | Superconducting Circuits, Trapped Ions
Parametric Drive | $\zeta^{*}a^{2}-\zeta a^{{\dagger}2}$ | Superconducting Circuits, Trapped Ions
Controlled Squeezing | $\sigma_{x}((a^{\dagger})^{2}+a^{2})$ | Trapped Ions, Neutral Atoms
Free Oscillator Evolution | $a^{\dagger}a$ | Trapped Ions, Superconducting Circuits
Controlled Beam-splitter | $\sigma_{x}(a^{\dagger}b+b^{\dagger}a)$ | Trapped Ions
Table A.2: Natively available interactions on various quantum computing
platforms.

###  A.1 Superconducting and Cavity QED

####  A.1.1 Fundamentals

Let us consider a single oscillator coupled to a single transmon qubit via the
Hamiltonian [54, 596, 28, 27]

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{0}+H_{1}+V_{0},$ (617)
$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle\bar{\omega}_{\mathrm{R}}a^{\dagger}a+\bar{\omega}_{\mathrm{Q}}b^{\dagger}b,$
(618) $\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle-\frac{\bar{K}}{2}b^{\dagger}b^{\dagger}bb,$ (619)
$\displaystyle V_{0}$ $\displaystyle=$ $\displaystyle
g(a^{\dagger}+a)(b^{\dagger}+b).$ (620)

Here, $H_{0}$ is the Hamiltonian of the resonator with ‘bare’ (i.e.,
uncoupled) frequency $\bar{\omega}_{\mathrm{R}}$ and the transmon with bare
frequency $\bar{\omega}_{\mathrm{Q}}$ and $H_{1}$ describes the (bare,
negative) anharmonicity $-\bar{K}$ of the transmon. The microscopic
Hamiltonian of the transmon involves a cosine of the phase difference across
the Josephson junction. $H_{1}$ is an approximation to the transmon
Hamiltonian that results from expanding the cosine to fourth order and making
a rotating wave approximation (RWA) [596]. Because of the RWA, $H_{1}$ is
diagonal in the Fock (excitation number) basis and the eigenstates of
$H_{0}+H_{1}$ are unchanged by $H_{1}$. We still have harmonic oscillator
eigenstates, but for the transmon, the energy level spacing is not uniform due
to the presence of the anharmonicity from $H_{1}$. The spectrum of the exact
Hamiltonian is of course bounded from below. The approximate Hamiltonian
containing $H_{1}$ is not bounded from below, implying that it must eventually
break down at large excitation numbers where the truncation of the cosine
potential of the Josephson junction ceases to be a valid approximation. We
will ignore this complication here.

$V_{0}$ is the coupling between the charge polarization (dipole moment) of the
transmon and the cavity electric field. The coupling constant $g$ is known in
the cavity- and circuit-QED literature as the vacuum Rabi coupling [54, 55,
28, 27].

If we assume that the transmon is always within the manifold spanned by its
ground state $|g\rangle$ and its first excited state $|e\rangle$, we can make
the two-level qubit approximation, mapping the transmon operators onto spin
operators. The mapping of the lowest two transmon states onto a qubit spin is
accomplished by

$\displaystyle b^{\dagger}b$ $\displaystyle\rightarrow$
$\displaystyle|e\rangle\langle e|=\frac{\hat{I}+\sigma_{z}}{2},$ (621)
$\displaystyle b$ $\displaystyle\rightarrow$
$\displaystyle\sigma^{-};\,\,\,b^{\dagger}\rightarrow\sigma^{+},$ (622)
$\displaystyle H_{1}$ $\displaystyle\rightarrow$ $\displaystyle 0,$ (623)
$\displaystyle H$ $\displaystyle\rightarrow$
$\displaystyle\omega_{\mathrm{R}}a^{\dagger}a+\frac{\omega_{\mathrm{Q}}}{2}\sigma_{z}+g\left(a+a^{\dagger}\right)\sigma_{x},$
(624)

where we have neglected an irrelevant constant in the energy in the last
expression. Eq. (624) describes the so-called Rabi Hamiltonian.

The Rabi Hamiltonian can be further simplified if we assume that the
cavity–qubit detuning
$\Delta=\bar{\omega}_{\mathrm{R}}-\bar{\omega}_{\mathrm{Q}}$ is small enough
relative to the cavity frequency $\bar{\omega}_{\mathrm{R}}$ that the rotating
wave approximation can be made to further simplify $V_{0}$ to

$V_{0}=g(a^{\dagger}b+ab^{\dagger}),$ (625)

which yields the celebrated Jaynes-Cummings Hamiltonian [54]

$H=\bar{\omega}_{\mathrm{R}}a^{\dagger}a+\frac{\bar{\omega}_{\mathrm{Q}}}{2}\sigma_{z}+g\left(a\sigma^{+}+a^{\dagger}\sigma^{-}\right).$
(626)

This ‘bare’ Jaynes-Cummings Hamiltonian has its parameters fixed at the time
of hardware fabrication and is not programmable, though in some circuit QED
architectures, the transmon frequency $\bar{\omega}_{\mathrm{Q}}$ can be tuned
by application of a magnetic field that threads flux through a SQUID loop.
Typical values for the anharmonicity of the transmon are $\bar{K}/(2\pi)\sim
50-200$ MHz and the vacuum Rabi coupling can be up to $g/(2\pi)\sim 150$ MHz.
Typical values for the operating frequencies of transmons and cavities are in
the range
$\bar{\omega}_{\mathrm{Q}}/(2\pi),\bar{\omega}_{\mathrm{R}}/(2\pi)\sim 5-10$
GHz.

An alternate approach to the Hamiltonian of a qubit coupled to a cavity that
is especially useful when higher levels of the transmon come into play, or
when the cavity and transmon are detuned from each other (‘dispersive regime’)
will be described below. In the dispersive regime, we will find that it is
possible to apply continuous microwave drives to the qubit and the cavity
which will yield an effective Jaynes-Cummings model with adjustable
parameters.

$H_{0}+V_{0}$ is the harmonic part of the total Hamiltonian and is of course
exactly soluble. We will use the strategy of expressing the anharmonic term,
$H_{1}$ based on exact normal modes of the harmonic part of the Hamiltonian.
For simplicity, we will make the rotating wave approximation for $V_{0}$ so
that the harmonic Hamiltonian can be written

$\displaystyle H_{0}+V_{0}$ $\displaystyle=$
$\displaystyle\begin{array}[]{c}\left(\begin{array}[]{cc}a^{\dagger}&b^{\dagger}\end{array}\right)\\\
\phantom{a}\end{array}\left(\begin{array}[]{cc}M_{11}&M_{12}\\\
M_{21}&M_{22}\end{array}\right)\left(\begin{array}[]{c}a\\\
b\end{array}\right),$ (634) (635)

where

$M=\left(\begin{array}[]{cc}\bar{\omega}_{\mathrm{R}}&g\\\
g&\bar{\omega}_{\mathrm{Q}}\end{array}\right)=\frac{\bar{\omega}_{\mathrm{R}}+\bar{\omega}_{\mathrm{Q}}}{2}\hat{I}+\left(\begin{array}[]{cc}\frac{\Delta}{2}&g\\\
g&-\frac{\Delta}{2}\end{array}\right).$ (636)

From this, we see that the two normal modes of the coupled oscillators have
frequencies

$\omega_{\pm}=\frac{\bar{\omega}_{\mathrm{R}}+\bar{\omega}_{\mathrm{Q}}}{2}\pm\sqrt{(\Delta/2)^{2}+g^{2}},$
(637)

and the normal mode operators are

$\displaystyle\left(\begin{array}[]{c}A\\\
B\end{array}\right)=\left(\begin{array}[]{cc}\cos\frac{\theta}{2}&+\sin\frac{\theta}{2}\\\
-\sin\frac{\theta}{2}&\cos\frac{\theta}{2}\end{array}\right)\left(\begin{array}[]{c}a\\\
b\end{array}\right),$ (644)

where $\tan\theta=\frac{2g}{\Delta}$, and

$\displaystyle\sin\frac{\theta}{2}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1-(\Delta/2)/\sqrt{g^{2}+(\Delta/2)^{2}}}{2}},$
(645) $\displaystyle\cos\frac{\theta}{2}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1+(\Delta/2)/\sqrt{g^{2}+(\Delta/2)^{2}}}{2}}.$
(646)

To express the anharmonic term $H_{1}$ in the normal-mode basis, it is
convenient to have the inverse transform

$\displaystyle\left(\begin{array}[]{c}a\\\
b\end{array}\right)=\left(\begin{array}[]{cc}\cos\frac{\theta}{2}&-\sin\frac{\theta}{2}\\\
+\sin\frac{\theta}{2}&\cos\frac{\theta}{2}\end{array}\right)\left(\begin{array}[]{c}A\\\
B\end{array}\right).$ (653)

Making the rotating wave approximation we have

$\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle-\frac{\bar{K}}{2}b^{\dagger}b^{\dagger}bb,$
$\displaystyle\approx$
$\displaystyle-\big{[}\chi_{BB}B^{\dagger}B^{\dagger}BB+\chi_{AB}B^{\dagger}BA^{\dagger}A$
$\displaystyle~{}~{}~{}~{}~{}~{}+\chi_{AA}A^{\dagger}A^{\dagger}AA\big{]},$
(655)

where the renormalized transmon self-Kerr (anharmonicity), the transmon-cavity
cross-Kerr (dispersive interaction), and the self-Kerr induced in the cavity
are respectively

$\displaystyle\chi_{BB}$ $\displaystyle=$
$\displaystyle\frac{K}{2}\cos^{4}(\theta/2),$ (656) $\displaystyle\chi_{AB}$
$\displaystyle=$ $\displaystyle 2K\cos^{2}(\theta/2)\sin^{2}(\theta/2),$ (657)
$\displaystyle\chi_{AA}$ $\displaystyle=$
$\displaystyle\frac{K}{2}\sin^{4}(\theta/2).$ (658)

In the so-called dispersive regime, $g\ll|\Delta|$, $\theta\ll 1$ and hence
$\chi_{BB}\gg\chi_{AB}\gg\chi_{AA}.$ In the so-called ‘strong-dispersive
regime,’ the cross-Kerr coupling is much larger than $\kappa$, the energy
dissipation rate of the cavity and $\gamma$, the energy relaxation rate of the
qubit $\chi_{AB}\gg\kappa,\gamma$ [597, 127].

If we now make the two-level qubit approximation for the transmon using the
analog of Eqs. (621-622), we obtain the Hamiltonian (up to an irrelevant
constant)

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{\mathrm{D}}+V,$ (659)
$\displaystyle H_{\mathrm{D}}$ $\displaystyle=$
$\displaystyle\omega_{\mathrm{R}}A^{\dagger}A+\frac{\omega_{\mathrm{Q}}}{2}\sigma_{z}+\frac{\chi}{2}\sigma_{z}A^{\dagger}A,$
(660) $\displaystyle V$ $\displaystyle=$
$\displaystyle-\frac{K_{A}}{2}A^{\dagger}A^{\dagger}AA,$ (661)

where $\chi\equiv\chi_{\mathrm{AB}}$ is the dispersive shift and $K_{A}\equiv
2\chi_{\mathrm{AA}}$ represents the weak negative anharmonicity inherited by
the cavity normal mode due to the coupling to the transmons. Here, the
renormalized cavity frequency is given by Eq. (637), with
$\omega_{\mathrm{R}}=\omega_{+}$ if the detuning $\Delta>0$, and
$\omega_{\mathrm{R}}=\omega_{-}$ if $\Delta<0$. The renormalized qubit
frequency is the remaining eigenvalue, $\omega_{\mathrm{Q}}=\omega_{-}$ if the
detuning $\Delta>0$, and $\omega_{\mathrm{Q}}=\omega_{+}$ if $\Delta<0$.

Because $K_{A}\ll\chi$, we will make the simplifying assumption that we can
neglect the induced cavity self-Kerr term $V$, and work only with the
dispersive Hamiltonian $H_{\mathrm{D}}$. The dispersive Hamiltonian $V_{0}$
can be interpreted two different ways. First, it tells us that the resonant
frequency of the cavity takes on different values,
$\omega_{\mathrm{R}}\pm\frac{\chi}{2}$, depending on the state of the qubit.
The frequency of the cavity is thus a proxy that can be used to make a high-
fidelity QND dispersive measurement of the state of the qubit [54, 55].

In the so-called ‘strong-dispersive regime,’ [597, 127] where $\chi$ greatly
exceeds the damping rate of the cavity and the transmon, the two resonance
peaks of the cavity are spectrally well-resolved. We will see shortly that
this will allow us to make operations on the cavity that are conditioned on
the qubit state, thus creating entanglement between the cavity and qubit.

The second interpretation of the dispersive Hamiltonian is that the qubit
frequency shifts by $\chi$ for each photon that is added to the cavity. We
will see shortly that this will allow us to make rotations of the qubit
conditioned on the precise number of photons in the cavity, giving us a second
way to entangle the qubit and cavity. In particular, it gives us the ability
to apply selective number-dependent arbitrary phase (SNAP, Sec. III.4.4) gates
[168] to each Fock state of the cavity.

Together, these two resources, qubit-conditioned cavity displacements and
cavity-conditioned qubit rotations, grant us full quantum control of the
combined system [44, 168, 233]. While some gates require only simple and
intuitive microwave control pulses, the execution of more general unitaries at
the fastest possible rates requires numerical determination of the control
pulses, using for example the GRAPE algorithm from optimal control theory
(OCT), which are difficult to interpret. For example, Heeres et al. [233]
carried out a highly non-Gaussian operation that executed a state transfer
from the $n=0$ Fock state of the cavity to the $n=6$ Fock state, using OCT
pulses applied simultaneously to the transmon and to the cavity. Efficiency
gains from bypassing qubit and gate abstractions to directly compile quantum
algorithm steps at the control pulse level have been discussed recently by Shi
et al. [598].

####  A.1.2 Conditional Cavity Gates

Conditional displacements. Qubit-controlled cavity displacement operations
have been previously discussed in [231, 49, 45, 599]. In this subsection, we
aim to sketch at a high level the implementation of the conditional
displacement gate $\mathrm{CD}(\alpha)$, an integral gate to the phase-space
ISA defined in Box III.7.

In Ref. [49, 45, 599], a clever echo sequence is used to make rapid
conditional displacements with a strong drive and simultaneously eliminate
unwanted Stark shift rotations of the ancilla. Adding a drive to the cavity
centered on frequency $\omega_{\mathrm{R}}$ and having an envelope
$\epsilon(t)$, the Hamiltonian in Eq. (660) in the frame doubly rotating at
the qubit and cavity frequencies becomes

$H=\frac{\chi}{2}\sigma_{z}A^{\dagger}A+i\epsilon(t)A^{\dagger}-i\epsilon^{*}(t)A.$
(662)

For the case where the complex drive phase is time-independent, we can remove
the drive using a time-dependent unitary transformation

$U(t)=e^{-\beta(t)A^{\dagger}+\beta^{*}(t)A}$ (663)

that shifts the oscillator to a displaced frame

$A\rightarrow A+\beta(t),$ (664)

where

$\beta(t)=\int_{0}^{t}d\tau\,\epsilon(\tau).$ (665)

The Hamiltonian in the new frame becomes

$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{\chi}{2}\sigma_{z}[A^{\dagger}+\beta^{*}(t)][A+\beta(t)]$
$\displaystyle=$
$\displaystyle\frac{\chi}{2}\sigma_{z}\left(\rule{0.0pt}{10.33327pt}\left[A^{\dagger}A\right]+\left[\beta^{*}(t)A+\beta(t)A^{\dagger}\right]+\left[|\beta(t)|^{2}\right]\right).$
(667)

The middle term is the generator of the desired conditional displacement. The
first and third terms are spurious but can be eliminated by a ‘double echo’
sequence in which we first drive $\beta$ to a large positive value, then flip
the ancilla state and drive $\beta$ to a large negative value. The details are
described in the Supplementary Material of Ref. [45], where the authors call
the sequence an _echoed-conditional-displacement (ECD) gate_. For the
appropriate choice of $\beta(t)$, this gate is equivalent (up to a bit-flip on
the ancilla) to the conditional displacement $\textrm{CD}(\alpha)$ defined in
Box III.7 and used throughout this work.

Conditional oscillator-oscillator entangling gates in the weak dispersive
regime. In Section V.3.1, we demonstrated a useful trick for compiling
conditional oscillator-oscillator entangling gates using controlled parity
gates (see the controlled beam-splitter gate in Eq. (395) as an example).
Importantly, controlled parity gates are realized by allowing the dispersive
interaction $e^{-i\chi ta^{\dagger}a}$ to act for a time $t=\pi/2\chi$. In the
low dispersive regime (where dispersive coupling $\chi$ is extremely small),
however, this strategy is undesirable as it results in a long gate duration.
Here, we show an alternative method to compile a conditional oscillator-
oscillator entangling gate – namely, a conditional two-mode squeezing – in the
weak dispersive regime. To achieve this, we leverage uncontrolled two-mode
squeezing operations, ${\rm TMS}(r,\phi)$, as defined in Table
LABEL:tab:gates-osc.

Using the definition of ${\rm TMS}(r,\phi)$, we first note that

$\displaystyle{\rm TMS}(\alpha,\pi)~{}a~{}{\rm
TMS}^{\dagger}(\alpha,\pi)=a\cosh{\alpha}+b^{\dagger}\sinh{\alpha}\,.$ (668)

This in turn implies

$\displaystyle{\rm TMS}(\alpha,\pi)e^{-i\chi ta^{\dagger}a\sigma_{z}}{\rm
TMS}^{\dagger}(\alpha,\pi)$
$\displaystyle=\mathrm{exp}\left[\rule{0.0pt}{10.33327pt}-i\chi
t\left(\rule{0.0pt}{10.33327pt}\cosh^{2}{(\alpha)}a^{\dagger}a+\sinh^{2}{(\alpha)}bb^{\dagger}\right.\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.\left.+\frac{1}{2}\sinh{(2\alpha)}(a^{\dagger}b^{\dagger}+ab)\right)\sigma_{z}\rule{0.0pt}{10.33327pt}\right]$
(669)

The next steps of the method are inspired by the construction of echoed-
conditional displacements (ECD) described above. Notice that the first two
terms in Eq. (669) will not change signs with $\alpha$, whereas the last two
will. Hence, running the pulse shown below yields an echoed two-mode
squeezing, in close analogy with the echoed displacement gate previously
discussed,

$\displaystyle{\rm TMS}(\alpha,\pi)e^{-i\chi ta^{\dagger}a\sigma_{Z}}{\rm
TMS}^{\dagger}(\alpha,\pi)\times\sigma_{x}$ $\displaystyle\times{\rm
TMS}(-\alpha,\pi)e^{-i\chi ta^{\dagger}a\sigma_{Z}}{\rm
TMS}^{\dagger}(-\alpha,\pi)$ $\displaystyle~{}~{}~{}~{}~{}=e^{-i\chi
t\sinh{2\alpha}(a^{\dagger}b^{\dagger}+ab)\sigma_{z}}$ (670)

The speed of the conditional two-mode squeezing gate, in this case, is decided
by $\chi t\sinh{2\alpha}$ instead of $\chi t$ as in the case where this gate
is synthesized with controlled parity gates. As $\sinh{2\alpha}$ is an
unbounded function, we can in principle increase it to extremely large values
by varying $\alpha$. Thus, in the weak dispersive regime, we can achieve fast
conditional oscillator-oscillator entangling gates by leveraging unconditional
two-mode squeezing with large $\alpha$ as a resource. Hamiltonian
amplification techniques [600] have been used in ion traps as well, see Ref.
[601, 602].

####  A.1.3 Cavity-controlled Qubit Gates

SQR Gates. A microwave electric field applied to the qubit couples to the
electric dipole moment of the qubit and causes transitions between the ground
and excited state. If the drive is applied at the qubit frequency it will
coherently rotate the qubit around an axis that is determined by the phase of
the microwave drive. Notice from the dispersive Hamiltonian $H_{\mathrm{D}}$
for the coupled qubit-cavity system in Eq. (660) that the qubit transition
frequency suffers a quantized light shift depending on the photon number in
the cavity

$\omega_{Q}+\chi\hat{n}.$ (671)

Hence we will be able to execute a qubit gate conditioned on the photon number
in the cavity. If we apply the qubit drive at frequency
$\omega_{Q}(m)=\omega_{Q}+m\chi$, it will be on resonance and therefore
effective at rotating the qubit on the Bloch sphere, if and only if, there are
exactly $m$ photons in the cavity (assuming $\chi$ is sufficiently large
relative to the drive amplitude). We will now derive this in detail.

Starting from the driven-dispersive Hamiltonian in the co-rotating frame of
the oscillator and qubit,

$\displaystyle H=\chi
a^{\dagger}a\sigma_{z}+\Omega^{*}(t)\sigma_{-}+\Omega(t)\sigma_{+},$ (672)

we perform a unitary frame transformation to eliminate the dynamics of the
dispersive term by using the unitary $U_{\chi}=\exp\left(-i\chi
ta^{\dagger}a\sigma_{z}/2\right)$. The transformed Hamiltonian
$\tilde{H}=U_{\chi}^{\dagger}HU_{\chi}+i\left(\partial_{t}U_{\chi}^{\dagger}\right)U_{\chi}$
is

$\displaystyle\tilde{H}=\sum_{n=0}^{\infty}\ket{n}\bra{n}\left(\Omega^{*}(t)\sigma_{-}e^{i\chi
nt}+\Omega(t)\sigma_{+}e^{-i\chi nt}\right).$ (673)

For selective qubit rotation gates up to Fock state $N$, we drive the qubit
resonant with each number-shifted transition. This drive takes the form
$\Omega(t)=\sum_{k=0}^{N}\Omega_{k}(t)e^{ik\chi t}$ where $\Omega_{k}(t)$ are
complex-valued envelopes. With this drive, the Hamiltonian becomes

$\displaystyle\tilde{H}$ $\displaystyle=\sum_{n=0}^{\infty}\ket{n}\bra{n}$
$\displaystyle\times\sum_{k=0}^{N}\left(\Omega^{*}_{k}(t)\sigma_{-}e^{i\chi\left(n-k\right)t}+\Omega_{k}(t)\sigma_{+}e^{-i\chi\left(n-k\right)t}\right).$
(674)

Under the condition of slow pulses (total time $T_{\text{gate}}>2\pi/\chi$)
$\Omega_{k}(t)$ will have a small bandwidth compared to $\chi$, and the first-
order rotating-wave-approximation (RWA) can be invoked. The implications of
this approximation, along with mitigation techniques for faster pulses, are
discussed in [449]. The RWA results in

$\displaystyle\tilde{H}_{\text{RWA}}=\sum_{n=0}^{N}\ket{n}\bra{n}\left(\Omega^{*}_{n}(t)\sigma_{-}+\Omega_{n}(t)\sigma_{+}\right).$
(675)

To realize an SQR gate, we pick $\Omega_{n}(t)=s(t)\theta_{n}e^{i\varphi_{n}}$
where $s(t)$ is a real-valued envelope function with
$\int_{0}^{T_{\text{gate}}}s(\tau)d\tau=1$. By transforming the operator back
to the initial frame, the total unitary after time $T_{\text{gate}}$ is

$\displaystyle
U=U_{\chi}\,\mathcal{T}\exp{\left(-i\int_{0}^{T_{\text{gate}}}\tilde{H}_{\text{RWA}}dt\right)}U_{\chi}^{\dagger}$
(676) $\displaystyle=\text{CR}\left(\chi
T_{\text{gate}}\right)\text{SQR}\left(\vec{\theta},\vec{\varphi}\right)\text{CR}\left(-\chi
T_{\text{gate}}\right).$ (677)

As indicated, the full gate includes additional conditional rotations, which
account for the cavity’s pull on the qubit during the gate. This CR can be
accounted for by either applying a counter CR gate before and after the
evolution, by constructing a suitable qubit echo sequence similar to the ECD
construction or, in cases where the initial and final qubit state are known
(such as when implementing the SNAP gate), can be tracked in software.

SNAP Gates. As described in Box III.10, two SQR gates can be combined to
produce the photon-number selective phase (SNAP) gate. One interesting choice
of phases is

$\varphi_{m}=m\varphi$ (678)

which simply rotates the cavity state in phase space by angle $-\varphi$.
Since the Berry phase is equal and opposite if the qubit starts and ends in
the excited state, one can achieve the qubit-dependent cavity rotation gate
defined in Box III.8

$\mathrm{CR}(\varphi)=e^{-i\varphi\hat{n}\sigma_{z}},$ (679)

which can be used to either enhance or compensate for the naturally occurring
time evolution under the dispersive Hamiltonian in Eq. (660).

For quantum simulation and control purposes, another interesting choice of
phases is

$\theta_{m}=\lambda m^{2},$ (680)

which yields

$U_{\mathrm{Kerr}}(\lambda)=e^{i\lambda\hat{n}^{2}\sigma_{z}},$ (681)

which can be used to dynamically cancel the undesired residual Kerr
interaction in the cavity found in Eq. (661). This gate can also be used to
enhance the cavity self-Kerr to provide a programmable interaction for
simulations of the boson Hubbard model as discussed in Sec. VI.4.5.

###  A.2 Trapped Ion Quantum Processors

In ion traps [15], one can produce the phase-space instructions set using the
interaction between ions and lasers or microwaves. As discussed in Fig. 1, the
mechanical oscillation modes of the ions form the bosonic modes. Unlike in the
superconducting case where the microwave resonators typically have GHz
frequencies, the mechanical modes in ion traps typically have frequencies in
the MHz range. Beyond the various gates we discuss below, multi-chromatic
laser drives that simultaneously couple many qubits to many mechanical modes
have been developed [603].

####  A.2.1 Qubit rotations

The qubits in ion-trap systems can be encoded in optical transitions of the
ions, but can also be encoded in hyperfine levels resulting from the coupling
between electron and nuclear spins in the ground state manifold of the ions.
The hyperfine energy level spacing is usually on the order of MHz and the
transitions can be driven using microwaves to achieve single-qubit rotations
much as in NMR experiments [604]. Given their long wavelength, microwaves can
be difficult to focus to address individual ions, but by engineering near-
field microwave gradients, individual addressability can be achieved [605]. An
alternative way to drive the qubit hyperfine transition is to use two detuned
laser beams which have a frequency difference that matches the hyperfine
splitting to drive the qubit transitions via a stimulated Raman process
(STIRAP). Because the stimulated Raman transition is 2nd-order (i.e., requires
one-photon absorption and one emission), it generally requires relatively
strong laser intensities.

####  A.2.2 Sideband gates

A sideband gate $\mathrm{JC}(\theta,\varphi)$ (see Table LABEL:tab:gates-
qubit-osc ) that flips one of the qubits and emits or absorbs one mechanical
quantum (phonon) is (for the case of a red sideband)

$\displaystyle\text{JC}(\theta,\varphi)=\exp\left[-i\theta\left(e^{i\varphi}\sigma_{-}a^{\dagger}+e^{-i\varphi}\sigma_{+}a\right)\right],$
(682)

which is universal when combined with an arbitrary single-qubit rotation
$R_{\varphi}(\theta)$ in the XY plane (see Table A.1) [47, 46, 74]. This gives
the ISA listed in Table A.3.

Similarly, a blue sideband gate $\text{AJC}(\theta,\varphi)$ (Table
LABEL:tab:gates-qubit-osc) could be used, and when combined with a red
sideband gate produces the controlled displacement gate $\mathrm{CD}(\alpha)$
of the phase-space ISA, as discussed below. We note that by use of a pulse
shaping technique, a uniform sideband gate has been realized wherein the Rabi
frequency between nearest neighbor Fock levels does not scale as the Fock
number but is a constant [248].

In the ion-trap setting, the finite momentum kick from the absorption or
emission of an optical photon allows for the possibility of higher-order
sideband transitions in which a spin flip is coupled to the creation or
destruction of two or more mechanical oscillation quanta (phonons). While such
higher sideband Hamiltonians are weaker than the lowest sideband coupling by
powers of the Lamb-Dicke parameter, they present interesting possibilities for
additional control schemes and instruction sets. We will however leave these
for future work and not discuss them further here.

We are not aware of a SNAP gate realization on trapped ion platforms because
of the lack of dispersive coupling between the qubits and the mechanical modes
of the system. Such a coupling can however be synthesized from sideband
transitions [74]. Beam-splitters can also be synthesized outside this sideband
gate set without using any lasers by parametrically modulating the electric
fields used to create the trapping potential seen by the ions [76, 74].
Alternatively, the second sidebands include beam-splitter interactions, which
are however weaker by an additional power of the Lamb-Dicke factor. Moreover,
arbitrary all-to-all connected beam-splitter Hamiltonians can be implemented
using multi-tone driving [57].

Trapped Ion ISA
---
Red sideband gate (Table LABEL:tab:gates-qubit-osc) | $\text{JC}(\theta,\varphi)$
Qubit rotation (Table A.1) | $R_{\varphi}(\theta)$
Table A.3: Minimal instruction set for a trapped ion platform. This is the
sideband ISA defined in Table III.4 with the important difference that any
qubit can have JC coupling with any motional mode in which it participates.
Application of the sideband gates to different qubits allows them to be
entangled through the exchange of quanta of one of the bosonic modes. A beam
splitter gate (see Box III.4) to couple different mechanical modes, $a$ and
$b$, can be synthesized from this gate set or can be created using RF control
of the trap parameters as noted in the caption of Fig. 1. Blue sideband gates
(AJC, Table LABEL:tab:gates-qubit-osc) are equally easy to perform and are a
convenient natively available gate that can be added to this minimal set.

####  A.2.3 Conditional Displacements

Recall that the red side-band produces coupling Hamiltonian in the form of
$\sigma^{+}a+\sigma^{-}a^{\dagger}$, whereas a blue side-band produces
coupling in the form of $\sigma^{-}a+\sigma^{+}a^{\dagger}$. If a bichromatic
laser drive [232] is applied so that both sideband drives are switched on
simultaneously and with equal strengths, we can realize the following
Hamiltonian as a sum of the two:

$\displaystyle(\sigma^{+}a+\sigma^{-}a^{\dagger})+(\sigma^{-}a+\sigma^{+}a^{\dagger})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle~{}~{}~{}=\sigma^{+}(a+a^{\dagger})+\sigma^{-}(a+a^{\dagger})~{}~{}\propto~{}~{}\sigma_{x}\hat{x}\,.$
(683)

This is nothing but a conditional displacement Hamiltonian, conditioned on the
eigenstates of the qubit Pauli-X operator. We may further recast this into the
familiar controlled displacement, $\textrm{CD}(\alpha)$ (Table
LABEL:tab:gates-qubit-osc), by conjugating the qubit using a Hadamard gate
since $\sigma_{z}=H\sigma_{x}H$.

Unconditional displacements can be achieved by using one of the qubits as an
ancilla that is not involved in the circuit. Initializing this spare qubit in
$\ket{+}$ and driving a conditional displacement realizes an unconditional
displacement.

###  A.3 Neutral Atom Quantum Processors

The motion of single atoms in optical tweezers [7, 606, 254] is also
approximately that of a quantum harmonic oscillator. In contrast to ions, each
atom is coupled to a maximum of three modes, where in practice one of those
modes is coupled to the atomic qubit, as illustrated in Fig. 1. This means
that there is no mode-crowding effect restricting the scaling to larger
systems, but this comes at the cost of reduced connectivity among the modes.
Moreover, due to the much weaker dipole force confining the atoms compared to
the electrostatic force confining ions, the motional mode frequency
$\omega_{r}$ is smaller, usually in the tens of kHz. Qubits may be encoded in
hyperfine levels (in alkali atoms such as rubidium) or optical clock qubits
(in strontium). By qubit-state-selective of the atoms into Rydberg states, the
Rydberg blockade effect can be used to create two-qubit entangling gates
[607]. Recently, combined qubit-boson operations have been demonstrated in a
tweezer array [72], showing coherence times of the modes over $100$ms and non-
classical motion of the atoms due to Fock-state preparation has been shown
[608].

Contrary to ions, the tweezer potential has a sizeable nonlinearity. The
maximum boson number occupation is given by the trap depth $V_{0}$, which in
current implementations is on the order of 1-10 MHz such that the maximum
accessible Fock state is $100$-$1000$. Another effect of the finite trap depth
is a finite Kerr nonlinearity of the potential, which is $0.1\%$-$1\%$ of the
mode frequency (coincidentally similar to superconducting transmon qubits).
The nonlinearity is in principle tunable by the power of the trapping laser.
However, despite this sizeable nonlinearity, qubit-boson gates are not
substantially altered since the qubit-boson coupling is so large.

Neutral atom ISA
---
Red sideband gate (Table LABEL:tab:gates-qubit-osc) | $\text{JC}(\theta,\varphi)$
Qubit rotation (Table A.1) | $R_{\varphi}(\theta)$
Qubit-qubit gate (Table A.1) | CZ
Table A.4: Minimal instruction set for neutral atom platforms. Sideband gates
couple a single qubit to a single mode. Beam splitter gates (see Box III.4) to
a couple of different mechanical modes, $a$ and $b$, can be synthesized via
the qubit-qubit Rydberg CZ gate.

####  A.3.1 Single qubit gates

Single-qubit gates are implemented just as in qubit-only applications of
neutral atoms by either Raman transitions [9] or microwaves combined with AC
Stark shifts [609].

####  A.3.2 Multi qubit gates

Multi-qubit gates are implemented by using a Rydberg level of the atoms. High
fidelity two-qubit CZ and even n-qubit multi-controlled Z gates have been
demonstrated [607] by using the blockade mechanism, which forbids two atoms
(within the blockade radius of each other) to be transferred to the Rydberg
state simultaneously [6]. When two atoms are simultaneously in the Rydberg
state, they produce forces that modify the bosonic state, as discussed in the
next subsection. However, during the Rydberg blockade entangling gates, the
atoms do not simultaneously occupy Rydberg states, thus preventing any impact
on the motional states.

####  A.3.3 Single mode gates

Phase space rotations are naturally implemented by the free evolution under
the tweezer potential. Moreover, the nonlinearity enables universal single-
mode control by modulating the tweezer potential [610].

####  A.3.4 Multi mode gates

Using two-qubit gates. Inter-mode gates such as beam-splitters can be compiled
using sideband gates and CZ gates. We employ the BCH formula

$e^{i\hat{A}\theta}e^{i\hat{B}\theta}e^{-i\hat{A}\theta}e^{-i\hat{B}\theta}=e^{-[\hat{A},\hat{B}]\theta^{2}+O(\theta^{3})}.$
(684)

First, we create a doubly-controlled displacement by choosing
$\hat{A}=\hat{Z}_{1}\hat{Z}_{2}$ and
$\hat{B}=\hat{X}_{1}(\hat{a}+\hat{a}^{\dagger})$. This way we get
$-[\hat{A},\hat{B}]\theta^{2}=2i\hat{Y}_{1}\hat{Z}_{2}(\hat{a}+\hat{a}^{\dagger})\theta^{2}\equiv
i\hat{C}\tilde{\theta}$ with
$\hat{C}=\hat{Y}_{1}\hat{Z}_{2}(\hat{a}+\hat{a}^{\dagger})$ and
$\tilde{\theta}=\sqrt{2}\theta$. Then, we can repeat the BCH trick,
conjugating $\hat{C}$ with $\hat{D}=\hat{X}_{2}(\hat{a}+\hat{a}^{\dagger})$.
From this concatenation, we obtain

$e^{-[\hat{C},D]\tilde{\theta}^{2}}=e^{i\hat{Y}_{1}\hat{Y}_{2}(\hat{a}+\hat{a}^{\dagger})(\hat{b}+\hat{b}^{\dagger})4\tilde{\theta}^{4}}.$
(685)

Hence, we have realized a doubly qubit-controlled beam-splitter/two-mode-
squeezing gate. The two-mode-squeezing or beam-splitter terms can be singled
out by applying the same sequence with different phases and employing the
Trotter formula.

Using dipolar forces in Rydberg state. Two neutral atoms excited to Rydberg
states exert strong dipolar forces on each other. By “dressing” the ground
state with the Rydberg state and modulating the tweezer potentials, beam-
splitters as well as two-mode squeezing can be realized [611].

Using tunneling. Another approach is to use tunneling between tweezers,
generalizing the approach discussed in the context of qubits encoded in the
lowest two levels of each tweezer [612].

####  A.3.5 Sideband gates

Qubit-boson coupling is implemented using sideband gates as in ions, which in
the Lamb-Dicke regime $\eta\sqrt{n}\ll 1$ is given by

$\hat{H}=\frac{1}{2}\Omega\eta\left(\hat{a}^{\dagger}\hat{\sigma}^{-}+\hat{a}\hat{\sigma}^{+}\right)$
(686)

with $\Omega$ the Rabi frequency, and the Lamb-Dicke factor $\eta$ is given in
terms of the momentum transfer $\delta k$ of the Raman beams and the atomic
mass $m$ by $\eta=\delta k/\sqrt{2m\omega_{r}}$. Due to the essentially
unlimited Rabi frequency in this context, the main restriction to the coupling
$\Omega\eta$ of these gates is given by the boson mode frequency - to resolve
the sideband transitions it needs to fulfill $\Omega\eta<\omega_{r}$.

####  A.3.6 Alternative architectures

Motional states in optical lattices. Alternatively, to the tweezer platform,
the motional modes in optical lattices can be similarly used. Lattice shaking
can then be used to excite oscillators, which may also be done in a qubit-
dependent way. Superlattices can be used to engineer tunneling interactions
for beam-splitters.

Bosonic atoms. A completely different approach is to use bosonic atoms in
optical lattices or tweezers, similar to the recent proposal for fermion-qubit
digital computation [36]. Recently, boson sampling has been realized this way
[613]. This approach is most naturally suited for problems conserving the
total number of bosons. However, displacements may be engineered by coupling
the tweezers to a Bose-Einstein-Condensate.

##  Appendix B Noise Models

We discuss a few common noise models for bosonic systems in different
experimental platforms in this section.

###  B.1 Kraus Operator Representation of Error Channel