Stroboscopic high-order nonlinearity for quantum optomechanics

TL;DR

This paper introduces a stroboscopic method to induce high-order quantum nonlinearity in levitated optomechanical systems, enabling the creation of nonclassical states with potential applications in quantum information processing.

Contribution

The authors propose a novel stroboscopic approach to generate high-order quantum nonlinearity in mechanical oscillators, overcoming decoherence and free evolution challenges.

Findings

Numerical analysis confirms feasibility with experimental parameters.

Method produces non-Gaussian states with negative Wigner functions.

Proposed setup enables transformation of ground states into nonclassical states.

Abstract

High-order quantum nonlinearity is an important prerequisite for the advanced quantum technology leading to universal quantum processing with large information capacity of continuous variables. Levitated optomechanics, a field where motion of dielectric particles is driven by precisely controlled tweezer beams, is capable of attaining the required nonlinearity via engineered potential landscapes of mechanical motion. Importantly, to achieve nonlinear quantum effects, the evolution caused by the free motion of mechanics and thermal decoherence have to be suppressed. For this purpose, we devise a method of stroboscopic application of a highly nonlinear potential to a mechanical oscillator that leads to the motional quantum non-Gaussian states exhibiting nonclassical negative Wigner function and squeezing of a nonlinear combination of mechanical quadratures. We test the method numerically by analysing highly instable cubic potential with relevant experimental parameters of the levitated optomechanics, prove its feasibility within reach, and propose an experimental test. The method paves a road for unique experiments instantaneously transforming a ground state of mechanical oscillators to applicable nonclassical states by nonlinear optical force.

Figures

Equations

$$\hat{H}={\hat{H}}_{\mathsf{HO}}+\alpha (t)V(\hat{x}),{\hat{H}}_{\mathsf{HO}}=\frac{1}{4}{\Omega }_{\mathsf{m}}({\hat{x}}^2+{\hat{p}}^2),$$

$$\hat{U}(t+M{T}_{\mathsf{m}},t)={\left[\hat{U}(t+T_{\mathsf{m}},t)\right]}^M\approx {\left\{\exp [-\mathrm{i}{\hat{H}}_{\mathsf{HO}}{T}_{\mathsf{m}}]\exp [-\mathrm{i}V(\hat{x})\tau ]\right\}}^M=\exp [-\mathrm{i}MV(\hat{x})\tau ].$$

$$\hat{U}(t+T_{\mathsf{m}},t)={\mathrm{e}}^{-\mathrm{i}{\hat{H}}_{\mathsf{HO}}{T}_{\mathsf{m}}/2}{\mathrm{e}}^{+\mathrm{i}V(\hat{x})\tau }{\mathrm{e}}^{-\mathrm{i}{\hat{H}}_{\mathsf{HO}}{T}_{\mathsf{m}}/2}{\mathrm{e}}^{-\mathrm{i}V(\hat{x})\tau }={\mathrm{e}}^{-2\mathrm{i}V(\hat{x})\tau },$$

$$\hat{U}(t+M{T}_{\mathsf{m}},t)=\exp [-\mathrm{i}\cdot 2M\cdot V(\hat{x})\tau ].$$

$$\dot{x}={\Omega }_{\mathsf{m}}p;\dot{p}=-{\Omega }_{\mathsf{m}}x-{\eta }_{\mathsf{m}}p+\sqrt{2{\eta }_{\mathsf{m}}}\xi ,$$

$$\rho (\tau )=\mathcal{N}\left[\rho (0)\right]=\mathcal{D}\left[\hat{U}(\tau ,0)\rho (0){\hat{U}}^{†}(\tau ,0)\right].$$

$$|{\gamma }_V⟩\propto \int \mathrm{d}x{\mathrm{e}}^{\mathrm{i}V(x)\tau }|x⟩,$$

$$\rho (V,r,n_0)={\mathrm{e}}^{\mathrm{i}V(\hat{x})\tau }{\rho }_0(r,n_0){\mathrm{e}}^{-\mathrm{i}V(\hat{x})\tau }.$$

$${\rho }_0(r,n_0)=\hat{S}(r){\rho }_0(0,n_0){\hat{S}}^{†}(r)$$

$${\sigma }_3^{\rho }(\lambda )\equiv \mathrm{Tr}(\rho (\hat{p}-\lambda {\hat{x}}^2{)}^2)-{\left(\mathrm{Tr}\rho (\hat{p}-\lambda {\hat{x}}^2)\right)}^2=\frac{v_{\mathsf{th}}}{s^2}+2(\lambda -\gamma {)}^2(s^2{v}_{\mathsf{th}}{)}^2,$$

$${\sigma }_3^{\text{vac}}(\lambda )=1+2{\lambda }^2.$$

$$W(x,p)=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}y{\mathrm{e}}^{-\mathrm{i}py}⟨x+y\left|\rho \right|x-y⟩.$$

$$W_{\mathsf{CPS}}(x,p)\propto \mathrm{Ai}\left[{\left(\frac{4}{\gamma }\right)}^{1/3}\left(\gamma x^2-p\right)\right],$$

$$\hat{H}={\hat{H}}_{\mathsf{HO}}+V(\hat{x})={\hat{H}}_p+{\hat{H}}_x,$$

$$\hat{U}(t+\tau,t)=\Big{[}\exp[-\mathrm{i}(\hat{H}_{\mathsf{HO}}+V(\hat{x}))\tfrac{\tau}{N}]\Big{]}^{N}=\Big{[}\hat{U}_{\mathsf{HO}}(\tfrac{\tau}{N})\hat{U}_{\mathsf{NL}}(\tfrac{\tau}{N})+O\left(\left(\tfrac{\tau}{N}\right)^{2}\right)\Big{]}^{N},$$

$$W_{\mathsf{th}}(x,p;n_0,s)=\frac{\exp \left(-\frac{1}{2}\left[\frac{(x/s{)}^2+(ps{)}^2}{2n_0+1}\right]\right)}{2\pi (2n_0+1)}.$$

$$W(x,p)=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}y{\mathrm{e}}^{-\mathrm{i}py}⟨x+y\left|\rho \right|x-y⟩,$$

$$⟨q|{\mathrm{e}}^{-\mathrm{i}{\hat{H}}_q(\hat{q})t}\rho {\mathrm{e}}^{\mathrm{i}{\hat{H}}_q(\hat{q})t}|q^{\prime }⟩=⟨q\left|\rho \right|q^{\prime }⟩{\mathrm{e}}^{-\mathrm{i}(H_q(q)-H_{q^{\prime }}(q^{\prime }))t}.$$

$$⟨x|{\hat{U}}_{\mathsf{NL}}(\delta t)\rho {\hat{U}}_{\mathsf{NL}}^{†}(\delta t)|x^{\prime }⟩=⟨x|\rho |x^{\prime }⟩{\mathrm{e}}^{-\mathrm{i}[V(x)-V(x^{\prime })]\delta t}.$$

$$W(x,p;t+\delta t)=W(x\cos \theta -p\sin \theta ,p\cos \theta +x\sin \theta ;t).$$

