Quantum limit-cycles and the Rayleigh and van der Pol oscillators

Self-oscillating systems, described in classical dynamics as limit cycles, are emerging as canonical models for driven dissipative nonequilibrium open quantum systems, and as key elements in quantum technology. We consider a family of models that interpolates between the classical textbook examples of the Rayleigh and the van der Pol oscillators, and follow their transition from the classical to the quantum domain, while properly formulating their corresponding quantum descriptions. We derive an exact analytical solution for the steady-state quantum dynamics of the simplest of these models, applicable to any bosonic system---whether mechanical, optical, or otherwise---that is coupled to its environment via single-boson and double-boson emission and absorption. Our solution is a generalization to arbitrary temperature of existing solutions for very-low, or zero, temperature, often misattributed to the quantum van der Pol oscillator. We closely explore the classical to quantum transition of the bifurcation to self-oscillations of this oscillator, while noting changes in the dynamics and identifying features that are uniquely quantum.


I. INTRODUCTION
Self-oscillating systems are ubiquitous-from humanmade clocks and transistors, through heart cells and neurons in the living body, to flashing fireflies and circadian rhythms-and are now emerging as canonical models for driven dissipative nonequilibrium open quantum systems, and as key elements in quantum technology. The dynamics of self oscillation are captured mathematically by the notion of a limit-cycle. Here we consider a family of models that interpolates between the Rayleigh [1] and the van der Pol (vdP) [2] oscillators, which are probably the most common textbook examples of limit-cycles in classical nonlinear dynamics. These models consist of a simple harmonic oscillator, driven by a time-independent energy pump in the form of "negative damping." When the pumping rate exceeds that of the normal damping rate, self-oscillations develop, which are then saturated by a nonlinear form of damping. The frequency of the oscillation is set by the physical parameters of the oscillator, while the magnitude of the oscillation is set by the ratio of the linear to the nonlinear damping rates. This provides a convenient knob with which to transition the oscillator from large-amplitude classical behavior to small-amplitude quantum behavior, which is our focus here.
Existing models for quantum limit cycles [3] consist of a harmonic, or possibly anharmonic, quantum oscillator, with linear as well as nonlinear coupling to the environment, which are expressed in terms of quantum Lindblad operators. These models are currently being used to study quantum entrainment [4], synchronization [5][6][7] and the phenomenon of "oscillation collapse" or "amplitude death" [8,9] in systems of coupled self-sustained oscillators, as well as the nonequilibrium spectral prop- * Corresponding author: ronlif@tau.ac.il erties [10], and the critical response to external drive [11], of single oscillators.
Our current focus is more basic. The classical Rayleigh and vdP oscillators are known for exhibiting a Hopf bifurcation, from a state of no motion at all to a state of self-oscillations at a fixed amplitude. We seek to characterize this bifurcation as the system transitions from the classical to the quantum domain. Our goal is to find answers to such questions as: How exactly should one model the Rayleigh and vdP oscillators in quantum mechanics? Can the quantum model analytically be solved, at least in its steady state? Is the quantum bifurcation different from the classical one? What experimentally observable indications are there to distinguish between quantum and classical behavior? What would be the first corrections to classical dynamics as one approaches the quantum domain?
Answers to these questions are relevant to a broad range of physical systems exhibiting quantum behavior, including lasers, or more generally photonic systems with nonlinear loss [12][13][14], as well as trapped ions [5,15] and electronic or superconducting circuits [16]. Particularly interesting is the attempt to observe such quantum behavior in nanotechnology-based human-made mechanical systems [17]. Indeed, modern nanomechanical resonators show exceptional behavior, as they routinely operate in the GHz range [18]. With nano-electromechanical systems (NEMS) [19] and nano-optomechanical systems (NOMS) [20] it is now possible to perform ultrasensitive measurements of physical quantities [21] such as single spins [22], minute charges [23], and tiny masses [24]. Relatively weak drive is needed in order for nonlinearity to be evident in the dynamics of nanomechanical systems [25,26], which is experimentally observed [27] and also exploited for applications [28]. Most importantly, at GHz frequencies, one need only cool to temperatures on the order of tens to even hundreds of mK for the thermal energy to become comparable to the quantum energylevel spacing of the mechanical resonator. This allows now to cool mechanical resonators down to their quantum ground state [29], and to start investigating fundamental physical questions on the borderline between the quantum and the classical worlds [30], as it applies to human-made macroscopic nonlinear mechanical objects. This, in turn, requires a well-based quantum theoretical framework.
We employ a phase-space approach to study the correspondence between classical and quantum limit-cycles. Since classical notions like a particle trajectory do not have a straightforward quantum analog, it is reasonable to compare quantum expectation values with classical statistical ensemble averages. We do so by solving the classical equations of motion for many different initial conditions (typically N = 10 4 ) taken from a Gaussian distribution, and keeping track of the different trajectories, thus representing a statistical distribution over phase space. The width of the initial distribution in phase space is taken to be the same as the quantum uncertainties ∆x and ∆p of an initial coherent-state wave function. In addition to expectation values, we also compare the full classical distribution with the quantum Wigner function W (x, p). The quantum dynamics are those of an open quantum system, and therefore described by a density matrix and its master equation, which dictates the steady state, and more generally, the dynamics of the quantum system.
We begin in section II with theoretical background for the classical dynamics of a family of models described by a generalized Rayleigh-van der Pol equation of motion (5), which interpolates continuously between the pure Rayleigh oscillator and the pure vdP oscillator. We provide a perturbative steady-state solution for limit cycles that are nearly-circular in phase space, obtained for weak driving just above the Hopf bifurcation to the oscillatory state. Moreover, we note that this solution is exact, and the limit cycles are always circular, for the model that lies exactly halfway between the pure Rayleigh and pure vdP oscillators, which we call the Rayleigh-van der Pol (RvdP) oscillator. In section III we introduce three quantum models, differing in the form of the nonlinear coupling of the oscillator to the environment. We discuss the basic features of these quantum models, and show that, for weak driving, their classical limits correspond to the RvdP oscillator (sec. III A), and to the pure vdP (sec. III B), and pure Rayleigh oscillators (sec. III C). In sec. III D we employ time correlation functions to elucidate some of the differences between these models. In section IV we derive an exact analytical solution for the steady-state dynamics of the quantum RvdP oscillator, which is a generalization to arbitrary temperature of existing solutions for very-low, or zero, temperature, often misattributed to the quantum vdP oscillator. In section V we consider in some detail the transition from classical to quantum dynamics of the RvdP oscillator, identifying dynamical behavior that is unique to the quantum domain. We conclude with a few summarizing remarks in section VI.

