Efficient quantum simulation of nonlinear interactions using SNAP and Rabi gates

Advances in the strong coupling regime of qubit-oscillator interactions in trapped ions and superconducting circuits [1–3] has opened new revenues for building qubit gates [4–6] and, mainly, controllable gates of the oscillators [7–9] and first nonlinear gates [10–12]. Trapped ions [13–15] and superconducting circuit QED with transmons [16–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 [16, 19]. 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 [12, 19, 25–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. Matsos et al [34] prepared various states from non-linear interactions using linear operations, optimising the pulse drive instead of the gate parameters, possibly applicable to general quantum simulations. 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 [35, 36] in nonlinear optomechanics. This interaction is fundamentally nonlinear, involving third power of ladder operators in the Hamiltonian [35]. 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 [37–40], will substantially advance current quantum technologies [41–48] and therefore, it is attractive for hybrid quantum simulations [49]. Nonlinear pressure-like phenomena are predicted to enable breakthrough applications in quantum information processing, metrology, communication, fundamental tests of quantum mechanics [50–52].

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 [35], removing the underlying nonlinear dynamics. Experimental strategies like amplification have also been employed to access the single-photon scale nonlinear optomechanical phenomena [53–59]. The membrane-in-the-middle configuration provides enhanced nonlinear optomechanical interactions [60], 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 [61]. There are several methods available for engineering continuous-variable gates [62, 63]. One such method involves simulating a conditional form of optomechanical interaction [64] using photonic ancillas and homodyne detectors [65]. Despite its potential, the achievable coupling strength is weak, and the implementation is probabilistic. Previous methods [23, 66, 67] 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 [66]. Recently, schemes utilizing the optical Fredkin gate and cross-Kerr interaction induced by a strong drive have been proposed [68–70]. 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 [63], efficient simulation in terms of the finite dimensions of Hilbert space has been explored, yet its applicability in infinite-dimensional spaces remains unclear.

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) [71, 72] 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–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.

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 (figure 1(a)).

Zoom In Zoom Out Reset image size

2.1. Dispersive-Rabi gate method

Our first approach, which serves a benchmark to a new method in section 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 figure 1(b). 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 [73], 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, 62, 66]. 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 [23, 67, 74]. The BCH formula can approximate the method in [66] 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]$ [1, 2] as in the following sequence makes a weak bosonic interaction:

$$\begin{align} \hat{O}_{1,1}\left[2\epsilon^2\right] = \hat{U}_x \hat{U}_z \hat{U}_x^\dagger \hat{U}_z^\dagger\approx \exp\left[-2\mathrm{i}\epsilon^2 \hat{\sigma}_y \hat{n}_1\hat{X}_2\right] \end{align} \tag{ 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\vert\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 [75]. Recent advancements in the experimental realization of these interactions have been made in cavity [76, 77] and circuit QED [25, 78–80] and trapped ion systems [7, 29]. However, this incremental approach has limited validity at a high strength, requiring a large number $M\gg1$ of repetitions to achieve a high strength target gate [23, 62, 73]. 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 [81]. 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 [23].

2.2. SNAP-Rabi gate method

To overcome the limitations of the incremental method, we introduce a non-incremental approach that leverages SNAP gates [71, 72] 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\vert 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\vert n\right\rangle\left\langle n\right\vert$:

$$\begin{align} \hat{S}[n] & \equiv \exp\left[-\mathrm{i} \frac{\pi}{2} \hat{\Pi}_n \hat{\sigma}_z\right]\exp\left[\mathrm{i} T_{1,1} n \frac{\hat{1}+\hat{\sigma}_x}{2}\hat{X}_2\right]\exp\left[\mathrm{i} \frac{\pi}{2} \hat{\Pi}_k \hat{\sigma}_z\right]\nonumber\\ & = \left(1-\hat{\Pi}_n\right)\left(\left\vert+\right\rangle\left\langle+\right\vert\exp[\mathrm{i} T_{1,1} n \hat{X}_2]]+\left\vert-\right\rangle\left\langle-\right\vert\right)+\hat{\Pi}_n \left(\left\vert+\right\rangle\left\langle+\right\vert+\left\vert-\right\rangle\left\langle-\right\vert\exp[\mathrm{i} T_{1,1} n \hat{X}_2]]\right), \end{align} \tag{ 2 }$$