$$W(x,p;t+\frac{\pi }{\Omega })=\underset{-\infty }{\overset{\infty }{\iint }}\mathrm{d}u\mathrm{d}v{W}_i(x-u,p-v;t)W_{\mathsf{B}}(u,v),$$

$$W_{\mathsf{B}}(u,v)=\frac{1}{2\pi {\sigma }_{\mathsf{th}}}\exp \left[-\frac{u^2+v^2}{2{\sigma }_{\mathsf{th}}}\right],$$

$$\rho (\tau )=\hat{U}(\tau ,0)\rho (0){\hat{U}}^{†}(\tau ,0),$$

$$\exp [-\mathrm{i}({\hat{H}}_1+{\hat{H}}_2)N\delta \tau ]\approx {\left(\exp [-\mathrm{i}{\hat{H}}_1\delta \tau ]\times \exp [-\mathrm{i}{\hat{H}}_2\delta \tau ]\right)}^N,$$

$$\text{Method 1:}{\hat{H}}_1={\hat{H}}_p,{\hat{H}}_2={\hat{H}}_x,$$

$$\text{Method 2:}{\hat{H}}_1={\hat{H}}_{\mathsf{HO}},{\hat{H}}_2=V(\hat{x}),$$

$$\hat{H}=\frac{\Omega_{\mathsf{m}}}{4}\Big{[}\hat{p}^{2}+(1-\alpha(t))\hat{x}^{2}\Big{]}+\alpha(t)V(\hat{x}).$$

$$\dot{x}={\Omega }_{\mathsf{m}}p;\dot{p}=-{\Omega }_{\mathsf{m}}x-{\eta }_{\mathsf{m}}p+\sqrt{2{\eta }_{\mathsf{m}}}\xi ,$$

$$x^{\prime }$$

$$p^{\prime }$$

$$\delta x({\tau }_{\mathsf{m}})$$

Full text

\setcapwidth

[c] \setcapindent0pt

\setkomafont

author

Stroboscopic high-order nonlinearity for quantum optomechanics

Andrey A. Rakhubovsky 111Corresponding author: andrey.rakhubovsky@gmail.com 1

  

Radim Filip222filip@optics.upol.cz 1

1Department of Optics, Palacký University, 17. Listopadu 12, 771 46 Olomouc, Czech Republic

{onecolabstract}

High-order quantum nonlinearity is an important prerequisite for the advanced quantum technology leading to universal quantum processing with large information capacity of continuous variables. Levitated optomechanics, a field where motion of dielectric particles is driven by precisely controlled tweezer beams, is capable of attaining the required nonlinearity via engineered potential landscapes of mechanical motion. Importantly, to achieve nonlinear quantum effects, the evolution caused by the free motion of mechanics and thermal decoherence have to be suppressed. For this purpose, we devise a method of stroboscopic application of a highly nonlinear potential to a mechanical oscillator that leads to the motional quantum non-Gaussian states exhibiting nonclassical negative Wigner function and squeezing of a nonlinear combination of mechanical quadratures. We test the method numerically by analysing highly instable cubic potential with relevant experimental parameters of the levitated optomechanics, prove its feasibility within reach, and propose an experimental test. The method paves a road for experiments instantaneously transforming a ground state of mechanical oscillators to applicable nonclassical states by nonlinear optical force.

1 Introduction

Quantum physics and technology with continuous variables (CVs) [1] has achieved noticeable progress recently. A potential advantage of CVs is the in principle unlimited energy and information capacity of single oscillator mode. In order to fully gain the benefits of CVs and to achieve universal quantum processing one requires an access to a nonlinear operation [2, 3], that is, at least a cubic potential. Additionally, the CV quantum information processing can be greatly simplified and stabilized if variable higher-order potentials are available [4]. Variability of nonlinear gates can also help to overcome limits for fault tolerance [5]. Nanomechanical systems profit from a straightforward feasible way to achieve the nonlinearity by inducing a controllable classical nonlinear force of electromagnetic nature on a linear mechanical oscillator  [6, 7, 8]. Such a nonlinear force needs to be fast, strong and controllable on demand to access different nonlinearities required for an efficient universal CV quantum processing. Therefore, the field of optomechanics [9, 10, 11, 12] is a promising candidate to provide the key element for the variable on-demand nonlinearity. Optomechanical systems have reached a truly quantum domain recently, demonstrating the effects ranging from the ground state cooling [13, 14] and squeezing [15, 16] of the mechanical motion to the entanglement of distant mechanical oscillators [17, 18]. Of particular interest are the levitated systems in which the potential landscape of the mechanical motion is provided by a highly developed device — an optical tweezer [19, 20, 21, 22]. Levitated systems have proved useful in force sensing [23, 24], studies of quantum thermodynamics [25, 26, 27], testing fundamental physics [28, 29, 30] and probing quantum gravity [31, 32]. From the technical point of view, the levitated systems have recently demonstrated noticeable progress in the controllability and engineering, particularly, cooling towards [33, 34, 35, 36] and eventually reaching the motional ground state [37]. Further theoretical studies of preparation of entangled states of levitated nanoparticles are underway [38, 39]. Besides the inherently nonlinear optomechanical interaction met in the standard bulk optomechanical systems the levitated ones possess the attractive possibility of engineering the nonlinear trapping potential [40, 41, 42, 43, 26, 44, 37]. Moreover, the trapping potentials can be made time-dependent and manipulated at rates exceeding the rate of mechanical decoherence and even the mechanical frequency [45]. In this manuscript, we assume a similar possibility to generate the non-linear potential for a mechanical motion and control it in a fast way (faster than the mechanical frequency). Our findings do not rely on the specific method of how the nonlinearity is created.

Here we propose a broadly applicable nonlinear stroboscopic method to achieve high-order nonlinearity in optomechanical systems with time- and space-variable external force. The method builds on the possibility to control the nonlinear part of the mechanical potential landscape and introduce it periodically, adjusted in time with the mechanical harmonic oscillations. Such periodic application inhibits the effect of the free motion and the restoring force terms in the Hamiltonian and allows approaching the state arising from the nonlinear potential only. This is achieved similarly to how a stroboscopic measurement enables a quantum non-demolition detection of displacement [46, 47]. To prove feasibility of the method, we theoretically investigate realistic dynamics of a levitated nanoparticle in presence of simultaneously a harmonic and a strong, stroboscopically applied, nonlinear potentials enabled by the engineering of the trapping beam. To run numerical simulations, we advance the theory of optomechanical systems beyond the covariance matrix formalism appropriate for Gaussian states. Using direct Fock-basis and Suzuki-Trotter [48] simulations we model the simultaneous action of the nonlinear potential and harmonic trap, and obtain the Wigner functions of the quantum motional states achievable in this system. We predict very nonclassical negative Wigner functions [49, 50] generated by highly nonlinear quantum-mechanical evolution for time shorter than one mechanical period. The oscillations of Wigner function reaching negative values, in accordance with estimates based on unitary dynamics, witness that the overall quantum state undergoes unitary transformation $\exp(\mathrm{i}V(\hat{x})\tau)$ sufficient for universal quantum processing [2]. To justify it, we prove a nonlinear combination of the canonical quadratures of the mechanical oscillator to be squeezed below the ground state variance which is an important prerequisite of this state being a resource for the measurement-based quantum computation [51, 52]. For numerical simulations, we focus our attention to realistic version of the key nonlinearity, namely the cubic one with $V(x)\propto x^{3}$, and find good agreement of the predictions based on experimentally feasible dynamics with the lossless and noiseless unitary approximation. The method allows straightforward extension to more complex nonlinear potentials which can be used for flexible generation of other resources for nonlinear gates and their applications [4, 52]. In comparison with simultaneously developed superconducting quantum circuits [53], an advantage of our approach stems from a much larger flexibility of nonlinear potentials. Stroboscopic driving of an optomechanical cavity in a linear regime was considered in [54] for the purpose of cooling and Gaussian squeezing of the mechanical mode.

