Laser Theory for Optomechanics: Limit Cycles in the Quantum Regime

Optomechanical systems can exhibit self-sustained limit cycles where the quantum state of the mechanical resonator possesses nonclassical characteristics such as a strongly negative Wigner density, as was shown recently in a numerical study by Qian et al. [Physical Review Letters, 109, 253601 (2012)]. Here we derive a Fokker-Planck equation describing mechanical limit cycles in the quantum regime which correctly reproduces the numerically observed nonclassical features. The derivation starts from the standard optomechanical master equation, and is based on techniques borrowed from the laser theory due to Haake's and Lewenstein. We compare our analytical model with numerical solutions of the master equation based on Monte-Carlo simulations, and find very good agreement over a wide and so far unexplored regime of system parameters. As one main conclusion, we predict negative Wigner functions to be observable even for surprisingly classical parameters, i.e. outside the single-photon strong coupling regime, for strong cavity drive, and rather large limit cycle amplitudes. The general approach taken here provides a natural starting point for further studies of quantum effects in optomechanics.


I. INTRODUCTION
Optomechanical systems provide a test bed to study a broad range of paradigmatic quantum optical processes at so far unexplored meso-and macroscopic mass and length scales [1][2][3]. That quantum effects can play an important and even dominating role in the dynamics of these systems has been shown in a number of recent experiments demonstrating cooling to the quantum ground state [4,5], ponderomotive squeezing of light [6,7], back action noise limited position sensing [8,9], coherent state transfer [10], and entanglement [11].
In the most elementary optomechanical setup a single cavity mode couples to a single mechanical oscillator through, e.g., radiation pressure or dipole gradient forces. The dynamics of the system depends crucially on the frequency of the external driving field applied to the cavity: For the purpose of position or force sensing as in [8,9] the driving field is chosen resonant, while for back action cooling or state transfer the field is tuned below (to the red side of) the cavity frequency [4,5,10]. For blue detuning the system exhibits a rather complex nonlinear behavior. When the driving field is swept from the red to the blue side the nonlinear dynamics sets in as a parametric amplification process where phonons and photons are created correlated in pairs [12]. This lies at the heart of the recently reported generation of optomechanical entanglement [11]. The amplification will finally go over into a regime of self-sustained limit cycles due to the nonlinearity inherent to the optomechanical coupling. The classical dynamics in this regime has been observed experimentally [13][14][15][16][17][18] and is well studied theoretically [19][20][21][22]. Motivated by the impressive progress towards quantum effects in optomechanical systems also the quantum regime of optomechanical limit cycles received significant attention in theoretical studies [23][24][25][26][27][28][29].
In particular, a recent numerical study of the full optomechanical master equation in the limit cycle regime showed that the Wigner function of the mechanical oscillator can become strongly negative [27]: Negativities of the Wigner function occur for driving fields at the blue sidebands and -more pronounced -also for resonant drive. Limit cycle states with negative Wigner density even exist in regions of red detuning where a (simple) classical model would not predict limit cycles at all. The numerical findings were independently confirmed in [28]. This reference predicts negative Wigner density even on higher sidebands and compares the extend of negativity found for different detunings in more detail. In view of these findings it is important to strive for a deeper understanding of these effects and the underlying mechanisms on the basis of an appropriate analytical model.
The transition from parametric amplification to optomechanical limit cycles can be understood in analogy to the threshold behaviour of a laser (or maser) cavity [30][31][32] where the roles of the laser cavity and the laser medium are played by, respectively, the mechanical oscillator and the optomechanical cavity [33]. Along this line a semiclassical rate equation model was derived in [21,33] for optomechanical systems. Rodrigues and Armour [25,26] developed a quantum mechanical treatment employing a truncated Wigner function approach to derive a Fokker-Planck equation (FPE) for the mechanical oscillator. The FPE predicted in particular a sub-Poissonian, or number-squeezed, phonon statistics in the limit cycle when the driving field is blue detuned from the cavity resonance by the mechanical oscillation frequency.
In the present article we apply the laser theory due to Haake and Lewenstein [31,32] to describe optomechanical limit cycles in the quantum regime. Our model correctly reproduces the characteristics of limit cycles mentioned above. It identi-arXiv:1310.1298v2 [quant-ph] 10 Feb 2014 fies general requirements on system parameters (such as coupling strength, driving power, sideband resolution, temperature etc.) for the occurrence of sub-Poissonian phonon statistics and negative Wigner functions, and establishes a tight connection between the two phenomena. We find that negative Wigner functions can be achieved also in rather classical parameter regimes where the coupling per single photon g 0 is smaller than the cavity line width, and where the cavity is driven strongly and limit cycle amplitudes are large. The associated small Fano factors are lower bounded by, and can reach, the sideband parameter κ/ω m (ratio of cavity line width to mechanical resonance frequency) for sufficiently strong optomechanical cooperativity.
Starting from the standard optomechanical master equation [1,2] an effective FPE is derived for the quasi-probability distribution (such as e.g. the Wigner-, P-or Q-function) of the mechanical oscillator under adiabatic elimination of the cavity mode. The nonlinearity of the optomechanical interaction gives rise to nonlinear drift and diffusion coefficients in the FPE which describe, respectively, the (classical) nonlinear physics of limit cycles [19,20] and the impact of quantum noise of the cavity. The approach taken here permits to work in a picture which interpolates between the dressed state picture introduced in [34,35] through a polaron transformation and the bare state picture of the standard master equation [1,2,25,28,36]. Remarkably, in analogy to the polaron picture, this intermediate picture explicitly separates the optical Kerr-nonlinearity inherent to the radiation pressure from the optomechanical interaction. In contrast to the polaron picture, the interaction term is not removed from the master equation and both the mechanical oscillator and the cavity remain separate systems as in the standard master equation picture. The entanglement of cavity and oscillator in the polaron picture would complicate the study of them as separate systems as required in the context of limit cycles. As we will show, the novel treatment of the optomechanical Kerr nonlinearity presented in this article can become essential to understand the physics of limit cycles.
The effective FPE derived here exactly reproduces the one of Rodriguez and Armour [25,26] when neglecting the different description of the Kerr nonlinearity of the cavity, which is treated in the standard master equation picture there. In comparison to [25,26] our approach does not require truncation of higher order derivatives, and gives a consistent and natural account of the Kerr nonlinearity.
The article is organized as follows: In Sec. II we give an executive summary of the main results, as far as they relate to the appearance of nonclassical mechanical states. In Sec. III we introduce the main idea of Haake and Lewenstein's laser theory in the context of optomechanics, and apply it to derive the effective FPE for the mechanical oscillator. Sec. IV we discuss the implications of the FPE equation for optomechanical limit cycles in the quantum regime. In principle each of these Sections can be read independently. Readers who are interested only in one particular aspect are encouraged to jump directly to the respective section.