II. THE CLASSICAL RAYLEIGH AND VAN DER POL OSCILLATORS
Consider the following classical equation of motion, describing a harmonic oscillator with effective mass m and natural frequency ω, where tildes denote physical parameters that are soon to be rescaled. The oscillator is driven by a velocitydependent force or "negative damping", with coefficient κ 1 ≥ 0, as described earlier. It also experiences normal linear damping, with coefficientγ 1 ≥ 0, which is unavoidable in most physical systems, as well as two types of nonlinear damping mechanisms: vdP damping with coefficientη ≥ 0, which is proportional to the velocity and the squared displacement of the oscillator, and Rayleigh damping with coefficientζ ≥ 0, which is proportional to the cubed velocity of the oscillator.
To obtain a dimensionless equation of motion we (a) measure mass in units of m, effectively setting m in Eq. (1) to unity; (b) measure inverse time in units of the oscillator frequency ω by defining which effectively sets ω to unity; (c) measure length in units of x 0 = /mω by setting in anticipation of the quantum treatment below; and consequently, (d) measure the pumping and damping rates with respect to the chosen units of mass and time, by defining (4) where γ 2 > 0 is an overall dimensionless nonlinear damping rate, and η and ζ are numerical factors, indicating the relative contributions of the two nonlinear damping mechanisms. Without loss of generality, one can set the larger of the two to unity, and the smaller to a number between 0 and 1.
Finally, we divide the original equation of motion (1) by the characteristic unit of force, mω 2 x 0 , yielding a scaled dimensionless equation of the form where = κ 1 − γ 1 , and dots denote derivatives with respect to the dimensionless time t. This generalized Rayleigh-van der Pol equation is usually studied in one of its following limiting cases: (1) the Rayleigh oscillator [1] with η = 0, ζ = 1; (2) the van der Pol (vdP) oscillator [2] with η = 1, ζ = 0; and (3) the Rayleigh-van der Pol (RvdP) oscillator with η = ζ = 1, which is sometimes refered to as the harmonic RvdP oscillator [31]. All these variants are known to generate steady-state limit cycles for positive , as shown in Fig. 1. In the weak-drive limit of small , with nearly circular orbits, one can use secular perturbation theory [25,26] to obtain an approximate solution for the generalized RvdP equation of motion (5), and determine the amplitude of limit-cycle oscillations. The solution is written as a slow modulation of harmonic oscillations at unit frequency, with -dependent corrections where T = t is a slow time scale, characteristic of the rate of relaxation toward the limit cycle, as opposed to the fast time scale t of the oscillations themselves. As usual, c.c. stands for the complex conjugate. The slow time variation of the complex amplitude A(T ) also provides the freedom to eliminate secular terms, and to ensure that the perturbative correction x 1 (t), as well as all higher-order corrections do not diverge. Substituting the solution (6) into the equation of motion (5) indeed generates such a secular term [26,Section 11.4], which when required to vanish leads to a first-order differential equation for determining the slowly varying amplitude, The parameter A c = /γ 2 sets the overall scale of the oscillations, but each variant has its own particular sat-urated oscillation-amplitude, depending on the relative contributions of the Rayleigh and van der Pol damping mechanisms. Steady-state oscillations are obtained when Eq. (7) is set to zero, and the amplitude satisfies RvdP. (8) Note that in the small-amplitude slow limit, without a particular model at hand, it is difficult to discern the nonlinear terms from one another, as they merely combine into a single effective coefficient η eff = η + 3ζ. However, in the large-amplitude strong-drive limit, with 1, as can be seen in Fig. 1, the limit cycles look qualitatively very different. In particular, the RvdP oscillator, with η = ζ, is unique in that it is invariant under phase-space rotations, producing circularly-symmetric limit cycles, or harmonic oscillations [31], for arbitrary drive strength . In fact, one can easily verify that the zeroth-order term of the expanded solution (6), gives the exact steady-state solution, x(t) = A c cos t, for the RvdP oscillator, with all higher-order corrections cancelling out. As we shall see below, the RvdP oscillator is also the simplest to treat quantum mechanically. Finally, as expected for an autonomous or timeindependent equation of motion (5), the complex amplitude equation (7) is independent of phase, which drops out of both sides. This implies that with purely deterministic dynamics the oscillator will maintain any initial arbitrary phase, but in the presence of thermal, or any other source of noise, the phase of the oscillator will diffuse over time. This is demonstrated numerically in Fig.  2 for the vdP oscillator with weak thermal noise, where an initial Gaussian-distributed ensemble of independent oscillators quickly relaxes to the expected amplitude 2A c , and eventually spreads over the whole limit cycle. The simplest quantum model of a limit cycle-which is often mistaken for "the quantum vdP oscillator"employs standard Lindblad formalism to describe the interaction of the oscillator with its environment, whereby the energy pump, or negative damping, is implemented in terms of single-phonon absorption, and the nonlinear damping is described as two-phonon emission ("phonon" should be replaced with "photon", "polaron", "magnon", or any other bosonic excitation, depending on the particular physical realization of the oscillator). The physical realization we have in mind follows the framework that was introduced by Dykman and Krivoglaz in the 1970's, whereby the nonlinear damping [32] appears as a result of nonlinear interaction of the oscillator with a continuum of bath oscillators, while energy injection [33] is introduced in the form of an off-resonance pump, detuned a frequency ∆ 1 away from the oscillator frequency ω. Within this realization, and as expected in most other alternative realizations, the coupling of the oscillator to the bath inevitably will induce normal linear damping with single-phonon emission, in addition to the two-phonon processes above.
Consequently, the master equation for the density matrix ρ of the oscillator-considered at T = 0, for the time being-contains three Lindblad operators of the form and is given bẏ where H 0 = ω(a † a + 1/2) is the Hamiltonian of the harmonic oscillator, and a is its annihilation operator. This master equation (10) differs conceptually from those that are commonly used in the literature [4][5][6][7][8][9][10]. Common models assume that as in the classical regime the effect of the pump, or negative linear damping, combines with the normal linear damping to give one physical process, with coefficient (κ 1 −γ 1 ) = mω 0 . Thus they either omit the first Lindblad operator below the threshold of self oscillations, whenκ 1 <γ 1 , or omit the second Lindblad operator above threshold, forκ 1 >γ 1 . Consequently, as will become evident below, even though they obtain limit cycles in the steady state, they miss important physical effects in the quantum regime, related to the fact that at zero temperature there are three rather than only two sources of spontaneous quantum processes that affect the quantum oscillator and its phase stability.
In order to facilitate the direct comparison between classical and quantum dynamics of limit cycles, we use the same scaling here for the quantum master equation (10) as we did earlier for the classical equation of motion (1). This, again, amounts to using the effective mass m of the oscillator as the unit of mass, and its inverse frequency 1/ω as the unit of time, thereby effectively setting both m and ω to unity. The choice of x 0 = /mω as the unit of length, and correspondingly p 0 = √ m ω as the unit of momentum, amounts to using as the unit of action with which phase-space area is measured, thereby effectively setting to unity [34]. With this choice of scaling, energy is measured in units of ω, the Hamiltonian becomes H = p 2 + x 2 /2 = a † a + 1/2, where the creation and annihilation operators are defined as and the commutator [x, p] = i. The resulting dimensionless zero-temperature master equatioṅ (12) can be used to study the dynamics of the density matrix itself, or any dynamical quantity that can be derived from it. For example, Fig. 3 shows the characteristic behavior of the time evolution of the Wigner function calculated numerically [35], for an oscillator initiated as a coherent state with α = 0.25(1 + i)A c . As in the classical case, shown in Fig. 2, one can see how the quantum oscillator first approaches the fixed-amplitude orbit of the limit cycle and only later loses its phase. Note that the amplitude of the quantum limit cycle is A c rather than 2A c , which according to Eq. (8) seems to indicate that this limit cycle may in fact be the quantum version of the RvdP oscillator, and not that of the vdP oscillator. One may use the master equation (12) to obtain the equation of motion for any expectation value, O = Tr{ρO}. In the Schrödinger picture, where operators are time-independent, one has FIG. 3: Time evolution the Wigner function of a quantum limit cycle, calculated numerically using the RvdP master equation (12), at T = 0 with κ 1 = 0.1, γ 1 = 0, γ 2 = 1/640, and therefore A c = 8, starting at t = 0 with an initial coherent state with α = (1 + i)A c /4. The Wigner function approaches the limit cycle, and then loses its initial phase over time. A square of area is shown in panel (a). The dashed red circles have radius A c .
Thus, for the annihilation operator a-using the fact that the trace of a product of operators is invariant under their cyclic permutations-the scaled master equation (12) gives the zero-temperature equation of motion We see that the nonlinear term is proportional to (a † a)a, or (x 2 + p 2 )a, again, as one would expect for the RvdP oscillator rather than the vdP oscillator. To see this more clearly, we take the semiclassical limit where a † a 1, and therefore a † a aa † . The semiclassical amplitude equation for α = a is then readily derived from Eq. (15) by replacing a † aa with |α| 2 α to give where A 2 c = /γ 2 , as defined earlier. In order to use an equivalent ansatz to the classical one in Eq. (6) we note that, according to the definition of the creation and annihilation operators in Eq. (11), α is a factor of √ 2 smaller than the complex amplitude of the oscillator. We therefore take with T = t as before, and find that the slow amplitude equation is given by which corresponds to the classical amplitude equation (8) as long as one takes η eff = 4, or η = ζ = 1, as expected for the RvdP oscillator, and in agreement with the amplitude of the limit cycle observed numerically in Fig. 3. Finally, using the definition of a in Eq. (11), we can take the real and imaginary parts of Eq. (16) to obtain the equations of motion for the expectation values of the position and momentum operators [36, Section 7.4], Differentiating Eq. (19a) with respect to time, and substituting Eq. (19b) for ṗ , yields a second-order equation of motion for x of the form Neglecting corrections of order 2 , this explicitly agrees with the classical equation of motion (5) for the Rayleighvan der Pol oscillator, with η = ζ = 1. We wish to emphasize the circular symmetry of the steady-state Wigner function in Fig. 3(d). In order for the steady-state Wigner function to lack such symmetry, the steady-state density matrix must contain off-diagonal elements that do not decay to zero. This can be seen by noting that the Wigner function (13) is a linear function of the density matrix, which can be expressed as a sum of its elements, and that the Wigner function of a diagonal element-given by [37,Section 4.4 where L n (x) is the n th Laguerre polynomial-is rotationally invariant.
As previous authors [10,12] have noted, the master equation (12) for the RvdP oscillator does not couple density-matrix elements that are not on the same diagonal. To see this, it is helpful to relabel the matrix elements ρ n,n+m = n|ρ|n + m according to their degree m of off-diagonality using a transformation, similar to the one used by Simaan and Loudon [12] (12), with γ 1 = 0, γ 2 = /16, and therefore A c = 4 for all different values of = κ 1 . All off-diagonal matrix elements decay to zero in the steady state, yielding the same circular limit cycle, independent of for constant A c . Compare with Figs. 6 and 7 below for the quantum vdP and Rayleigh oscillators.
The rate equations for the transformed matrix elements are theṅ where evidently, matrix elements are coupled only if they have the same degree m of off-diagonally. Thus, each diagonal can be considered as a separate "block" of the density matrix, evolving independently of all the others, allowing the off-diagonal elements to decay to zero, as one expects, independent of the principal diagonal elements, which are the only ones to survive in the steady state. This is confirmed numerically in Fig. 4.

B. The Quantum van der Pol Oscillator
One can obtain a master equation whose classical limit gives the vdP oscillator, at least to first order in , and is capable of producing quantum limit cycles that are non-circular in phase space. This is done by changing the Lindblad operator for the nonlinear damping term in Eq. (12) from γ 2 D[a 2 ] to γ 2 D[xa/ √ 2], breaking the rotational symmetry in phase space. The zero-temperature master equation then becomeṡ where we recall that x = (a + a † )/ √ 2. Consequently, the nonlinear term in Eq. (15) for the dynamics of a becomes −γ 2 x 2 a /4, which in the semiclassical limit, where x 2 a x 2 α, yields in place of Eq. (16). Finally, by taking the real and imaginary parts of Eq. (25), and as in Eqs. (19), differentiating Eq. (26a) with respect to time, and substituting Eq. (26b) for ṗ , we obtain a second-order equation of motion for x of the form which to within corrections of O 2 , is indeed the classical equation of motion (5), for the van der Pol oscillator, with η = 1 and ζ = 0. Figure 5 shows the steady-state Wigner functions that are obtained numerically from the vdP master equation (24) for different values of at T = 0. A comparison with the phase-space distributions of 10 4 classical van der Pol oscillators at T = 0.1, confirms that for small values of the quantum and classical models agree very well. For large values of , the quantum master equation clearly deviates from the classical vdP behavior, as expected, yet it retains the non-circular limit cycles associated with the relaxation-oscillation behavior of largeamplitude vdP oscillators.
The rate equations for the transformed density-matrix elements (22), obtained from the vdP master equa-  (24), with A c = 2 and γ 1 = 0. All the odd diagonals are free to decay to zero, while the even diagonals, which are coupled to the principal diagonal, are not. Compare with Fig. 4 above for the quantum RvdP oscillator.  Fig. 6, but for the quantum Rayleigh master equation (29), with A c = 2 √ 3 and γ 1 = 0. Compare with Fig. 4 above for the quantum RvdP oscillator. tion (24), take the forṁ n,m = κ 1 (n + m) n−1,m − 1 2 (2n + m + 2) n,m One can see that matrix elements on the m th diagonal are now coupled to elements from the m ± 2 diagonals, thus coupling the even diagonals to each other, and the odd diagonals to each other. Given the fact that the principle m = 0 diagonal cannot decay to zero, the rate equations (28) feed the even diagonals that are coupled to it, generically hindering their decay in the steady state. This is demonstrated numerically in Fig. 6, where we plot the Wigner functions and the absolute values of the density-matrix elements for different values of , while keeping the ratio between and γ 2 , and therefore A c , constant. For small values of the coupling between the off-diagonals is relatively weak, making the density matrix nearly diagonal and the Wigner function nearly circular. Increasing γ 2 increases the coupling between the even diagonals, which become non-zero. Note that the odd diagonals, which are not coupled to the principle diagonal do vanish in the steady state. Compare with the corresponding Fig. 4 for the RvdP master equation (12), where the limit cycles remain circular and the density matrix remains diagonal, even for large γ 2 .

C. The Quantum Rayleigh Oscillator
For completeness, let us present a quantum model whose classical limit at T = 0 yields the classical Rayleigh oscillator of Eq. (5), with η = 0 and ζ = 1, to within corrections of O 2 . To do so, we change the Lindblad operator for the nonlinear damping term in Eq. (12) from . The Rayleigh master equation is then given bẏ After making this change to the nonlinear term in the master equation, the nonlinear term in Eq. (15) for the dynamics of a becomes −γ 2 x 2 a + 2 p 2 a /4, which in the classical limit yields an amplitude equation of the form in place of Eq. (16). Finally, by taking the real and imaginary parts of (30), and as in Eqs. (19), differentiating Eq. (31a) with respect to time, and substituting Eq. (31b), we obtain a secondorder equation of motion for x of the form which up to corrections of O 2 , is the classical Rayleigh equation, given by Eq. (5) with η = 0 and ζ = 1. Wigner functions for the quantum Rayleigh oscillator are plotted in Fig. 7 alongside their density matrices, and comparisons between quantum and classical limit cycles are shown in Fig. 8.

D. Correlations and spectral distributions
It is convenient to consider time correlation functions of various operators, along with their Fourier spectral distributions [32], in order to characterize the different quantum limit-cycle models. It is not our intention to provide a thorough analysis of these quantities here, but only to demonstrate that the models do differ in their dynamics. To compare the models side by side we use parameters that generate limit cycles with equal amplitudes A, maintaining the same κ 1 and γ 1 , and varying γ 2 accordingly, thus setting γ vdP 2 = 4γ RvdP 2 and γ Rayleigh 2 = γ vdP 2 /3. In all the examples shown here we initiate the dynamics with the steady state density matrices, thus following the decay of correlations, while the oscillators are already in their steady state. Recall that we are still operating at T = 0, thus the decay of correlations, which results from noise-induced phase diffusion, is caused by quantum rather than thermal fluctuations. Fig. 9 shows the displacement correlation function x(t)x(0) , along with its spectral distribution, for limit cycles of moderate amplitude A = √ 10 and different driving strengths = κ 1 , with γ 1 = 0. We see that for very small , where the steady-state limit cycles are all circular, the relaxation dynamics are also very similar, with the correlations for the RvdP oscillator decaying only slightly slower than for the other two oscillators. Correspondingly, the RvdP spectral peak at ω = 1 is slightly  sharper. Recall that the RvdP oscillator is the only one that performs exact simple harmonic motion at frequency ω = 1, for any value of . As increases, as shown in Figs. 5 and 8, the vdP and Rayleigh limit cycles deviate from perfect circles, and the differences between the three spectral peaks become more evident.
Because the steady-state density matrices of the vdP and Rayleigh oscillators contain non-zero elements in x 2 (t)x 2 (0) (bottom panels), and their corresponding spectral distributions S x 2 x 2 (ω) (top panels), calculated for the same parameters as in Fig. 9 for the three quantum models. FIG. 11: Bottom panels: Two-phonon correlation function a † 2 (t)a 2 (0) . Top panels: The corresponding spectral distributions S a 2 a 2 (ω). Calculated for the same parameters as in Fig. 9 for the three quantum models.
their even off-diagonals, as shown in figures 6 and 7, it is interesting to examine the squared-displacement correlation function x 2 (t)x 2 (0) , shown in Fig. 10, whose calculation involves simultaneous creation or annihilation of pairs of phonons. Again, at very small , the behavior is quite similar in all three models, showing two spectral peaks, at ω = 0 and ω = 2, as expected from the squaring of x(t). As increases, deviations between the models quickly become noticeable. In particular, note the rather large shift to higher frequencies of the spectral peak of the vdP oscillator, as the limit cycle becomes less and less circular. Also note the different asymptotic values of the correlation functions, which tend to x 2 (0) 2 , which is greatest for the Rayleigh oscillator owing to its larger r.m.s. displacement (see the classical limit-cycle shapes in Fig. 1). Compare the squared-displacement correlation function with the two-phonon correlation a † 2 (t)a 2 (0) , shown in Fig. 11. This is only one of the terms appearing in the calculation of the squared-displacement correlation function, annihilating two phonons at time equal 0, and recreating them at a later time t. It directly probes the m = 2 off-diagonal of the density matrix [10], which are non-zero for the vdP and Rayleigh models. Indeed, for these models the two phonon correlations decay to a nonzero asymptotic value a † 2 (0) a 2 (0) , which is greater for the vdP oscillator. As a final example, we wish to demonstrate that both κ 1 and γ 1 , and not only their difference affect the dynamics of the oscillators as independent sources of quantum noise, especially as one approaches the quantum regime. For this purpose we return to the displacement correlation function, and consider limit cycles of a smaller amplitude A = 0.1, weakly driven with = 0.01, at γ RvdP 2 = 1 . The correlation functions and their spectral distributions are shown in Fig. 12, for three different values of the ratio r = γ 1 /κ 1 between the linear damping and the pumping rates. One sees very clearly that the RvdP oscillator exhibits a much slower decay of its displacement correlation function, as well as a higher sensitivity to the value of r, than the other two oscillators. To see why this is so, consider the rate equations (28) for the density matrix elements of the vdP oscillator, and notice that all the off-diagonal (m = 0) terms have some negative coupling due to the nonlinear damping γ 2 , which causes these terms quickly to decay for the large values of γ 2 in the quantum limit. The same holds for the Rayleigh oscillator. On the other hand, inspection of the rate equations (23) for the the RvdP density matrix, reveals that it has a single off-diagonal element ρ 0,1 that does not have a negative coupling term proportional to γ 2 . Assuming that the large nonlinear damping rate quickly depletes all other off-diagonal matrix elements, one remains with this last element, whose decay rate, which is governed by the much smaller rates γ 1 and κ 1 , is indifferent to the nonlinear damping rate, and increases as r approaches 1. Also note that the contributions of noise in the energy pump and noise in the linear damping mechanism to the decay rate are additive. In particular, the decay rate at the bifurcation, where = 0, tends to 2κ 1 = 2γ 1 . Thus, the oscillator experiences critical slowing down as it crosses the bifurcation only if κ 1 and γ 1 are both zero, which may be difficult to arrange experimentally.

IV. ANALYTICAL SOLUTION FOR THE STEADY-STATE DENSITY MATRIX OF THE RAYLEIGH-VAN DER POL OSCILLATOR
An analytical solution for the steady state of the T = 0 quantum RvdP oscillator can be found in previous work [13,33,38,39], along with approximate solutions for T ≥ 0 in the limit of k B T ω [13,33]. Here we provide a general analytical solution for arbitrary temperature. In doing so, we consider a slightly more general physical system than the one described by our master equation (12) above, by adding to the model a process of two phonon absorption at rate κ 2 . This additional process, while only recently demonstrated in a micromechanical system [40], might be quite relevant for other physical systems, such as optical ones, where two-photon absorption might be as likely as two-photon emission.
The revised temperature-dependent master equation is then written aṡ where the last line is responsible for two phonon absorption, andn(ω) = (e ω/kBT − 1) −1 is the Bose-Einstein distribution through which the temperature T is introduced. We start by defining four temperature-dependent effective rates, that reduce back to the original rates in the limit of T → 0, Using these, we rewrite the revised RvdP master equation (34) more compactly as, As discussed earlier, the off-diagonal elements of the RvdP density matrix decay to zero in the steady state, as they are decoupled from the principal diagonal. The remaining rate equations for the diagonal elements P n ≡ n,0 = ρ nn = n|ρ|n are given by where we have rescaled all the rates by the nonlinear damping coefficientΓ 2 , In the steady state, withṖ n = 0, the set of equations (37) provide recurrence relations for the Fock-state probabilities P n , giving the steady-state value of each level in terms of the four levels preceding it. Dykman [33] and others [13,38,39] solve these recurrence relations for special limiting cases, by using the method of generating functions, which yields a second-order differential equation for the generating function. We use the same method here, but before doing so we note that when summing consecutive rate equations, one obtains a telescopic sum in which many terms cancel out. Thus by summing the first n+1 equations (37), from 0 to n, and dividing by an overall factor of (n + 1), we obtain a simpler equation to solve, where the maximum power of n is 1 rather than 2, reducing the corresponding differential equation from second to first order. We solve these simplified recurrence relations using the generating function By multiplying Eq. (39) by x n+1 , and summing from n = 0 to ∞, we replace the infinite set of recurrence relations with a single differential equation with respect to the auxiliary variable x, This nonhomogeneous first-order differential equation can be solved in a standard manner, using an integrating factor. It should be noted, though, that the apparent constant term on the right-hand side of the equation depends linearly on the solution itself, with P 0 = A(0) and P 1 = A (0) [41]. Therefore, the solution of the associated homogeneous equation as well as any particular solution of the full nonhomogeneous equation are both determined only to within a multiplicative factor. As a consequence, the space of solutions is a 2-dimensional vector space, and we still require two constraints, or boundary conditions, to pin down the physically relevant solution. We shall use the fact that the coefficients P n in the expansion (40) of A(x) are probabilities. As such, their values are constrained to be between 0 and 1; they are normalized such that their sum ∞ n=0 P n = A(1) = 1; (42) and their alternating sum lies a distance not greater than unity away from the origin, ∞ n=0 (−1) n P n = |A(−1)| ≤ 1.
Before solving Eq. (41) we perform the substitution with where for convenience we set a = √ K 2 . After some algebra, we obtain a differential equation for f (z) of the form which after defining b = aΓ 1 + K 1 2a (1 − a) , and c = becomes and we remember that C 1 is a constant to be determined through the boundary conditions. The general solution to this equation is given by where C 2 is a constant of integration multiplying the solution of the associated homogeneous equation, and is the hypergeometric function, where (x) n is the socalled Pochhammer symbol, denoting the rising factorial, Because (1) n = n!, the expansion (50) in our case reduces to and the solution (49) can equivalently be expressed as is the incomplete beta function. Although the power series (52) diverges for |z(x)| ≥ 1, we need only to evaluate its derivatives at x = 0, and the condition |z(0)| < 1 is fulfilled as long as the nonnegative parameter a < 1. Recall that a 2 = K 2 is the ratio between nonlinear absorption and emission, thus the physical interpretation of a < 1 is that there is no steady-state solution when the nonlinear gain is stronger than the nonlinear damping, which is indeed the case.
In terms of the original variable, the solution (49) for the generating function becomes where the new constants D 1 and D 2 still need to be determined. Clearly, for c > 1, the solution of the associated homogeneous equation has a singularity at x = −1, in contradiction to the condition of Eq. (43) that the alternating sum be bounded, requiring us to set D 2 = 0. The normalization condition (42) then yields the final form of the generating function where the normalization constant With the generating function at hand, we can calculate the probabilities where f (k) (x) denotes the k th derivative of f (x). Finally, the derivatives of 2 F 1 can be evaluated using the relation [42, see their equation (5.2.2)] to give In Fig. 13 we plot the analytical solution given by Eq. (60), alongside numerical calculations of the steadystate solutions of the temperature-dependent RvdP master equation (34) for different parameter values, showing perfect agreement.
As noted above, previous authors [13,33,38,39] used Eq. (37) directly, without the telescopic sum (39). Instead of our first-order nonhomogeneous equation (41), they obtained a homogeneous second-order differential equation of the form Differentiation of the first-order equation (41) yields this second-order equation (61), thus solutions to the firstorder equation solve the second-order equation as well.
One can obtain the solution (63), for the special case of zero temperature and no two-phonon absorption, from our general solution (56) by taking the limit of a → 0, where In this limit, the hypergeometric function 2 F 1 (a, b; c; z) reduces to the confluent hypergeometric function such that

V. CLASSICAL TO QUANTUM TRANSITION OF THE RAYLEIGH-VAN DER POL OSCILLATOR
We saw in section III A that when taken to its classical limit, with A c 1, the quantum RvdP limit-cycle forms self-oscillations at an amplitude given by A c . However, as A c is reduced toward unity, in terms of the quantum unit of length x 0 , the oscillator approaches zero-point motion, and one expects quantum effects to take place. In this section we examine how the limit cycle behaves as the oscillator transitions into this quantum regime. Figures 14 and 15 show the steady-state Wigner functions and Fock-state distributions of the quantum RvdP limit-cycle of Eq. (12), for different values of A c . Figure 14 shows that when coming from the classical regime, by reducing A c from 8 down to 2, the radius of the limit cycle is approximately A c , and relatively many Fock states are populated. On the other hand, Figure 15 shows that when entering the quantum regime, as A c is lowered further from 1 down to 0.1, only a few Fock states are populated, and the radius of the limit cycle does not get much smaller than x zp = x 0 / √ 2. To see this more quantitatively, we follow Steiner [44] and sum all the even, or alternatively all the odd, rate equations (37) for the Fock-state probabilities, thereby telescopically eliminating all the two-phonon transitions, and finding that in the steady state (69) In the quantum limit of large nonlinear damping γ 2 , or small A c , and low temperature k B T ω, both Γ 1 and K 1 , defined in Eq. (38), tend to zero with corrections of order γ −1 2 , while K 2 tends to zero with corrections of order γ −1 2 or exp{−2 ω/k B T }. An inspection of the first few rate equations (37) then shows that all P n , with n > 1, are smaller than P 0 and P 1 at least by an order of γ −1 2 or exp{−2 ω/k B T }. Neglecting all these higher states, with n > 1, in Eq. (69) then yields a relation between the occupation probabilities of the remaining two lowest states, given by Thus, with Tr{ρ} = 1, in the low temperature quantum limit, with γ 2 → ∞, we find that the density matrix becomes where the temperature-dependent ratio tends to the bare ratio r ≡ γ 1 /κ 1 of the linear damping rate to the pumping rate, when T → 0. Also note that R approaches 1 as T increases, it is equal to 1 if and only if r = 1 at arbitrary T , and it approaches exp{− ∆ 1 /k B T } for fixed T as r tends to zero.
As was previously understood, in this limit only the |0 and |1 states are occupied, because all phonons in any other state are immediately annihilated by the infinitely strong nonlinear damping. But, contrary to the zero-temperature result of previous authors [4][5][6][7][8][9][10], who take γ 1 = 0 above the bifurcation, we find that the actual occupation depends on the ratio r = γ 1 /κ 1 , and is not universal. This is demonstrated numerically in Fig. 16, where we compare the steady-state zero-temperature Fock-state distributions of the RvdP oscillator with thermal distributions for the same average phonon occupation, which according to Eq. (71) is given by Cross sections through the corresponding Wigner functions for the same parameter values are shown in Fig. 17, where one can observe the onset of the bifurcation at N = 1/4 as the mean phonon number N gradually increases from 0 to 1/3. Note the quantitative differences between the Wigner functions that appear even below the bifurcation.
In the case of the quantum RvdP oscillator, we choose to associate the amplitude A of the limit-cycle oscillations with the maxima |α| max of its circular Wigner function, which we evaluate either numerically or using Eq. (21), while recalling the factor of √ 2 which arises from the definition of Eq. (11). For the extreme quantum-limit steady-state density matrix of Eq. (71) this yields W q (α, α * ) = 1 π (74) whose maximum determines the limit-cycle amplitude where in the zero-temperature limit, the ratio R = Γ 1 /K 1 appearing in Eqs. (74) and (75) is replaced by r = γ 1 /κ 1 . Note that the bifurcation occurs at R = 1, which according to Eq. (72) happens if and only if r = 1 regardless of the temperature. In the case of finite γ 2 , we expect these expressions to have corrections of O γ −1 2 , as higher Fock states become populated.
In the zero-temperature quantum limit, the Wigner functions still exhibit a clear bifurcation to selfoscillations with an amplitude that grows continuously from zero, as the κ 1 = γ 1 threshold is crossed. Nevertheless, the nature of this bifurcation is quite different from the classical Hopf bifurcation. In the classical regime, one expects the amplitude of steady-state oscillations to scale as the square root of the reduced pumping, A c = /γ 2 , where = κ 1 − γ 1 , and therefore for the oscillations to die out for infinite nonlinear damping (unless the pumping rate κ 1 is infinite as well). This is shown by a straight black line in Fig. 18. However, in the quantum regime, the |1 state is protected from nonlinear damping, which enables the oscillator to undergo a  : Amplitude A of the RvdP limit cycle as a function of temperature, for different values of r in the quantum limit, with κ 1 = 1 and γ 2 = 10 5 , and for pump detunings of (a) ∆ 1 = 0.1, and (b) ∆ 1 = 1. Numerical values (scattered points), obtained by solving the steady-state master equation (12), are compared with the approximate expression of Eq. (75) (solid lines), showing good agreement at low temperatures, particularly for small detuning. As the temperature increases, and R approaches 1, the amplitude decreases to zero.
bifurcation into self-oscillations, at an amplitude given by Eq. (75), even when the nonlinear damping is infinitely strong. The linear pumping rate κ 1 need only be large compared to the linear damping rate γ 1 . This is purely a quantum effect. Accordingly, as we noticed earlier in Fig. 15, as γ 2 tends to infinity rather than decaying to zero as /γ 2 , the zero-temperature steady-state amplitude saturates at (1 − r)/2 = /2κ 1 . This is demonstrated numerically by the colored curves in Fig. 18 for a few values of the ratio r.
This quantum effect is somewhat smeared out when temperature is turned on and the amplitude saturates at (1 − R)/2, rather than (1 − r)/2, decreasing with FIG. 20: Amplitude A of the RvdP limit cycle with γ 1 = 1 as a function of κ 1 = 1/r for different temperatures, in (a) the quantum limit with γ 2 = 10 5 , and (b) the classical limit with γ 2 = 1. The temperature seems to have no effect on the overall shape of the curves in the classical limit, whereas in the quantum limit it causes the amplitude to saturate at lower values as κ 1 increases or r decreases. Numerical values, obtained by solving the steady-state master equation (12), are compared in panel (a) to solid lines showing the infinite γ 2 , low temperature, approximate solution of Eq. (75).
temperature towards zero, as R increases from r towards 1. This is confirmed numerically in Fig. 19, showing the oscillation amplitude in the quantum limit decaying to zero as the temperature increases. As expected, the approximate expression of Eq. (75) holds better at low temperatures and for small pump detuning ∆ 1 .
In the limit of r → 0, as κ 1 increases or γ 1 decreases, the infinite-γ 2 oscillation amplitude tends to with an exponential dependence on temperature. This is demonstrated in Fig. 20(a) for γ 2 = 10 5 , while Fig. 20(b) shows essentially no temperature dependence of the amplitude in the classical limit with γ 2 = 1. A closer inspection of this exponential temperature dependence for r = 0.1 is shown in Fig. 21, where we plot the Fockstate distributions and Wigner-function cross sections, for T ≤ 0.5. One can see how the increase in temperature gradually smears out the limit cycle. On one hand, as can be infered form Eq. (71), the increase in R causes an increase of the occupation probability P 0 of the |0 state, while at the same time increasing the neglected corrections of O(exp{−2 ω/k B T }) in the form of nonzero occupation probabilities of the |2 and |3 states.

VI. CONCLUSIONS
We have studied a collection of master equations that yield quantum limit cycles in their steady-state dynamics. They all describe a simple harmonic oscillator, inter- acting with the environment through a combination of Lindblad operators, responsible for linear and nonlinear damping and energy injection, or pumping, in the form of single-phonon or double-phonon emission and absorption processes. We have established the correct correspondence between these quantum master equations and their classical counterparts, noting that the commonly used quantum model-which is symmetric under phasespace rotations and therefore always yields circular limitcycles-is often mistaken to be the "van der Pol (vdP) oscillator", even though it actually corresponds to the classical "Rayleigh-van der Pol (RvdP) oscillator". We have also noted that, in all cases, the correspondence holds only for oscillations just above the bifurcation, namely, only to first order in the bifurcation parameter . We have analyzed a generalized version of the quantum RvdP limit cycle, applicable to a broad range of physical systems, such as nanomechanical oscillators, optical oscillators or lasers, electronic or superconducting oscillating circuits, and cold ions. We have obtained an exact analytical solution to the master equation in its steady state for arbitrary temperature, and considered its small-amplitude quantum limit-obtained by increasing the nonlinear damping rate-in some detail. A number of features emerge in this quantum regime, some of which were previously overlooked. Most important is the fact that, at T = 0, the |1 state of the quantum oscilla-tor is protected from nonlinear damping. One therefore still obtains limit-cycle oscillations, even with an infinite nonlinear damping rate, yet these quantum limit cycles are strongly affected by both the linear damping and the pumping rates, and are not universal as previously believed. We show that whereas in the classical regime it is only the difference between the linear pumping and the linear damping rates that affects the zero-temperature dynamics, in the quantum regime the ratio of the two rates plays a significant role as well, as they each contribute an independent source of spontaneous quantum processes. We have also described the effect of temperature in smearing out these nonclassical bifurcations.
We have performed a numerical comparison between classical and quantum dynamics of the different models, showing perfect correspondence-where expectedbetween the quantum Wigner functions and the corresponding classical phase-space distributions. The agreement holds not only for the steady-state limit cycle dynamics, but for the transients as well, whereby an initial oscillating coherent state first quickly relaxes, or drifts, to the expected amplitude, and only then slowly diffuses around the limit cycle losing its initial phase. Deviations between the two occur in the quantum regime, as just mentioned above, where rather than decaying to zero as nonlinear damping increases, the quantum limit-cycle is protected, with its amplitude saturating at around zeropoint motion, at a value that depends on the ratio of the linear pumping and damping rates. Deviations also occur far above the bifurcation, where the quantum and classical models no longer agree with each other. It should be emphasized that the Wigner functions that describe all the limit cycles are "essentially classical", developing no negative regions for any choice of parameters. This is a well-known property of the simple harmonic oscillator, which persists in these open systems, as long as the oscillator is linear [45] and is uncoupled to additional oscillators or other degrees of freedom.
Our results should provide a firmer theoretical basis for ongoing studies of physical phenomena such as quantum entrainment and synchronization, and more generally, nonequilibrium nonlinear quantum dynamics involving self-sustained oscillators. We hope that our analytical results could be tested experimentally in the near future, where they should provide better tools with which to analyze the measured data.