# Efficient quantum simulation of nonlinear interactions using SNAP and Rabi
gates

Kimin Park, Petr Marek and Radim Filip  Department of Optics, Palacky
University, 77146 Olomouc, Czech Republic

###### Abstract

Quantum simulations provide means to probe challenging problems within
controllable quantum systems. However, implementing or simulating deep-strong
nonlinear couplings between bosonic oscillators on physical platforms remains
a challenge. We present a deterministic simulation technique that efficiently
and accurately models nonlinear bosonic dynamics. This technique alternates
between tunable Rabi and SNAP gates, both of which are available on
experimental platforms such as trapped ions and superconducting circuits. Our
proposed simulation method facilitates high-fidelity modeling of phenomena
that emerge from higher-order bosonic interactions, with an exponential
reduction in resource usage compared to other techniques. We demonstrate the
potential of our technique by accurately reproducing key phenomena and other
distinctive characteristics of ideal nonlinear optomechanical systems. Our
technique serves as a valuable tool for simulating complex quantum
interactions, simultaneously paving the way for new capabilities in quantum
computing through the use of hybrid qubit-oscillator systems.

##  I Introduction

Advances in the strong coupling regime of qubit-oscillator interactions in
trapped ions and superconducting circuits [1, 2, 3] has opened new revenues
for building qubit gates [4, 5, 6] and, mainly, controllable gates of the
oscillators [7, 8, 9] and first nonlinear gates [10, 11, 12]. Trapped ions
[13, 14, 15] and superconducting circuit QED with transmons [16, 17, 18] are
highly promising platforms for qubit-oscillator quantum information
processing. The fundamental qubit-oscillator interaction is the resonant Rabi
coupling [1, 2], while dispersive interaction appears far from resonance [19,
16]. Strong Rabi couplings have been experimentally demonstrated in cavity QED
systems, including trapped ions [20] and superconducting circuits [21]. These
enable diverse non-Gaussian operations such as nonlinear phase gates [22, 23],
which are essential for CV quantum information processing [24]. Both trapped
ion and superconducting circuits can dynamically switch between interaction
types. In particular, superconducting circuits show promise for implementing
tunable dispersive and resonant qubit-oscillator couplings [19, 12, 25, 26,
27, 28], whereas achieving dispersive coupling remains challenging in trapped
ion systems [29]. In this regime beyond the rotating-wave approximation,
theory suggests that we can achieve universal control of the oscillator using
sequences of Rabi gates [23] and dispersive gates [30]. Such methodology
offers quantum simulations of nonlinear effects with general bosonic systems
[31, 32]. Recently, Kudra et al. [33] experimentally generated high-fidelity
Wigner-negative states, including Schrödinger-cat states and Gottesman-Kitaev-
Preskill states, using novel SNAP gates in 3D microwave cavities. It is
intriguing to consider whether such a sequence of pulsed gates can construct a
variety of nonlinear interactions between oscillators and efficiently simulate
more complex hybrid nonlinear phenomena in quantum physics.

A prototypical example is the strong nonlinear radiation pressure interaction
between a movable mirror in a Fabry-Perot cavity and the intracavity field
[34, 35] in nonlinear optomechanics. This interaction is fundamentally
nonlinear, involving third power of ladder operators in the Hamiltonian [34].
At the quantum level, such a resulting pressure-like interaction enables, in
principle, challenging single-photon-scale displacements of the mirror.
Exploring this highly nonlinear coupling, which can generate nonclassical
superposition states from phase-insensitive noise [36, 37, 38, 39], will
substantially advance current quantum technologies [40, 41, 42, 43, 44, 45,
46, 47] and therefore, it is attractive for hybrid quantum simulations [48].
Nonlinear pressure-like phenomena are predicted to enable breakthrough
applications in quantum information processing, metrology, communication,
fundamental tests of quantum mechanics [49, 50, 51].

However, the inherent weakness of nonlinear radiation pressure interactions
presents a challenge in directly observing quantum nonlinear phenomena without
the introduction of additional noise. Although strong optical pumping can
enhance effective coupling strengths, it inevitably leads to linearization of
the coupling [34], removing the underlying nonlinear dynamics. Experimental
strategies like amplification have also been employed to access the single-
photon scale nonlinear optomechanical phenomena [52, 53, 54, 55, 56, 57, 58].
The membrane-in-the-middle configuration provides enhanced nonlinear
optomechanical interactions [59], making it an attractive resource for quantum
simulation of the hybrid system. Nevertheless, proof-of-principle
demonstrations and analysis of these higher-order nonlinearities remain
challenging, and finding new routes to harness these interactions has been an
attractive prospect for quantum technologies.

An effective strategy for exploring deep-strong and nonlinear quantum
phenomena involves the use of quantum simulation based on engineered
interactions [60]. There are several methods available for engineering
continuous-variable gates [61]. One such method involves simulating a
conditional form of optomechanical interaction [62] using photonic ancillas
and homodyne detectors [63]. Despite its potential, the achievable coupling
strength is weak, and the implementation is probabilistic. Previous methods
[64, 65, 23] to engineer operations using an auxiliary qubit mainly focused on
higher-order single-mode Hamiltonians. However, extending these to two-mode
simulations presented significant complexities, particularly the necessity of
an exotic three-body interaction for [64]. Recently, schemes utilizing the
optical Fredkin gate and cross-Kerr interaction induced by a strong drive have
been proposed [66, 67, 68]. However, these methods involve linearizing a
naturally weak, higher-order optical interaction into a lower-order
interaction, leading to decreased accuracy and resource efficiency. Extending
these methods to simulate various nonlinear couplings, especially between
hybrid modes at the different platforms, is challenging.

In this work, we introduce a highly efficient constructive approach to
simulate the dynamics of deep-strong quantum multimode bosonic coupling. We
show how alternating a selective dispersive interaction conditioned on the
individual Fock state in a cavity mode (SNAP gate) [69, 70] with Rabi gates
can effectively reproduce arbitrary-order radiation-pressure interaction. Our
approach does not rely on existing nonlinearities between oscillators and
offers an exponential enhancement in the resource usage over other methods,
reproducing distinctive effects of radiation pressure and membrane-in-the-
middle optomechanical systems. Superconducting circuits are a promising
platform for realizing this approach. Recent advancements have demonstrated
the feasibility of coupling a transmon qubit to multiple cavity modes and
achieving switchable dispersive and resonant interactions [25, 26, 27, 28].
Our work not only facilitates the exploration of experimentally challenging
nonlinear coupling regimes but also significantly enhances the capabilities of
quantum simulation, opening up new capabilities for hybrid quantum information
processing using qubits and oscillators and tests of new quantum nonlinear
phenomena.

##  II Method

Our objective is to implement general deep-strong cavity optomechanical
Hamiltonians,  denoted as
$\hat{H}_{k,l}=\Omega_{k,l}\hat{n}_{L}^{k}\hat{X}_{M}^{l}$. Here
$\Omega_{k,l}$ represents the coupling strength, and $k,l\in\mathbb{Z}$. Here,
the number operator $\hat{n}_{L}=\hat{a}^{\dagger}\hat{a}$ and a quadrature
position operator $\hat{X}_{M}=\frac{\hat{b}+\hat{b}^{\dagger}}{\sqrt{2}}$
where $\hat{a}$ ($\hat{b}$) is the annihilation operator for the optical
(mechanical) mode. The mode index $L,M$ will be set as $1,2$ in the
simulation. The unitary evolution operator by this Hamiltonian is given as
$\hat{U}_{k,l}[T_{k,l}]=\exp[\mathrm{i}T_{k,l}\hat{n}_{L}^{k}\hat{X}_{M}^{l}]$
where $T_{k,l}$ is the dimensionless strength. The challenge lies in both
nonlinear regimes at the level of a few energy quanta and the deep-strong
coupling that largely dominates over free oscillations of both optical and
mechanical modes. We focus on the lowest orders to demonstrate feasibility
with $k=1$ and $l=1$ representing deep-strong coupling in movable cavity
mirror systems causing radiation-pressure-like coupling linear in the
mechanical operator and the $l=2$ representing the membrane-in-the-middle
optomechanics (Fig. 1a).

###  II.1 Dispersive-Rabi gate method

Figure 1:  Quantum simulation circuits of optomechanical interactions
switching between two types of qubit-oscillator gates. a Target nonlinear
interactions involving radiation pressure coupling in movable cavity mirror
$\hat{H}_{1,1}$ or membrane-in-the-middle coupling $\hat{H}_{1,2}$ between an
optical cavity mode and a mechanical mode. We aim to simulate these systems
using physical experimental setups such as trapped ion or circuit QED by
mapping cavity and mechanical modes to two harmonic oscillators. b Circuit
diagram illustrating the incremental dispersive-Rabi gate approach serving as
a benchmark, alternating dispersive (off-resonant) and Rabi (on-resonant)
qubit-oscillator interactions to induce an effective radiation-pressure
coupling between the modes. This sequence should be repeated multiple times
$M\gg 1$ to simulate a strong target coupling accurately. c Circuit diagram
showing the non-incremental SNAP-Rabi gate approach using $N$ photon-number
selective dispersive gates (SNAP gates) and Rabi gates to accurately reproduce
the movable cavity mirror Hamiltonians. In both methods, the membrane-in-the-
middle coupling can be simulated by substituting $\hat{X}$ by $\hat{X}^{2}$.

Our first approach, which serves a benchmark to a new method in Sec. II. B,
considers a qubit that interacts with two oscillators in the different regimes
of resonance alternatingly between dispersive and Rabi qubit-oscillator gates,
mediating a virtual effective interactions between the oscillators
deterministically without the need for any measurement, as in Fig. 1b. The
basic premise is that by switching between these two different coupling
regimes and leveraging the non-commuting Pauli operators $\hat{\sigma}_{j}$
with $j=\mathrm{x,y,z}$, the resulting operator transformations can
progressively build up higher-order nonlinear interactions. This circuit takes
advantage of the Baker–Campbell–Hausdorff (BCH) formula [71], which allows the
multiplication of non-commuting operators to be approximated as a series of
commutators. This formula possesses the power to systematically generate
higher order nonlinearity when applied recursively [23, 61, 64]. It first
builds a tripartite interaction, where a qubit non-demolitionally mediates the
indirect effective interactions between the two oscillators. A previous
version of this incremental approach exploiting a mediating qubit has been
employed for various purposes [65, 72, 73]. The BCH formula can approximate
the method in [64] for small parameters, reducing it to the dispersive-Rabi
gate method.

For example, alternating two different kinds of unitary interactions such as
dispersive gates
$\hat{U}_{z}=\exp[\mathrm{i}\epsilon\hat{\sigma}_{z}\hat{n}_{1}]$ [25] and
Rabi gates $\hat{U}_{x}=\exp[\mathrm{i}\epsilon\hat{\sigma}_{x}\hat{X}_{2}]$
[2, 1] as in the following sequence makes a weak bosonic interaction:

$\displaystyle\hat{O}_{1,1}[2\epsilon^{2}]=\hat{U}_{x}\hat{U}_{z}\hat{U}_{x}^{\dagger}\hat{U}_{z}^{\dagger}\approx\exp[-2\mathrm{i}\epsilon^{2}\hat{\sigma}_{y}\hat{n}_{1}\hat{X}_{2}]$
(1)