II. PREVIEW OF THE MAIN RESULTS
The aim of this section is to give a preview of our most important results, and to indicate how these results could be derived in a relatively simple (quantum noise) approach. The main idea is to find the width of the mechanical limit cycles in phase space and to deduce from that the spread in phonon numbers. For simplicity, we will assume here that the optomechanical interaction dominates (i.e. formally zero mechanical damping). The full optomechanical laser-theory analysis will go significantly beyond this, but it will reproduce the features discussed here.
In the following, we will find it convenient to characterize the optomechanical coupling in several ways: as the cavity frequency shift per displacement G, via the single-photon coupling strength g 0 = Gx ZPF , and via the dimensionless ratio η = 2g 0 /ω m . x ZPF = √ /2mω m is the zero-point amplitude of the mechanical oscillator with mass m and frequency ω m . We start by assuming mechanical oscillations at a fixed amplitude r such that x(t) = x ZPF Re[re −iω m t ]. At each instant of time, the radiation pressure force F = Ga † a (a is the photon annihilation operator) will feed energy into the mechanical oscillations at a rate (power) P = F(t)ẋ(t). Following the classical approach [19], we can predict the slow drift of the mechanical oscillation amplitude by calculating the average power input, P(t) . Here · denotes the quantum expectation value, while the overbar averages over a time-window comprising several oscillation periods. We note that the power balance equation is analogous to the loss-gain equations in a laser and that the laser analogy will be heavily used throughout the manuscript.
In addition to this drift, however, there will be diffusion of the mechanical oscillator's energy, due to the fundamental radiation pressure shot noise fluctuations. The energy diffusion constant is given by D E = 1 2 +∞ −∞ dτ δP(t + τ)δP(t) , where δP(t) = P(t) − P(t) denotes the fluctuations. In order to discuss the quantum dynamics of optomechanical limit cycles, it turns out to be crucial to study the behaviour of this diffusion constant as a function of cycle amplitude. In particular, we will show that the appearance of nonclassical mechanical states can only be understood by a rather subtle cancellation of some term that would usually dominate, leaving the diffusion constant small and leading to a narrowing of the phonon distribution by the sideband ratio κ/ω m (where κ denotes the cavity amplitude decay rate).
Our task of calculating this diffusion constant is complicated by the fact that we are dealing with shot noise inside an optical cavity whose resonance frequency oscillates sinusoidally. We thus have to solve the equation for the light field inside such a mechanically driven cavity, i.e. da/dt = 2κa in (t), where ∆ = ω L − ω c is the detuning between the laser at frequency ω L and the bare cavity resonance at ω c . The solution for a(t) can be expressed via the extra phase θ(t) = ηIm[re −iω m t ] accumulated in the cavity field due to the oscillations. It reads where is the standard cavity filter function (and θ = 0 recovers the usual case).
The light intensity oscillates at harmonics of the mechanical motion. α(t) = e iθ(t) a(t) = n α n e inω m t is the average cavity amplitude (modulo the phase), which can be obtained by evaluating Eq. (1). For a constant laser drive, with an amplitude √ 2κ a in ≡ E, we obtain α n = Ee −inφ J −n (ηr)/h n , where h n = κ + i(nω m − ∆). These are the Bessel amplitudes that also determine the appearance of multiple stable attractors in the classical analysis of the optomechanical instability [19]. These attractors can be found by noting that the drift of the amplitude r is governed by the power input P(t) , as the energy of the mechanical oscillator is given by mω 2 m r 2 /2. In the regime of interest here, this drift can be approximated aṡ The limit cycle amplitude is thus fixed at the zeroes of the Bessel function, in the absence of additional mechanical damping. This will be crucial further below.
In addition, there are the electromagnetic vacuum fluctuations δa in (t) = a in (t) − a in entering the cavity. In order to evaluate the mechanical energy diffusion constant that is governed by those fluctuations, we need the force-force correlator F(t)F(t ) , i.e. ultimately the shot-noise (irreducible) part of the photon number correlator. By using the vacuum noise correlator δa in (t)δa † in (0) = δ(t), we find directly For a constant α, this reduces to the shot noise correlator employed in the quantum noise approach to optomechanical cooling [37]. Now, we can proceed to evaluate the energy diffusion constant D E introduced above. The resulting slightly lengthy expression (Eq. (B24) in Appendix B) can be simplified in the regime of interest here to Here we have introduced D W as the diffusion constant for the amplitude r of the limit cycle. This amplitude is connected to the energy via E = ω m r 2 , such that one obtains the relation between D E and D W shown here. It is now crucial to observe that the diffusion constant has a minimum right at the first limit cycle. This is because the first contribution in Eq. (3), which dominates at smaller amplitudes, is completely suppressed at the limit cycle, where J 0 (ηr) = 0, see Eq. (2). Thus, only the second term survives, which is suppressed by a factor κ 2 /ω 2 m , i.e. the sideband ratio squared. We show in the main text that this suppression is caused by squeezing terms that exactly cancel the corresponding incoherent diffusion terms in leading order. Now we can combine these results to discuss the width σ 2 of the distribution in the amplitude r. The compromise between the diffusion at rate D W and the restoring force that drives r back to the limit cycle results in a width σ 2 = −D W /µ . For a fixed limit cycle amplitude, both diffusion and drift scale as g 2 0 E 2 , such that the laser power and the optomechanical coupling drop out of this expression. This will change in the presence of mechanical damping and thermal fluctuations, but it still correctly describes the behavior once the optomechanical damping rate overwhelms the thermal fluctuations.
In order to estimate when the limit cycle may turn into a nonclassical mechanical quantum state, we will now look at the variance of the phonon number Var(n). Since r is already measured in terms of the zero-point amplitude x ZPF , we have r 2 = n. Thus Var(n) = 4 n σ 2 . This can be minimized by choosing an optimal detuning (∆ = κ), where we find Var(n) = n (κ/ω m ). In other words, in the resolved sideband regime (κ ω m ), one can get close to a mechanical Fock state, Var(n) < 1, as long as the limit cycle is sufficiently small, n < ω m /κ. Note that the optomechanical coupling strength g 0 enters indirectly here, since (in the absence of mechanical damping) the limit cycle amplitude is determined by J 0 (ηr) = 0, with η = 2g 0 /ω m . Taking this into account, Var(n) < 1 is equivalent to g 2 0 /ω m κ > 1.4. However, it turns out that it is easier to produce a nonclassical state, i.e. one where the Wigner density has negative components. For the type of states relevant here, we numerically find that it is sufficient to have Var(n) < 0.6 n 0.7 for this purpose. Thus, the condition for nonclassicality reads approximately which is less stringent than the condition for achieving a Fock state, since one could still admit g 0 /κ < 1 if the sideband ratio ω m /κ is sufficiently large. In the simplified description given here, we have neglected several factors which will be discussed in our full analysis. This includes the effects of the mechanical damping, which will decrease the limit cycle amplitude (shifting away from the point of minimum diffusion constant). In addition, thermal fluctuations will add to the diffusion. Nevertheless, this effect can be overcome if the scale of the optically induced damping rate, γ opt µ (r) ∝ g 2 0 E 2 /ω 3 m , dominates the influx of thermal phonons wheren is the thermal phonon number of the bath, and γ is the mechanical damping rate. This is equivalent to the condition for ground state cooling, but here applied for the instable regime. It does not involve the coupling per single photon g 0 , but only the linearized coupling g ∝ g 0 E, such that Eq. (5) essentially represents a condition on the strength of the driving field. Another important aspect neglected here is the shift of the cavity resonance by the Kerr effect. This leads to an effective detuning ∆ eff that will enter all expressions instead of ∆. The impact of this change is especially large near ∆ ≈ 0, which is precisely the regime which we find to be optimal for nonclassical states.
The heuristic reasoning applied here and the resulting conditions (4) and (5) for achieving nonclassical mechanical states will receive a rigorous justification in Sec. IV on the basis of the Fokker-Planck equation derived in the next Section.