2 Results

2.1 The nonlinear stroboscopic protocol

To implement the stroboscopic method, it is possible to versatilely use a levitated nanoparticle [45] with optical [6], electric [7] or magnetic [8] trapping. It is also possible to use a mirror equipped with a fully-optical spring [55], or a membrane with electrodes allowing its nonlinear actuation and driving [56]. In any of such systems, the mechanical mode can be posed into nonlinear potential $V(x)$, particularly, the cubic potential $V_{3}(x)\propto x^{3}$ for the pioneering test. In this manuscript we focus on the experimental parameters peculiar to the levitated nanoparticles [36, 37], although the principal results remain valid for the other systems as well. We also focus here solely on the evolution of the mechanical mode of the optomechanical cavity, assuming that the coupling to the optical cavity mode (blue in Fig. 1) is switched off.

The mechanical mode is a harmonic oscillator of eigenfrequency $\Omega_{\mathsf{m}}$, described by position and momentum quadratures, respectively, $\hat{x}$ and $\hat{p}$, normalized such that $\left[\hat{x},\hat{p}\right]=2\mathrm{i}$. The oscillator is coupled to a thermal bath at rate $\eta_{\mathsf{m}}$. We also assume fast stroboscopic application of an external nonlinear potential $\alpha(t)V(x)$ with a piecewise constant $\alpha(t)$ illustrating the possibility to periodically switch the nonlinear potential on and off as depicted at Fig. 1(a). The Hamiltonian of the system, therefore, reads ($\hbar=1$)

[TABLE]

To illustrate the key idea behind the stroboscopic method, we first examine the regime of absent mechanical damping and decoherence. In this case, the unitary evolution of the oscillator is given by $\rho(t)=\hat{U}(t,t_{0})\rho(t_{0})\hat{U}^{\dagger}(t,t_{0})$, with $\hat{U}(t+\delta t,t)=\exp[-\mathrm{i}\hat{H}\delta t]$. When the nonlinearity is switched on permanently ($\alpha(t)=\alpha_{0}$), the free evolution dictated by $H_{\mathsf{HO}}$ mixes the quadratures of the oscillator which prohibits the resulting state from possessing properties of the target nonlinear quantum state arising purely from $V(x)$ regardless of the nonlinearity strength (see Supplementary Note 4 for more details). Willing to obtain a unitary transformation as close to $\exp[-\mathrm{i}V(\hat{x})\tau]$ as possible despite a constant presence of $\hat{H}_{\mathsf{HO}}$, we assume that the nonlinearity is repeatedly switched on for infinitesimally short intervals of duration $\tau$. If the duration $\tau$ is sufficiently short for the harmonic evolution to be negligible, the resulting evolution during $\tau$ is approximately purely caused by the $V(x)$ part of the Hamiltonian. To enhance the magnitude of the effect of the nonlinear potential, $V(\hat{x})$ we can apply it every $2\pi/\Omega_{\mathsf{m}}$ for short enough intervals as shown in Fig. 1 (b,c). This allows to establish an effective rotating frame within which the nonlinearity is protected from the effect of harmonic evolution. Realistically, the stroboscopic application corresponds to $\alpha(t)=\sum_{k}\delta_{\tau}(t-2\pi k/\Omega_{\mathsf{m}}),$ with $k\in\mathbb{Z}$, where $\delta_{\tau}$ is a physical approximation of Dirac delta function with width $\tau$ much shorter than the period of mechanical oscillations: $\tau\Omega_{\mathsf{m}}\ll 1$. Then we can consider the evolution over a number $M$ of harmonic oscillations as consisting of subsequent either purely harmonic or purely nonlinear steps, and the evolution operator can be approximately written as

[TABLE]

For the last equality we used the fact that the unitary harmonic evolution through a single period of oscillations is an identity map: $\hat{U}_{\mathsf{HO}}(t+T_{\mathsf{m}},t)\equiv\exp[-\mathrm{i}\hat{H}_{\mathsf{HO}}T_{\mathsf{m}}]~{}=~{}\hat{\mathds{1}}$. Motion of a real mechanical oscillator can be approximated by a harmonic unitary evolution with good precision because of very high quality of mechanical modes of optomechanical systems [57, 35, 36]. Eq. 2 shows that the effect of sufficiently short pulses of the strong nonlinear potential timed to be turned on precisely once per a period of mechanical oscillations $M$ times, is equivalent to an $M$-fold increase of the nonlinearity.

A further improvement is possible by noting that undamped harmonic evolution over a half period simply flips the sign of the two quadratures $(\hat{x},\hat{p})\mapsto-(\hat{x},\hat{p})$. Therefore, it is possible to similarly apply the potential twice per period, switching its sign each second time. This can be formalized as setting $\alpha(t)=\sum_{k}(-1)^{k}\delta_{\tau}(t-\pi k/\Omega_{\mathsf{m}}),\ k~{}\in~{}\mathbb{Z}$. In this case,

[TABLE]

and therefore, after $M$ periods

[TABLE]

The idealised scheme proposed above in reality faces two potentially deteriorative factors: the finiteness of the duration of the nonlinearity $\tau$, and the mechanical decoherence caused by the thermal environment. We take a proper account of these two factors by considering the evolution as consisting of two kinds: (i) unitary undamped dynamics in a sum of the quadratic and nonlinear potentials, (ii) damped harmonic evolution between those. These two kinds of evolution are subsequently repeated, as shown in Fig. 1 (b,c). We develop an advanced method based on Suzuki-Trotter simulation (STS) to simulate the quantum state of a realistic optomechanical system after the application of our proposed protocol. It is worth noting that the STS is typically used to approximately achieve a novel evolution $\hat{U}$ from experimentally available exact unitaries [58, 59]. In our work, we use STS to simulate and justify that the experimentally available compound evolution $\hat{U}$ can approximate one of its building blocks $\hat{U}_{\mathsf{NL}}(\delta t)\equiv\exp(-\mathrm{i}V(x)\delta t)$. To verify the convergence of STS we also perform simulations in the Fock-state basis that allow direct computation of the propagator corresponding to the Hamiltonian (1). Fock-state-basis simulations unfortunately do not grant access to phase-space distributions such as Wigner function [60] which makes use of STS the primary strategy. Excellent agreement between the very distinct methods (STS and Fock-state-basis simulations) indicates that our results are correct. The details of the simulation methods are presented in Supplementary Notes 1 and 2.