for a very small $\epsilon\ll 1$. The qubit degree of freedom can be ignored
if we initially prepare the eigenstate of the engineered Pauli operator
$\hat{\sigma}_{y}$, such as $\left|\pm_{\mathrm{i}}\right\rangle$. The
dispersive interaction can arise from the Rabi coupling far from resonance
[19], and can be linearized to Rabi intereaction by applying strong
displacement operations [74]. Recent advancements in the experimental
realization of these interactions have been made in cavity [75, 76] and
circuit QED [77, 78, 79, 25] and trapped ion systems [7, 29]. However, this
incremental approach has limited validity at a high strength, requiring a
large number $M\gg 1$ of repetitions to achieve a high strength target gate
[23, 71, 61]. The total resource usage in terms of total coupling strength is
counted as $4M\epsilon$. We can improve the approximation using concatenated
operations as in Appendix A, and online squeezing [80]. A more complex circuit
using concatenated BCH operations can create $\hat{O}_{k,l}[T_{k,l}]$ that
simulates ideal operations $\hat{U}_{k,l}[T_{k,l}]$ for higher $k$ and $l$ as
theoretically shown in [72].

###  II.2 SNAP-Rabi gate method

To overcome the limitations of the incremental method, we introduce a non-
incremental approach that leverages SNAP gates [69, 70] to address each Fock
state independently. We utilize a circuit that bears resemblance with [33].
For approximate simulation, we use SNAP gates that flip the qubit selectively
for a specific Fock state $\left|n\right\rangle$. Combined with Rabi gates,
the following sequence constructs a selective optomechanical coupling on
$n$-th Fock subspace $\hat{\Pi}_{n}\equiv\left|n\right\rangle\left\langle
n\right|$:

$\displaystyle\hat{S}[n]\equiv\exp[-\mathrm{i}\frac{\pi}{2}\hat{\Pi}_{n}\hat{\sigma}_{z}]\exp[\mathrm{i}T_{1,1}n\frac{\hat{1}+\hat{\sigma}_{x}}{2}\hat{X}_{2}]\exp[\mathrm{i}\frac{\pi}{2}\hat{\Pi}_{k}\hat{\sigma}_{z}]$
$\displaystyle=(1-\hat{\Pi}_{n})(\left|+\right\rangle\left\langle+\right|\exp[\mathrm{i}T_{1,1}n\hat{X}_{2}]]+\left|-\right\rangle\left\langle-\right|)+\hat{\Pi}_{n}(\left|+\right\rangle\left\langle+\right|+\left|-\right\rangle\left\langle-\right|\exp[\mathrm{i}T_{1,1}n\hat{X}_{2}]]),$
(2)

where $\left|\pm\right\rangle$ are $\hat{\sigma}_{x}$ eigenstates. Repeating
this sequence for all $n$ upto a finite number $N$ yields an $N$-th order
approximation of the optomechanical interaction without incremental
limitations:

$\displaystyle\hat{\mathbb{S}}_{1,1}[T_{1,1}]$
$\displaystyle\equiv\prod_{n=1}^{N}\hat{S}[n].$ (3)

In the qubit $\left|-\right\rangle$ subspace, the erroneous operation is
eliminated by orthogonality.

An improved SNAP-Rabi gate sequence in terms of simulation accuracy and
resource efficiency is given as (Fig. 1c):

$\displaystyle\mathbb{S}_{1,1}[T_{1,1}]=\exp[-\mathrm{i}\frac{\pi}{2}\sum_{n=0}^{N-1}\hat{\Pi}_{n}\sigma_{z}]\prod_{n=0}^{N-1}\overbrace{\exp[\mathrm{i}T_{1,1}\frac{\hat{1}-\hat{\sigma}_{x}}{2}\hat{X}_{2}]\exp[\mathrm{i}\frac{\pi}{2}\hat{\Pi}_{n}\sigma_{z}]}^{\hat{S}[n]}$
$\displaystyle=\left|-\right\rangle\left\langle-\right|(\hat{1}-\sum_{n=0}^{N}\hat{\Pi}_{n})\exp[\mathrm{i}NT_{1,1}\hat{X}_{2}]+\left|+\right\rangle\left\langle+\right|(\hat{1}-\sum_{n=0}^{N}\hat{\Pi}_{n})+\left|-\right\rangle\left\langle-\right|\sum_{n=0}^{N}\exp[\mathrm{i}nT_{1,1}\hat{X}_{2}]\hat{\Pi}_{n}+\left|+\right\rangle\left\langle+\right|\sum_{n=0}^{N}\exp[\mathrm{i}(N-n)T\hat{X}_{2}]\hat{\Pi}_{n}$
(4)

where the term
$\exp[-\mathrm{i}\frac{\pi}{2}\sum_{n=0}^{N-1}\hat{\Pi}_{n}\sigma_{z}]$ can be
omitted for resource efficiency by sacrificing the accuracy. A visual
depiction of how the SNAP-Rabi gate sequence constructs the simulated
optomechanical interaction in a progressive, selective manner is provided in
Appendix C.  This gate uses only total strength $T_{1,1}N/2$ of the Rabi
gates, as opposed to the $T_{1,1}N(N+1)/4$ in the previous approach in (3),
and thus is even more resource-efficient for a large $N$. The total resource
usage of total coupling strength is given as $T_{1,1}N/2+\pi N$. It also
yields better fidelities under qubit errors, due to the absence of an large
erroneous displacement depending on $N(N+1)$ in (3) and a better behavior in
the high Fock subspace beyond $\hat{\Pi}_{N}$.

It is noteworthy that by using a Rabi interaction of the form
$\exp[\mathrm{i}tF(n)\frac{\hat{1}+\hat{\sigma}_{x}}{2}\hat{X}_{2}]$ with a
modulated strength function $F(n)$, for example $F(n)=n^{k>1}$, instead of a
standard Rabi interaction in (2), we can simulate an arbitrary-order
optomechanical interaction. If we modulate the Rabi gate strength to implement
a second order optomechanical interaction, we obtain:

$\displaystyle\mathbb{S}_{2,1}[T_{2,1}]=\exp[-\mathrm{i}\frac{\pi}{2}\sum_{n=0}^{N-1}\hat{\Pi}_{n}\sigma_{z}]\prod_{n=0}^{N-1}\exp[\mathrm{i}(2n+1)T_{2,1}\frac{\hat{1}-\hat{\sigma}_{x}}{2}\hat{X}_{2}]\exp[\mathrm{i}\frac{\pi}{2}\hat{\Pi}_{n}\sigma_{z}].$
(5)

The total Rabi gate usage is given as $T_{2,1}N^{2}/2$, in contrast to
$O[N^{3}]$ order total strength in the previous method. In general,
implementation of $\mathbb{S}_{k,1}[T_{k,l}]$ uses Rabi gates with total
strength $O[N^{k}]$, in contrast to the previous method that uses
$O[N^{k+1}]$. The relative resource scaling advantages of the SNAP-Rabi
approach compared to other methods is analyzed in more depth in Appendix E.

To engineer a higher optomechanical interaction of the form
$\hat{\mathbb{S}}_{1,2}[T]$ in a similar manner, we require access to a second
order Rabi interaction of the form
$\exp[\mathrm{i}t\hat{\sigma}_{x}\hat{X}^{2}]$ with a varied strength $t$. The
standard Rabi gate can be transformed into the second order Rabi gate of the
form $\hat{\sigma}_{x}\hat{X}^{2}$ at large detuning, with a noise-controlled
squeezing occurring concurrently, as observed in [81]. Alternatively, the
second and higher order Rabi interaction can be efficiently engineered using a
decomposition method [23]. By employing higher-order Rabi interactions
engineered via decomposition methods or derived from dispersive interactions
at large detuning, our technique can simulate the membrane-in-the-middle
interaction $\hat{U}_{\mathrm{1,2}}$ and more generally
$\hat{U}_{\mathrm{k,l}}$, combining such interactions with SNAP gates as
outlined in (2). This SNAP-Rabi technique can be extended further to
multipartite gates such as
$\exp[\mathrm{i}\chi\hat{n}_{1}^{k_{1}}\hat{n}_{2}^{k_{2}}\hat{X}_{3}^{k_{3}}\hat{X}_{4}^{k_{4}}..]$,
providing a versatile toolbox for multimode bosonic simulations. The following
sections compare these simulated dynamics by engineered optomechanical
interactions with the benchmark Kerr-based method. The novel SNAP-Rabi gate
method provides exponential improvements in both accuracy and resource usage
beyond those offered by benchmarks. This method can be also adapted to
engineer a cross-Kerr gate by simultaneously using SNAP gates on both modes.

##  III Results and discussion

The discussion in this section focuses on assessing the non-Gaussian
properties of the simulated nonlinear dynamics that are not accessible by the
linearized Gaussian one. Additionally, we compare our methods with an approach
based on naturally existing higher-order Kerr nonlinear interactions, which is
briefly described in Appendix. B. This Kerr-based method also needs many
repetitions by $M$ times, due to its incremental limitation.

In our system, the implemented optomechanical gates $\hat{O}_{k,l}[T_{k,l}]$
from incremental dispersive-Rabi method (1) and
$\hat{\mathbb{S}}_{k,l}[T_{k,l}]$ from non-incremental SNAP-Rabi methods (3)
are not identical to the ideal optomechanical interactions
$\hat{U}_{k,l}[T_{k,l}]$, and therefore, it is crucial to investigate the
extent to which the observable effects can reproduce the target nonlinear
dynamics for various input states, especially at deep-strong coupling regime
$T_{k,l}\approx 1$. Nonlinear optomechanical processes allow the phase-
insensitive generation of coherent phase-sensitive effects (displacement and
squeezing) on the mechanical mode, which is not possible with linearized
versions. Therefore, phase-insensitive coherent states
$\rho_{\mathrm{prc}}[\beta]=\int_{\phi=0}^{2\pi}d\phi\left|\beta\mathrm{e}^{\mathrm{i}\phi}\right\rangle_{1}\left\langle\beta\mathrm{e}^{\mathrm{i}\phi}\right|=\mathrm{e}^{-|\beta|^{2}}\frac{\beta^{2\hat{n}}}{\hat{n}!}$
with Poissonian statistics and thermal states
$\rho_{\mathrm{th}}[\bar{n}]=\left(\frac{\bar{n}}{\bar{n}+1}\right)^{\hat{n}}/(\bar{n}+1)$
with Bose-Einstein statistics are chosen as input states, whose mean photon
number $\bar{n}$ can be experimentally adjusted. A thermal state represents
the maximum entropy average of arbitrary pure states, subject to the
constraint of a given average number of quanta.  These classical states can be
easily prepared for any oscillators, providing experimentally straightforward
initial states to test quantum effects from the simulations. For entanglement
analysis of simulated interaction, we use phase-sensitive coherent states to
understand the impact of amplitude and phase noise. The state in the on-
resonant mode simulating mechanical mode is chosen to be in a ground state or
in thermal states with the same average quanta from realistic experimental
considerations.

Figure 2:  Simulation infidelity versus resource usage of various methods.
Logarithmic infidelity vs. simulation order $N$ proportional to the resource
used in the SNAP-Rabi method and logarithmic repetition number $\log M$ used
in the other methods  for a,b Poissonian noise and c,d Bose-Einstein noise
both with $\bar{n}=1$ at a,c movable cavity mirror dynamics at strength
$T_{1,1}=1$ and b,d simulation of membrane-in-the-middle dynamics at
$T_{1,2}=1$. Different curves show the impact of qubit loss (dashed) and
dephasing errors (dot-dashed) during the Rabi gates.