A. Haake-Lewenstein Laser Theory Ansatz in Optomechanics
Master Equation -The standard master equation of an optomechanical system is [1,2] where The three Liouvillians L m , L c , and L int refer to the mechanical oscillator, the cavity, and their interaction respectively. a and b denote the annihilation operators of the cavity and the mechanical oscillator. The frequency of the mechanical oscillator is ω m , its amplitude damping rate is γ = ω m /Q m , and its mean phonon number in thermal equilibriumn. We use the notation D[A]ρ = 2AρA † − A † Aρ − ρA † A for Lindblad operators. κ is the cavity amplitude decay rate, ∆ = ω L − ω c is the detuning from cavity resonance at ω c of the driving field E = √ 2κP L / ω L with power P L and frequency ω L . The master equation is written in a frame rotating at the frequency ω L of the driving field. The optomechanical coupling per single photon is denoted by g 0 , and essentially determines the dispersive shift of the cavity frequency with the displacement of the oscillator in units of the mechanical zero-point amplitude [38].
Note that in contrast to e.g. [19,23,27] the definitions of γ and κ used here refer to the decay rate of the amplitude and will be used for all analytical results, in order to make the equations more readable. The corresponding decay rates for the energy κ E = 2κ and γ E = 2γ are the standard convention from [1]. For comparison to most experimental and numerical studies, we provide also the energy decay rates in the numerical results.
Our primary aim is to derive an effective equation of motion for the mirror based on the assumption that the dynamics of the cavity adiabatically follows the mechanical oscillator. This will be strictly the case when the cavity decay rate κ is larger than the characteristic coupling strength of the oscillator and the cavity mode (i.e. g 0 or the linear coupling g = g 0 α enhanced by the mean cavity field α at the position of the limit cycle). As we will see, the resulting effective equation of motion for the mechanical oscillator gives good results for the stationary state also when this condition is fulfilled barely, and even when it is mildly violated.
Quasiprobability distribution -Most importantly, we will not assume the usual linearization of the optomechanical coupling when we perform the adiabatic elimination. This is achieved by means of an Ansatz inspired by laser theory [31,32], which allows us to use a different adiabatic reference state of the cavity field for each point in phase space of the mechanical oscillator. The idea is to switch to a phase-space representation for the mechanical degree of freedom. In principle any quasi-probability distribution (e.g. P-distribution, Wigner function etc.) can be used, but we will in the following mostly focus on the (Husimi) Q function which yields the simplest formulas for the calculation presented below. In this formalism the density operator ρ is replaced by where |β is a coherent state of the mechanical oscillator. In the Appendix we provide an extension and comparison of the present approach based on the Q function to a general (sparameterized) quasi-probability distribution including the Pdistribution and Wigner function as special cases. σ(β, β * ) is a density operator for the cavity field and a quasi-probability distribution for the oscillator over the complex phase space variables (β, β * ). The reduced density operator for the cavity is obtained by integrating over phase space, and the quasi-probability distribution (Q function) for the oscillator follows on taking the trace over the cavity, σ(β, β * ) itself still contains all information about the state of both systems, and is fully equivalent to the density operator ρ.
For the Q function the replacement rules [32] b and their adjoints can be applied to the master equation (6) in order to arrive at an equivalent description in phase space of the oscillator. We use the notation ∂ β to denote the partial derivative with respect to a variable β. The translated equation of motion is with The Liouvillian L m affects only the mechanical oscillator, and is just the Fokker-Planck version of Eq. (7). A crucial point in this formalism is that the nonlinear optomechanical interaction L int from Eq. (9) makes a contribution to both, the new Liouvillian for the cavity L c and the new interaction L int . Parts of the interaction can thus formally be treated as a shift of the detuning by 2g 0 Re(β) which depends on the phase space variables (β, β * ). Note that Eq. (12) is still exactly equivalent to (6).
A Semi-Polaron-Transformation -The parametric dependence of the cavity detuning on the phase space variables can be transformed into one of the driving field E by means of a transformatioñ When transforming the equation of motion (12) care has to be taken on commuting the unitary operators in (15) with derivatives with respect to (β, β * ) in L m and L int due to the β-dependence of θ. Details are given in App. C. The resulting equation of motion forσ(β, β * , t) can be written again in the form of Eq. (12), where L m and L int remain unchanged as in (13) and (14), and the Liouvillian operator for the cavity becomes In this picture the phase of the driving field is different for each point in phase space (via θ(β, β * )), and the cavity acquires an effective Kerr nonlinearity of strength We point out that the effective Kerr nonlinearity of the optomechanical interaction gives rise to ponderomotive squeezing of light, as was recently observed in [6,7]. The equation of motion forσ(β, β * , t), Eq. (16), is an approximation. In principle it contains further terms which are of order Q −1 m and whose explicit form is given in App. A 1. For high quality oscillators these terms provide only small corrections and, therefore, will be dropped in the following. Apart from this approximation Eq. (16) still contains the full nonlinear dynamics of the system, while the aspect of the optical Kerr nonlinearity is explicitly separated from the nonlinearity in the optomechanical interaction. It is also important to note that the quasiprobability distribution for the reduced state of the oscillator still follows from the transformed statẽ σ(β, β * , t) in Eq. (15) by taking the partial trace over the cavity Semi-Polaron-versus Polaron-Transformation -The transformation in Eq. (15) has many parallels with the polaron transformation [39] which has been applied fruitfully to optomechanical systems in order to describe single-photon strong coupling effects [34,35]. The polaron transformation is effected by a unitary transformation of the density operator which should be compared to the transformation in Eq. (15). Instead of (6) the transformed stateρ fulfills a transformed master equatioṅ where L m is given in Eq. (7). This equation is again correct up to terms of order Q −1 m . It is instructive to compare the master equation in the polaron picture (20) to the equation of motion (16) attained in our "semi-polaron transformation". In both equations of motion the Liouvillians for the cavity, Eqs. (17) and (21) respectively, exhibit a Kerr nonlinearity and contain a driving field whose phase depends on the momentum of the oscillator. Crucially, the polaron transformation changes the jump operator describing cavity decay from a to e iη(b−b † )/2 a, and entirely removes the interaction term (9). Moreover, since the polaron picture corresponds to a transformation into dressed states of the optomechanical system the partial trace ofρ pol over the (dressed) cavity mode does not give the reduced state of the mechanical oscillator, cf. Eq. (19). In contrast, the semi-polaron transformation introduced here retains a nonlinear interaction L int , Eq. (14), leaves the jump operator for cavity decay unchanged, and conserves the important relation (18). These properties are crucial in order to perform second order perturbation theory in L int , and to derive an effective equation of motion for the mechanical oscillator as a separate system. For further comments on the semi-polaron transformation we refer to Appendix C.