We omit damping and thermal decoherence during the fast unitary action of the combined potential. Such applications are assumed to happen every half of a mechanical oscillation and have durations much shorter than the mechanical period, therefore, due to high quality of state-of-the-art mechanical oscillators, such omission is justified. The joint action is simulated using STS and verified by Fock-state-basis simulation. Between the applications of the nonlinear potential $V(x)$ the mechanical oscillator experiences damped harmonic evolution described by the linear Heisenberg-Langevin equations

[TABLE]

where $\xi$ is the quantum Langevin force, obeying $\left[\xi(t),x(t)\right]=\mathrm{i}\sqrt{2\eta_{\mathsf{m}}}$ and $\frac{1}{2}\left\langle\left\{\xi(t),\xi(t^{\prime})\right\}\right\rangle=(2n_{\mathsf{th}}+1)\delta(t-t^{\prime})$ with $n_{\mathsf{th}}$ being the mean occupation of the bath. The experimentally accessible value of the heating rate $H_{\mathsf{m}}$ is given by $H_{\mathsf{m}}=\eta_{\mathsf{m}}n_{\mathsf{th}}$. The density matrix $\rho(T_{\mathsf{m}}/2)$ of the particle after half of a period of oscillations including the action of the nonlinear potential and subsequent damping can be formally written as

[TABLE]

where $\hat{U}(\tau,0)$ describes the particle’s unitary dynamics in combined potential, and $\mathcal{D}[\ \bullet\ ]$ indicates the mapping performed by the damping. Generalization of (6) to include the second half of the period, and the subsequent generalization to multiple periods, is straightforward.

Using these advanced numerical tools, further elaborated in Section 4 we evaluate the quantum state of the mechanical oscillator after the protocol and explore the limits of the achievable nonlinearities in the optomechanical systems that are accessible now or are within reach.

2.2 Application to the cubic nonlinearity

Motivated by the role of cubic nonlinearity for the universal continuous-variable quantum processing, we illustrate the devised method by numerically evaluating the evolution of a levitated particle under a stroboscopic application of a cubic potential $V(x)\propto x^{3}$. The nonlinear phase gate $\mathrm{e}^{-\mathrm{i}V(x)}$ is a limit case of motion, it modifies only momentum of the object without any change of its position. This nondemolition aspect is crucial for use in universal quantum processing. A nonlinear phase state (particularly, the cubic phase state introduced in [3]) as the outcome of evolution of the momentum eigenstate $\left|p=0\right\rangle$ in a nonlinear potential $V(x)$ is defined as

[TABLE]

where $V(x)$, is a highly nonlinear potential, and $\left|x\right\rangle$ the position eigenstate $\hat{x}\left|x\right\rangle=x$. The state (7) requires an infinite squeezing possessed by the ideal momentum eigenstate before the nonlinear potential is applied. More physical is an approximation of this state obtained from a finitely squeezed thermal state $\rho_{0}(r,n_{0})$, ideally, vacuum, by the application of $V$:

[TABLE]

The initial state $\rho_{0}$ is the result of squeezing a thermal state with initial occupation $n_{0}$

[TABLE]

where $\hat{S}(r)=\exp[\tfrac{1}{2}r^{*}(\hat{a})^{2}-\tfrac{1}{2}r(\hat{a}^{\dagger})^{2}]$ is a squeezing operator ($\hat{a}=(\hat{x}+\mathrm{i}\hat{p})/2$), and the initial state $\rho_{0}$ is thermal with mean occupation $n_{0}$. Phase of the squeezing parameter $r=|r|\mathrm{e}^{\mathrm{i}\theta}$ determines the squeezing direction. When $n_{0}=0$, $|r|\to+\infty$, and $\theta=\pi$, the initial state is infinitely squeezed in $p$, and  Eq. 8 approaches the ideal cubic state (7).

The quantum state obtained as a result of the considered sequence of interactions approximates the ideal state given by Eq. 7. The quality of the approximation can be assessed by evaluating the variance of a nonlinear quadrature $p-\lambda x^{2}$, or the cuts of the Wigner functions of the states. A reduction in nonlinear quadrature variance below the vacuum is a necessary condition for application of these states in nonlinear circuits [51, 52]. On the other hand, the phase-space interference fringes of the Wigner function reaching negative values are a very sensitive witness of quantum non-Gaussianity of the states used in the recent experiments [61, 62, 63, 64]. Fidelity happens to be an improper measure of the success of the preparation of the quantum resource state [65] because it does not predict either applicability of these states as resources or their highly nonclassical aspects.

A noise reduction in the cubic phase gate can be a relevant first experimental test of the quality of our method. The approximate cubic state obtained from a squeezed thermal state (that is, the state (8) with $V(x)=\gamma x^{3}/(6\tau)$) should possess arbitrary high squeezing in the variable $p-\lambda x^{2}$ for $n_{0}=0$ given sufficient squeezing of the initial mechanical state. The state (8) obtained from a squeezed in momentum state, has the following variance of the nonlinear quadrature

[TABLE]

where $v_{\mathsf{th}}=2n_{0}+1$ is the variance of each canonical quadrature in the initial thermal state before squeezing, and $s=\mathrm{e}^{r}$ is the magnitude of squeezing. An important threshold is the variance of the nonlinear quadrature attained at the vacuum state ($n_{0}=0$, $s=1$)

[TABLE]

Application of a unitary cubic evolution to the initial vacuum state displaces this curve along the $\lambda$ axis by $\gamma$. Further, as follows from Eq. 10, squeezing the initial state allows reducing the minimal value of $\sigma_{3}$, and increase of the initial occupation $n_{0}$ causes also increase of $\sigma_{3}$. Suppression of fluctuations in the nonlinear quadrature is a convenient figure of merit because it is a direct witness of the applicability of the quantum state as a resource for measurement-based quantum information processing [51, 52] and a witness of non-classicallity [66]. It can be evaluated in an optomechanical systems with feasible parameters using pulsed optomechanical interaction [66] without a full quantum state tomography.

In Fig. 2 (a) we show the variance $\sigma_{3}$ at the instants $t=M_{\mathsf{T}}T_{\mathsf{m}}/2$. Each of the curves takes into account also the interaction with the thermal environment lasting $T_{\mathsf{m}}/2$ after each application of the nonlinearity. The heating rate parameter $H_{0}=4\pi H_{\mathsf{m}}/\Omega_{\mathsf{m}}$ assumes value $H_{0}=2\times 10^{-3}$. For an oscillator of eigenfrequency $\Omega_{\mathsf{m}}=2\pi\times$100\text{,}\mathrm{kHz}$$ and $Q=10^{6}$ is equivalent to occupation of the environment equal to $n_{\mathsf{th}}\approx 10^{7}$ phonons. This is the equilibrium occupation of such an oscillator at the temperature of $50\text{\,}\mathrm{K}$. A recent experiment of ground state cooling of a levitated nanoparticle [37] reported the heating rate corresponding to $H\approx 10^{2}H_{0}$. A proof of robustness of the method against such heating is in the Supplementary Note 3.

