Stroboscopic high-order nonlinearity for quantum optomechanics

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


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 a 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. In addition, 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-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 toward [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 26,37,[40][41][42][43][44] .
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 nonlinear 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 nondemolition detection of displacement 46,47 . To prove feasibility of the method, we theoretically investigate realistic dynamics of a levitated nanoparticle in the 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 (WFs) of the quantum motional states achievable in this system. We predict very nonclassical negative WFs 49,50 generated by highly nonlinear quantum-mechanical evolution for time shorter than one mechanical period. The oscillations of WF reaching negative values, in accordance with estimates based on unitary dynamics, witness that the overall quantum state undergoes unitary transformation expðiVðxÞτÞ 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 measurementbased 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) ∝ 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 ref. 54 for the purpose of cooling and Gaussian squeezing of the mechanical mode.

RESULTS
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) ∝ 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 Ω m , described by position and momentum quadratures, respectively,x andp, normalized such thatx;p ½ ¼ 2i. The oscillator is coupled to a thermal bath at rate η m . We also assume fast stroboscopic application of an external nonlinear potential α(t) V(x) with a piecewise constant α(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 (ћ = 1) 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 ρðtÞ ¼Ûðt; t 0 Þρðt 0 ÞÛ y ðt; t 0 Þ, withÛðt þ δt; tÞ ¼ exp½ÀiĤδt. When the nonlinearity is switched on permanently (α(t) = α 0 ), the free evolution dictated by H 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½ÀiVðxÞτ as possible despite a constant presence ofĤ HO , we assume that the nonlinearity is repeatedly switched on for infinitesimally short intervals of duration τ. If the duration τ is sufficiently short for the harmonic evolution to be negligible, the resulting evolution during τ is approximately purely caused by the V (x) part of the Hamiltonian. To enhance the magnitude of the effect of the nonlinear potential, VðxÞ we can apply it every 2π/Ω m for short enough intervals as shown in Fig. 1b, 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 where δ τ is a physical approximation of Dirac delta function with width τ much shorter than the period of mechanical oscillations: τΩ m ≪ 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 aŝ For the last equality we used the fact that the unitary harmonic evolution through a single period of oscillations is an identity map: U HO ðt þ T m ; tÞ exp½ÀiĤ HO T m ¼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 35,36,57 . Equation (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 ðx;pÞ 7 ! Àðx;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 αðtÞ ¼ P k ðÀ1Þ k δ τ ðt À πk=Ω m Þ; k 2 Z. In this case, Uðt þ T m ; tÞ ¼ e ÀiĤHOTm=2 e þiVðxÞτ e ÀiĤHOTm=2 e ÀiVðxÞτ ¼ e À2iVðxÞτ ; ( and therefore, after M periodŝ The idealized scheme proposed above in reality faces two potentially deteriorative factors: the finiteness of the duration of the nonlinearity τ, and the mechanical decoherence caused by the thermal environment. We take a proper account of these two Fig. 1 Scheme of the proposed stroboscopic method. a A levitated optomechanical system as an illustration of mechanical oscillator in a nonlinear potential. A dielectric subwavelength particle (P) is trapped by a tweezer (not shown). The particle feels a total potential U(x) = Ω m x 2 /4 + α(t)V(x) that is a sum of the quadratic (green) and the nonlinear (orange, here: cubic) parts, both provided by the trapping beam. The particle can be placed inside a high-Q cavity and probed by the laser light. b, c Stroboscopic application of nonlinear potential. The nonlinear part of the potential is switched on for only a short fraction of the mechanical period (orange segments). The quadratic trapping potential (green segments) is present throughout all the evolution. d Suzuki-Trotter simulation of stroboscopic evolution of the mechanical mode. In the figure, orange segments represent action of the nonlinear potential, empty and filled green segments correspond, respectively, to unitary and damped harmonic evolution.
A.A. Rakhubovsky and R. Filip 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. 1b, 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Û from experimentally available exact unitaries 58,59 . In our work, we use STS to simulate and justify that the experimentally available compound evolutionÛ can approximate one of its building blockŝ U NL ðδtÞ expðÀiVðxÞδ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 WF 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 where ξ is the quantum Langevin force, obeying ξðtÞ; xðtÞ ½ ¼ i ffiffiffiffiffiffiffiffi 2η m p and 1 2 ξðtÞ; ξðt 0 Þ f g h i¼ ð2n th þ 1Þδðt À t 0 Þ with n th being the mean occupation of the bath. The experimentally accessible value of the heating rate H m is given by H m = η m n th . The density matrix ρ(T 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 whereÛðτ; 0Þ describes the particle's unitary dynamics in combined potential, and D½ 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 "Methods" 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.
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) ∝ x 3 . The nonlinear phase gate e −iV(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 p ¼ 0 j iin a nonlinear potential V(x) is defined as where V(x), is a highly nonlinear potential, and x j i the position eigenstatex x j i ¼ 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 ρ 0 (r, n 0 ), ideally, vacuum, by the application of V: ρðV; r; n 0 Þ ¼ e iVðxÞτ ρ 0 ðr; n 0 Þe ÀiVðxÞτ : The initial state ρ 0 is the result of squeezing a thermal state with initial occupation n 0 ρ 0 ðr; n 0 Þ ¼ŜðrÞρ 0 ð0; n 0 ÞŜ y ðrÞ whereŜðrÞ ¼ exp½ 1 2 r Ã ðâÞ 2 À 1 2 rðâ y Þ 2 is a squeezing operator (â ¼ ðx þ ipÞ=2), and the initial state ρ 0 is thermal with mean occupation n 0 . Phase of the squeezing parameter r = |r|e iθ determines the squeezing direction. When n 0 = 0, |r| → +∞, and θ = π, 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 − λ x 2 , or the cuts of the WFs 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 WF 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) = γx 3 /(6τ)) should possess arbitrary high squeezing in the variable p − λ 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 where v th = 2n 0 + 1 is the variance of each canonical quadrature in the initial thermal state before squeezing, and s = 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) Application of a unitary cubic evolution to the initial vacuum state displaces this curve along the λ axis by γ. Further, as follows from Eq. (10), squeezing the initial state allows reducing the minimal value of σ 3 , and increase of the initial occupation n 0 causes also increase of σ 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 system with feasible parameters using pulsed optomechanical interaction 66 without a full quantum state tomography.
In Fig. 2a, we show the variance σ 3 at the instants t = M T T m /2. Each of the curves takes into account also the interaction with the thermal environment lasting T m /2 after each application of the nonlinearity. The heating rate parameter H 0 = 4πH m /Ω m assumes value H 0 = 2 × 10 −3 . For an oscillator of eigenfrequency Ω m = 2π × 100 kHz and Q = 10 6 is equivalent to occupation of the environment equal to n th ≈ 10 7 phonons. This is the equilibrium occupation of such an oscillator at the temperature of 50 K. A recent experiment of ground state cooling of a levitated nanoparticle 37 reported the heating rate corresponding to H ≈ 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 γ. 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 σ 3 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 e ÀiĤnτ ρe iĤnτ 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 (/p 2 ) and harmonic (/x 2 ) terms in the Hamiltonian.
The non-Gaussian character of the prepared quantum state can be witnessed via its WF W(x, p) which for a quantum state ρ reads 67 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. 2b are WFs W(0, p) of the mechanical oscillator computed for the same approximate states as in the panel (a). The WF of an ideal cubic phase state, i.e. the state given by Eq. (7) for V(x) = x 3 γ/(6τ), reads 3,50 where 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 γ > 0. The WFs 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 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 Tr½ρ red ρ cyan =Tr½ρ 2 cyan ¼ 0:9877. Despite an excellent overlap of the ideal and approximate WFs, the important resource, squeezing σ 3 shows a noticeable deviation from the ideal scenario. It is for this reason that we choose the squeezing σ 3 as the main figure of merit for the protocol.

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 nonclassical states by a continuous evolution in presence of nonlinear terms in the Hamiltonian 68-70 , our method allows approaching the states that approximate evolution according to the unitary e iV(x) where V (x) is the nonlinear potential profile. We tested our method on a cubic nonlinearity /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 m can be increased approximately 100-fold until the curve σ 3 (λ) fails to overcome the vacuum boundary 1 + 2λ 2 . We are able to prove that when the duration of the application τ is increased tenfold such that the corresponding product V(x)τ 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 σ 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 73 , including electromechanical systems 74 , 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 realization of the potential V( is the dimensional displacement of the oscillator, m is its mass, k B is the Boltzmann constant and T temperature. From this value we can make a very approximate estimate for the nonlinear gain γ = μ(ћ/(2mΩ m )) 3/2 τ ≈ 1.2 × 10 −3 assuming duration τ = π/(50Ω m ), temperature T = 300 K, frequency Ω m = 2π × 1 kHz, and a mass m = 4 × 10 −15 g of a silica nanoparticle of 70 nm radius.
Experimental implementation of the proposed method can guarantee preparation of a strongly non-Gaussian 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 75 . An analysis of the estimation of the nonlinear mechanical quadrature variance via pulsed optomechanical interaction can be found in ref. 66 . The optical read out can be improved using squeezed states of light 76 . This experimental step will open applications of the proposed method to other nonlinear potentials relevant for quantum computation 4,51,52,77 , quantum thermodynamics 78,79 and quantum force sensing 80,81 .
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 76 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 analyze other higher-order mechanical nonlinearities such as V(x) ∝ x 4 or tilted double-well potentials required for tests of recently disputed quantum Landauer principle 82 , counter-intuitive Fock state preparation 83 and approaching macroscopic superpositions [84][85][86][87] .

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 WF 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 readŝ whereĤ p ¼ Ωm Our important simulation tool is the STS (see ref. 48 ) forÛ whereÛ HO ðδtÞ expðÀiH HO δtÞ,Û NL ðδtÞ expðÀiVðxÞδtÞ, N is called the Trotter number. The accuracy of the approximation is, thereby, increasing with decreasing τ/N. Despite τ being now sufficiently large to take into account noticeable free rotation through an angle Ω m τ in the phase space, we still assume that τ is much shorter than the mechanical decoherence timescale, set by the heating rate H 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) ∝ 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 Fockstate basis. Equation (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 ρ(0), which has a representation by the WF in the phase space W th ðx; p; n 0 ; sÞ ¼ The WF corresponding to a quantum state ρ is defined 67 as Wðx; pÞ ¼ 1 2π dy e Àipy x þ yjρjx À y h i ; (17) and the corresponding density matrix element can be obtained from the WF 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 ρ under action of a Hamiltonian proportional to a quadratureq can be straightforwardly computed in the basis of this quadrature, where it amounts to multiplication of density matrix elements with c-numbers: q h je ÀiĤqðqÞt ρe iĤqðqÞt q 0 j i ¼ qjρjq 0 h ie ÀiðHqðqÞÀH q 0 ðq 0 ÞÞt : In particular, the nonlinear evolution reads x h jÛ NL ðδtÞρÛ y NL ðδtÞ x 0 j i ¼ x h jρ x 0 j ie Ài½VðxÞÀVðx 0 Þδt : The undamped harmonic evolution driven byĤ HO can be represented by the rotation of WF in the phase space. A unitary rotation through an angle θ = Ω m δt in the phase space maps the initial WF W(x, p; t) onto the final W(x, p; t + δt) as Wðx; p; t þ δtÞ ¼ Wðx cos θ À p sin θ; p cos θ þ x sin θ; tÞ: 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 WF W(x, p; t) with a thermal kernel du dv W i ðx À u; p À v; tÞW B ðu; vÞ; where the expression for the kernel reads with σ th = (2n th + 1)2πη m /Ω m , where n th ≈ k B T/(ћΩ m ) is the thermal occupation of the bath set by its temperature T. In terms of the heating rate Hm, σ th = 4πH m /Ω m . Equation (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 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 A.A. Rakhubovsky and R. Filip explore the limits of the achievable nonlinearities in optomechanical systems that are accessible now or are within reach.