B. Fokker-Planck Equation for the Mechanical Oscillator
Interaction picture -Our goal is now to adiabatically eliminate the cavity field from the dynamics, similar to the analysis of sideband-cooling [36]. This requires that the cavity dynamics, governed byL c in (17) with dominant characteristic time scale κ, is fast as compared to all other time scales in L m and L int . Since we aim to cover in particular also the resolved sideband regime, ω m > κ, we move to an interaction picture with respect to the free harmonic motion of the mirror. The equation of motion is still given by Eq. (16) where L m describes thermal decay only, andL c and L int become explicitly time-dependent, The phase of the driving field is θ(β, β * , t) = η Im β e −iω m t .
In the adiabatic elimination it is assumed that the cavity essentially remains in the (quasi) stationary state of its undisturbed (by L int ) dynamics,ρ c =L c ρ c (23) withL c given by (22). This Liouvillian describes the dynamics of a Kerr nonlinear cavity driven by an amplitude-and phase-modulated field, where J n are Bessel functions. Note that the partial amplitudes depend on the mechanical phase space variable β = |β|e iφ . We do not attempt to solve Eq. (23) exactly. While in fact an exact solution for the stationary state of a Kerr nonlinear cavity exists [40] for the case of a constant driving field (i.e. θ ≡ const.), no such state can be expected for the present situation. Due to the periodic modulation of the driving field the cavity will not settle into a strictly stationary state, but rather to quasi-stationary state with a periodic time-dependence. If the Kerr nonlinearity is neglected an exact solution for this quasi-stationary state can be constructed by means of a Floquet series Ansatz [41]. However, in the present case both aspects, modulated drive and Kerr nonlinearity, are important and shall be taken into account. In order to arrive at an approximate solution of Eq. (23) which can serve as a (β-dependent) reference state for the adiabatic elimination of the cavity we will follow two complementary approaches in the paragraphs below: In the first one we assume the cavity is driven to a state of large mean amplitude, which we determine self-consistently from an essentially classical nonlinear dynamics. The fluctuations around this mean field will be treated in a linearized model as Gaussian noise. The second approach concerns the case of a weak driving fields for which the cavity essentially stays close to its ground (vacuum) state, which corresponds to the regime considered in [25,26]. In this case the master equation Eq. (23) can be expanded and directly solved on the low lying Fock states.
In both cases we aim to retain a nonlinear dynamics for the mean cavity amplitude, and use a linearized description for fluctuations. Formally this is done by switching to a displaced frame, definingσ = D(α(t))σD † (α(t)) where D(α) = exp(αa † − αa). We choose α(t) ∈ C such that the terms of dominant order in α are canceled from the transformed master equation forσ. In the case of |α(t)| 1 (|α(t)| 1) we cancel the terms of third (up to first) order in α and then neglect the terms up to first (of third) order in α. The full equation of this transformation may be found in Eq. (A13) of the Appendix, we proceed here with its most important features: Displaced frame for the limit |α(t)| 1 -In the limit |α(t)| 1 we identify α(t) with the long time solution oḟ which formally follows from the requirement that terms of third order in α in the resulting master equation for the displaced stateσ are canceled. Due to the Kerr nonlinearity the dynamics described by this equation of motion can be bistable. On assuming a single stable solution we preclude bistable regimes from our description. For a constant phase θ bistability occurs only for driving fields which are red detuned with respect to the cavity resonance for detunings ∆ < − √ 3κ, see [40]. For the present case of a modulated drive no such simple condition can be given. However, it is reasonable to expect that bistability will become an issue only when the driving field has sufficient spectral weight for frequencies with a detuning below − √ 3κ. In the following we are mainly concerned with the cases of resonant or blue detuned drive, for which it turns out that bistability is not an issue [42][43][44].
From Eq. (25) and (24) we can expect that in the long time limit the cavity amplitude will be of the form Inserting this expression into (25) one sees that the effective detuning experienced by the cavity will be dominantly given by the DC component of |α(t)| 2 , such that it is useful to define an effective detuning Eq. (27) has to be read as a non-linear equation for ∆ eff . In regimes where more than one solution exists, the system will be bi-or multistable, and we have to expect large photon number fluctuations. The validity of our approach will thus be limited to regions where only a single stable solution for ∆ eff exists, as discussed above. We seek an approximate solution to (25) by assuming a fixed effective detuning ∆ eff , such that where we follow the notation of [25,26]. In total α(t) in (26) depends on the mechanical phase space variable β through both ∆ eff (β, β * ) and the β-dependent driving field Ee iθ(β,β * ,t) .
We will see that the β-dependence in ∆ eff (β, β * ) is a crucial effect for the case of resonant cavity-drive (for which ∆ eff κ). The Liouvillians after the transformation with D(α(t)) are The Liouvillian for the mechanical oscillator, L m , acquires an additional drift term (second line in (30)) with a nonlinear nonlinear drift coefficient ∝ e iω m t |α(β, β * , t)| 2 which contains in particular the nonlinear DC force and dynamic back action effects (i.e. optical damping and frequency shifts), as will be discussed below. In the Liouvillian for the cavity,L c , terms of order α(t) and lower have been dropped. The leading terms of order α 2 describe squeezing dynamics and an effective detuning. Finally, in L int only the term of linear order in α has been kept. Note also that when moving to the displaced frame commutators of the (β-dependent) displacement operators and derivatives with respect to β have been neglected. They would add corrections to the Liouvillians of higher order in g 0 . We have now removed the driving field from the dynamics of the cavity. The remaining Liouvillian (31) describe the Gaussian evolution of fluctuations: The ponderomotive squeezing of the light field is naturally contained in the α 2 a † 2 -term and its hermitian conjugate. While in this article we will study parameters for which this squeezing is negligible, the effect of ponderomotive squeezing back on the mirror after the adiabatic elimination of the cavity is an interesting perspective for future applications of our new formalism: In the adiabatic elimination of the cavity one would have to use a squeezed reference state, which can introduce possibly additional diffusion terms in the motion of the mirror. Applied to the situation of limit cycles, this may cause the state of the oscillator to become more classical.
Curiously, the Kerr nonlinearity induces a different effective detuning for the mean field α than for the fluctuations (compare Eqs. (25) and (31)). This is consistent with results for a Kerr nonlinear cavity [40]. We therefore definẽ The fast decay rate to the vacuum is still given by κ from the original master Eq. (8).
This can be used in order to adiabatically eliminate the cavity taking into account second order effects in the optomechanical interaction ∝ g 0 , Eq. (32), very much in the spirit of laser cooling theory [36]. Details of the calculation can be found in Appendix B. The result is an effective equation of motion for the mechanical oscillator in the form of a Fokker-Planck equatioṅ for the Q-function of the mechanical oscillator. In analogy to h n we defineh n = κ + i(nω m −∆ eff ) with∆ eff given in (33). The drift and diffusion coefficients in the Fokker-Planck equation (34) do not depend on the phase of β as a consequence of the rotating wave approximation involved in its derivation. We therefore transform the Fokker-Planck equation to polar coordinates β = re iφ , and focus on the time evolution of the oscillator amplitude r by integrating out the phase variable φ. The time evolution for r is then a one dimensional Fokker-Planck equation (on a half space), with drift µ(r) and diffusion coefficient D(r) The details of the transformation may again be found in Appendix B. Equation (35) admits a potential solution in steady state which is given by (up to normalization) This solution is valid for any value of ∆, such that it covers both the regime of optomechanical cooling and the regime of self-induced oscillations. In [36] an effective equation of motion for the oscillator was derived in order to study the limits of sideband cooling under linearization of the dynamics and adiabatic elimination of the cavity using a coherent state as a reference state. The present approach generalizes this calculation to the nonlinear regime by using a different reference state for each phase space point of the oscillator. The nonlinear quantum dynamics has been described analytically using a method based on the classical theory for limit cycles [45] and by means of Langevin equations [25], and has been applied in great detail to limit cycles, but also to the cooling regime [26]. The results of our calculation reproduce these results in the regime of a negligible Kerr-Term and provide suitable extensions in those cases where the Kerr nonlinearity of the cavity becomes a dominant effect. In the next section we will compare the analytical expression for the steady state of the mechanical to numerical solutions of the exact master equation (6) to study the limit cycle regime. We conclude this section by briefly stating the corresponding results for the limit of small intracavity field amplitude, followed by a comparison of limit cycles studied in different laser setups. Displaced frame and adiabatic elimination for the limit |α(t)| 1 -In the case of |α(t)| 1 all steps can be performed in analogy. The difference is that we need to cancel the terms up to first order in α and then neglect the terms of third order in α. The effective detuning now is given by i.e. the bare detuning is just shifted by the constant Kerr term in this extreme regime. No distinction between ∆ eff and∆ eff needs to be made. The adequate choice for the displacement amplitude is the long term solution oḟ The result of the adiabatic elimination is structurally the same, α n and h n are given as in equations (28) and (29), but with the effective detuning now given as in (39). In Eq. (37) theh n are simply replaced by h n .
Quantum limit cycles in lasers -It seems natural to base a model of optomechanical limit cycles on theory used in the context of lasers [30][31][32], where quantum limit cycles have been extensively studied most prominently. The standard laser system consists of a reservoir of many atoms which forms a bath for the cavity mode. The pumped atoms will drive the laser mode to a high amplitude limit cycle, where it settles into a coherent state with random phase.
A setup that can be driven to highly sub-Poissonian states is the regularly pumped laser [46,47], where excited atoms fly through a cavity. The mechanism works in the situation, where at each time approximately only one atom interacts with the light mode and the interaction is a swapping of excitations. In the case of more regular than Poissonian statistics of the pump, the fluctuation of the transmitted energy decreases and the light mode will have sub-Poissonian phonon statistics. This setup is sometimes also referred to as one atom laser/maser or micro maser, because the events when more than one atom interacts with the field can be neglected.
The one atom laser is different from the 'one-and-the-same' atom laser [48], where a single atom is trapped inside a cavity and drives the laser mode. Also in this setup a sub-Poissonian steady state can be reached and the explanation again relies on counting the number of interactions exchanged between the atom and the cavity [49]. In our optomechanical system a single laser mode is the bath driving the mechanical oscillator. The bath consisting of only a single mode is in analogy to some extend to the micro maser, as stressed in [28], and even more similar to the 'one-and-the-same' atom laser. Even though we also describe sub-Poissonian boson statistics, the analytical techniques developed e.g. in [49] cannot be readily applied to our situation, because they crucially rely one the preservation of total excitations by the interaction, which is not given in the optomechanical setup. Our analytical model III is based on [31], which was first developed for the standard setup without sub-Poissonian statistics.
For the creation of non-Gaussian states a nonlinearity is required. In the optomechanical setup the nonlinearity stems from the interaction, while in the 'one-and-the-same' atom laser it stems from the two level nature of the bath, which is equivalent to a highly nonlinear cavity.