We now assess the accuracy of the simulated optomechanical interactions based
on three important criteria: i) fidelity between output states of ideal and
simulated interactions, ii) ability to reproduce optomechanical phenomena like
noise-driven displacement and squeezing, iii) accuracy in generating target
entanglement. The rough evaluation is given by the fidelity, defined as
$F=\mathrm{Tr}[\sqrt{\sqrt{\rho_{\mathrm{id}}}\rho_{\mathrm{re}}\sqrt{\rho_{\mathrm{id}}}}]^{2}$,
between the output states by the target ideal interactions
$\rho_{\mathrm{id}}$ and the approximate interactions $\rho_{\mathrm{re}}$.
This quantifies how close the output states from the simulated dynamics are to
the ideal target dynamics and serves as the average measure of closeness of
the states. Higher fidelity indicates the simulation more accurately
reproduces the desired nonlinear effects. In Fig. 2, the SNAP-Rabi gate method
on Poissonian and Bose-Einstein noise shows a linear scaling in the resource
usage for the infidelity drop, corresponding to an exponential improvement
over the incremental dispersive-Rabi and Kerr-based methods that show
exponential scaling in the resource usage for the infidelity drop. For
example, at $N=10$ of the SNAP-Rabi gate, the total resource gate strength
consumed is about $0.5\%$ of that by the dispersive-Rabi with similar
fidelity. This exponential enhancement in resource usage is maintained even
for simulating the higher-order membrane-in-the-middle interaction. The
results demonstrate the significantly higher accuracy and resource-efficiency
of the SNAP-Rabi gate approach in reproducing the target nonlinear dynamics.
The fidelity scaling trends and analysis of how the different simulation
methods react to increased resources are further detailed in Appendix E. In
our stability analysis, we considered the effects of losses and noise. Each
Rabi and SNAP gate was subjected to a constant level of qubit dephasing
($dq$), boson loss ($d\eta$), and qubit loss ($dr$), although in practical
scenarios these values could depend on gate strengths. These losses led to a
saturation in accuracy for large $N$, a consequence of the interplay between
enhanced approximations and accumulating errors. Notably, these losses had a
more substantial impact on the benchmark methods (dispersive-Rabi method and
Kerr-based method) due to their significantly higher resource usage compared
to the SNAP-Rabi gate method. A detailed summary of the effect of boson loss
is provided in Appendix G.

Secondly, we investigate optical noise-induced mechanical quantum coherence,
which is a basic nonlinear effect at a low number of quanta. The nonlinear
optomechanical interaction can induce displacement in the movable cavity
mirror. It is quantified by the signal-to-noise ratio (SNR) of displacement
comparing the mean induced displacement to the initial uncertainty quantifying
how well the simulated interaction displaces the mechanical state beyond
Gaussian noise (see Appendix D for its definition). The signal-to-noise ratio
(SNR) compares the induced displacement of the mechanical state to its initial
uncertainty. A higher SNR indicates the nonlinear optomechanical interaction
can displace the mechanical state beyond the original noise levels. Squeezing
defined by the smallest eigenvalue of a covariance matrix can be induced in
the mechanical mode by membrane-in-the-middle interaction from phase-
insensitive states in the optical mode depending on the phase-insensitive
photon number. The mechanical quantum Gaussian coherence arising from the
noise in the optical mode serves as a key effective witness of the genuine
nonlinear optomechanical interaction at the quantum level. The noise-induced
squeezing of the membrane-in-the middle also appears qualitatively different
from many other methods used to induce squeezing [82]. The accuracy of the SNR
and induced squeezing has a similar resource scaling as the fidelity. The
details about these metrics is explained in Appendix D.

In Fig 3 a), we demonstrate that the SNRs for the Poissonian and Bose-Einstein
noise with different $\bar{n}$ values in the simulated optical mode, obtained
from the simulated movable cavity mirror interactions $\hat{O}_{1,1}[T_{1,1}]$
and $\hat{\mathbb{S}}_{1,1}[T_{1,1}]$, approach those from the ideal
optomechanical interaction $\hat{U}_{1,1}[T_{1,1}]$ for various $T_{1,1}$,
especially for a low $\bar{n}$. In Fig. 3 b), we illustrate the minimum
variance of quadratures by various processes simulating
$\hat{U}_{1,2}[T_{1,2}]$. These results demonstrate the squeezing generated
beyond a Gaussian approximation of linearized interaction, as indicated by the
minimum variance falling below shot-noise limit $\Delta X_{p}^{2}<1/2$, which
implies the existence of a nonclassical phase-sensitive state. For both the
Poissonian and Bose-Einstein noises in the simulated optical mode with a
ground state in the simulated mechanical mode at the initial time, the general
trends of the SNRs and induced squeezing by the SNAP-Rabi gate are highly
overlapping with that by the ideal processes even at a deep-strong strength
beyond incremental regime, suggesting that our approach effectively achieves
the desired nonlinear effect with a higher accuracy than the Kerr-based
approach.

Figure 3:  Signal-to-noise ratio (SNR) of displacement and squeezing in the
simulated mechanical mode induced by the ideal and simulated dynamics. a SNR
of the displacement in the movable cavity mirror. The incremental dispersive-
Rabi gate was made of $M=50$ repetitions of the unit setups, Kerr-based gate
was made of $M=10$ repetitions, and SNAP-Rabi gate approach utilized $N=10$
approximation order. The target strength was set as $T_{1,1}=1$. b Minimum
quadrature variance $\Delta X_{p}^{2}$ induced by the ideal and approximate
membrane-in-the-middle dynamics at $T_{1,2}=1$. Repetition of $M=50$ was used
for incremental dispersive-Rabi method, while Kerr-based approach was without
repetition. Variance below the shot noise limit (0.5) indicates the generation
of squeezing and nonclassical mechanical states. In all cases, the SNAP-Rabi
gate method shows excellent agreement with ideal results, validating its
ability to accurately reproduce the target dynamics.  Figure 4:  Entanglement
negativity and Gaussian negativities generated by the ideal and simulated
dynamics. The simulation parameters are the same as in Fig. 3. Entanglement
above 1 ebit ($N_{c}=1/2$) and Gaussian entanglement (gray shade) indicates
presence of non-qubit and non-Gaussian entanglement. In b, the output state
has purely non-Gaussian entanglements due to the inherent non-Gaussian nature
of the interaction, as indicated by $0$ Gaussian negativities. In all cases
the output state has a Schmidt rank $\infty$ demonstrating the continuous
variable nature of the entanglement. The SNAP-Rabi gate method closely
reproduces both the Gaussian and non-Gaussian entanglement of the ideal
interaction. This demonstrates its ability to simulate the complex hybrid
continuous-discrete variable entanglement dynamics.

The third criterion characterizes the non-local aspect of various dynamics via
entanglement negativity. We quantify the entanglement of a bipartite state
density operator $\rho$ of simulated optical and mechanical modes by the
entanglement negativity
$\mathbb{N}[\rho]=\frac{\mathrm{Tr}[|\rho^{\mathrm{PT}}|]-1}{2}$, where
$\rho^{\mathrm{PT}}$ is the partial transpose of $\rho$ and
$\mathrm{Tr}[|\cdot|]$ denotes the trace norm [83, 84]. Entanglement
negativity quantifies the amount of entanglement between the optical and
mechanical modes. Higher negativity indicates greater quantum entanglement.
Note that bipartite qubit entanglement of the form
$\left|\psi_{\mathrm{qub}}\right\rangle\equiv
2^{-1/2}(\left|1\right\rangle_{1}\left|0\right\rangle_{2}\pm\left|0\right\rangle_{1}\left|1\right\rangle_{2})$
has a maximal negativity of
$\mathbb{N}[\left|\psi_{\mathrm{qub}}\right\rangle\left\langle\psi_{\mathrm{qub}}\right|]=1/2$
or 1 ebit. This exceeds the entanglement predicted by the covariance matrix
under Gaussian approximation. On the other hand, continuous variable states
such as Gaussian states can have more than 1 ebit of entanglement, but they
lack quantum non-Gaussian properties. Entanglement exceeding 1 ebit with non-
Gaussian characteristics beyond the covariance matrix belongs to a completely
distinct class. An ideal optomechanical interaction has a specific entangling
nature due to the product of two different operators, $\hat{n}_{L}$ and
$\hat{X}_{M}$, which have discrete and continuous spectra of eigenvalues.
Generating such hybrid entanglement is a key signature of nonlinear quantum
effects.

In Fig. 4, the entanglement negativity is compared to Gaussian entanglement
negativity, which refers to the entanglement predicted by the covariance
matrix under Gaussian approximation, produced by the ideal and simulated
dynamics. Both the simulated interactions and ideal optomechanical coupling
generate entanglement that exceeds the Gaussian entanglement negativity,
indicating the presence of authentic non-Gaussian entanglement beyond what is
predicted within the Gaussian approximation framework. States diagonal in the
Fock basis, such as a Bose-Einstein and a Poissonian noise in simulated
optical mode, cannot generate entanglement and only result in classical
correlations. In the simulated optical mode, therefore, we investigate input
coherent states with various average photon number $\bar{n}$ that lacks
quantum nonclassicality and non-Gaussianity. Both target ideal and simulated
dynamics generate entanglement exceeding $1$ ebit. This entanglement surpasses
the entanglement of any Gaussian states with the same covariance matrix [85],
demonstrating both authentic non-Gaussian, beyond-two-qubit entanglement. The
output state by the membrane-in-the-middle dynamics has non-Gaussian
entanglements without any entanglement predicted by the covariance matrix
within Gaussian approximation. The agreement of the SNAP-Rabi gate with the
target ideal coupling is much closer than both the dispersive-Rabi gate and
the Kerr-based simulation. The exponential enhancement in resource-scaling of
the accuracy in the simulation of the entanglement negativity is maintained.
This result demonstrates that SNAP-Rabi gate method reproduces the complex
continuous-discrete variable entanglement arising from the nonlinear coupling
more precisely than the Kerr-based method and the dispersive-Rabi gate method.
The intensity-intensity correlation between the optical and mechanical modes
is analyzed for the different simulation methods and summarized in Appendix F,
providing further validation of the approaches.

Our work primarily operates within the interaction picture and omits the
effective frequency. In current experimental optomechanics setups, the
effective mechanical frequency ($\omega_{m}$) typically falls within the MHz
to GHz range. The actual experimental optomechanical interaction strength
($g_{m}$), however, is typically much smaller, with $g_{m}/\omega_{m}$ on the
order of $10^{-3}$ to $10^{-6}$. In our simulation, we can achieve an
ultrastrong coupling where the resource coupling (denoted as $T/2$) is on the
order of $1$, such that $T_{1,1}/\omega_{m}$ reaching a value of approximately
$2$. This results in a simulated optomechanical interaction strength that is
enhanced by several orders of magnitude relative to current state-of-the-art
experiments. This significant enhancement of the interaction strength
underscores the potential of our proposed method to advance research in this
field. The SNAP-Rabi gate approach demonstrates a substantial efficiency
advantage over conventional techniques such as the dispersive-Rabi and Kerr-
based methods. This considerable gain in resource efficiency allows the SNAP-
Rabi method to accurately simulate deep-strong nonlinear optomechanical
interactions that remain difficult to directly probe and investigate in
existing experimental optomechanical systems. It allows for the exploration of
deep-strong nonlinear couplings between bosonic oscillators, a domain that
remains challenging to access with current technology.

##  IV Conclusions