Thin lines with markers show the analytic curves defined by Eq. 10 for the corresponding values of $\gamma$. The good quantitative correspondence between the approximate states resulting from the realistic stroboscopic application of nonlinear potential and the analytic curves again proves the validity of the stroboscopic method. Importantly, each of the curves has an area where it lies below the corresponding ground-state level $\sigma_{3}^{\text{vac}}$. This means each of the corresponding states gives advantage over vacuum if used as ancilla for the cubic gate. The dashed lines show the simulated quantum states obtained from the same initial state by longer unitary evolution according to the full Hamiltonian from Eq. 1, that is $\mathrm{e}^{-\mathrm{i}\hat{H}n\tau}\rho\mathrm{e}^{\mathrm{i}\hat{H}n\tau}$ where $n=1,2,3$ for blue, yellow and red correspondingly. Further divergence from ideal than that of the ones corresponding to the stroboscopic method witnesses an advantage of the latter in generation of the resource for the measurement-based computation.

The stroboscopic application of a fixed limited gain nonlinearity, therefore, indeed allows amplification of nonlinearity in accordance with Eq. 4. Importantly, even despite requiring longer evolution in a noisy environment, the stroboscopic method allows better amplification than a unitary longer application of the nonlinearity in presence of the free evolution ($\propto\hat{p}^{2}$) and harmonic ($\propto\hat{x}^{2}$) terms in the Hamiltonian.

The non-Gaussian character of the prepared quantum state can be witnessed via its Wigner function $W(x,p)$ which for a quantum state $\rho$ reads [67]

[TABLE]

Wigner function (WF) shows a quasiprobability distribution over the phase space spanned by position $x$ and momentum $p$ and its negativity is a prerequisite of the non-classicallity of a quantum state. In Fig. 2 (b) are Wigner functions $W(0,p)$ of the mechanical oscillator computed for the same approximate states as in the panel (a). The Wigner function of an ideal cubic phase state, i.e. the state given by Eq. 7 for $V(x)=x^{3}\gamma/(6\tau)$, reads [3, 50]

[TABLE]

where $\operatorname*{Ai}(x)$ is the Airy function. This function with apparent non-Gaussian shape in the phase space exhibits fast oscillations reaching negative values in the positive momentum for any $\gamma>0$. The Wigner functions of the states obtained by application of the stroboscopic protocol approach the one of Eq. 8. Each of them exhibits areas with negative values which proves quantum non-Gaussian character of the resulting state. Moreover, with increased number of stroboscopic applications involved, the resulting approximate state corresponds to stronger nonlinearity. For the last curve where $M_{\mathsf{T}}=3$ we also show the result of an idealized instantaneous unitary application of an equal total nonlinearity by a cyan line with markers which is indistinguishable from the line for the approximate evolution. For this pair of curves we also have an estimate for the overlap $\operatorname*{\operatorname{Tr}}[\rho_{\mathsf{red}}\rho_{\mathsf{cyan}}]/\operatorname*{\operatorname{Tr}}[\rho_{\mathsf{cyan}}^{2}]=0.9877$. Despite an excellent overlap of the ideal and approximate Wigner functions, the important resource, squeezing $\sigma_{3}$ shows a noticeable deviation from the ideal scenario. It is for this reason that we choose the squeezing $\sigma_{3}$ as the main figure of merit for the protocol.

3 Discussion

In this article we have proposed and theoretically analyzed a protocol to create a nonlinear motional state of the mechanical mode of an optomechanical system with controllable nonlinear mechanical potential. The method uses the possibility to apply the nonlinear potential to the motion of mechanical object in a stroboscopic way, twice per a period of oscillations. This way of application allows reducing the deteriorative effect of the free oscillations and approach the effect of pure action of nonlinear potential. In contrast to other methods of creating non-classical states by a continuous evolution in presence of nonlinear terms in the Hamiltonian [68, 69, 70], our method allows approaching the states that approximate evolution according to the unitary $\mathrm{e}^{\mathrm{i}V(x)}$ where $V(x)$ is the nonlinear potential profile. We tested our method on a cubic nonlinearity $\propto\hat{x}^{3}$ though the method is applicable to a wide variety of nonlinearities. Our simulations prove that application of the protocol allows one to obtain the squeezing in the nonlinear quadrature below the shot-noise level even if the initial state of the particle is not pure. The nonlinear state created in a stroboscopic protocol clearly outperforms as a resource the vacuum for which the bound (11) holds. Moreover, the stroboscopic states approximate the one defined by Eq. 8 obtained in absence of the free rotation and thermal decoherence. We also verify that the stroboscopic application of the cubic nonlinear potential generates nonclassical states under conditions that are further from optimal than the ones of Fig. 2. In particular, we find that the heating rate $H_{\mathsf{m}}$ can be increased approximately 100-fold until the curve $\sigma_{3}(\lambda)$ fails to overcome the vacuum boundary $1+2\lambda^{2}$. We are able to prove that when the duration of the application $\tau$ is increased tenfold such that the corresponding product $V(x)\tau$ remains constant, (that is, a less stiff potential is applied stroboscopically for proportionally longer time intervals), the resulting nonlinear state as well shows squeezing in $\sigma_{3}$. This proves robustness of the method to the two major imperfections.

We have shown the method to work for the parameters inspired by recent results demonstrated by the levitated optomechanical systems [71, 72]. The optical trap with a cubic potential has been already used in the experiments [43, 44]. Levitated systems [22], including electromechanical systems [73], have recently shown significant progress in the motional state cooling [35, 36, 37] and feedback-enhanced operation [34] which lays solid groundwork of the success of the proposed protocol. The authors of Ref. [44] report the experimental realisation of the potential $V(X)=\mu X^{3}/6$ with $\mu\approx 8k_{\mathsf{B}}T$\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-3}, where $X=x\sqrt{\hbar/(2m\Omega_{\mathsf{m}})}$ is the dimensional displacement of the oscillator, $m$ is its mass, $k_{\mathsf{B}}$ is the Boltzmann constant and $T$ temperature. From this value we can make a very approximate estimate for the nonlinear gain $\gamma=\mu(\hbar/(2m\Omega_{\mathsf{m}}))^{3/2}\tau\approx 1.2\times 10^{-3}$ assuming duration $\tau=\pi/(50\Omega_{\mathsf{m}})$, temperature $T=$300\text{\,}\mathrm{K}, frequency $\Omega_{\mathsf{m}}=2\pi\times$1\text{,}\mathrm{kHz}, and a mass $m=$4\text{\times}{10}^{-15}\text{\,}\mathrm{g} of a silica nanoparticle of $70\text{\,}\mathrm{nm}$ radius.

Experimental implementation of the proposed method can guarantee preparation of a strongly nongaussian quantum motional state. Further analysis of such a state will require either a full state tomography or better suited well-tailored methods to prove the nonclassicality [74]. An analysis of the estimation of the non-linear mechanical quadrature variance via pulsed optomechanical interaction can be found in [66]. The optical read out can be improved using squeezed states of light [75]. This experimental step will open applications of the proposed method to other nonlinear potentials relevant for quantum computation [4, 76, 51, 52], quantum thermodynamics [77, 78] and quantum force sensing [79, 80].