A. Introduction
As an introduction to our study, we sum up some known results on limit cycles that the rest of the article refers to. First, we introduce the theory for the amplitude of classical limit cycles as developed in [19,23], and then we recapitulate the numerical results on nonclassical states of quantum limit cycles as reported in [27]. When comparing these findings to our analytical treatment we will be mainly concerned with the special case of close to resonant driving field, ∆ eff κ ω m . Therefore, we start out by stating some approximate expres- sions for this case.
In the sideband-resolved regime, the equation for the effective detuning, Eq. (27), becomes a third order poynomial in ∆ eff and in the limit ∆ eff κ it even simplifies to a simple and explicit expression ∆ eff (r) ∆ + 2 KE 2 κ 2 J 2 0 (ηr).
Following Eq. (41) the combined intrinsic and optically induced damping of the oscillator γ eff (r) close to resonance is then given as the sum of the intrinsic mechanical damping γ and the amplitude-dependent optical damping Note that the sign of the optically induced damping at r = 0 coincides with the sign of ∆ eff . For negligible intrinsic damping, γ γ opt , one can then expect limit cycles to always start for ∆ eff > 0 (whereas the dynamics will be stable for ∆ eff < 0). The possible amplitudes r 0 for limit cycles are given by the conditions γ eff (r 0 ) = 0 and γ eff (r 0 ) > 0. The first condition is equivalent to The left hand side of this equation has infinitely many roots as the Bessel functions oscillate at a constant amplitude, cf. Fig. 1 a). The envelope is given by the r −1 decay. As illustrated in Fig. 1 a) the exact position of the limit cycle and the number of possible amplitudes is then determined by the right hand side of Eq. (44).

B. Outline
In the following subsections we will explain two features of limit cycles on resonance that can be heavily influenced by the Kerr term: First, in Sec. IV C we show that in the strong driving limit |α| 2 1 the phase transition between optomechanical cooling and self-induced oscillations is crucially determined by the dynamical dependence of the effective detuning on the intracavity amplitude and its corresponding nonlinear dependence on the cycle amplitude, cf. Eq. (42). This behavior can also be explained in a classical picture.
We then develop an explanation of the interesting numerical result for limit cycles in the quantum regime reported in [27,28]: For approximately resonant driving fields, ∆ 0, and at the blue detuned sideband resonance, ∆ ω m , the steady state of the mechanical oscillator can have a Wigner function with a negative area. The requirement on the strength of the optomechanical coupling g 0 is more stringent at the sideband than on resonance where non-classical limit cycles appear already for weaker coupling.Curiously, on resonance the numerical solution to the master equation predicts (nonclassical) limit cycles also for parameters where classically the effective detuning ∆ eff < 0, and one would expect a stable cooling dynamics. Fig. 2 shows the steady-state Wigner function of the mechanical oscillator for such parameters.
We will use the analytical description of limit cycles with the Fokker-Planck equation to explain the features displayed in Fig. 2, and to predict general requirements on system parameters to achieve a non-positive Wigner function. In section IV E we show that the occurrence of negative Wigner functions in turn is intimately linked to achieving a small variance of the phonon statistics, as characterized by a small Fano factor F = ∆n 2 / n , along with a small cycle amplitude r 0 . We analyze the variance of the phonon number in section IV D and find that the conditions for small Fano factor are favorable at the ∆ = 0-resonance.
In section IV F we describe the numerical method used to check the analytical predictions. It allows for the first time to numerically study quantum features of optomechanical limit cycles in the regime of large mechanical amplitudes and strong laser drive, populating many states of the cavity. We find that the analytical model can still be applied and even for g 0 < κ negativity of the Wigner function can be observed.