Quantum simulation of nonlinear bosonic interactions in challenging deep-
strong coupling can allow proof-of-principle tests of many new phenomena.  We
introduced an accurate and efficient method for simulating challenging
nonlinear dynamics by alternating qubit-oscillator coupling regimes. Our
technique does not rely on existing bosonic nonlinearities and provides an
exponential improvement in accuracy resource efficiency over previous methods.
Our simulations demonstrate that the simulated dynamics exhibit the essential
nature of a new nonlinear generator of coherence and nonclassical states and
hybrid entangler, validating the accuracy of our scheme in reproducing the
target dynamics in a previously unexplored regime. The signal-to-noise ratio,
and entanglement negativity analysis confirm the ability to reproduce the
complex dynamics of ideal optomechanical systems beyond the accuracy of
existing Kerr nonlinearity-based approaches. The SNAP-Rabi gate method based
on SNAP gates [69, 70], in particular, shows excellent agreement with ideal
optomechanical interactions. The proposed approach can be extended to simulate
higher-order nonlinearities by utilizing higher-order Rabi gates or dispersive
gates, broadening the class of dynamics that can be simulated.

These first proposals focused on specific targets provide broad and  valuable
tools for hybrid continuous-discrete variable quantum information processing
and simulations of nonlinear bosonic systems. However, implementing this
method needs further experimental development, including switching a single
qubit between resonant and dispersive coupling with multiple cavity modes, a
capability never tested in current superconducting circuit QED and ion trap
systems limited to fixed interaction regimes. While advancements in achieving
switchable qubit-oscillator couplings on platforms like superconducting
circuits suggest that an experimental demonstration may be feasible with
continued progress [19, 12, 25, 26, 28], more research is required to fully
explore their capabilities and overcome limitations in current systems
regarding fixed interaction regimes. Our analysis serves as an important step
toward accessing intriguing dynamics for future quantum technologies.

## Acknowledment

K.P. and R.F. acknowledge the grant No. 22-27431S. P.M. acknowledges the
European Union’s HORIZON Research and Innovation Actions under Grant Agreement
no. 101080173 (CLUSTEC). We also acknowledge EU H2020-WIDESPREAD-2020-5
project NONGAUSS (951737) under the CSA - Coordination and Support Action.

##  Appendix A Incremental dispersive-Rabi gate method

The incremental approach utilizes the commutation relations of the Pauli
operators to systematically generate higher-order nonlinear interactions
through operator transformations. By alternating unitary operations containing
Pauli operators with opposite signs, undesired terms can be suppressed via
cancellations. However, the validity of the incremental approximations is
limited for high interaction strengths, necessitating many repetitions of weak
gate operations.  The basic construction of operators for the incremental
method can be achieved by a complex arrangement of two types of interactions:

$\displaystyle\exp[-\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}]\exp[2\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}]\exp[-\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]$
$\displaystyle=\exp\left[-\mathrm{i}\left(2t_{1}\hat{A}+\frac{\pi}{2}\right)\sigma_{y}\exp[2\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{x}]-\mathrm{i}\frac{\pi}{2}\sigma_{y}\right]$
$\displaystyle\approx\exp[\mathrm{i}\sin[2t_{2}\hat{B}]\cos[2t_{1}\hat{A}]\hat{\sigma}_{z}+\mathrm{i}\sin[4t_{1}A]\sin[Bt_{2}]^{2}\hat{\sigma}_{y}]$
$\displaystyle\approx\exp[\mathrm{i}2t_{2}\hat{B}\cos[2t_{1}\hat{A}]\hat{\sigma}_{z}],$
(6)

and using a sequence with different signs as

$\displaystyle\exp[-\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}]\exp[2\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]\exp[-\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}]\exp[-\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]$
$\displaystyle=\exp\left[2\mathrm{i}\hat{A}t_{1}\sigma_{y}\left(\cos[2\hat{B}t_{2}]\exp[-2\mathrm{i}\hat{A}t_{1}\sigma_{y}]-\mathrm{i}\sigma_{x}\sin[\hat{B}t_{2}]\right)\right]$
$\displaystyle=\exp\left[2\mathrm{i}\hat{A}t_{1}\sigma_{y}\exp[-2\mathrm{i}\hat{A}t_{1}\sigma_{y}]\left(\cos[2\hat{B}t_{2}]-\mathrm{i}\exp[2\mathrm{i}\hat{A}t_{1}\sigma_{y}]\sigma_{x}\sin[\hat{B}t_{2}]\right)\right]$
$\displaystyle=\exp\left[2\mathrm{i}\hat{A}t_{1}\sigma_{y}\exp[-2\mathrm{i}\hat{A}t_{1}\sigma_{y}]\exp[-\mathrm{i}\hat{B}t_{2}\exp[2\mathrm{i}\hat{A}t_{1}\sigma_{y}]\sigma_{x}]\right]$
$\displaystyle\approx\exp[\mathrm{i}\sin[2t_{2}\hat{B}]\sin[2t_{1}\hat{A}]\hat{\sigma}_{x}-\mathrm{i}\sin[4t_{1}A]\sin[t_{2}B]^{2}\hat{\sigma}_{y}]$
$\displaystyle\approx\exp[\mathrm{i}2t_{2}\hat{B}\sin[2t_{1}\hat{A}]\hat{\sigma}_{x}].$
(7)

where the last line was obtained in a weak strength limit $t_{1},t_{2}\ll 1$.
Here, a product form of $\hat{B}f(\hat{A})$ with an infinite order polynomial
$f$ has been achieved in the exponent, which can become the sources of high-
order nonlinear interactions.

This approximation can be further improved by suppressing the errors by
concatenating with similar operations with opposite signs by which the
unwanted terms are cancelled out. More rigorously, we can use the Fourier
series expansion to obtain the target interaction efficiently. Using
$\exp[-\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]\hat{\sigma}_{z}\exp[\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]=\cos[2t_{1}\hat{A}]\hat{\sigma}_{z}+\sin[2t_{1}\hat{A}]\hat{\sigma}_{x}=\hat{\sigma}_{z}\exp[\mathrm{i}\hat{\sigma}_{y}2t_{1}\hat{A}]$
and
$\exp[-\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]\hat{B}\exp[\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]=\hat{B}$,
we have:

$\displaystyle\exp[-\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}]\exp[\mathrm{i}t_{1}\hat{A}\hat{\sigma}_{y}]=\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}\exp[\mathrm{i}2t_{1}\hat{A}\hat{\sigma}_{y}]].$
(8)

Note that this is an exact equation, the validity of which is not limited by
the parameters of $t_{1},t_{2}$. Now these operations can be combined to make
the target function $\exp[\mathrm{i}T\hat{\sigma}_{\theta}\hat{A}\hat{B}]$ for
certain strength $T$ and a certain Pauli matrix $\hat{\sigma}_{\theta}$. Let
us assume for simplicity that $2t_{1}$ takes discrete values, such that
$2t_{1}=k\tau$ with $k=...,-2,-1,0,1,2,...$ and $t_{2}$ can take arbitrary
values. Then we are given the operations as
$\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}\exp[\mathrm{i}k\tau\hat{A}\hat{\sigma}_{y}]]$.
Here, we use the Fourier series expansion
$x=\mathrm{i}\sum_{k=1,2,...}\frac{(-1)^{k}}{k}\exp[\mathrm{i}kx]-\frac{(-1)^{k}}{k}\exp[-\mathrm{i}kx]$,
to have

$\displaystyle\prod_{k=1,2,...}^{k_{\mathrm{max}}=\infty}\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}\frac{(-1)^{k}}{k}\exp[\mathrm{i}k\tau\hat{A}\hat{\sigma}_{y}]]\exp[-\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}\frac{(-1)^{k}}{k}\exp[-\mathrm{i}k\tau\hat{A}\hat{\sigma}_{y}]]$
$\displaystyle\approx\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}\left(\sum_{k=1,2,...}^{k_{\mathrm{max}}}\frac{(-1)^{k}}{k}\exp[\mathrm{i}k\tau\hat{A}\hat{\sigma}_{y}]-\frac{(-1)^{k}}{k}\exp[-\mathrm{i}k\tau\hat{A}\hat{\sigma}_{y}]\right)]=\exp[-\mathrm{i}t_{2}\hat{\sigma}_{x}\tau\hat{A}\hat{B}].$
(9)

For a finite order approximation with $k_{\mathrm{max}}$, we can get a fairly
good approximation as the strength required is sufficiently weak for a high
order contributions.

A better approximation can be obtained using the following approximation:
$AB\approx
0.6(-\mathrm{e}^{\mathrm{i}(-A-B)}+\mathrm{e}^{\mathrm{i}(A-B)}+\mathrm{e}^{-\mathrm{i}(A-B)}-\mathrm{e}^{\mathrm{i}(A+B)})$.
Therefore, we obtain

$\displaystyle\exp[-\mathrm{i}t_{2}\hat{\sigma}_{y}\tau\hat{A}\hat{B}]\approx\exp[-\mathrm{i}t_{2}\hat{\sigma}_{y}(-\mathrm{e}^{\mathrm{i}\hat{\sigma}_{x}(-A-B)}+\mathrm{e}^{\mathrm{i}\hat{\sigma}_{x}(A-B)}+\mathrm{e}^{-\mathrm{i}\hat{\sigma}_{x}(A-B)}-\mathrm{e}^{\mathrm{i}\hat{\sigma}_{x}(A+B)})]$
$\displaystyle=\exp[\mathrm{i}t_{2}\hat{\sigma}_{y}\mathrm{e}^{\mathrm{i}\hat{\sigma}_{x}(-A-B)}]\exp[-\mathrm{i}t_{2}\hat{\sigma}_{y}\mathrm{e}^{\mathrm{i}\hat{\sigma}_{x}(A-B)}]\exp[-it_{2}\hat{\sigma}_{y}\mathrm{e}^{-i\hat{\sigma}_{x}(A-B)}]\exp[it_{2}\hat{\sigma}_{y}\mathrm{e}^{i\hat{\sigma}_{x}(A+B)}].$
(10)

Here we used $\hat{\sigma}_{x}^{2}=1$. Here, each term can be engineered from
(8) with substitution $A\rightarrow A\pm B,B\rightarrow 1$.

We can go a step further by using Fourier transform corresponding to a
continuous limit of Fourier series expansions. We note that for any quantity
$x$, $\int\mathrm{d}xx\exp[\mathrm{i}xp]=-\mathrm{i}\delta^{\prime}(p)$ by the
Fourier transform with a integration variable $p$ where $\delta^{\prime}(p)$
is the derivative of Dirac delta function. Now noticing the approximation
$\delta(p)=\lim_{a\rightarrow 0}\frac{\exp[-(p/a)^{2}]}{\sqrt{\pi}a}$, we can
approximate
$\delta^{\prime}(p)\approx-\frac{2p\mathrm{e}^{-\frac{p^{2}}{a^{2}}}}{\sqrt{\pi}a^{3}}$.
This implies the following approximation
$x\approx\mathrm{i}\int\mathrm{d}p\frac{2p\mathrm{e}^{-\frac{p^{2}}{a^{2}}}}{\sqrt{\pi}a^{3}}\exp[-\mathrm{i}xp]$.
Now arbitrarily choosing small $a$ and substituting integration by summation
by setting $\mathrm{d}p=a$, we have
$x\approx\mathrm{i}\sum_{k}\frac{2k\mathrm{e}^{-k^{2}}}{\sqrt{\pi}a}\exp[-\mathrm{i}kax]$.
Using this equation, we will get

$\displaystyle\exp[-\mathrm{i}t_{2}\hat{\sigma}_{x}\tau\hat{A}\hat{B}]=\exp[-\mathrm{i}t_{2}\hat{\sigma}_{z}\hat{B}i\int\mathrm{d}p\frac{2p\mathrm{e}^{-\frac{p^{2}}{a^{2}}}}{\sqrt{\pi}a^{3}}\exp[-\mathrm{i}\tau\hat{\sigma}_{y}p\hat{A}]]\approx\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}\left(\mathrm{i}\sum_{k}\frac{2k\mathrm{e}^{-k^{2}}}{\sqrt{\pi}a}\exp[-\mathrm{i}\hat{\sigma}_{y}ka\tau\hat{A}]\right)]$
$\displaystyle\approx\prod_{k}\exp[\mathrm{i}t_{2}\hat{B}\hat{\sigma}_{z}\mathrm{i}\frac{2k\mathrm{e}^{-k^{2}}}{\sqrt{\pi}a}\exp[-\mathrm{i}\hat{\sigma}_{y}ka\tau\hat{A}]].$
(11)