where $\left\vert\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:

$$\begin{align} \hat{\mathbb{S}}_{1,1}\left[T_{1,1}\right]&\equiv\prod_{n = 1}^N \hat{S}\left[n\right]. \end{align} \tag{ 3 }$$

In the qubit $\left\vert-\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 (figure 1(c)):

$$\begin{align} \mathbb{S}_{1,1}\left[T_{1,1}\right] & = \exp\left[-\mathrm{i} \frac{\pi}{2}\sum_{n = 0}^{N-1}\hat{\Pi}_{n} \sigma_z\right]\prod_{n = 0}^{N-1} \overbrace{\exp\left[\mathrm{i} T_{1,1} \frac{\hat{1}-\hat{\sigma}_x}{2}\hat{X}_2\right]\exp\left[\mathrm{i} \frac{\pi}{2}\hat{\Pi}_{n} \sigma_z\right]}^{\hat{S}\left[n\right]}\nonumber\\ & = \left\vert-\right\rangle\left\langle-\right\vert\left(\hat{1}-\sum_{n = 0}^N \hat{\Pi}_n\right)\exp\left[\mathrm{i} N T_{1,1} \hat{X}_2\right]+\left\vert+\right\rangle\left\langle+\right\vert\left(\hat{1}-\sum_{n = 0}^N \hat{\Pi}_n\right)\nonumber\\ & \quad + \left\vert-\right\rangle\left\langle-\right\vert\sum_{n = 0}^N \exp\left[\mathrm{i} n T_{1,1} \hat{X}_2\right]\hat{\Pi}_n+ \left\vert+\right\rangle\left\langle+\right\vert\sum_{n = 0}^N \exp\left[\mathrm{i} \left(N-n\right) T \hat{X}_2\right]\hat{\Pi}_n \end{align} \tag{ 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} t F(n)\frac{\hat{1}+\hat{\sigma}_x}{2}\hat{X}_2]$ with a modulated strength function F(n), for example $F(n) = n^{k\gt1}$, 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:

$$\begin{align} &\mathbb{S}_{2,1}\left[T_{2,1}\right] = \exp\left[-\mathrm{i} \frac{\pi}{2}\sum_{n = 0}^{N-1}\hat{\Pi}_{n} \sigma_z\right]\prod_{n = 0}^{N-1} \exp\left[\mathrm{i} \left(2n+1\right) T_{2,1} \frac{\hat{1}-\hat{\sigma}_x}{2}\hat{X}_2\right]\exp\left[\mathrm{i} \frac{\pi}{2}\hat{\Pi}_{n} \sigma_z\right]. \end{align} \tag{ 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 [82]. 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.

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} \mathrm{d}\phi \left\vert\beta \mathrm{e}^{\mathrm{i}\phi}\right\rangle_1\left\langle\beta \mathrm{e}^{\mathrm{i}\phi}\right\vert = \mathrm{e}^{-|\beta|^2}\frac{\beta^{2 \hat{n}}}{\hat{n}!}$ with Poissonian statistics and thermal states $\rho_\mathrm{th}[\bar{n}] =$ $(\frac{\bar{n}}{\bar{n}+1})^{\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.

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 figure 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.

Zoom In Zoom Out Reset image size

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 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 [83]. 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 figure 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 figure 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\lt1/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.

Zoom In Zoom Out Reset image size

The third criterion characterizes the non-local aspect of various dynamics via entanglement negativity. We quantify the entanglement of a bipartite state density operator ρ 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 ρ and $\mathrm{Tr}[|\cdot|]$ denotes the trace norm [84, 85]. 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\vert\psi_\mathrm{qub}\right\rangle\equiv2^{-1/2}(\left\vert1\right\rangle_1\left\vert0\right\rangle_2\pm \left\vert0\right\rangle_1\left\vert1\right\rangle_2)$ has a maximal negativity of $\mathbb{N}[\left\vert\psi_\mathrm{qub}\right\rangle\left\langle\psi_\mathrm{qub}\right\vert] = 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 figure 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 [86], 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.

Zoom In Zoom Out Reset image size

Our work primarily operates within the interaction picture and omits the effective frequency. In current experimental optomechanics setups, the effective mechanical frequency (ω_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⁻³ to 10⁻⁶. 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.

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 [71, 72], 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 [12, 19, 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.

K P and R F acknowledge the Grant No. 22-27431S. K P acknowledges the Grant number LUAUS23012, funded by the Ministry of Education, Youth and Sports. 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.

All data that support the findings of this study are included within the article (and any supplementary files).

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:

$$\begin{align} &\exp\left[-\mathrm{i} t_1\hat{A} \hat{\sigma}_y\right] \exp\left[\mathrm{i} t_2\hat{B} \hat{\sigma}_z\right]\exp\left[2\mathrm{i} t_1\hat{A} \hat{\sigma}_y\right]\exp\left[\mathrm{i} t_2\hat{B} \hat{\sigma}_z\right]\exp\left[-\mathrm{i} t_1\hat{A} \hat{\sigma}_y\right]\nonumber\\ & \quad = \exp \left[-\mathrm{i} \left(2 t_1\hat{A} +\frac{\pi }{2}\right) \sigma _y \exp\left[2 \mathrm{i} t_2\hat{B} \hat{\sigma}_x\right]-\mathrm{i}\frac{\pi}{2} \sigma _y\right]\nonumber\\ & \quad \approx\exp\left[\mathrm{i} \sin\left[2t_2\hat{B}\right]\cos\left[2 t_1\hat{A}\right]\hat{\sigma}_z+\mathrm{i}\sin\left[4t_1A\right]\sin\left[B t_2\right]^2\hat{\sigma}_y\right]\nonumber\\ & \quad \approx\exp\left[\mathrm{i} 2t_2\hat{B}\cos\left[2 t_1\hat{A}\right]\hat{\sigma}_z\right], \end{align} \tag{ A1 }$$

and using a sequence with different signs as

$$\begin{align} &\exp\left[-\mathrm{i} t_1\hat{A} \hat{\sigma}_y\right] \exp\left[\mathrm{i} t_2\hat{B} \hat{\sigma}_z\right]\exp\left[2\mathrm{i} t_1\hat{A} \hat{\sigma}_y\right]\exp\left[-\mathrm{i} t_2\hat{B} \hat{\sigma}_z\right]\exp\left[-\mathrm{i} t_1\hat{A} \hat{\sigma}_y\right]\nonumber\\ & \quad = \exp \left[2 \mathrm{i} \hat{A} t_1 \sigma _y \left(\cos\left[2 \hat{B} t_2\right] \exp\left[-2 \mathrm{i} \hat{A} t_1 \sigma _y\right]-\mathrm{i} \sigma _x \sin\left[\hat{B} t_2\right]\right)\right]\nonumber\\ & \quad = \exp \left[2 \mathrm{i} \hat{A} t_1 \sigma _y \exp\left[-2 \mathrm{i} \hat{A} t_1 \sigma _y\right]\left(\cos\left[2 \hat{B} t_2\right] -\mathrm{i} \exp\left[2 \mathrm{i} \hat{A} t_1 \sigma _y\right] \sigma _x \sin\left[\hat{B} t_2\right]\right)\right]\nonumber\\ & \quad = \exp \left[2 \mathrm{i} \hat{A} t_1 \sigma _y \exp\left[-2 \mathrm{i} \hat{A} t_1 \sigma _y\right]\exp\left[-\mathrm{i}\hat{B} t_2\exp\left[2 \mathrm{i} \hat{A} t_1 \sigma _y\right] \sigma _x\right]\right]\nonumber\\ & \quad \approx\exp\left[\mathrm{i} \sin\left[2t_2\hat{B}\right]\sin\left[2 t_1\hat{A}\right]\hat{\sigma}_x-\mathrm{i} \sin\left[4 t_1 A \right] \sin\left[t_2B \right]^2\hat{\sigma}_y\right]\nonumber\\ & \quad \approx\exp\left[\mathrm{i} 2 t_2\hat{B}\sin\left[2 t_1\hat{A}\right]\hat{\sigma}_x\right]. \end{align} \tag{ A2 }$$

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[2 t_1\hat{A}]\hat{\sigma}_z+\sin[2 t_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:

$$\begin{align} &\exp\left[-\mathrm{i} t_1\hat{A} \hat{\sigma}_y\right]\exp\left[\mathrm{i} t_2\hat{B}\hat{\sigma}_z \right]\exp\left[\mathrm{i} t_1\hat{A} \hat{\sigma}_y\right] = \exp\left[\mathrm{i} t_2\hat{B}\hat{\sigma}_z\exp\left[\mathrm{i} 2t_1 \hat{A}\hat{\sigma}_y\right]\right]. \end{align} \tag{ A3 }$$

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 = {\ldots},-2,-1,0,1,2,{\ldots}$ and t₂ 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,{\ldots}}\frac{(-1)^k}{k}\exp[\mathrm{i} k x]-$ $\frac{(-1)^k}{k}\exp[-\mathrm{i} k x]$, to have

$$\begin{align} &\prod_{k = 1,2,{\ldots}}^{k_\mathrm{max} = \infty}\exp\left[\mathrm{i} t_2\hat{B}\hat{\sigma}_z\frac{\left(-1\right)^k}{k}\exp\left[\mathrm{i} k \tau \hat{A}\hat{\sigma}_y\right]\right]\exp\left[-\mathrm{i} t_2\hat{B}\hat{\sigma}_z\frac{\left(-1\right)^k}{k}\exp\left[-\mathrm{i} k \tau \hat{A}\hat{\sigma}_y\right]\right]\nonumber\\ & \quad \approx \exp\left[\mathrm{i} t_2\hat{B}\hat{\sigma}_z\left(\sum_{k = 1,2,{\ldots}}^{k_\mathrm{max}}\frac{\left(-1\right)^k}{k}\exp\left[\mathrm{i} k \tau \hat{A}\hat{\sigma}_y\right]-\frac{\left(-1\right)^k}{k}\exp\left[-\mathrm{i} k \tau \hat{A}\hat{\sigma}_y\right]\right)\right] = \exp\left[-\mathrm{i} t_2\hat{\sigma}_x \tau \hat{A}\hat{B}\right].\end{align} \tag{ A4 }$$

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: $A B\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

$$\begin{align} &\exp\left[-\mathrm{i} t_2\hat{\sigma}_y \tau \hat{A}\hat{B}\right]\approx \exp\left[-\mathrm{i} t_2\hat{\sigma}_y \left(-\mathrm{e}^{\mathrm{i} \hat{\sigma}_x \left(-A-B\right)}+\mathrm{e}^{\mathrm{i} \hat{\sigma}_x \left(A-B\right)}+\mathrm{e}^{-\mathrm{i} \hat{\sigma}_x \left(A-B\right)}-\mathrm{e}^{\mathrm{i} \hat{\sigma}_x\left(A+B\right)}\right)\right]\nonumber\\ & \quad = \exp\left[\mathrm{i} t_2\hat{\sigma}_y \mathrm{e}^{\mathrm{i}\hat{\sigma}_x \left(-A-B\right)}\right]\exp\left[-\mathrm{i} t_2\hat{\sigma}_y \mathrm{e}^{\mathrm{i} \hat{\sigma}_x \left(A-B\right)}\right]\exp\left[-i t_2\hat{\sigma}_y \mathrm{e}^{-i \hat{\sigma}_x\left(A-B\right)}\right]\exp\left[i t_2\hat{\sigma}_y \mathrm{e}^{i \hat{\sigma}_x\left(A+B\right)}\right]. \end{align} \tag{ A5 }$$

Here we used $\hat{\sigma}_x^2 = 1$. Here, each term can be engineered from (A3) 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}x x \exp[\mathrm{i} x p] = -\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{2 p \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{2 p \mathrm{e}^{-\frac{p^2}{a^2}}}{\sqrt{\pi } a^3} \exp[-\mathrm{i} x p]$. 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{2 k \mathrm{e}^{-k^2}}{\sqrt{\pi } a} \exp[-\mathrm{i} k a x ]$. Using this equation, we will get

$$\begin{align} \exp\left[-\mathrm{i} t_2\hat{\sigma}_x \tau \hat{A}\hat{B}\right] & = \exp\left[-\mathrm{i} t_2\hat{\sigma}_z \hat{B} i\int \mathrm{d}p \frac{2 p \mathrm{e}^{-\frac{p^2}{a^2}}}{\sqrt{\pi } a^3} \exp\left[-\mathrm{i} \tau\hat{\sigma}_y p \hat{A} \right]\right]\nonumber\\ & \approx\exp\left[\mathrm{i} t_2\hat{B}\hat{\sigma}_z\left( \mathrm{i} \sum_{k} \frac{2 k \mathrm{e}^{-k^2}}{\sqrt{\pi } a} \exp\left[-\mathrm{i} \hat{\sigma}_y k a \tau \hat{A} \right]\right)\right]\nonumber\\ & \approx\prod_{k} \exp\left[\mathrm{i} t_2\hat{B}\hat{\sigma}_z \mathrm{i} \frac{2 k \mathrm{e}^{-k^2}}{\sqrt{\pi } a} \exp\left[-\mathrm{i} \hat{\sigma}_y k a \tau \hat{A} \right]\right].\end{align} \tag{ A6 }$$

Alternatively, we can enhance the strength of the achieved gates by applying the additional squeezings [87], 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.

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 [70, 88]. Specifically, we can use the following sequence of interactions:

$$\begin{align} \exp\left[-\mathrm{i} \frac{\alpha^2}{2}t \hat{n}_1\right] \exp\left[\mathrm{i} \alpha \hat{P}_2\right]\exp\left[\mathrm{i} t \hat{n}_1\hat{n}_2\right]\exp\left[-\mathrm{i} \alpha \hat{P}_2\right]\exp\left[-\mathrm{i} t \hat{n}_1\hat{n}_2\right]\approx\exp\left[\mathrm{i} t\alpha \hat{n}_1\hat{X}_2\right]. \end{align} \tag{ B1 }$$

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 [68, 70]. Another approach to achieve $\hat{U}_{1,1}$ is to transform a Fredkin gate by a strong drive as described in [69], the process of which is simply described as :

$$\begin{align} \exp\left[\mathrm{i} t \hat{a}^\dagger \hat{a}\left(\hat{b}^\dagger \hat{c}+\hat{c}^\dagger \hat{b}\right)\right]\left\vert\xi\right\rangle_c\approx \exp\left[\mathrm{i} t \hat{a}^\dagger \hat{a}\left(\hat{b}^\dagger \xi+\xi^* \hat{b}\right)\right]\left\vert\xi\right\rangle_c. \end{align} \tag{ B2 }$$

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:

$$\begin{align} \exp\left[- \mathrm{i} \frac{t}{2} \left(\mathrm{e}^{2r}-1\right)\hat{n}_1\right]\exp\left[-\mathrm{i} t \mathrm{e}^{2r}\hat{n}_1\hat{n}_2\right]S_2\left[r\right]\exp\left[\mathrm{i} t \hat{n}_1\hat{n}_2\right]S_2\left[-r\right]\approx \exp\left[\mathrm{i} t \sinh\left[2r\right]\hat{n}_1\hat{X}_2^2\right] \end{align} \tag{ B3 }$$

where the squeezing operator on the mechanical mode $S_2[r] = \exp[-\frac{r}{2} \hat{b}^{\dagger2}+\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 [89].

Figure 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.

Zoom In Zoom Out Reset image size

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

$$\begin{align} &\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. \end{align} \tag{ D1 }$$

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

$$\begin{align} &\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. \end{align} \tag{ D2 }$$

The average momentum shift in (D1) can be non-zero for an input Fock state in the optical mode 1, while that in (D2) 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 ${\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 ${\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{\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:

$$\begin{align} &\hat{P}_2^{^{\prime}} = \hat{P}_2+2T\hat{n}_1\hat{X}_2. \end{align} \tag{ D3 }$$

These relations allow us to obtain the momentum variance:

$$\begin{align} \Delta P_2^2 & = \langle P_2^2\rangle_T-\langle P_2\rangle_T^2\nonumber\\ & = \langle P_2^2\rangle-\langle P_2\rangle^2+4T \langle\hat{n}_1\rangle\left(\left\langle\frac{\hat{X}_2\hat{P}_2-\hat{P}_2\hat{X}_2}{2}\right\rangle-\langle X_2\rangle\langle P_2\rangle\right)\nonumber\\ & \quad +4T^2\langle\hat{n}_1^2\rangle\left(\langle\hat{X}_2\rangle^2-\langle\hat{X}_2\rangle^2\right). \end{align} \tag{ D4 }$$

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.

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 figure 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 δ_R is termed as the infidelity decay rate representing how simulation accuracy reacts the resource usage. In figure 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.

Zoom In Zoom Out Reset image size

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\vert n\right\rangle$ can be sampled with a constraint $-1\unicode{x2A7D} \mathrm{Re}[c_n],\mathrm{Im}[c_n]\unicode{x2A7D} 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.

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

$$\begin{align} g^{\left(2\right)}\left(0\right) = \frac{\langle\hat{n}_1\hat{n}_2\rangle}{\langle\hat{n}_1\rangle\langle\hat{n}_2\rangle}, \end{align} \tag{ F1 }$$

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

$$\begin{align} \exp\left[-i t \hat{n}_1\hat{X}_2\right]\hat{n}_1\exp\left[i t \hat{n}_1\hat{X}_2\right] & = \hat{n}_1, \end{align} \tag{ F2 }$$

$$\begin{align} \exp\left[-i t \hat{n}_1\hat{X}_2\right]\hat{n}_2\exp\left[i t \hat{n}_1\hat{X}_2\right] & = \hat{n}_2+t \hat{n}_1 \hat{P}_2+\frac{t^2}{2}\hat{n}_1^2. \end{align} \tag{ F3 }$$

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:

$$\begin{align} \langle\hat{n}_1\rangle = \bar{n}, ~\langle\hat{n}_2\rangle = \frac{t^2}{2}\left(\bar{n}^2+\bar{n}\right), ~\langle\hat{n}_1\hat{n}_2\rangle = \frac{t^2}{2}\left(\bar{n}^3+3\bar{n}^2+\bar{n}\right), \end{align} \tag{ F4 }$$

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

$$\begin{align} g^{\left(2\right)}_\mathrm{Poisonnian}\left(0\right) = \frac{\bar{n}^2+3\bar{n}+1}{\bar{n}^2+\bar{n}}, \end{align} \tag{ F5 }$$

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

$$\begin{align} \langle\hat{n}_1\rangle = \bar{n}, ~\langle\hat{n}_2\rangle = \frac{t^2}{2}\left(2\bar{n}^2+\bar{n}\right), ~\langle\hat{n}_1\hat{n}_2\rangle = \frac{t^2}{2}\left(6\bar{n}^3+6\bar{n}^2+\bar{n}\right), \end{align} \tag{ F6 }$$

and

$$\begin{align} g^{\left(2\right)}_\mathrm{Bose-Einstein}\left(0\right) = \frac{6\bar{n}^2+6\bar{n}+1}{2\bar{n}^2+\bar{n}}, \end{align} \tag{ F7 }$$

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.

Zoom In Zoom Out Reset image size

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 γ. This loss process can be described by the master equation in the Lindblad form:

$$\begin{align} \dot{\rho} = \sum_{j = 1}^{2}\gamma \left(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}\right), \end{align} \tag{ G1 }$$

where $\hat{a}_{j}$ ($\hat{a}_{j}^{\dagger}$) is the annihilation (creation) operator for the j-th bosonic mode. Any single mode density matrix ρ 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 G.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⁻⁴ to 10⁻¹. For both systems, the infidelity increases with loss rate and decreases with more resources used.

Zoom In Zoom Out Reset image size

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 figure 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.