C. Drift and dynamical detuning
In this section we study in detail the time evolution of the mean amplituder, which is determined by the drift µ(r) in (41). In particular we show how the dynamical dependence of ∆ eff (r) on r gives new results which are not observed in any model based on a static detuning (like the one we used above). We focus on the regime wherer is larger than its standard deviation ∆r, such that we can derive the time evolution ofr viaṙ = µ(r) directly from (41) aṡ These assumptions are fulfilled for small η = 2g 0 ω m , because ηr is the argument of the Bessel functions and hencer ∝ 1 η . With the oscillator initially in the ground state, it is the sign of ∆ eff (0) that determines if the limit cycle starts at all: For ∆ eff (0) < 0 the optical damping is initially positive and no oscillation starts, but for ∆ eff (0) > 0 it is negative and may be larger than the intrinsic damping γ, so that a self-induced oscillations can start. The oscillator arrives at its steady state, whenṙ = 0. Neglecting the small corrections due to γ, this is equivalent to the condition ∆ eff (r)J 0 (ηr)J 1 (ηr) = 0. If the effective detuning ∆ eff is independent of r, the smallest root of this product is always the first root of J 0 . This corresponds to the standard situation (as discussed above) valid for a negligible Kerr parameter or in the weak driving limit, cf. Eq. (39).
In the converse case, for large amplitudes |α| 2 1 and non-negligible Kerr parameter, the dynamic nature of the effective detuning can become important: The smallest root of the product ∆ eff (r)J 0 (ηr)J 1 (ηr) is then determined either by J 0 or ∆ eff , depending on the sign of ∆. If the bare detuning is on the blue (heating) side, ∆ 0, the condition for the limit cycle is still J 0 (ηr 0 ) = 0 as in the case of a static detuning. However, if the bare detuning is on the red (cooling) side ∆ < 0 the effective detuning for a small cycle amplitude can still be positive as ∆ eff (0) = ∆ + 2KE 2 /κ 2 , cf. Eq. (42). This is the case in particular for a driving field E larger than a critical value of E crit = κ √ 2g 0 √ |∆|ω M . The sign of ∆ eff (r) will then depend on, and ultimately change with, the increasing amplitude r of the oscillation since ∆ eff = ∆ < 0 at the roots of J 2 0 (ηr). With increasing oscillator amplitude r the DC-component of the cavity occupation and hence (via the Kerr nonlinearity) also the shift of the detuning drops. The steady state amplitude r 0 of the limit cycle is reached when ∆ eff (r 0 ) = 0. Using again approximation (42) the condition ∆ eff (r 0 ) = 0 is equivalent to J 0 (ηr 0 ) = κ √ 2g 0 E √ |∆|ω M . Thus, the Kerr nonlinearity smoothens the transition from cooling to amplification. This is in contrast to models with a static detuning where a sharp transitions occurs at ∆ eff = 0.
We numerically check the dynamical nature of the detuning by integrating the equations of motioṅ which are the classical analogue to the master equation (6). Fig. 3 illustrates the time dependence of the detuning with an example of a time evolution where the bare detuning ∆ < 0, so that the limit cycle amplitude r 0 in steady state is determined by the condition ∆ eff (r) = 0. Fig. 4 shows that this condition gives a good prediction for r 0 as a function of ∆.
An approximation similar to equation (45) for the case of a laser drive close to the first blue sideband, ∆ ≈ ω m , shows that there the position of the limit cycle does not depend on the exact value of ∆. It is approximately given by the first root of J 1 (ηr). Thus the limit cycle amplitude is generally smaller on resonance than on the sideband. We will use this observation in section IV E, where we will see that a small limit cycle amplitude is favorable for the occurrence of a negative area in the Wigner function.   (46), with initial condition r = 0 for ∆ 0 but ∆ eff (r = 0) > 0. Effective detuning ∆ eff (t) (a) with scale on the left (blue) axis, oscillator amplitude r(t) (b) and DC-shift in position (c) with scale on the right (black) axis. A positive effective detuning at r ≈ 0 ensures that the limit cycle starts. With increasing oscillator amplitude the intra-cavity photon number n |α n | 2 from Eq. (27) drops and hence also ∆ eff . As µ ∝ ∆ eff , see Eq. (41), the oscillator settles in steady state as soon as this drop reaches ∆ eff = 0. The parameters in this plot are (E, g 0 , κ E = 2κ, γ E = 2γ) = (4.0, 0.05, 0.3, 2 · 10 −5 ) × ω m .