Alternatively, we can enhance the strength of the achieved gates by applying
the additional squeezings [86], which transforms a quadrature operator as
$S[-r]\hat{X}S[r]=\hat{X}\mathrm{e}^{r}$. This transformation can also be
applied to operators containing quadrature operators $\exp[f(\hat{X})]$ with
an arbitrary function $f(\cdot)$. In our analysis, we focus on the simplest
case of
$S_{2}[-r]\exp[\mathrm{i}t\hat{n}_{1}\hat{X}_{2}]S_{2}[r]=\exp[\mathrm{i}t\mathrm{e}^{r}\hat{n}_{1}\hat{X}_{2}]$,
which shows that squeezing enhances the gate strength from $t$ to
$t\mathrm{e}^{r}$. On the other side, we can also enhance the strength of the
number operator as $S[-r]\hat{n}S[r]\approx\hat{n}\cosh[r]^{2}$, thus
obtaining a similar enhancement by squeezing
$S_{1}[-r]\exp[\mathrm{i}t\hat{n}_{1}\hat{X}_{2}]S_{1}[r]\approx\exp[\mathrm{i}t\cosh[r]^{2}\hat{n}_{1}\hat{X}_{2}]$.
This technique can increase the feasibility of experimental schemes and
overcome the limit of weak strengths of individual operations, although it
cannot enhance the non-Gaussian nature of the interaction.

##  Appendix B Benchmarks based on higher order interactions

An alternative approach to achieve a high strength radiation-pressure
interaction $\hat{U}_{1,1}$ is to linearize a naturally arising higher-order
cross-Kerr interaction by a large displacement, as demonstrated in [87, 68].
Specifically, we can use the following sequence of interactions:

$\displaystyle\exp[-\mathrm{i}\frac{\alpha^{2}}{2}t\hat{n}_{1}]\exp[\mathrm{i}\alpha\hat{P}_{2}]\exp[\mathrm{i}t\hat{n}_{1}\hat{n}_{2}]\exp[-\mathrm{i}\alpha\hat{P}_{2}]\exp[-\mathrm{i}t\hat{n}_{1}\hat{n}_{2}]\approx\exp[\mathrm{i}t\alpha\hat{n}_{1}\hat{X}_{2}].$
(12)

This approximation linearizes one interacting mode, and holds well when $t\ll
1$, and $\alpha\gg 1$. Again, repetition of this sequence by $M$ times can
increase the strength of the simulated gate or the accuracy of the simulation.
Here, the Kerr interaction can also be engineered by setting $A=\hat{n}_{1}$
and $B=\hat{n}_{2}$ from two dispersive interactions as in (1), or we can use
the naturally existing Kerr interaction existing in QED systems [66, 68].
Another approach to achieve $\hat{U}_{1,1}$ is to transform a Fredkin gate by
a strong drive as described in [67], the process of which is simply described
as :

$\displaystyle\exp[\mathrm{i}t\hat{a}^{\dagger}\hat{a}(\hat{b}^{\dagger}\hat{c}+\hat{c}^{\dagger}\hat{b})]\left|\xi\right\rangle_{c}\approx\exp[\mathrm{i}t\hat{a}^{\dagger}\hat{a}(\hat{b}^{\dagger}\xi+\xi^{*}\hat{b})]\left|\xi\right\rangle_{c}.$
(13)

However, we do not delve into this approach, as it requires an even less
accessible three-mode interaction than the cross-Kerr interaction. We can also
simulate the membrane-in-the-middle interaction $\hat{U}_{1,2}$ using Kerr
interactions and squeezing. This can be achieved using the following
transformation:

$\displaystyle\exp[-\mathrm{i}\frac{t}{2}(\mathrm{e}^{2r}-1)\hat{n}_{1}]\exp[-\mathrm{i}t\mathrm{e}^{2r}\hat{n}_{1}\hat{n}_{2}]S_{2}[r]\exp[\mathrm{i}t\hat{n}_{1}\hat{n}_{2}]S_{2}[-r]\approx\exp[\mathrm{i}t\sinh[2r]\hat{n}_{1}\hat{X}_{2}^{2}]$
(14)

where the squeezing operator on the mechanical mode
$S_{2}[r]=\exp[-\frac{r}{2}\hat{b}^{\dagger 2}+\frac{r}{2}\hat{b}^{2}]$. This
approximation holds well when $t\ll 1$ and $|r|\gtrsim 1$. These Kerr-
interaction-based methods are inefficient in terms of resource usage due to
the necessity of inverse interactions, but are the main targets for comparison
with our protocols. A similar technique was used to enhance the strengths of
Kerr gate in [88].

##  Appendix C SNAP-Rabi gate sequence

Figure C.1: A diagram visualizing the SNAP-Rabi gate sequence on joint qubit-
oscillator states. Each horizontal level shows the effect of the gate set
$S[n]$. It flips the qubit in X-basis for $n$-th Fock space, and displaces the
mechanical mode conditioned by the higher Fock spaces. The upward arrows
represent conditional displacement $D[\mathrm{i}T/\sqrt{2}]$ of the mechanical
oscillator. This induces selective displacements that build up the simulated
nonlinear optomechanical interaction in each Fock subspace.

Fig. C.1 describes intuitively how SNAP-Rabi gate builds up the simulated
movable cavity mirror dynamics. By combining conditional qubit flips mediated
by the SNAP gates and conditioned displacements mediated by the Rabi gates,
the sequence builds up the simulated nonlinear interaction in a selective,
progressive manner within each Fock subspace. Overall, the total displacement
applied to the mechanical mode is proportional to the Fock number $n$ in the
optical mode, as required to accurately simulate the optomechanical dynamics.

##  Appendix D Signal to noise ratio and induced mechanical squeezing

In addition to fidelity, other metrics can also characterize how well the
simulated interactions reproduce the target optomechanical dynamics. The
first-order optomechanical interaction $\hat{U}_{1,1}$ possesses phase
sensitivity in the mechanical mode, which results in displacement beyond the
noise that is proportional to the photon number in the optical mode. In
contrast, the linearized interaction cannot generate net mechanical
displacement for phase-insensitive states, such as PRC or thermal states.
Therefore, the successful implementation of the optomechanical interaction
$\hat{U}_{1,1}[T]$ can be demonstrated by producing a mechanical displacement
beyond the linearized regime, which generates displacement beyond shot noise
and results in a mechanical state with non-zero off-diagonal elements in the
Fock state basis in the mechanical mode. These off-diagonal elements indicate
the quantum coherence induced by the optical noise.

To evaluate the effectiveness of implementing the optomechanical interaction
$\hat{U}_{1,1}[T]$, we can use the SNR to measures the displacement generated
relative to the initial Gaussian uncertainty located at the phase space
origin. In the Heisenberg picture, any operator $\hat{O}$ under a unitary
evolution $\hat{U}$ is described as $\hat{U}^{\dagger}\hat{O}\hat{U}$. The
momentum operator shifted by the adjoint action of an ideal $\hat{U}_{1,1}$
interaction is described by

$\displaystyle\hat{P}^{\prime}_{2}=\hat{P}_{2}-t\hat{n}_{1},~{}~{}\delta\hat{P}_{2}=\hat{P}^{\prime}_{2}-\hat{P}_{2}=-t\hat{n}_{1}.$
(15)

In contrast, the momentum operator shifted by the linearized interaction is
described by:

$\displaystyle\hat{P}^{\prime}_{2}=\hat{P}_{2}-t\hat{X}_{1},~{}~{}\delta\hat{P}_{2}=\hat{P}^{\prime}_{2}-\hat{P}_{2}=-t\hat{X}_{1}.$
(16)

The average momentum shift in (15) can be non-zero for an input Fock state in
the optical mode $1$, while that in (16) is zero due to the phase
insensitivity of the Fock states, or the states diagonal in the Fock basis
used in the simulation. Therefore, a non-zero average momentum shift can serve
as evidence of the nonlinear optomechanical interaction, but a sufficient
strength interaction is required to have the shift beyond shot-noise limit.

We can calculate the SNR using the formula
${\color[rgb]{0,0,0}\mathrm{SNR}=}\langle\delta\hat{P}_{2}\rangle/\sqrt{\mathrm{Var}(\delta\hat{P}_{2})}$
for various average photon number $\bar{n}$ in the optical mode, where the
average $\langle\cdot\rangle$ and the variance $\mathrm{Var}(\cdot)$ is
calculated for the initial states. Initially, the simulated mechanical mode is
assumed to be in a vacuum state without quantum coherence in any basis,
including the Fock state basis. For PRC in the optical mode, the SNR is
calculated to be equal to
${\color[rgb]{0,0,0}\mathrm{SNR}^{\mathrm{(PRC)}}=}\sqrt{\bar{n}}=\alpha$ for
the ideal interaction $\hat{U}_{1,1}[T]$, which is the coherent amplitude of
the state, and is independent of the strength $T$. For a thermal state in the
optical mode, the SNR increases approximately as
$\mathrm{SNR}^{\mathrm{(th)}}=\sqrt{\frac{\bar{n}}{\bar{n}+1}}\stackrel{{\scriptstyle\bar{n}\ll
1}}{{\to}}\sqrt{\bar{n}}$, and asymptotically saturates to the value $1$ when
$\bar{n}\gg 1$. Again, the SNR is independent of the strength $T$.

On the other hand, a higher order interaction such as
$\hat{U}_{\mathrm{1,2}}[T]=\exp[\mathrm{i}T\hat{n}_{1}\hat{X}_{2}^{2}]$ can
induce both the quantum coherence and the nonclassical squeezing effects. The
presence of a square term in $\hat{X}$ in the Hamiltonian provides a mechanism
for quantum squeezing in the mechanical mode in general. By using this second
order ideal optomechanical interaction, we can derive operator relations in
Heisenberg picture as follows:

$\displaystyle\hat{P}_{2}^{\prime}=\hat{P}_{2}+2T\hat{n}_{1}\hat{X}_{2}.$ (17)

These relations allow us to obtain the momentum variance:

$\displaystyle\Delta P_{2}^{2}=\langle P_{2}^{2}\rangle_{T}-\langle
P_{2}\rangle_{T}^{2}$ $\displaystyle=\langle P_{2}^{2}\rangle-\langle
P_{2}\rangle^{2}+4T\langle\hat{n}_{1}\rangle\Big{(}\langle\frac{\hat{X}_{2}\hat{P}_{2}-\hat{P}_{2}\hat{X}_{2}}{2}\rangle-\langle
X_{2}\rangle\langle P_{2}\rangle\Big{)}$
$\displaystyle+4T^{2}\langle\hat{n}_{1}^{2}\rangle(\langle\hat{X}_{2}\rangle^{2}-\langle\hat{X}_{2}\rangle^{2}).$
(18)

We observe that the uncertainty $\Delta P_{2}^{2}$ monotonically increases
with $T$ due to the product of covariances of mechanical mode $2$ and photon
number moments of optical mode $1$. This increased uncertainty in $P$ implies
the existence of squeezing along some direction of quadrature in phase space.

##  Appendix E Resource Scaling Analysis of Simulation Methods