In our simulations we focused solely on the dynamics of the mechanical mode and assumed the conventional optomechanical interaction absent. This interaction, well developed in recent years, provides a sufficiently rich toolbox that allows incorporation of the mechanical mode into the optical circuits of choice [9]. As an option, a prepared nonlinear state can be transferred to traveling light pulse [75] using optomechanical driving. In a more complicated scenario, one can add optomechanical interaction to the stroboscopic evolution to obtain even richer dynamics. A complete investigation of such dynamics, however, goes beyond the scope of the present research focused on the preparation of nonlinear motional states.

In parallel with the experimental verification, the stroboscopic method can be used to analyse other higher-order mechanical nonlinearities such as $V(x)\propto x^{4}$ or tilted double-well potentials required for tests of recently disputed quantum Landauer principle [81], counter-intuitive Fock state preparation [82] and approaching macroscopic superpositions [83].

4 Methods

4.1 Tools of numerical simulations

The Eq. 6 approximates the damped evolution of an oscillator in a nonlinear potential by a sequence of individual stages of harmonic, nonlinear and damped harmonic evolution (see Fig. 1). Below we elaborate on how it is possible to simulate such dynamics using Wigner function in the phase space, density matrix in position, momentum and Fock basis.

First, we evaluate the action of the nonlinearity during the stroboscopic pulse. While the nonlinearity is on, the Hamiltonian of the system reads

[TABLE]

where $\hat{H}_{p}=\frac{\Omega_{\mathsf{m}}}{4}\hat{p}^{2},$ and $\hat{H}_{x}=\frac{\Omega_{\mathsf{m}}}{4}\hat{x}^{2}+\frac{\gamma}{6}\hat{x}^{3}.$

Our important simulation tool is the Suzuki-Trotter simulation (STS, see Ref. [48]) for $\hat{U}$

[TABLE]

where $\hat{U}_{\mathsf{HO}}(\delta t)\equiv\exp(-\mathrm{i}H_{\mathsf{HO}}\delta t)$, $\hat{U}_{\mathsf{NL}}(\delta t)\equiv\exp(-\mathrm{i}V(x)\delta t)$, $N$ is called the Trotter number. The accuracy of the approximation is, thereby, increasing with decreasing $\tau/N$. Despite $\tau$ being now sufficiently large to take into account noticeable free rotation through an angle $\Omega_{\mathsf{m}}\tau$ in the phase space, we still assume that $\tau$ is much shorter than the mechanical decoherence timescale, set by the heating rate $H_{\mathsf{m}}$. This is well justified for the current experiments [33, 36, 37], also see Supplementary Note 2. The STS requires the summands forming the Hamiltonian to be self-adjoint which is not always the case of $V(x)$, in particular $V(x)\propto x^{3}$. We take the necessary precautions by considering such nonlinearities over only short time in a finite region of the phase space. Thus we cautiously take care of the quantum motional state being limited to this finite region. To further verify the correctness of numerics via STS, we cross-check it using numerical simulations in Fock-state basis.

Eq. 14 shows the two possibilities to split the full Hamiltonian into summands to use the STS. We use these two possibilities to compute independently the mechanical state in order to verify the correctness of the STS in Supplementary Note 1.

First, we start from a squeezed thermal state $\rho(0)$, which has a representation by the Wigner function in the phase space

[TABLE]

The Wigner function corresponding to a quantum state $\rho$ is defined [67] as

[TABLE]

and the corresponding density matrix element can be obtained from the Wigner function by an inverse Fourier transform. It is therefore possible to extend this approach to any $W(x,p)$ beyond the Gaussian states.

The evolution of a state $\rho$ under action of a Hamiltonian proportional to a quadrature $\hat{q}$ can be straightforwardly computed in the basis of this quadrature, where it amounts to multiplication of density matrix elements with $c$-numbers:

[TABLE]

In particular, the nonlinear evolution reads

[TABLE]

The undamped harmonic evolution driven by $\hat{H}_{\mathsf{HO}}$ can be represented by the rotation of Wigner function (WF) in the phase space. A unitary rotation through an angle $\theta=\Omega_{\mathsf{m}}\delta t$ in the phase space maps the initial WF $W(x,p;t)$ onto the final $W(x,p;t+\delta t)$ as

[TABLE]

The unitary transformation of the density matrix can be as well computed in the Fock state basis.

Finally, damped harmonic evolution of a high-Q harmonic oscillator over one half of an oscillation can also be evaluated in the phase space as a convolution of the initial Wigner function $W(x,p;t)$ with a thermal kernel

[TABLE]

where the expression for the kernel reads

[TABLE]

with $\sigma_{\mathsf{th}}=(2n_{\mathsf{th}}+1){2\pi\eta_{\mathsf{m}}}/{\Omega_{\mathsf{m}}}$, where $n_{\mathsf{th}}\approx k_{\mathsf{B}}T/(\hbar\Omega_{\mathsf{m}})$ is the thermal occupation of the bath set by its temperature $T$. In terms of the heating rate $H_{\mathsf{m}}$, $\sigma_{\mathsf{th}}=4\pi H_{\mathsf{m}}/\Omega_{\mathsf{m}}$. Eq. 21 is obtained by solving the joint dynamics of the oscillator and bath followed by tracing out the latter. The detailed derivation of Eq. 21 can be found in Supplementary Note 2.

Using these techniques, one can evaluate the action of the map $\mathcal{N}$ defined by Eq. 6 on the state of the quantum oscillator. This yields the quantum state of the particle after one half of a mechanical oscillation. Repeatedly applying the same operations, one can obtain the state after multiple periods of the mechanical oscillations. Our purpose is then to explore the limits of the achievable nonlinearities in optomechanical systems that are accessible now or are within reach.

Data Availability

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

Acknowledgments

A.A.R. acknowledges the support of the project 20-16577S of the Czech Science Foundation and national funding from the MEYS under grant agreement No. 731473 and from the QUANTERA ERA-NET cofund in quantum technologies implemented within the European Union’s Horizon 2020 Programme (project TheBlinQC). R.F. acknowledges the project 21-13265X of the Czech Science Foundation.

Author Contributions

R.F. conceived the theoretical idea. A.A.R. performed the numerical computations. Both authors participated in writing the manuscript.

Competing Interests

The authors declare no competing interests.

Supplementary Note \fpeval5-4 Convergence of the Suzuki-Trotter approximation

In this section we numerically demonstrate that the Suzuki-Trotter approximation that we use to simulate the combined dynamics, indeed converges. For this purpose we provide the plots of the variance of nonlinear variable $\operatorname*{\operatorname{Var}}(p-\lambda x^{2})$ computed in the approximate nonlinear state

[TABLE]

where $\hat{U}(\tau,0)=\exp[-\mathrm{i}\hat{H}\tau]$ using different methods.

For the first two, we use the Suzuki-Trotter expansion of the unitary operator $\hat{U}$

[TABLE]

where, depending on the method, we choose either

[TABLE]

or

[TABLE]