D. Diffusion and Fano factor
Having discussed the conditions for a limit cycle to start and having derived the mean amplitude in steady state for the ∆ 0 resonance, we now consider the fluctuations caused by the diffusion D around this mean value to derive a prediction for the Fano factor F = ( n 2 − n 2 )/ n = ∆n 2 / n , which is a measure for number squeezing: For a coherent state the phonon distribution is Poissonian so that ∆n 2 = n and F = 1. A state with sub-Poissonian phonon variance can hence be characterized by F < 1.
We will use the term Fano factor in the context of limit cycles as follows: For generic parameters an optomechanical system can exhibit several limit cycles, such that the Fano factor of the full density matrix typically is larger than one. The oscillations at each of these attractors are metastable, such that it is possible to consider the phonon statistics at a particular limit cycle. Especially in the relatively classical regime where g 0 /ω m is not too large the cycles will be well separated. When we refer to Fano factor, we will implicitly always mean the Fano factor of one particular attractor.
We obtain the mean and variance of the phonon number n via [32] a r (a † ) s where W(α, α * ) is the Wigner function. We use here the Wigner function because it gives better agreement with the numerical analysis for the the statistics of the phonon number  Figure 4. Amplitude r 0 for the first stable limit cycle versus bare detuning ∆. In the limit of an amplitude-independent effective detuning (red) the values for large amplitudes are predicted correctly. It is known from [19] that for small amplitudes at the onset of limit cycles, the amplitude follows a square root (red). With inclusion of the dynamical effective detuning ∆ eff (r) (blue) the limit cycle amplitude r 0 follows J 0 (ηr 0 ) = κ √ 2g 0 E √ |∆|ω M , both limit cases are reproduced, and the whole transition between the regimes of damping and antidamping can be described. In this figure  than other quasi-probability distributions. Drift and diffusion coefficients for the Wigner function are calculated in App. A along the same lines as shown above for the Q-function. In particular, close to resonance the radial diffusion coefficient as relevant to the Wigner function is where we applied to equation (B24) the same approximations as in Sec. IV A for the drift coefficient. For most amplitudes the J 2 0 -term is dominant, as it is enhanced by at least (ω m /κ) 2 over the J 2 1 -term. For parameters where the optical anti-damping is much stronger than the intrinsic mechanical decay, a curious cancellation of the diffusion occurs in steady state: The limit cycle will then settle exactly at the first root of J 0 as discussed in section IV C. There the term proportional to J 2 1 , which is suppresed by (κ/ω m ) 2 , becomes the only relevant term in the diffusion. This suppression is illustrated in figure 1 b) and can be intuitively explained: The last two terms in equation (B11) (or equivalently (B24)) are the (coherent) squeezing terms. For n = 1 they exactly cancel the corresponding (incoherent) diffusion terms ∝ ∂ β * ∂ β * in leading order and only the higher order terms in κ 2 /ω 2 m remain. Because of this suppression of diffusion in the sideband-resolved regime one can obtain a very small Fano factor of the mechanical oscillator, as we show below.
The phase space distribution in steady state is given by Eq. (38). In the limit of small g 0 /ω m , where ∆n n , and for the case of only a single stable limit cycle centered around a position r 0 with µ(r 0 ) = 0, we linearize µ(r) µ(r 0 ) + µ (r 0 )(r − r 0 ) around this r 0 and set D(r) D(r 0 ) so that the corresponding solution for W is approximately with σ 2 = −D(r 0 )/µ (r 0 ). One can then derive the approximate expression F 4σ 2 for the limit ω m /g 0 > σ. In the sideband-resolved regime and with the limit cycle position at the first root of J 0 this gives where ζ 0.27 is the numerical value of J 2 1 at the position of the limit cycle. The Fano factor is minimal at an effective detuning ∆ eff (r 0 ) = κ where it takes on the value Note first that Eq. (51) implies that the Fano factor is lower bounded by the sideband resolution and that this bound is achieved for sufficiently large driving field E = √ 2κP L / ω L (laser power P L ). Furthermore Eq. (51) implies that the condition for sub-Poissonian statistics 1 > F is exactly equivalent to . This can be interpreted as a condition for the driving power which for small κ/ω m becomes It is instructive to express this also in terms of the (thermal, linearized) cooperativity parameter where we used that the relevant average intracavity amplitude at the optomechanical limit cycles is α = α 1 E/ω m , cf. Eq. (28). Condition (53) then takes the form (in the limitn 1) Note that this is essentially a requirement on the linearized optomechanical coupling (g ∝ g 0 E), and not on the coupling per single photon g 0 . The condition in Eqs. (53) and (55), and the lower bound in Eq. (52) are the main result regarding sub-Poissonian phonon statistics. The possibility of a sub-Poissonian number distribution was discussed in [25,26] for the resonance at the first (and higher) blue sidebands. The prediction of the analytical model is especially good for the regime with small g 0 that results in larger limit cycle amplitudes. In Figure 5, which compares the Fano factors as derived from our analytical model and from solving the master equation, the good agreement can be seen. For larger g 0 (not depicted in Figure 5) the condition neccesary for adiabatic elimination is less satisfied and also the linear approximation (49) gets worse, because ∆n ≈ n . Thus the quantitative agreement gets worse. Still the resonances for F at ∆ ≈ 0, ω m are qualitatively reproduced.
In [25,26] the Fano Factor has been calculated with a derivation using the truncated Wigner function approximation and solving the resulting Langevin equation. If we use the Wigner function as the phase space distribution, our calculation, which does not rely on this truncation, gives the same result in the regime where the Kerr parameter K is negligible.
For limit cycles with the cavity close to its ground state, different approaches to treat the Kerr effect have been taken in the literature: [28] uses the classical part of the Kerr effect, as derived with the standard master equation approach, to introduce a renormalized detuning with a shift proportional to the cavity occupation. An additional constant (independent of the cavity occupation) shift of the detuning by K = g 2 0 /ω m , was numerically observed in [26] and then introduced by hand, to match the numerical data. It is one of the main results of the the semi-polaron approach, that the separate Kerr term for the cavity is naturally derived for limit cycles. It causes exactly the additional quantum shift of ∆ observed in [26], which is most striking in the |α| 1 limit, cf. Eq. (39).

E. Nonpositive Wigner Function
Finally we use the Fano factor to predict the occurence of a negative area in the Wigner function. For a Fock state the Fano factor F is of course zero and, except for the vacuum, all Fock states have a pronounced negativity of the Wigner function. Both F and the Wigner function are continuous functions of the state ρ. Hence, for a given mean phonon number n 0 there is a critical value F c , such that for a state with F < F c the Wigner function has a negative area. For simple set of Ansatz states given by a density matrix diagonal in Fock basis with Gaussian probability distribution we numerically determined the corresponding critical Fano factor F c . The result is illustrated in Fig. 6. We use this particular Ansatz, because the typical steady state density matrix of our problem is approximately of this form when g 0 /ω m is not too large. In [26] the steady state as a Gaussian distribution in Fock states is derived in more detail. Figure 6 shows that this threshold F c is smaller for larger amplitude r 0 . We infer that in order to see negativity of the Wigner function in steady state, small limit cycle amplitudes with small Fano factors are favorable. Applied to the results of [27,28] this explains the more favorable condition for negativity at the ∆ 0-resonance as compared to the ∆ ω mresonance, because the limit cycle there has a smaller amplitude (given by the first root of J 0 (ηr) as compared to J 1 (ηr), as discussed in section IV C). Independent of ∆, the amplitude scales with the inverse of g 0 /ω m , such that for a large ratio  (56). From this plot one can read of, how small the Fano factor needs to be for a given r 0 , to see a negative value in the Wigner function. Implicitely this is also a requirement on g 0 because r 0 ∝ ωm g 0 , see Section IV C.
g 0 /ω m a non-positive Wigner function is achieved already for larger Fano factors. More precisely, we can conclude from Fig. 6 that where the constant ξ depends on how negative the Wigner function should be. In order to achieve a ratio of minimal to maximal value of the Wigner function of e.g. −0.1 this constant is found to be ξ 0.6. As a comparison, this negativity ratio can reach (approximately) -2.5 for odd Fock states and -0.4 for even Fock states.
Since the amplitude of the first limit cycle is r 0 ω m /g 0 the condition F < F c , together with Eqs. (51) and (57), is equivalent to (ζ 0.27) Thus, one necessary condition for negative Wigner function is that the square bracket on the left side is positive. This is a condition on the single photon optomechanical coupling g 0 , that can be written equivalently as both Note that this condition for the occurrence of a quantum state is weaker than the condition g 0 /κ > 1 which one would have expected naively. Assuming this condition to be well fulfilled we can drop the second terms on both left and right hand side of (58) and get the power requirement Note that this is stronger than the requirement (53) for sub-Poissonian statistics, as one would expect. In terms of the cooperativity (for anyn) this becomes Note also that even for zero temperature,n → 0, there is now a threshold for the power (cooperativity) in contrast to the condition for sub-Poissonian statistics. Condition (59) on the strength of the optomechanical coupling per single photon, and condition (60) (or (61)), which reproduce the heuristically derived conditions (4) and (5) from Sec. II, are the main results regarding negative Wigner functions.