Figure E.1: Scaling behavior against resource usage over various mean photon
number of Poissonian and Bose-Einstein noise models. a Simulation of movable
cavity mirror dynamics. b Simulation of Membrane-in-the-middle dynamics.

The metrics for various method of simulation can be compared in their scaling
against the resources used and the mean photon number for a complete picture.
This can be quantified by the linear slope of the curves in Fig. 2. The
fidelities generally improve exponentially with more resource usage, modeled
empirically as $1-F[R]=(1-F[0])e^{\delta_{R}R}$, where $R$ represents the
resource used, and $\delta_{R}$ is termed as the infidelity decay rate
representing how simulation accuracy reacts the resource usage. In Fig. E.1,
the infidelity decay rate for movable cavity mirror and membrane-in-the-middle
dynamics are shown, with Poissonian and Bose-Einstein noise models. The
resource usage was counted as $N$ for SNAP-Rabi gate method and $M$ for the
other methods. The SNAP-Rabi gate method shows a much larger decay rate,
meaning the input resource rapidly improves the simulation fidelity.
Noticeably the decay rate of SNAP-Rabi gate method depends more on the mean
photon number of the noise model due to its finite dimension nature in Fock
subspace. Again, the Kerr-based method has a very small decay rates, meaning
its repetition is not very effective in improving the simulation accuracy.

Superposition states of a few Fock states can exhibit nonlinear optomechanical
effects that differ from statistical mixtures, highlighting the role of
quantum coherence. Such superpositions may be useful for generating highly
non-Gaussian quantum states, although preparing and characterizing specific
Fock state superpositions presents additional experimental challenges. Here,
we test the performance of the simulated first-order optomechanical operation
on randomly sampled states on the off-resonant mode $1$, and a vacuum state on
the resonant mode $2$. The complex coefficients of $c_{n}$ random sampled
states $\sum_{n}\mathrm{e}^{-n/2}c_{n}\left|n\right\rangle$ can be sampled
with a constraint $-1\leq\mathrm{Re}[c_{n}],\mathrm{Im}[c_{n}]\leq 1$.

The SNAP-Rabi gate method yields advantageously much higher fidelities, with
an infidelity as low as $1-F=10^{-8}$ using $10$ SNAP gates for these states
with average excitation number $1$ at $T=1$. We also observe high fidelities
around $1-F\approx 10^{-4.5}$ for randomly sampled states. Fock states in the
off-resonant modes give the optimal fidelity for the simulation performed
here.

##  Appendix F Other non-local characteristics

Intensity-intensity correlation Intensity-intensity correlation quantifies the
correlation between the photon number fluctuations in the optical and
mechanical modes. This metric directly captures the nonlocal energy transfer
caused by the nonlinear coupling. It is defined as

$\displaystyle
g^{(2)}(0)=\frac{\langle\hat{n}_{1}\hat{n}_{2}\rangle}{\langle\hat{n}_{1}\rangle\langle\hat{n}_{2}\rangle},$
(19)

where $\hat{n}_{1,2}$ are the photon number operators for the optical and
mechanical modes, respectively. For ideal optomechanical interaction, we can
use the following equations

$\displaystyle\exp[-it\hat{n}_{1}\hat{X}_{2}]\hat{n}_{1}\exp[it\hat{n}_{1}\hat{X}_{2}]=\hat{n}_{1},$
(20)
$\displaystyle\exp[-it\hat{n}_{1}\hat{X}_{2}]\hat{n}_{2}\exp[it\hat{n}_{1}\hat{X}_{2}]=\hat{n}_{2}+t\hat{n}_{1}\hat{P}_{2}+\frac{t^{2}}{2}\hat{n}_{1}^{2}.$
(21)

For an ideal optomechanical interaction on the coherent states or equivalently
on the phase randomized coherent states and a vacuum state in the mechanical
mode, the radiation pressure coupling results in:

$\displaystyle\langle\hat{n}_{1}\rangle=\bar{n},~{}\langle\hat{n}_{2}\rangle=\frac{t^{2}}{2}(\bar{n}^{2}+\bar{n}),~{}\langle\hat{n}_{1}\hat{n}_{2}\rangle=\frac{t^{2}}{2}(\bar{n}^{3}+3\bar{n}^{2}+\bar{n}),$
(22)

where $\bar{n}$ is the average photon number in the optical mode. Therefore,
the intensity-intensity correlation is given as

$\displaystyle
g^{(2)}_{\mathrm{Poisonnian}}(0)=\frac{\bar{n}^{2}+3\bar{n}+1}{\bar{n}^{2}+\bar{n}},$
(23)

a value independent from the strength of the interactions, and diverging at
$\bar{n}\to 0$ and saturating as $\lim_{\bar{n}\to\infty}=1$. For an initial
thermal state in the optical mode, the ideal interaction gives

$\displaystyle\langle\hat{n}_{1}\rangle=\bar{n},~{}\langle\hat{n}_{2}\rangle=\frac{t^{2}}{2}(2\bar{n}^{2}+\bar{n}),~{}\langle\hat{n}_{1}\hat{n}_{2}\rangle=\frac{t^{2}}{2}(6\bar{n}^{3}+6\bar{n}^{2}+\bar{n}),$
(24)

and

$\displaystyle g^{(2)}_{\mathrm{Bose-
Einstein}}(0)=\frac{6\bar{n}^{2}+6\bar{n}+1}{2\bar{n}^{2}+\bar{n}},$ (25)

also diverging at $\bar{n}\to 0$ and saturating as
$\lim_{\bar{n}\to\infty}=3$.

For separable states, this quantity will give $g^{(2)}(0)=1$, and deviation
from this value will witness the existence of entanglement. Deviations from
these ideal values can diagnose the accuracy of simulated optomechanical
interactions. Calculating and comparing $g^{(2)}(0)$ for the output states of
both simulated and ideal interactions could provide a good metric of how well
the nonlinear coupling is reproduced.

Figure F.1 shows the intensity-intensity correlation $g^{(2)}(0)$ between the
optical and mechanical modes as a function of the mean photon number in the
optical mode. The correlation obtained from the ideal optomechanical
interaction is compared against the dispersive-Rabi gate, SNAP-Rabi gate, and
Kerr-based simulation methods. For a coherent state input, all three
simulation methods can reproduce the ideal correlation relatively well across
the range of $\bar{n}$, with the SNAP gate approach providing the closest
match. This provides further evidence that the simulated dynamics mimic the
desired nonlocal effects of the target nonlinear coupling. However, for a
thermal state input, simulating the ideal correlation becomes more
challenging, especially in the large $\bar{n}$ limit, due to the long tail of
the thermal photon number distribution.

Figure F.1: Intensity-intensity correlation $g^{(2)}(0)$ between optical mode
and mechanical mode vs the average photon number $\bar{n}$ in the optical mode
for different simulation methods. Initial states are  (a) an initial thermal
state and (b) an initial phase-randomized coherent state in the optical mode.
The mechanical mode is assumed to be in the vacuum state. The SNAP gate method
provides better agreement with the ideal optomechanical correlation,
particularly for larger $\bar{n}$.

##  Appendix G Boson loss

Figure G.1:  Logarithmic infidelity versus resource usage by SNAP-Rabi gate
method for simulating bosonic Gaussian operations using a movable cavity
mirror and membrane-in-the-middle. Different lines correspond to various loss
rates per gate, from $1-d\eta=10^{-4}$ to $10^{-1}$ where $d\eta$ is the loss
parameter occurring at each gate. As resource usage increases, infidelity
eventually saturates due to the constant loss rate per gate.

Bosonic loss can significantly affect the efficiency of quantum simulation
protocols. To account for these losses, we consider a loss model where each
bosonic mode undergoes amplitude damping at a rate $\gamma$. This loss process
can be described by the master equation in the Lindblad form:

$\dot{\rho}=\sum_{j=1}^{2}\gamma(2\hat{a}_{j}\rho\hat{a}_{j}^{{\dagger}}-\hat{a}_{j}^{{\dagger}}\hat{a}_{j}\rho-\rho\hat{a}_{j}^{{\dagger}}\hat{a}_{j}),$
(26)

where $\hat{a}_{j}$ ($\hat{a}_{j}^{{\dagger}}$) is the annihilation (creation)
operator for the $j$-th bosonic mode. Any single mode density matrix $\rho$
undergoing bosonic loss channel $\Gamma_{\eta}$ evolves as
$\Gamma_{\eta}[\rho]=\sum_{l=0}^{\infty}\frac{(1-\eta)^{l}}{l!}\eta^{\hat{n}/2}\hat{a}^{l}\rho\hat{a}^{\dagger
l}\eta^{\hat{n}/2}$, where the loss parameter is defined as
$\eta=\mathrm{e}^{-\gamma t}\in(0,1]$ and $\gamma t$ is the dimensionless
damping parameter, equivalent to Kraus operator notation.

To evaluate the impact of bosonic losses on the accuracy of our simulation
schemes, we calculate the infidelity of our simulation protocols under
different loss rates per gate. Figure H.1 shows the logarithmic infidelity
versus resource usage for the SNAP-Rabi gate method for simulating bosonic
Gaussian operations using a movable cavity mirror and membrane-in-the-middle.
Different lines correspond to various loss rates per gate, from $10^{-4}$ to
$10^{-1}$. For both systems, the infidelity increases with loss rate and
decreases with more resources used.

As for the question about dissipation in the referee report, we would like to
point out that in realistic scenarios, dissipation can indeed impact the
performance of our simulation protocols. However, our SNAP-Rabi gate method
shows a certain level of resilience against bosonic losses. As can be seen
from Fig. G.1, even under moderate loss rates, the method still maintains a
relatively high fidelity, indicating that it could be robust against
dissipation in practical implementations. Moreover, the rapid decrease of
infidelity with increased resource usage suggests that the detrimental effects
of the mild losses can be mitigated by using more resources, even though for a
larger resource usage, the infidelity saturates due to the fixed per gate loss
rate.