with notations of Eq. 14. The action of individual unitary operators in Method 1 is computed in the corresponding quadrature basis according to Eq. 18. In Method 2 the action of $\hat{H}_{1}$ is simulated by rotation of the Wigner function in the phase space according to Eq. 20, and the action of $\hat{H}_{2}$ in the position basis as in Eq. 18. In the main text, Method 1 is used to produce the results.

For the Method 3, we write the Hamiltonian $\hat{H}$ in the Fock state basis. This allows to directly compute matrix elements of $\hat{U}(\tau,0)$ in the Fock basis as well, without the need to use STS, and directly obtain the approximate nonlinear state. This method, unfortunately, does not allow to obtain Wigner function directly.

Comparison of the nonlinear variances from different methods is in Supplementary Fig. S\fpeval3-2. We show that as the Trotter number $N$ increases, the simulation using Method 1 converges. The similarly converging simulation using Method 2 is not shown to avoid cluttering of the figure. Both methods approximately converge to the simulations in the Fock basis. Method 2 shows larger deviation since the operation of rotation of the Wigner function in the phase space involves interpolation and is therefore less accurate. This is the reason to choose Method 1 for the main text. Good numerical coincidence of the results of all the three different methods proves convergence of the STS.

As an experiment, we also simulate the regime of operation that is accessible only to the levitated nanoparticles. In such systems, one has in principle full control over the potential that the particle is exposed to, and can, therefore alternate between the quadratic trapping potential and the higher-order nonlinear one. Formally this corresponds to having a Hamiltonian that, instead of Eq. 1, reads

[TABLE]

We also simulate the action of the nonlinear stage of this Hamiltonian, that is the action of $\hat{H}=\Omega_{\mathsf{m}}\hat{p}^{2}/4+\gamma\hat{x}^{3}/6$, using Methods 1 and 3. That is we simulate this regime using STS and Fock-state expansion. The results of simulation show good agreement what again proves the STS convergence. The nonlinear quadrature is squeezed stronger in this regime, which suggests that the strategy to switch between the potentials can be advantageous for the levitated nanoparticles, or other systems where a full control over the potential is available.

Supplementary Note \fpeval6-4 Impact of the thermal noise

Owing to recent progress in design and manufacturing of the nanomechanical devices, the thermal noise is strongly suppressed. In this section we check the robustness of our scheme to the thermal noise from the environment and show that for the state-of-the art systems it is not hampering the quantum performance.

The thermal noise can be included by writing the standard Langevin equations:

[TABLE]

where $\xi$ is the quantum Langevin force, obeying the commutation relation $\left[\xi(t),x(t)\right]=\mathrm{i}\sqrt{2\eta_{\mathsf{m}}}$ and markovian autocorrelation $\frac{1}{2}\left\langle\left\{\xi(t),\xi(t^{\prime})\right\}\right\rangle=2n_{\mathsf{th}}+1$.

The solution of these equations corresponding to one half of a period of the mechanical oscillations ($\tau_{\mathsf{m}}=\pi/\Omega_{\mathsf{m}}$) for the experimentally relevant regime of high-Q mechanical oscillator ($Q\equiv\Omega_{\mathsf{m}}/\eta_{\mathsf{m}}\gg 1$) reads

[TABLE]

with $\zeta=\Omega_{\mathsf{m}}\tau_{\mathsf{m}}\sigma$, $\sigma=\sqrt{1-(\eta_{\mathsf{m}}/2\Omega_{\mathsf{m}})}$ and $\epsilon=\eta_{\mathsf{m}}/(2\Omega_{\mathsf{m}}\sigma)$. The definitions for the noise operators $\delta x,\delta p$ read

[TABLE]

These are the canonical quadratures of a mode in a thermal state with variance $\sigma_{\mathsf{th}}=2n_{\mathsf{th}}+1$.

The transformation of the Wigner function can be found as follows. Consider that at instant $t=0$ the mechanical oscillator has WF $W_{\mathsf{m}}(x(0),p(0))$. The mode of the bath is in a thermal state with WF

[TABLE]

Assuming the joint evolution of the mechanical oscillator and bath to be unitary, one can write the WF of the composite system (mechanical oscillator+bath)

[TABLE]

Here the first argument of $W_{\mathsf{m}}$, $x[\dots]$, means the solution of Eqs. S\fpeval29-22 and S\fpeval30-22 for $x$ etc. Since in the contemporary experiments the quality-factors of the mechanical oscillators can exceed $Q=10^{6}$ (see e.g. [57, 35, 36]), one can approximate $\sigma\approx 1$, $\epsilon=0$, $\zeta=\pi$ with accuracy $o(Q^{-1})$. Moreover, $\exp[-\eta_{\mathsf{m}}\tau_{\mathsf{m}}/2]\approx 1$. Therefore

[TABLE]

To obtain the WF of the mechanical oscillator after this evolution, one has to trace out the degrees of freedom of the environment

[TABLE]

Making a substitution $(u,v)=\theta\cdot(\delta x,\delta p)$ we arrive to the simple expression

[TABLE]

where $W^{\prime}_{\mathsf{B}}(u,v)=W_{\mathsf{B}}(u,v,\sigma_{\mathsf{th}}\theta^{2})$. The Eq. S\fpeval37-22 describes a convolution of the initial Wigner function $W_{\mathsf{m}}$ with a WF of a thermal state, whose variance is reduced by the mechanical $Q$-factor. Due to this rescaling, $W^{\prime}_{\mathsf{B}}$ is a very narrow Gaussian distribution, with variance much below $1$. In this manuscript we use primarily the value $\sigma_{\mathsf{th}}\theta^{2}=4\pi\frac{H_{\mathsf{m}}}{\Omega_{\mathsf{m}}}\equiv H_{0}=0.002$ which, for an oscillator of eigenfrequency $\Omega_{\mathsf{m}}=2\pi\times$100\text{,}\mathrm{kHz}$$ and $Q=10^{6}$ is equivalent to occupation of the environment equal to $n_{\mathsf{th}}\approx 10^{7}$ phonons. This is the equilibrium occupation of such an oscillator at the temperature of $50\text{\,}\mathrm{K}$.

Supplementary Note \fpeval7-4 Robustness of the stroboscopic method

In this section, to prove the feasibility of the stroboscopic method, we evaluate the robustness of our method against the finiteness (non-instantaneousness) of the duration of the nonlinearity and the impact of the thermal noise. It is convenient to measure the temporal extent $\tau$ of the nonlinear potential application in units of the mechanical motion period $\Theta=\Omega_{\mathsf{m}}\tau$, so that $\Theta=2\pi$ corresponds to a full mechanical oscillation. In Supplementary Fig. S\fpeval4-2 (a) we plot the variance of nonlinear quadrature as a function of the free parameter for different values of the phase rotation $\Theta$. The plots illustrate that as the phase gain $\Theta$ from the quadratic evolution decreases, the result of the combined dynamics approaches the action of the pure nonlinearity. In particular, for the value in the figure, difference in the minimal value of the nonlinear variance is less than 10%. As the strength of the applied nonlinear potential increases, the nonlinear squeezing becomes more sensitive to the free rotation. The reason is apparently because the state with stronger nonlinearity is further from the classical domain and hence more fragile to perturbations.