F. Numerical Analysis
In this section we compare the predictions from the sections above with the numerical result for the master equation of Eq. (6). To do the calculation for large Hilbert space dimension, we applied the Monte-Carlo wave function method from [50][51][52] as implemented in QuTiP [53,54], the quantum toolbox for python. The advantage is that one needs to simulate only wave functions and not density matrices, so that the Hilbert space dimension required for the simulation scales only with the number of possible pure states N instead of N 2 . In this method the individual trajectory of an initially pure state is calculated, conditioned on the history of fictive photon and phonon counters measuring the particles leaking out of the system. With this knowledge of the environment an initially pure state stays pure. The density matrix is then retrieved by averaging over a large ensemble of such conditional states. The ensemble average can be replaced by the time average for calculating a steady state density matrix.
Our implementation was done with an adaptive Hilbert space, where the Fock states are not only limited from above, but also from below and after each mechanical oscillation the Hilbert space is updated so that it is centered around the current state. To make sure that not too much of the Hilbert space is truncated, the number of states to be used is scaled with the standard deviation in energy of the state in the previous step. This flexibility of the Hilbert space during the calculation allows to run the simulation without much a priori knowledge of the steady state and even fewer basis states are required.
The solution is obtained in the following steps: For speed up of the calculation, the initial state is chosen to be a coherent state with an amplitude close to the expected steady state. It is then evolved for some period until at a time t 0 the conditional state's amplitude and Fano factor stop to drift and only fluctuate. We then make use of the fact, that in steady state the time average corresponds to the ensemble average, and calculate the steady state of the oscillator as where |ψ t is the conditional state at time t, and T spans many mechanical oscillations. This procedure is performed many times in parallel on a cluster and the resulting matrices ρ M are averaged. The deviation of the individual ρ M provides an error estimate for the method. As a further benchmark and control, we also calculated the steady state with the biconjugate gradient steadystate solver from scipy [55], which is however limited to a comparably small Hilbert space dimension.
The algorithm described above allows for the first time to numerically study optomechanical limit cycles in the experimentally relevant regime of large amplitudes of the mechanical oscillator (as caused by a relatively small g 0 /ω m ) and with more than only a few photons in the cavity. In previous studies the question was posed, whether the analytical theory can be applied to this regime [26] and if the nonclassical features survive [28] for more than one photon in the cavity. We answer this question affirmative: Fig. 7 shows an example of a Wigner function in this regime with small Fano factor and some negative density.
Strictly speaking the steady state calculated here is only metastable if γ is so small that there is more than one attractor for the limit cycle, cf. Fig. 1. The timescale for switching between different attractors is much longer than the time to relax in a given metastable steady state. Thus it is not considered in this article. In order to choose the metastable attractor for the numerical simulation, we choose an initial state in the vicinity of our preferred attractor, in this case the limit cycle with lowest possible amplitude. Also in the analytical expressions for the Fano factor, we always treat possibly metastable states as steady states. For very large g 0 /ω m different metastable attractors start to merge and the analysis becomes more involved. This merging of attractors and its effect on nonclassical features was studied in detail by [28].

V. CONCLUSIONS
We studied the quantum regime of optomechanical limit cycles. Based on the Laser theory of Haake and Lewenstein [31] we derived an effective Fokker-Planck equation for an optomechanical system. The analytical prediction for the oscillator's steady state is in agreement with the work of Rodrigues and Armour [25,26] for driving fields on the first blue sideband.
Our treatment naturally includes also the Kerr effect, which becomes important for large g 2 0 /ω m . One consequence important for the quantum theory of limit cycles is the shift of the detuning of equation (39), which occurs even without photons in the cavity, and had to be introduced phenomenologically in [26]. This shift explains the possibility of limit cycles on the blue sideband in [27] or for the parameters of Fig. 2.
The effective cavity detuning is usually approximated as a static variable. Within our framework one can describe its dynamical nature, which is a classical phenomenon scaling proportional to the Kerr parameter. Figures 4 and 3 show how this smoothens the phase transition between optomechanical cooling and self-induced oscillations.
We studied the quantum limit cycles on resonance and found the simple analytical expression (50), that predicts the possibility of very small values for the Fano factor F of the mechanical oscillator. We found that in the sideband resolved regime a large value of m γn , i.e. a large linearized optomechanical coupling, is required to minimize F.
We then established a relation between sub-Poissonian phonon statistics and negativity of the Wigner density for typical parameters of limit cycles: The oscillator's steady state has an approximately Gaussian number distribution at each metastable limit cycle. For these states the requirement on F to see negativity of the Wigner function is given by the function of figure 6.
Using a Monte-Carlo method with an adaptive Hilbert space, we numerically checked this scaling even for limit cycles with very large amplitude and many photons in the cavity, where an ordinary steady state solver cannot be applied. The numerical simulation depicted in figure 5 show that indeed the criterion of a small Fano factor can predict the negativity of the Wigner function. For currently more feasible experimental parameters with even smaller g 0 /ω m , the negativity disappears according to Fig. 6, but the very small Fano factors remain.
We believe that the present approach provides a suitable starting point for further studies of optomechanical system in the limit of strong couplings. We point out once more that in the "semi-polaron picture" introduced here the Kerr nonlinearity and the optomechanical interaction occur as independent terms. This enables in principle to take into account the squeezed noise of the cavity when deriving effective equations of motion of the mechanical oscillator. While for the parameters considered in this article we could neglect this effect, additional diffusion for the mechanical oscillator is to be expected for very strong laser drive. This would apply to the case of limit cycles, but could also become important in the cooling regime.
In the main text we introduced the semi-polaron tranformation only for the special case of the Q-function, to make the equations more readable. Here we drop this restriction and assume the more general case of an s-parameterized phase space distribution P s with s ∈ [−1, 1]. For the convenience of the calculation we define p = s+1 2 ∈ [0, 1] and q = 1 − p. Note that for q = 0 this corresponds to the Glauber-Sudarshan P-representation, for q = 1 2 to the Wigner-representation, and for q = 1 to the Husimi Q-representation.
Depending on wether one wants to study the regime |α| 1 or |α| 1, either the terms with low or high order in α can be neglected at this point and a different choice of α(t) is required to cancel all displacement-like terms.
3. Displaced frame for |α| 1 If we restrict the analysis to only the lowest two Fock states, the operators consisting of three creation/annihlationoperators resulting from transformation (A13) can be approximated with just one operator, e.g. aa † a ≈ a. This time we neglect the terms proportional to K of third order in α. By imposing that α(t) solves this timė we can cancel the remaining displacement-like terms. The Liouvillians are now L c σ = − i −∆a † a − K a † a 2 − K α 2 (a † ) 2 + h.c. , σ + L c σ (A19) L int σ = − ig 0 (q∂ β t − p∂ β * t )σ(α * a + αa † + a † a) + ig 0 (q∂ β * t − p∂ β t )(α * a + αa † + a † a)σ (A20) When written in the second form we can apply the replacement rules (11) to write the last equation in operator representation In the second line we expressed the generator for the semi polaron transformation in terms of a commutator with a Hamiltonian and three Lindblad terms. The semi polaron transformation in operator representation is thus ρ = exp ηL s−pol ρ.
It becomes equivalent to the polaron transformation if the Lindblad terms in the generator L s−pol are dropped. Thus, the semi polaron transformation is non-unitary. In the context of adiabatic elimination of a cavity mode in the bad cavity limit a similar transformation to a "dissipation picture" was employed in [56,57].