## References

  * Forn-Díaz  _et al._ [2019] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Reviews of Modern Physics 91, 25005 (2019).
  * Frisk Kockum  _et al._ [2019] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, F. Nori, A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Nature Reviews Physics 1, 19 (2019).
  * Vijay  _et al._ [2012] R. Vijay, C. Macklin, D. H. Slichter, S. J. Weber, K. W. Murch, R. Naik, A. N. Korotkov, and I. Siddiqi, Nature 490, 77 (2012).
  * Devoret and Schoelkopf [2013] M. H. Devoret and R. J. Schoelkopf, Science 339, 1169 (2013).
  * Tan  _et al._ [2015] T. R. Tan, J. P. Gaebler, Y. Lin, Y. Wan, R. Bowler, D. Leibfried, and D. J. Wineland, Nature 528, 380 (2015).
  * Milne  _et al._ [2020] A. R. Milne, C. L. Edmunds, C. Hempel, F. Roy, S. Mavadia, and M. J. Biercuk, Physical Review Applied 13, 024022 (2020).
  * Flühmann  _et al._ [2018] C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. Home, Nature 566, 513 (2018).
  * Puri  _et al._ [2020] S. Puri, L. St-Jean, J. A. Gross, A. Grimm, N. E. Frattini, P. S. Iyer, S. T. A. Krishna, L. Jiang, A. Blais, S. T. Flammia, and S. M. Girvin, Science advances 6, eaay5901 (2020).
  * Grimm  _et al._ [2020] A. Grimm, N. E. Frattini, S. Puri, S. O. Mundhada, S. Touzard, M. Mirrahimi, S. M. Girvin, S. Shankar, and M. H. Devoret, Nature 584, 205 (2020).
  * Yoshihara  _et al._ [2016] F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Nature Physics 13, 44 (2016).
  * Heeres  _et al._ [2017] R. W. Heeres, P. Reinhold, N. Ofek, L. Frunzio, L. Jiang, M. H. Devoret, and R. J. Schoelkopf, Nature Communications 8, 94 (2017).
  * Gao  _et al._ [2018] Y. Y. Gao, B. J. Lester, Y. Zhang, C. Wang, S. Rosenblum, L. Frunzio, L. Jiang, S. M. Girvin, and R. J. Schoelkopf, Physical Review X 8, 021073 (2018).
  * Bruzewicz  _et al._ [2019] C. D. Bruzewicz, J. Chiaverini, R. McConnell, and J. M. Sage, Applied Physics Reviews 6, 021314 (2019).
  * Wolf  _et al._ [2019] F. Wolf, C. Shi, J. C. Heip, M. Gessner, L. Pezzè, A. Smerzi, M. Schulte, K. Hammerer, and P. O. Schmidt, Nature Communications 10, 2929 (2019).
  * Katz and Monroe [2023] O. Katz and C. Monroe, Phys. Rev. Lett. 131, 033604 (2023).
  * Blais  _et al._ [2021] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Reviews of Modern Physics 93, 025005 (2021).
  * Gu  _et al._ [2017] X. Gu, A. F. Kockum, A. Miranowicz, Y. xi Liu, and F. Nori, Physics Reports 718-719, 1 (2017).
  * Ofek  _et al._ [2016] N. Ofek, A. Petrenko, R. Heeres, P. Reinhold, Z. Leghtas, B. Vlastakis, Y. Liu, L. Frunzio, S. M. Girvin, L. Jiang, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, Nature 536, 1 (2016).
  * Zueco  _et al._ [2009] D. Zueco, G. M. Reuther, S. Kohler, and P. Hänggi, Physical Review A - Atomic, Molecular, and Optical Physics 80, 033846 (2009).
  * Lv  _et al._ [2018] D. Lv, S. An, Z. Liu, J.-N. Zhang, J. S. Pedernales, L. Lamata, E. Solano, and K. Kim, Phys. Rev. X 8, 021027 (2018).
  * Marković  _et al._ [2018] D. Marković, S. Jezouin, Q. Ficheux, S. Fedortchenko, S. Felicetti, T. Coudreau, P. Milman, Z. Leghtas, and B. Huard, Physical Review Letters 121, 040505 (2018).
  * Park  _et al._ [2016] K. Park, P. Marek, and R. Filip, Physical Review A 94, 062308 (2016).
  * Park  _et al._ [2018a] K. Park, P. Marek, and R. Filip, New Journal of Physics 20, 053022 (2018a).
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys. 84, 621 (2012).
  * Eickbusch  _et al._ [2022] A. Eickbusch, V. Sivak, A. Z. Ding, S. S. Elder, S. R. Jha, J. Venkatraman, B. Royer, S. M. Girvin, R. J. Schoelkopf, and M. H. Devoret, Nature Physics 18, 1464 (2022).
  * Campagne-Ibarcq  _et al._ [2020] P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zalys-Geller, N. E. Frattini, V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf, L. Frunzio, M. Mirrahimi, and M. H. Devoret, Nature 584, 368 (2020).
  * Stehlik  _et al._ [2021] J. Stehlik, D. M. Zajac, D. L. Underwood, T. Phung, J. Blair, S. Carnevale, D. Klaus, G. A. Keefe, A. Carniol, M. Kumph, M. Steffen, and O. E. Dial, Phys. Rev. Lett. 127, 080505 (2021).
  * Pan  _et al._ [2023] X. Pan, J. Schwinger, N.-N. Huang, P. Song, W. Chua, F. Hanamura, A. Joshi, F. Valadares, R. Filip, and Y. Y. Gao, Physical Review X 13, 021004 (2023).
  * Flühmann and Home [2020] C. Flühmann and J. P. Home, Phys. Rev. Lett. 125, 043602 (2020).
  * Krastanov  _et al._ [2015] S. Krastanov, V. V. Albert, C. Shen, C. L. Zou, R. W. Heeres, B. Vlastakis, R. J. Schoelkopf, and L. Jiang, Physical Review A - Atomic, Molecular, and Optical Physics 92, 1 (2015).
  * Li  _et al._ [2018] X. Li, Y. Ma, J. Han, T. Chen, Y. Xu, W. Cai, H. Wang, Y. Song, Z.-Y. Xue, Z.-q. Yin, and L. Sun, Phys. Rev. Appl. 10, 054009 (2018).
  * Langford  _et al._ [2017] N. K. Langford, R. Sagastizabal, M. Kounalakis, C. Dickel, A. Bruno, F. Luthi, D. J. Thoen, A. Endo, and L. Dicarlo, Nature Communications 8, 1715 (2017).
  * Kudra  _et al._ [2022] M. Kudra, M. Kervinen, I. Strandberg, S. Ahmed, M. Scigliuzzo, A. Osman, D. P. Lozano, M. O. Tholén, R. Borgani, D. B. Haviland, G. Ferrini, J. Bylander, A. F. Kockum, F. Quijandría, P. Delsing, and S. Gasparinetti, PRX Quantum 3, 030301 (2022).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Reviews of Modern Physics 86, 1391 (2014).
  * Millen  _et al._ [2019] J. Millen, T. S. Monteiro, R. Pettit, and A. N. Vamivakas, Reports on Progress in Physics 83, 026401 (2019).
  * O’Connell  _et al._ [2010] A. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bialczak, M. Lenander, E. Lucero, M. Neeley, D. Sank, H. Wang, M. Weides, _et al._ , Nature 464, 697 (2010).
  * Wollman  _et al._ [2015] E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwald, F. Marquardt, A. A. Clerk, and K. C. Schwab, Science 349, 952 (2015).
  * Ockeloen-Korppi  _et al._ [2018] C. Ockeloen-Korppi, E. Damskägg, J. Pirkkalainen, T. Heikkilä, F. Massel, and M. Sillanpää, Nature 556, 478–482 (2018).
  * Clark  _et al._ [2016] J. B. Clark, F. Lecocq, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Nature Physics 12, 683 (2016).
  * Shkarin  _et al._ [2019] A. B. Shkarin, A. D. Kashkanova, C. D. Brown, S. Garcia, K. Ott, J. Reichel, and J. G. E. Harris, Phys. Rev. Lett. 122, 153601 (2019).
  * Wiederhecker  _et al._ [2019] G. S. Wiederhecker, P. Dainese, and T. P. Mayer Alegre, APL Photonics 4, 071101 (2019).
  * Piergentili  _et al._ [2018] P. Piergentili, L. Catalini, M. Bawaj, S. Zippilli, N. Malossi, R. Natali, D. Vitali, and G. D. Giuseppe, New Journal of Physics 20, 083024 (2018).
  * Tavernarakis  _et al._ [2018] A. Tavernarakis, A. Stavrinadis, A. Nowak, I. Tsioutsios, A. Bachtold, and P. Verlot, Nature Communications 9, 662 (2018).
  * Lombardi  _et al._ [2018] A. Lombardi, M. K. Schmidt, L. Weller, W. M. Deacon, F. Benz, B. de Nijs, J. Aizpurua, and J. J. Baumberg, Phys. Rev. X 8, 011016 (2018).
  * Renninger  _et al._ [2018] W. H. Renninger, P. Kharel, R. O. Behunin, and P. T. Rakich, Nature Physics 14, 601–607 (2018).
  * Delić  _et al._ [2020] U. Delić, D. Grass, M. Reisenbauer, T. Damm, M. Weitz, N. Kiesel, and M. Aspelmeyer, Quantum Science and Technology 5, 025006 (2020).
  * Johnsson  _et al._ [2021] M. T. Johnsson, P. M. Poggi, M. A. Rodriguez, R. N. Alexander, and J. Twamley, Physical Review Research 3, 1 (2021), 2101.12364 .
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153 (2014).
  * Stannigel  _et al._ [2010] K. Stannigel, P. Rabl, A. S. Sørensen, P. Zoller, and M. D. Lukin, Physical Review Letters 105, 220501 (2010).
  * Marinković  _et al._ [2018] I. Marinković, A. Wallucks, R. Riedinger, S. Hong, M. Aspelmeyer, and S. Gröblacher, Phys. Rev. Lett. 121, 220404 (2018).
  * Pikovski  _et al._ [2012] I. Pikovski, M. R. Vanner, M. Aspelmeyer, M. S. Kim, and Č. Brukner, Nature Physics 8, 393 (2012).
  * Xuereb  _et al._ [2012] A. Xuereb, C. Genes, and A. Dantan, Phys. Rev. Lett. 109, 1 (2012).
  * Rimberg  _et al._ [2014] A. J. Rimberg, M. P. Blencowe, A. D. Armour, and A. P. D. Nation, New J. Phys. 16, 055008 (2014).
  * Heikkilä  _et al._ [2014] T. T. Heikkilä, F. Massel, J. Tuorila, R. Khan, and M. A. Sillanpää, Phys. Rev. Lett. 112, 3 (2014).
  * Pirkkalainen  _et al._ [2015] J.-M. Pirkkalainen, S. U. Cho, F. Massel, J. Tuorila, T. T. Heikkilä, P. J. Hakonen, and M. A. Sillanpää, Nat. Commun. 6, 6981 (2015).
  * Lü  _et al._ [2015] X.-Y. Lü, Y. Wu, J. R. Johansson, H. Jing, J. Zhang, and F. Nori, Phys. Rev. Lett. 114, 093602 (2015).
  * Lemonde  _et al._ [2016] M.-A. Lemonde, N. Didier, and A. A. Clerk, Nat. Commun. 7, 11338 (2016).
  * Wang and Safavi-Naeini [2017] Z. Wang and A. H. Safavi-Naeini, Nat. Commun. 8, 15886 (2017).
  * Dumont  _et al._ [2019] V. Dumont, S. Bernard, C. Reinhardt, A. Kato, M. Ruf, and J. C. Sankey, Optics Express 27, 25731 (2019).
  * Hartmann [2016] M. J. Hartmann, Journal of Optics 18, 104005 (2016).
  * Kalajdzievski and Quesada [2021] T. Kalajdzievski and N. Quesada, Quantum 5, 1 (2021).
  * Park  _et al._ [2015] K. Park, P. Marek, and R. Filip, Physical Review A - Atomic, Molecular, and Optical Physics 92, 033813 (2015).
  * Park  _et al._ [2014] K. Park, P. Marek, and R. Filip, Physical Review A - Atomic, Molecular, and Optical Physics 90, 013804 (2014).
  * Jacobs [2007] K. Jacobs, Physical Review Letters 99, 117203 (2007).
  * Park  _et al._ [2017] K. Park, P. Marek, and R. Filip, Scientific reports 7, 11536 (2017).
  * Liao  _et al._ [2020] J. Q. Liao, J. F. Huang, L. Tian, L. M. Kuang, C. P. Sun, and C. P. Sun, Physical Review A 101, 063802 (2020).
  * Yin  _et al._ [2022] X. L. Yin, Y. H. Zhou, and J. Q. Liao, Physical Review A 105, 013504 (2022).
  * Zhou  _et al._ [2021] Y.-H. Zhou, X.-L. Yin, and J.-Q. Liao, Optics Express 29, 28202 (2021).
  * Schuster  _et al._ [2007] D. I. Schuster, A. A. Houck, J. A. Schreier, A. Wallraff, J. M. Gambetta, A. Blais, L. Frunzio, J. Majer, B. Johnson, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Nature 445, 515 (2007).
  * Heeres  _et al._ [2015] R. W. Heeres, B. Vlastakis, E. Holland, S. Krastanov, V. V. Albert, L. Frunzio, L. Jiang, and R. J. Schoelkopf, Physical Review Letters 115, 137002 (2015).
  * Van-Brunt and Visser [2015] A. Van-Brunt and M. Visser, Journal of Physics A: Mathematical and Theoretical 48, 1 (2015).
  * Park  _et al._ [2018b] K. Park, P. Marek, and R. Filip, New Journal of Physics 20, 053022 (2018b).
  * Tufarelli  _et al._ [2009] T. Tufarelli, M. S. Kim, and S. Bose, Physical Review A - Atomic, Molecular, and Optical Physics 80, 062317 (2009).
  * van Loock  _et al._ [2008] P. van Loock, W. J. Munro, K. Nemoto, T. P. Spiller, T. D. Ladd, S. L. Braunstein, and G. J. Milburn, Physical Review A 78, 22303 (2008).
  * Auffeves  _et al._ [2003] A. Auffeves, P. Maioli, T. Meunier, S. Gleyzes, G. Nogues, M. Brune, J. M. Raimond, and S. Haroche, Phys. Rev. Lett. 91, 230405 (2003).
  * Friebe [2018] K. Friebe,  _On dispersive interactions between a trapped ion and a cavity field_ , Ph.D. thesis (2018).
  * Premaratne  _et al._ [2017] S. P. Premaratne, F. C. Wellstood, and B. S. Palmer, Nature Communications 8, 1 (2017).
  * Peronnin  _et al._ [2019] T. Peronnin, D. Marković, Q. Ficheux, and B. Huard, Phys. Rev. Lett. 124, 180502 (2019).
  * Touzard  _et al._ [2019] S. Touzard, A. Kou, N. E. Frattini, V. V. Sivak, S. Puri, A. Grimm, L. Frunzio, S. Shankar, and M. H. Devoret, Phys. Rev. Lett. 122, 080502 (2019).
  * Hastrup  _et al._ [2020] J. Hastrup, K. Park, R. Filip, and U. L. Andersen, Physical Review Letters 126, 153602 (2020).
  * Sandu [2015] T. Sandu, Romanian Journal of Physics 60, 711 (2015).
  * Sete and H. [2014] E. A. Sete and E. H., Physical Review A 29, 013841 (2014).
  * Vidal and Werner [2002] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
  * Lee  _et al._ [2000] J. Lee, M. S. Kim, Y. J. Park, and S. Lee, J.Mod.Opt 47, 2151 (2000).
  * Laurat  _et al._ [2005] J. Laurat, G. Keller, J. A. Oliveira-Huguenin, C. Fabre, T. Coudreau, A. Serafini, G. Adesso, and F. Illuminati, Journal of Optics B: Quantum and Semiclassical Optics 7, S577 (2005).
  * Albarelli  _et al._ [2018] F. Albarelli, M. G. Genoni, M. G. A. Paris, and A. Ferraro, Physical Review A 98, 52350 (2018).
  * Yanagimoto  _et al._ [2020] R. Yanagimoto, T. Onodera, E. Ng, L. G. Wright, P. L. Mcmahon, and H. Mabuchi, Physical Review Letters 124, 240503 (2020).
  * Bartkowiak  _et al._ [2014] M. Bartkowiak, L. A. Wu, and A. Miranowicz, Journal of Physics B: Atomic, Molecular and Optical Physics 47, 145501 (2014).