In Supplementary Fig. S\fpeval4-2 (b) we show the nonlinear variance from Supplementary Fig. S\fpeval4-2 for the value $\Theta=\pi/20$ after one half of mechanical period in the thermal environment. We can see that in the setup of the present experiment with levitated nanoparticles [37], for which $H=10^{2}H_{0}$ most of the nonlinear squeezing is lost, however for a narrow region of the parameter values the mechanical state can still be used as a resource of the nonlinear squeezing. In case of the heating rate being suppressed tenfold, the nonlinear squeezing persists in the thermal environment. It has to be noted that in electromechanical or optomechanical experiments with bulk oscillators, typical heating rates are much smaller. For instance, in electromechanical experiment [85], $H_{\mathsf{m}}=1.4\times 10^{-3}\Omega_{\mathsf{m}}$, and in the experiment with optomechanical crystal reported in Ref. [17], $H_{\mathsf{m}}<0.1\times 10^{-3}\Omega_{\mathsf{m}}$.

Supplementary Note \fpeval8-4 On duration of the nonlinear potential application

One possibility to use the nonlinear potential that seems straightforward, is to apply it simultaneously with the quadratic one and investigate the steady state. This strategy, unfortunately does not allow seeing coherent effects of nonlinearity. In this section we explore application of nonlinear potential for durations comparable to the mechanical period and show that the stroboscopic application is optimal.

First, we have to note that a mechanical oscillator in a sum of quadratic and an odd-order potentials will have a steady state only in the case of weak cubic nonlinearity when there exists a local minimum. Then the oscillator is going to be trapped in its vicinity where the effective potential is a displaced quadratic one. Therefore, we only have to consider finite-time action of the cubic potential and show that increasing the duration of nonlinear potential’s action decreases its effect.

The unitary dynamics of an oscillator in presence of both quadratic and polynomial potential is described by the operator which in our notations for the nonlinear gain $\gamma$ can be written as

[TABLE]

The resulting action of this unitary is determined by two parameters: total phase rotation $\Omega\tau$ and nonlinear gain $\gamma$. In the experiment, however, the parameters that are defined by the external control (e.g., optical driving) are the stiffnesses of quadratic and polynomial potentials, respectively, $\Omega$ and $\varkappa=\gamma/\tau$.

In this light, there are two different ways to compare longer interactions:

• •

Assuming constant stiffnesses but longer duration

[TABLE]

One could expect improvement in nonclassicality generation from this strategy as longer durations here cause larger gains of nonlinearity. Surprisingly, after a certain value of the gain, stronger nonlinearities are cancelled by the quadratures mixing caused by free rotation. Moreover, in an experiment, too high gain can cause the partice to escape the trap, or damage to the setup in case of a bulk system. Indeed, the original stroboscopic application assumed very short strong nonlinearity that, if applied for significantly longer durations, can exceed reasonable bounds.

• •

Assuming longer interaction but with proportionally reduced nonlinear stiffness, that causes equal nonlinear gain

[TABLE]

This way, equal nonlinearity faces larger rotation in the phase space, so no improvement is expected.

These findings are summarized at Supplementary Fig. S\fpeval5-2 where squeezing of the nonlinear quadrature is plotted against the parameter $\lambda$. The nonlinear quadrature is computed with respect to the state obtained from the initial squeezed thermal state, same as in the main text, by application one of the unitaries given by (S\fpeval39-22) or (S\fpeval40-22) assuming cubic nonlinearity ($k=3$).

Rather expectedly, if the duration of the nonlinearity is increased while preserving the total gain, increase of duration does not increase the nonlinear squeezing. Eventually the nonclassical effects as negative Wigner function and nonlinear squeezing from nonlinearity become insignificant.

If the nonlinearity of constant stiffness is applied for longer times, one could, in principle, expect stronger nonclassical effects. Unfortunately, this is not the case, and longer joint application of quadratic and cubic potential results in the mechanical mode being driven to the thermal state with increased occupation. Importantly, increasing the stiffness of the nonlinearity keeping the duration intact only increases the occupation of the resuting noisy state. That is, the dephasing action of the free rotation causes the nonlinearity to produce only noise in mechanics even for stiff nonlinear potentials. In conclusion, our simulations show that the presence of even a weak quadratic term in the potential makes the effect of the cubic nonlinearity vanish entirely.

Supplementary Note \fpeval9-4 Expectations of quadrature products $\left\langle\hat{x}^{n}\hat{p}^{m}\right\rangle$

In this section we describe a way to compute the expectation of an operator $\hat{x}^{n}\hat{p}^{m}$. This material can as well be found in textbooks such as [67, 84].

Following Weyl quantization rules [86], one can connect the classical expectation $\left\langle x^{n}y^{m}\right\rangle_{\mathsf{cl}}$ with the expectation of a combination of the operators $\hat{x}^{r}\hat{p}^{m}\hat{x}^{n-r}$:

[TABLE]

In this equation the classical and quantum expectations $\left\langle\bullet\right\rangle_{\mathsf{cl}}$ and $\left\langle\bullet\right\rangle$, respectively, are computed as

[TABLE]

where $W$ is the Wigner function corresponding to the state $\rho$.

Eq. S\fpeval41-22 is not particularly useful when one requires to compute a product of the form $\left\langle\hat{x}^{n}\hat{p}^{m}\right\rangle$, e.g., $\left\langle\hat{x}^{2}\hat{p}\right\rangle$. The efficient computation can be done recurrently.

First, note that

[TABLE]

therefore

[TABLE]

Importantly, all the terms with commutators in this expression simplify to expressions whose total power is smaller than $n+m$. This allows to write a recurrent relation

[TABLE]

The upper limit of the sum is decreased because for $r=n$, $\left[\hat{p}^{m},\hat{x}^{n-r}\right]=\left[\hat{p}^{m},\mathds{1}\right]=0$.

The equation above connects the expectation of the quantum operator $\hat{x}^{n}\hat{p}^{m}$ with a classical value $\left\langle x^{n}p^{m}\right\rangle_{\mathsf{cl}}$ and a combination of expectations of operators of lower total power. Substituting Eq. S\fpeval46-22 into the right part of itself one at the end arrives to an expression that consists entirely of classical expectations.

To illustrate this method, let us compute the variance of the nonlinear quadrature $\hat{p}-\lambda\hat{x}^{2}$ in the quantum state $\rho$, which involves computing the quantum expectation of $\hat{x}^{2}\hat{p}$.

First, we write the expression for the variance of the nonlinear quadrature

[TABLE]

The expectations of the powers of individual quadratures directly map to the classical expectations

[TABLE]

For the cross-terms, on the one hand, one can write a simple chain of equations

[TABLE]

For the purpose of illustrating Eq. S\fpeval46-22 we first reorder the cross-term to have the powers of $\hat{x}$ on the left

[TABLE]

The first term can be evaluated as

[TABLE]

Combining Eqs. S\fpeval50-22 and S\fpeval51-22 yields

[TABLE]

Supplementary Figures