# On Machine Learning DoS Attack Identification from Cloud Computing Telemetry

João Henrique Corrêa, Patrick Marques Ciarelli,
Moises R. N. Ribeiro, Rodolfo da Silva Villaca  0000-0002-8124-8985 Federal
University of Espírito SantoAv. Fernando Ferrari, 514 - Goiabeiras - Vitória -
Espírito Santo - BrasilVitória - Espírito Santo - Brasil<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

The detection of Denial of Service (DoS) attacks remains a challenge for the
cloud environment, affecting a massive number of services and applications
hosted by such virtualized infrastructures. Typically, in the literature, the
detection of DoS attacks is performed solely by analyzing the traffic of
packets in the network. This work advocates for the use of telemetry from the
cloud to detect DoS attacks using Machine Learning algorithms. Our hypothesis
is based on richness of such native data collection services, with metrics
from both physical and virtual hosts. Our preliminary results demonstrate that
DoS can be identified accurately with k-Nearest Neighbors (kNN) and decision
tree (CART).

††copyright: none††doi: ††isbn: ††price:

##  1\. Introduction

Denial of Service (DoS) attacks are a pressing problem in computer networks,
and it becomes an even more complex problem in cloud computing environments.
According to the Verisign report of 2018, these kinds of attacks directly
impact business, financial, information technology and telecommunications
services (Verisign, 2018).

The migration of traditional services and applications to centralized cloud
environments amplifies their vulnerability to disruption in relation to the
impact of DoS attacks performed over conventional architectures. Thus, there
is an urgent need to perform the identification and characterization of DoS
traffic in cloud environments. In general, there are three ways to perform
traffic classification: i) identification using TCP/UDP port numbers; ii)
inspecting the contents of the network packets; and, finally, iii) applying
Machine Learning (ML) techniques to network data source (Amaral et al., 2016).
This paper advances the latter, and proposes the use of ML algorithms in order
to distinguish malicious traffic from legitimate traffic from clients in a
cloud environment. ML algorithms have been already used for traffic
identification and classification (Fadlullah et al., 2017), including
malicious ones, to spot the occurrence of a DoS attack.

According to a recent survey (Boutaba et al., 2018), the traditional data
source used by ML algorithms are traces of traffic in the network (packets).
It is well-known that DoS attacks systemically affects the usage of cloud
computing resources. Different from traditional approached based on traffic
traces, this work proposed the use of the telemetry from the cloud (such as
resources usage from physical and virtual hosts) as data source for ML
algorithms.

Large scale monitoring traffic in conventional networks usually involves
costly and complex architectures, probe packets and other artifices. In
contrast, clouds have native telemetry, i.e., data collection services.
Metrics from both physical and virtual hosts. Moreover, such rich data set can
be used with little or no overhead to the cloud. Therefore, the research
hypothesis of this paper is that it is possible to improve identification and
accuracy of DoS attacks using the information from the usage of computing
resources from physical and virtual hosts allocated to the applications and
services.

##  2\. Evaluation

As an initial scenario of this proposal, we monitor the usage of computing
resources of a cloud under the TCP SYN Flood attack. This attack has been
chosen because it is simple to perform and reproduce. An experimental setup
has been created with an OpenStack cloud (Rocky version), containing a web
server using Apache2, a virtual machine simulating legitimate web clients with
Siege 111https://www.joedog.org/siege-home/. Finally, another host performs
the role of malicious users and generates the SYN Flood Attack with
hping3222https://linux.die.net/man/8/hping3.

(a) CPU (b) Memory (c) HD Request read (d) HD Request write

(e) Interface Bytes in (f) Interface Bytes out (g) Interface packets in (h)
Interface packets out

Figure 1. Comparison between the metrics in both scenarios.

Two test scenarios have been performed. In the first test there is only
legitimate clients, generating traffic at from users so that the web server
can hold it. In the second test, however, along with legitimate clients, SYN
Flood attack is executed. During both tests, the following metrics are
monitored by the OpenStack cloud: CPU and Memory usage, number of disk read
requests, number of disk write requests, rate of incoming and outgoing Bytes
in each network interface, rate of incoming and outgoing packets in each
network interface.

In both tests, the total duration is 30 minutes. There are attack phases with
different intensities, which last 10 minutes each. There are 10s gaps at
beginning and at ending parts of our trace with no legitimate clients’ nor
attacker’s traffic. Metrics are collected every 5 seconds performing a total
of 360 samples. The attack is divided into three phases of 10 minutes each. In
the first one, the attack happens at a small rate i.e., at every 300
milliseconds. In the second phase, 10 minutes after, the malicious attacks
happens every 250 milliseconds. Finally, in the last phase, the attack is
generated at the maximum rate that the attacking host can perform.

Given the measurements of the two experiments, a Principal Component Analysis
(PCA) was used to verify which metrics are relevant in the case of a TCP SYN
Flood attack. A sample comparison of these metrics is shown in Figure 1. In
all the results, the blue lines represent the tests performed only in the
client scenario and the red lines in the client and attacker scenario.
According to these results the following metrics were chosen: CPU usage,
number of disk write requests, rate of incoming and outgoing Bytes, rate of
incoming and outgoing packets.

After the selection of the most important characteristics, the traffic
labeling in the data is performed. These samples are used as a training
dataset for two different ML algorithms: k-Nearest Neighbors (kNN) and a
decision tree (CART). These simple algorithms were chosen for the sake of
simplicity to obtain preliminary results in order to verify our hypothesis
regarding ML exploiting could telemetry. The data from this dataset has been
labeled in two classes: ”attack” and ”no attack”.

In addition, another experiment was carried out, with a duration of 120
minutes, alternating between moments of traffic with only legitimate clients,
only attacker and a combination between both legitimate clients and attackers.
With the collected metrics of this experiment, the trained kNN and CART tries
to identify attack periods.

The Table 1 presents the results (in percentage) of the Accuracy, Precision,
Recall and F1-Score metrics obtained in the dataset using the kNN and CART
algorithms, being the first line for the kNN and the last line for the
decision tree algorithm (CART). In these preliminary tests a high accuracy was
verified, with low incidence of false positives and false negatives, mainly in
the kNN algorithm.

Table 1. Result of ML algorithms ML Algorithms | Accuracy | Precision | Recall | F1-Score
---|---|---|---|---
kNN | 99.30 | 99.33 | 99.31 | 99.31
CART | 91.07 | 95.53 | 91.08 | 92.38

##  3\. Conclusion

Given the promising results, the proposal to use data from the resources usage
from cloud computing environment is a feasible hypothesis.

As future work, we aim at performing the characterization of other types of
DoS attacks, as well as combinations of these types. We intend not only to
detect that an attack is occurring (anomaly detection), but also, through ML
algorithms, to perform (real-time) the identification of what kind of attack
is going on. Another future work is to explore different ML algorithms,
including the use of deep learning.

## References

  * (1)
  * Amaral et al. (2016) P. Amaral, J. Dinis, P. Pinto, L. Bernardo, J. Tavares, and H. S. Mamede. 2016\.  Machine Learning in Software Defined Networks: Data collection and traffic classification. In 2016 IEEE 24th International Conference on Network Protocols (ICNP).
  * Boutaba et al. (2018) Raouf Boutaba, Mohammad A. Salahuddin, Noura Limam, Sara Ayoubi, Nashid Shahriar, Felipe Estrada-Solano, and Oscar M. Caicedo. 2018.  A comprehensive survey on machine learning for networking: evolution, applications and research opportunities.  Journal of Internet Services and Applications 9, 1 (21 Jun 2018), 16.
  * Fadlullah et al. (2017) Z. Fadlullah, F. Tang, B. Mao, N. Kato, O. Akashi, T. Inoue, and K. Mizutani. 2017\.  State-of-the-Art Deep Learning: Evolving Machine Intelligence Toward Tomorrow’s Intelligent Network Traffic Control Systems.  IEEE Communications Surveys Tutorials PP, 99 (2017), 1–1.
  * Verisign (2018) Verisign. 2018.  What Is a DDoS Attack?  https://www.verisign.com/en_US/security-services/ddos-protection/what-is-a-ddos-attack/index.xhtml. (2018).  Accessed: 2019-03-11.