## Critical Response of a Quantum van der Pol Oscillator

Shovan Dutta ∗ and Nigel R. Cooper †

T.C.M. Group, Cavendish Laboratory, University of Cambridge, JJ Thomson Avenue, Cambridge CB3 0HE, United Kingdom

(Dated: December 23, 2019)

Classical dynamical systems close to a critical point are known to act as efficient sensors due to a strongly nonlinear response. We explore such systems in the quantum regime by modeling a driven van der Pol oscillator. We find the classical response survives down to one excitation quantum. At very weak drives, genuine quantum features arise, including diverging and negative susceptibilities. Further, the linear response is greatly enhanced by using a strong incoherent pump. These results are largely generic and can be probed in current experimental platforms suited for quantum sensing.

A key insight from the theory of phase transitions is that systems close to a critical point are highly sensitive to perturbations [1]. This is evident, for example, in a diverging compressibility near a liquid-gas phase transition. Such divergences in susceptibility also occur in classical dynamical systems close to a bifurcation [2]. A particularly important prototype is an oscillator with nonlinear damping, which can transition from a dormant to a limitcycle phase across a Hopf bifurcation [3]. Generically, the response of such a van der Pol (vdP) oscillator [4] to a resonant drive Ω grows as Ω 1 / 3 at the critical point [5]. This nonlinearity is what enables the ear and other biological sensors to detect very weak stimuli and maximally process environmental inputs [5-9].

of energy scales and should generalize to other systems. After briefly reviewing the iconic classical vdP oscillator, we present results for the quantum response, concluding with a discussion of experimental realizations.

In this Letter, we ask if this increased sensitivity persists into the quantum regime. This could have important applications in the growing frontier of quantum sensing [10, 11]. Besides, understanding the dynamics of a quantum vdP oscillator is also of fundamental interest, as it represents an easy-to-implement nonequilibrium setting where coherent drive competes with incoherent dissipation. Such open quantum systems are at the forefront of modern physics [12, 13], led by advances in synthetic experimental platforms with unprecedented control and tunability [14-18].

At the semiclassical level, the vdP oscillator is realized in the physics of lasers [19, 20]. Recent theoretical studies have also examined a quantum vdP oscillator, largely in the context of synchronization [21-27]. However, its critical properties and response are as yet unexplored. Here, we fully characterize a driven quantum vdP oscillator by a standard master equation, finding surprising features which can be probed in experiments.

We find the classical response persists all the way down to one excitation quantum. At weaker drives, the quantum oscillator exhibits both divergent as well as negative susceptibilities. Further, in the limit-cycle phase, the linear response is only limited by two-body loss, providing a strong sensitivity enhancement over a passive system. Several of these features originate from a competition

∗ E-mail: sd843@cam.ac.uk

† E-mail: nrc25@cam.ac.uk

The classical vdP model describes a harmonic oscillator with nonlinear damping, ¨ = x -ω x 2 +(2 γ 1 -8 γ x 2 2 ) ˙ . x Here, ω is the natural frequency, γ 1 is a negative damping which may arise from an energy source, and γ 2 is the leading nonlinear term which damps the system at large x . For | γ 1 | glyph[lessmuch] ω , the small-amplitude oscillations are sinusoidal with a slowly-varying amplitude, x = Re[ α t e ( ) i ωt ], where ˙ α = γ α 1 -γ 2 | α | 2 α [4]. Thus, if γ 1 &lt; 0, the system is fully damped and all oscillations die out, whereas for γ 1 &gt; 0, they reach a stable amplitude with | α | = √ γ /γ 1 2 . When subjected to a resonant drive F cos( ωt + φ ), the equation of motion gains a forcing term,

$$\dot { \alpha } = \gamma _ { 1 } \alpha - \gamma _ { 2 } | \alpha | ^ { 2 } \alpha + \Omega$$

where Ω = -i e i φ F/ (2 ω ). We choose φ = π/ 2 so that Ω is real and nonnegative. From Eq. (1), the steady-state response at the critical point γ 1 = 0 is purely nonlinear, α = (Ω /γ 2 ) 1 / 3 . Hence, the susceptibility χ cl ≡ d α/ dΩ is divergent at zero drive, χ cl ∝ Ω -2 / 3 . Away from criticality, the response is of the form α = ˜ Ω 1 / 3 f (˜ γ / 1 ˜ Ω 2 / 3 ) where ˜ γ 1 ≡ γ /γ 1 2 and Ω ˜ ≡ Ω /γ 2 (details in Supplement [28]). In particular, the limiting behavior is such that

$$\Gamma _ { \ell }, \quad \alpha = \begin{cases} \sqrt { \gamma _ { 1 } / \gamma _ { 2 } } + \Omega / ( 2 \gamma _ { 1 } ) & \text{for } \tilde { \gamma } _ { 1 } \gg \tilde { \Omega } ^ { 2 / 3 } \\ \Omega / | \gamma _ { 1 } | & \text{for } \tilde { \gamma } _ { 1 } \ll - \tilde { \Omega } ^ { 2 / 3 } \end{cases}, \quad ( 2 )$$

which are consistent with a diverging susceptibility, χ cl ∝ 1 / γ | 1 | , at the critical point γ 1 → 0.

We study a quantum version of the vdP model, where a quantum harmonic oscillator is subjected to several dissipative processes, as in Refs. [21-27]. These are (i) oneparticle loss with rate γ -1 , (ii) one-particle gain with rate γ + 1 , and (iii) two-particle loss with rate γ 2 . As we will discuss later, such processes can be engineered in current experimental setups by coupling the oscillator to a suitable environment. Additionally, the oscillator is subject to a resonant drive of amplitude Ω, given by the Hamiltonian ˆ = H ωa a ˆ ˆ + Ω(i † e i ωt ˆ a † +H.c.). Here, ˆ annihilates a para ticle in the oscillator mode and glyph[planckover2pi1] = 1. It is convenient to

FIG. 1. Susceptibility χ and response 〈 ˆ a 〉 of a resonantly driven van der Pol oscillator at the critical point γ ± 1 = 0. Here, Ω is the drive amplitude and γ + 1 , γ -1 , and γ 2 are the rates of one-particle gain, one-particle loss, and two-particle loss, respectively. The classical and quantum responses match down to 〈 ˆ a 〉 ∼ 1, below which quantum fluctuations cut off the classical divergence, causing χ to saturate.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-f7f83cb2-4310-4f36-8143-dac349691d58.jpg)

work in the rotating frame where the Hamiltonian reads ˆ H = iΩ(ˆ a † -ˆ). a In general, the full dynamics of such a driven-dissipative system are governed by a master equation for the density matrix ˆ, obtained by tracing out the ρ environment degrees of freedom [29]. However, for typical atomic/photonic setups, the environment relaxes to equilibrium on optical timescales, several orders of magnitude faster than the system dynamics [30]. Under such a routine Markov approximation [31, 32], the master equation reduces to the Lindblad form

$$\dot { \hat { \rho } } = - \text{i} [ \hat { H }, \hat { \rho } ] + \gamma _ { 1 } ^ { + } \mathcal { D } [ \hat { a } ^ { \dagger } ] \hat { \rho } + \gamma _ { 1 } ^ { - } \mathcal { D } [ \hat { a } ] \hat { \rho } + \gamma _ { 2 } \mathcal { D } [ \hat { a } ^ { 2 } ] \hat { \rho } \,, \ \ ( 3 ) \quad \text{from} \, \mathfrak { c }$$

where D [ˆ]ˆ x ρ ≡ ˆˆˆ xρx † -{ ˆ ˆ x x, ρ † ˆ } / 2 [33, 34]. Using Eq. (3), or operator equations of motion [35], the expectation value 〈 ˆ a 〉 can be shown to obey

$$\langle \dot { \hat { a } } \rangle = [ ( \gamma _ { 1 } ^ { + } - \gamma _ { 1 } ^ { - } ) / 2 ] \langle \hat { a } \rangle - \gamma _ { 2 } \langle \hat { a } ^ { \dagger } \hat { a } \hat { a } \rangle + \Omega \,. \quad ( 4 ) \quad \text{ of er}$$

Replacing ˆ with a complex number a α in Eq. (4) leads to the classical limit in Eq. (1) with γ 1 = ( γ + 1 -γ -1 ) / 2. Hence, the critical case corresponds to γ + 1 = γ -1 . Note that although Eq. (1) was derived in the small-amplitude limit, we will refer to Eq. (3) more generally as a quantum vdP oscillator, similar to Refs. [21-27].

We find the response by numerically solving Eq. (3) for the steady-state density matrix, then computing 〈 ˆ a 〉 = Tr(ˆˆ). It is most straightforward to use the Fock basis, aρ ˆ = ρ ∑ nn ′ ρ n,n ′ | n 〉〈 n ′ | , where n = 0 1 2 , , , . . . denotes the number of particles in the oscillator mode. The response is given by 〈 ˆ a 〉 = ∑ n √ nρ n,n -1 . Past studies have found closed-form expressions for ρ n,n ′ in special cases [36-38].

First we consider the simplest critical case, γ + 1 = γ -1 = 0. Figure 1 shows the susceptibility χ ≡ d ˆ 〈 a / 〉 dΩ as a function of drive. Note that it coincides with the classical result for 〈 ˆ a 〉 glyph[greaterorsimilar] 1. However, at weaker drives, the classical divergence is cut off and χ saturates, producing a linear response. We can understand this low-energy cutoff

FIG. 2. Zero-drive susceptibility as a function of damping γ -1 (with γ + 1 = 0, red axis) and antidamping γ + 1 (with γ -1 = 0, blue axis). Both classical and quantum susceptibilities diverge as 1 /γ ± 1 , although at γ ± 1 = 0, the quantum susceptibility is finite (black circle, also see Fig. 1). The inset shows transitions among the lowest three oscillator levels for γ + 1 , Ω glyph[lessmuch] γ 2 and γ -1 = 0. Since any two-particle excitation decays very rapidly, the dynamics are confined to this low-energy subspace.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-5eda03df-30b4-4a31-89bf-f4449b5e390d.jpg)

as follows. At Ω = 0, both | 0 〉 and 1 | 〉 are steady states as there is only two-particle decay. A nonzero drive couples these neighboring Fock states, yielding ρ 10 , ρ 21 ∼ Ω and ρ 22 ∼ Ω . Hence, to linear order, one can restrict the dy2 namics to the lowest three levels, which gives 〈 ˆ a 〉 = 2Ω /γ 2 (see Supplement [28]). Stronger drives inject more particles, inducing a crossover to the classical limit.

We now consider departures from γ ± 1 = 0 by allowing either γ + 1 or γ -1 to be nonzero, spanning the transition from quiescent ( γ + 1 = 0, γ -1 &gt; 0) to the limit-cycle regime ( γ + 1 &gt; 0, γ -1 = 0). In Fig. 2 we plot the zero-drive susceptibility across the transition, which, remarkably, shows a divergence as 1 /γ ± 1 at small damping or antidamping. This divergence can be understood from a competition of energy scales. For γ + 1 = 0 and γ -1 &gt; 0, the system is fully damped and the undriven steady state is | 0 . 〉 For weak drives Ω &lt; γ -1 , the occupation remains small, so the two-particle decay is irrelevant. Thus, we can discard the nonlinear term in Eq. (4), which then reduces to a driven damped oscillator with steady state 〈 ˆ a 〉 = 2Ω /γ -1 , matching the classical result. The response for γ + 1 &gt; 0 and γ -1 = 0 is more involved, as the oscillator can have large occupations even at zero drive. In general, the linear response is given by Liouvillian perturbation theory [39-41]. However, for γ + 1 glyph[lessmuch] γ 2 , the dynamics are confined to the lowest three levels, as shown in Fig. 2 inset, which yields 〈 ˆ a 〉 = 2Ω (9 / γ + 1 ) (see Supplement [28]).

At first sight, the results in Figs. 1 and 2 might appear to be contradictory: at γ + 1 = γ -1 = 0, the zero-drive susceptibility is finite (black circle in Fig. 1); however Fig. 2 shows that it diverges as the critical point is approached from either side. The resolution is that the susceptibility χ (Ω → 0 , γ ± 1 → 0) depends on the order of the two limits. This is better understood by plotting the full response as a function of Ω for small γ ± 1 , as shown in Fig. 3. The re-

FIG. 3. Top: response 〈 ˆ a 〉 as a function of drive Ω for γ + 1 = 0 and γ -1 /γ 2 = 0 02. . The quantum response is linear for Ω glyph[lessorsimilar] γ -1 and shows negative susceptibility for γ -1 glyph[lessorsimilar] Ω glyph[lessorsimilar] ( γ -1 γ 2 ) 1 / 2 . For Ω glyph[greaterorsimilar] γ 2 , it agrees with the classical response, as in Fig. 1. Bottom: rescaled Wigner functions W α ( ) describing a quasiprobability distribution in phase space [28]. The response 〈 ˆ a 〉 is given by the center of mass of W α ( ), shown by black dots.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-9df28c32-67bf-4e8b-a8ea-ffd2e1c00ab2.jpg)

sponse is strongly nonlinear for Ω glyph[greaterorsimilar] γ ± 1 . Thus, although the linear susceptibility diverges, the linear region itself shrinks to zero as γ ± 1 → 0.

The full response exhibits four different regimes: (i) a linear response for Ω glyph[lessorsimilar] γ ± 1 , (ii) negative susceptibility for γ ± 1 glyph[lessorsimilar] Ω glyph[lessorsimilar] ( γ ± 1 γ 2 ) 1 / 2 , (iii) an extended quantum response for ( γ ± 1 γ 2 ) 1 / 2 glyph[lessorsimilar] Ω glyph[lessorsimilar] γ 2 , and (iv) classical response for Ω glyph[greaterorsimilar] γ 2 . The first three regions are well reproduced by a three-level approximation, which yields, to lowest order in γ ± 1 /γ 2 and Ω /γ 2 (see Supplement [28] for a derivation),

$$\langle \hat { a } \rangle \approx \frac { 2 \Omega } { \gamma _ { 2 } } \left [ \frac { ( \gamma _ { 1 } ^ { + } + \gamma _ { 1 } ^ { - } ) \gamma _ { 2 } + 8 \Omega ^ { 2 } } { ( 3 \gamma _ { 1 } ^ { + } + \gamma _ { 1 } ^ { - } ) ^ { 2 } + 8 \Omega ^ { 2 } } \right ]. \quad \ \ ( 5 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Hence, in the linear region, χ ≈ 2( γ + 1 + γ -1 ) / (3 γ + 1 + γ -1 ) 2 , which diverges for γ ± → 0, whereas for Ω ∼ ( γ ± 1 γ 2 ) 1 / 2 , 〈 ˆ a 〉 ≈ 2Ω /γ 2 + ( γ + 1 + γ -1 ) / (4Ω), which exhibits a local minimum. This nonmonotonic response is also evident in the Wigner functions shown in the lower panel of Fig. 3, which describe a quasiprobability distribution in phase space [34]. The Wigner function is widely used as an integral representation of the density matrix [28].

Physically, the separate regions originate from an interplay between drive and dissipation within a few-level manifold. Consider the case γ + 1 = 0 and 0 &lt; γ -1 glyph[lessmuch] γ 2 . A weak drive couples the steady state | 0 〉 to | 1 , produc-〉 ing a coherence ρ 10 ∼ Ω /γ -1 . However, ρ 11 also grows with Ω and saturates at 1 / 2. This saturation acts as a negative feedback for ρ 10 , as seen from the equation of

FIG. 4. Response for the critical condition γ + 1 = γ -1 = Γ 1 with Γ 1 /γ 2 = 5 (blue), 50 (red), and 1000 (green), showing the approach to classical limit (dotted black). Inset shows the sensitivity gain over a passive oscillator, G 0 ∼ √ Γ 1 /γ 2 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-01f039b0-c3a9-4036-8bfe-d058df1985fe.jpg)

motion ˙ ρ 11 = 2Ω ρ 10 -γ -1 ρ 11 . For ρ 11 ≈ 1 / 2, ρ 10 falls off as γ -1 / (4Ω), giving rise to negative susceptibility. At larger drives, higher-energy modes become accessible and coherences can grow again. Then the response is similar to the case γ ± 1 = 0, with a purely quantum regime for 〈 ˆ a 〉 glyph[lessorsimilar] 1 and a classical regime for 〈 ˆ a 〉 glyph[greaterorsimilar] 1 (cf. Fig. 1). A similar variation with drive is found for 0 &lt; γ + 1 glyph[lessmuch] γ 2 .

Note that the concomitant divergent and negative susceptibilities are genuine quantum features which result from a very general set of conditions involving a competition between coherent and incoherent energy scales. Thus, we expect them to show up generically in drivendissipative systems. For the vdP oscillator, we find they are robust to anharmonicity and detuning.

We have shown the quantum system exhibits a diverging linear susceptibility close to a critical point (Fig. 2), qualitatively similar to the classical system. Can this be advantageous in quantum sensing applications? More precisely, can one enhance the sensitivity of a damped quantum oscillator by introducing an incoherent pump? For such a passive oscillator, the equation of motion is given by ˆ = Ω ˙ a -( γ -1 / 2)ˆ, a so the steady state is 〈 ˆ a 〉 = 2Ω /γ -1 and the passive susceptibility is χ p = 2 /γ -1 . To discuss whether pumping enhances sensitivity, we define the sensitivity gain at weak drives, G 0 ≡ χ/χ p | Ω → 0 .

One regime where there is a gain, G &gt; 0 1, is the classical limit. Then the response is only a function of γ + 1 -γ -1 , and one can simply set γ + 1 = γ -1 to make the system infinitely sensitive to weak signals for any damping. Quantum mechanically, however, the response also depends on the average Γ 1 ≡ ( γ + 1 + γ -1 ) / 2. Figure 4 shows how the classical limit emerges with increasing Γ 1 for the critical case, γ + 1 = γ -1 . We see that G 0 scales as √ Γ 1 /γ 2 and one can have χ &gt; χ p if γ -1 glyph[greaterorsimilar] 3 γ 2 , which holds for typical experimental systems [42]. The relatively slow growth of G 0 with Γ 1 stems from large number fluctuations at criticality (which also limits the sensing efficiency, see below). Since γ + 1 = γ -1 , the undriven system would reach an infinite-temperature state in the limit γ 2 → 0. A

FIG. 5. (a)-(c) Classical trajectories and rescaled Wigner functions at increasing drive for γ + 1 /γ 2 = 50 and γ -1 /γ 2 = 20. The quantum and classical responses are shown by white and red dots, respectively. As the quantum system crosses over from a symmetry-broken limit cycle to the classical steady state, its susceptibility is only limited by two-body loss.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-a72e592b-871a-47bc-813c-ff6112a21064.jpg)

nonzero γ 2 results in a Gaussian number distribution of width √ Γ 1 /γ 2 [37], leading to the linear susceptibility χ ≈ 2 / √ π Γ 1 γ 2 (derivation in Supplement [28]).

The sensitivity can be enhanced much further by operating deep in the limit-cycle phase, γ + 1 -γ -1 glyph[greatermuch] γ 2 , as in Fig. 5. Here, a vanishingly small drive breaks the phaserotation symmetry, yielding a nonzero classical response. For the quantum oscillator, the coherence 〈 ˆ a 〉 builds up at a finite rate (see Supplement for a derivation [28])

$$\chi \approx \frac { 2 } { 3 } \, \gamma _ { 2 } ^ { - 1 } \left [ 1 - 2 \gamma _ { 1 } ^ { - } / ( 3 \gamma _ { 1 } ^ { + } ) \right ] \,. \quad \quad ( 6 ) \quad \text{inv} \quad.$$

Hence, the susceptibility is only limited by γ 2 and one can obtain a large enhancement. In particular, for γ + 1 glyph[greatermuch] γ -1 , G 0 ≈ γ -1 / (3 γ 2 ) glyph[greatermuch] 1. In other words, a sufficiently strong incoherent pump can negate the linear damping and yield a response bounded only by two-body loss. This heightened sensitivity persists until the oscillator switches over to the classical limit at Ω glyph[greaterorsimilar] ( γ + 1 γ 2 ) 1 / 2 , as shown in Fig. 5. Such a behavior is reminiscent of sensing protocols based on switching in optically bistable systems [43].

Apart from the 'signal' 〈 ˆ , the efficiency of sensors is a 〉 limited by the signal-to-noise ratio (SNR) [10, 44]. Estimating the noise from the spread of the Wigner function, we find one can have G &gt; 0 1 and SNR glyph[greaterorsimilar] 1 in the quantum domain, 〈 ˆ a 〉 ∼ 1 (see Supplement [28]). Further, one can measure weak signals in fewer shots in the limit-cycle regime than in the classical limit of the critical oscillator. This is because particle-number fluctuations diverge in this limit, causing the SNR to vanish.

The quantum vdP model can be probed in several experimental setups. In particular, as proposed in [21, 22], one can engineer the dissipation via sideband transitions which either add or remove energy quanta. For instance, one can laser excite a harmonically trapped ion to its blue or red motional sidebands to implement one-phonon gain or two-phonon loss, as illustrated in Fig. 6. In such a setup, the environment relaxation time is of order ∆ E -1 , where ∆ E is the level spacing between | g 〉 and | e 〉 , which dominates all other energy scales [30]. Hence, the Markov approximation holds. Additionally, one must have resolved sidebands and suppress off-resonant excitations.

FIG. 6. Experimental scheme for realizing a quantum vdP oscillator with a trapped ion, adapted from Ref. [21]. Motional modes | n 〉 act as the energy levels of the oscillator, while | g 〉 , | e 〉 , and | e ′ 〉 represent internal states. Sideband transitions, followed by spontaneous decay, effectively inject or remove energy quanta (phonons), realizing gain and loss processes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-70779e80-c275-4455-b7cb-e05cb05e2620.jpg)

As shown in Ref. [21], one can satisfy these constraints for several tens of low-energy modes. The response to a resonant drive can be measured with well-established techniques, including tomography [45] and more direct mapping of the Wigner function [46-50].

Similar techniques can be used for phonon modes of an optomechanical membrane [22] or photon modes in microwave resonators [51]. One can also realize an incoherent pump through spontaneous emission from a set of inverted qubits [52, 53]. Strong two-body loss can result from Josephson junctions in superconducting circuits [42] or polariton blockade in an optical cavity [54].

Finally, note that Eq. (3) is not the only quantum model that reproduces the classical limit of the vdP oscillator. In particular, one can harness energy-dependent one-body loss instead of two-body loss, as in Refs. [5558]. Numerics show that the key features are unaffected. There may also exist alternative descriptions of the vdP oscillator in terms of nonstandard Hamiltonians with possible quantization [59-61].

In summary, we have characterized the response of an open quantum system across a dynamical critical point by modeling a prototypical self-sustained oscillator. Like its classical counterpart, we find an increased sensitivity to external drives, which is ideal for sensing applications. Genuine quantum features arise at very weak drives, including a power-law diverging linear susceptibility and a concomitant negative susceptibility. To further explore the prospects of sensing, we compare the response with that of a passive system, where the sensitivity is limited by one-body loss. We show this linear damping is negated by a strong incoherent pump, yielding a response only bounded by weak nonlinear effects, e.g., two-body decay. Our findings are largely generic and can be probed in present-day experimental platforms, contributing to a broader understanding of dynamical criticality [9] and quantum-to-classical crossover [62, 63].

We thank Andreas Nunnenkamp and Chris Parmee for helpful discussions. This work was supported by EPSRC Grant EP/P009565/1 and a Simons Investigator Award.

- [1] M. Kardar, Statistical Physics of Fields (Cambridge University Press, Cambridge, UK, 2007).
- [2] S. H. Strogatz, Nonlinear Dynamics and Chaos (Westview Press, Cambridge, MA, 1994).
- [3] J. E. Marsden and M. McCracken, The Hopf Bifurcation and its Applications (Springer-Verlag, New York, 1976).
- [4] B. van der Pol, Jun., 'LXXXVIII. On 'relaxationoscillations',' Philos. Mag. 2 , 978 (1926).
- [5] T. A. J. Duke and F. J¨licher, 'Critical Oscillators as Acu tive Elements in Hearing,' in Active Processes and Otoacoustic Emissions in Hearing (Springer, 2008) p. 63.
- [6] V. M. Egu´ ıluz, M. Ospeck, Y. Choe, A. J. Hudspeth, and M. O. Magnasco, 'Essential Nonlinearities in Hearing,' Phys. Rev. Lett. 84 , 5232 (2000).
- [7] T. Mora and W. Bialek, 'Are Biological Systems Poised at Criticality?' J. Stat. Phys. 144 , 268 (2011).
- [8] M. A. Mun˜ oz, 'Colloquium: Criticality and dynamical scaling in living systems,' Rev. Mod. Phys. 90 , 031001 (2018).
- [9] A. Roli, M. Villani, A. Filisetti, and R. Serra, 'Dynamical Criticality: Overview and Open Questions,' J. Syst. Sci. Complex. 31 , 647 (2018).
- [10] C. L. Degen, F. Reinhard, and P. Cappellaro, 'Quantum sensing,' Rev. Mod. Phys. 89 , 035002 (2017).
- [11] D. Braun, G. Adesso, F. Benatti, R. Floreanini, U. Marzolino, M. W. Mitchell, and S. Pirandola, 'Quantumenhanced measurements without entanglement,' Rev. Mod. Phys. 90 , 035006 (2018).
- [12] L. M. Sieberer, M. Buchhold, and S. Diehl, 'Keldysh field theory for driven open quantum systems,' Rep. Prog. Phys. 79 , 096001 (2016).
- [13] I. Rotter and J. P. Bird, 'A review of progress in the physics of open quantum systems: theory and experiment,' Rep. Prog. Phys. 78 , 114001 (2015).
- [14] M. M¨ uller, S. Diehl, G. Pupillo, and P. Zoller, 'Engineered open systems and quantum simulations with atoms and ions,' Adv. At. Mol. Opt. Phys. 61 , 1 (2012).
- [15] A. A. Houck, H. E. T¨ ureci, and J. Koch, 'On-chip quantum simulation with superconducting circuits,' Nat. Phys. 8 , 292 (2012).
- [16] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, 'Cavity optomechanics,' Rev. Mod. Phys. 86 , 1391 (2014).
- [17] C. Noh and D. G. Angelakis, 'Quantum simulations and many-body physics with light,' Rep. Prog. Phys. 80 , 016401 (2017).
- [18] I. M. Georgescu, S. Ashhab, and F. Nori, 'Quantum simulation,' Rev. Mod. Phys. 86 , 153 (2014).
- [19] K. Vahala, M. Herrmann, S. Kn¨ unz, V. Batteiger, G. Saathoff, T. W. H¨nsch, a and Th. Udem, 'A phonon laser,' Nat. Phys. 5 , 682 (2009).
- [20] M. Sargent III, M. O. Scully, and W. E. Lamb Jr., Laser Physics (Addison-Wesley, Reading, MA, 1974).
- [21] T. E. Lee and H. R. Sadeghpour, 'Quantum synchronization of quantum van der Pol oscillators with trapped ions,' Phys. Rev. Lett. 111 , 234101 (2013).
- [22] S. Walter, A. Nunnenkamp, and C. Bruder, 'Quan-

tum synchronization of a driven self-sustained oscillator,' Phys. Rev. Lett. 112 , 094102 (2014).

- [23] N. L¨ orch, E. Amitai, A. Nunnenkamp, and C. Bruder, 'Genuine quantum signatures in synchronization of anharmonic self-oscillators,' Phys. Rev. Lett. 117 , 073601 (2016).
- [24] T. E. Lee, C.-K. Chan, and S. Wang, 'Entanglement tongue and quantum synchronization of disordered oscillators,' Phys. Rev. E 89 , 022913 (2014).
- [25] S. Walter, A. Nunnenkamp, and C. Bruder, 'Quantum synchronization of two Van der Pol oscillators,' Ann. Phys. 527 , 131 (2015).
- [26] K. Ishibashi and R. Kanamoto, 'Oscillation collapse in coupled quantum van der Pol oscillators,' Phys. Rev. E 96 , 052210 (2017).
- [27] O. Scarlatella, A. A. Clerk, and M. Schiro, 'Spectral functions and negative density of states of a drivendissipative nonlinear quantum resonator,' New J. Phys. 21 , 043040 (2019).
- [28] See Supplemental Material for exact classical solution, analytic results for the quantum response at weak drives, mapping between density matrix and the Wigner function, and signal-to-noise ratio estimates.
- [29] R.R. Puri and S.V. Lawande, 'On exact master equation for an open system,' Phys. Lett. A 62 , 143 (1977).
- [30] A. J. Daley, 'Quantum trajectories and open many-body quantum systems,' Adv. Phys. 63 , 77 (2014).
- [31] ´ . Rivas, A. D. K. Plato, S. F. Huelga, A and M. B. Plenio, 'Markovian master equations: a critical study,' New J. Phys. 12 , 113032 (2010).
- [32] ´ . A Rivas, S. F. Huelga, and M. B. Plenio, 'Quantum non-Markovianity: characterization, quantification and detection,' Rep. Prog. Phys. 77 , 094001 (2014).
- [33] H.-P. Breuer and F. Petruccione, The theory of open quantum systems (Oxford University Press, Oxford, UK, 2002).
- [34] C. Gardiner and P. Zoller, Quantum noise (Springer, New York, 2004).
- [35] S. Sachdev, 'Atom in a damped cavity,' Phys. Rev. A 29 , 2627 (1984).
- [36] H. D. Simaan and R. Loudon, 'Quantum statistics of single-beam two-photon absorption,' J. Phys. A 8 , 539 (1975).
- [37] V. V. Dodonov and S. S. Mizrahi, 'Exact stationary photon distributions due to competition between one- and two-photon absorption and emission,' J. Phys. A 30 , 5657 (1997).
- [38] K. V. Kheruntsyan, 'Wigner function for a driven anharmonic oscillator,' J. Opt. B 1 , 225 (1999).
- [39] A. C. Y. Li, F. Petruccione, and J. Koch, 'Perturbative approach to Markovian open quantum systems,' Sci. Rep. 4 , 4887 (2014).
- [40] V. V. Albert, B. Bradlyn, M. Fraas, and L. Jiang, 'Geometry and response of Lindbladians,' Phys. Rev. X 6 , 041031 (2016).
- [41] M. Konopik and E. Lutz, 'Quantum response theory for nonequilibrium steady states,' Phys. Rev. Research 1 , 033156 (2019).
- [42] Z. Leghtas et al., 'Confining the state of light to a quan-

tum manifold by engineered two-photon loss,' Science 347 , 853 (2015).

- [43] S. Aldana, C. Bruder, and A. Nunnenkamp, 'Detection of weak forces based on noise-activated switching in bistable optomechanical systems,' Phys. Rev. A 90 , 063810 (2014).
- [44] H.-K. Lau and A. A. Clerk, 'Fundamental limits and non-reciprocal approaches in non-Hermitian quantum sensing,' Nat. Commun. 9 , 4320 (2018).
- [45] A. I. Lvovsky and M. G. Raymer, 'Continuous-variable optical quantum-state tomography,' Rev. Mod. Phys. 81 , 299 (2009).
- [46] K. Banaszek and K. W´ odkiewicz, 'Direct probing of quantum phase space by photon counting,' Phys. Rev. Lett. 76 , 4344 (1996).
- [47] K. Banaszek, C. Radzewicz, K. W´dkiewicz, o and J. S. Krasi´ski, 'Direct measurement of the Wigner function n by photon counting,' Phys. Rev. A 60 , 674 (1999).
- [48] D. Leibfried, D. M. Meekhof, B. E. King, C. Monroe, W. M. Itano, and D. J. Wineland, 'Experimental determination of the motional quantum state of a trapped atom,' Phys. Rev. Lett. 77 , 4281 (1996).
- [49] L. G. Lutterbach and L. Davidovich, 'Method for direct measurement of the Wigner function in cavity QED and ion traps,' Phys. Rev. Lett. 78 , 2547 (1997).
- [50] P. Bertet, A. Auffeves, P. Maioli, S. Osnaghi, T. Meunier, M. Brune, J. M. Raimond, and S. Haroche, 'Direct measurement of the Wigner function of a one-photon Fock state in a cavity,' Phys. Rev. Lett. 89 , 200402 (2002).
- [51] J. R. Johansson, G. Johansson, and F. Nori, 'Optomechanical-like coupling between superconducting resonators,' Phys. Rev. A 90 , 053833 (2014).
- [52] J. Lebreuilly, M. Wouters, and I. Carusotto, 'Towards strongly correlated photons in arrays of dissipative nonlinear cavities under a frequency-dependent incoherent pumping,' C. R. Phys. 17 , 836 (2016).
- [53] R. Ma, C. Owens, A. Houck, D. I. Schuster, and J. Si-
- mon, 'Autonomous stabilizer for incompressible photon fluids and solids,' Phys. Rev. A 95 , 043811 (2017).
- [54] N. Jia, N. Schine, A. Georgakopoulos, A. Ryou, L. W. Clark, A. Sommer, and J. Simon, 'A strongly interacting polaritonic quantum dot,' Nat. Phys. 14 , 550 (2018).
- [55] N. L¨rch, S. E. Nigg, A. Nunnenkamp, R. P. Tiwari, o and C. Bruder, 'Quantum synchronization blockade: Energy quantization hinders synchronization of identical oscillators,' Phys. Rev. Lett. 118 , 243602 (2017).
- [56] S. E. Nigg, 'Observing quantum synchronization blockade in circuit quantum electrodynamics,' Phys. Rev. A 97 , 013811 (2018).
- [57] J. J. Viennot, X. Ma, and K. W. Lehnert, 'Phononnumber-sensitive electromechanics,' Phys. Rev. Lett. 121 , 183601 (2018).
- [58] S. Rips, M. Kiffner, I. Wilson-Rae, and M. J. Hartmann, 'Steady-state negative Wigner functions of nonlinear nanomechanical oscillators,' New J. Phys. 14 , 023042 (2012).
- [59] V. K. Chandrasekar, M. Senthilvelan, and M. Lakshmanan, 'On the Lagrangian and Hamiltonian description of the damped linear harmonic oscillator,' J. Math. Phys. 48 , 032701 (2007).
- [60] C. M. Bender, M. Gianfreda, N. Hassanpour, and H. F. Jones, 'Comment on 'On the Lagrangian and Hamiltonian description of the damped linear harmonic oscillator' [J. Math. Phys. 48, 032701 (2007)],' J. Math. Phys. 57 , 084101 (2016).
- [61] V. Chithiika Ruby, M. Senthilvelan, and M. Lakshmanan, 'Exact quantization of a PT-symmetric (reversible) Li´nard-type nonlinear oscillator,' J. Phys. A e 45 , 382002 (2012).
- [62] L. Davidovich, 'From quantum to classical: Schrdinger cats, entanglement, and decoherence,' Phys. Scr. 91 , 063013 (2016).
- [63] M. A. Schlosshauer, 'Quantum decoherence,' Phys. Rep. 831 , 1 (2019).

## Supplemental Material for 'Critical response of a quantum van der Pol oscillator'

## CONTENTS

SI. Steady state of a classical van der Pol oscillator

SII. Response of a quantum van der Pol oscillator at weak drives

- A. Linear response for γ + 1 = γ -1 = 0
- B. Nonmonotonic response for γ ± 1 glyph[lessmuch] γ 2
- C. Linear response for γ + 1 = γ -1 glyph[greatermuch] γ 2
- D. Linear response for γ + 1 -γ -1 glyph[greatermuch] γ 2

SIII. Mapping between density matrix and Wigner function

SIV. Signal-to-noise ratio

References

1

1

2

3

4

6

7

8

8

## SI. STEADY STATE OF A CLASSICAL VAN DER POL OSCILLATOR

Here we derive a closed-form expression for the steady-state response of a classical van der Pol (vdP) oscillator. As described in the main text, the equation of motion for such an oscillator is given by

$$\dot { \alpha } = \gamma _ { 1 } \alpha - \gamma _ { 2 } | \alpha | ^ { 2 } \alpha + \Omega$$

where α is the complex amplitude, γ 1 is the negative damping, γ 2 is the nonlinear damping, and Ω is the drive. The phase of the drive is chosen such that Ω is real and positive. Hence, α is real in steady state. To solve for the response, we introduce the rescaled quantities ˜ γ 1 ≡ γ /γ 1 2 , ˜ Ω ≡ Ω /γ 2 , and ˜ α ≡ α/ ˜ Ω 1 / 3 . In steady state,

$$( \tilde { \gamma } _ { 1 } / \tilde { \Omega } ^ { 2 / 3 } ) \tilde { \alpha } - \tilde { \alpha } ^ { 3 } + 1 = 0 \.$$

The only stable solution to Eq. (S2) is given by ˜ = α f (˜ γ / 1 ˜ Ω 2 / 3 ), where

$$f ( x ) \equiv \left ( \frac { 1 } { 2 } + \sqrt { \frac { 1 } { 4 } - \frac { x ^ { 3 } } { 2 7 } } \, \right ) ^ { 1 / 3 } + \text{sgn} ( x ) ^ { 2 / 3 } \left ( \frac { 1 } { 2 } - \sqrt { \frac { 1 } { 4 } - \frac { x ^ { 3 } } { 2 7 } } \, \right ) ^ { 1 / 3 }.$$

Here, sgn denotes the sign function and the powers are calculated using principal values of the arguments. Figure S1 shows the variation of f ( x ). Note that f (0) = 1, i.e., α = (Ω /γ 2 ) 1 / 3 for γ 1 = 0. This is the characteristic nonlinear response utilized in biological sensors [S1]. Conversely, for large x , f ( x ) can be approximated as

$$f ( x ) = \begin{cases} \sqrt { x } + 1 / ( 2 x ) + O ( x ^ { - 5 / 2 } ) & \text{for } x \gg 1 \\ 1 / | x | + O ( x ^ { - 4 } ) & \text{for } x \ll - 1 \,. \end{cases},$$

The former limit is realized for γ 1 &gt; 0 and Ω /γ 1 glyph[lessmuch] √ γ /γ 1 2 , where the steady state corresponds to a perturbed limit cycle, α = √ γ /γ 1 2 +Ω (2 / γ 1 ) + O [ (Ω /γ 1 ) 2 √ γ /γ 2 1 ] . The latter limit is realized for γ 1 &lt; 0 and Ω / γ | 1 | glyph[lessmuch] √ | γ 1 | /γ 2 . Here the system is fully damped and the response grows linearly with drive, α = Ω / γ | 1 | + O [(Ω /γ 1 ) 3 ( γ /γ 2 1 )].

## SII. RESPONSE OF A QUANTUM VAN DER POL OSCILLATOR AT WEAK DRIVES

In this section, we find analytic expressions describing the response of a weakly driven quantum vdP oscillator in special cases of interest. We start from the full dynamics in rotating frame (see main text),

$$\dot { \hat { \rho } } = [ \Omega ( \hat { a } ^ { \dagger } - \hat { a } ), \hat { \rho } ] + \gamma _ { 1 } ^ { + } \mathcal { D } [ \hat { a } ^ { \dagger } ] \hat { \rho } + \gamma _ { 1 } ^ { - } \mathcal { D } [ \hat { a } ] \hat { \rho } + \gamma _ { 2 } \mathcal { D } [ \hat { a } ^ { 2 } ] \hat { \rho } \,,$$

FIG. S1. Classical steady-state response as a function of rescaled parameters. The function f ( x ) is defined in Eq. (S3).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-ee5e9347-54dd-4a52-9620-5fed2e65b17d.jpg)

where ˆ is the density matrix, ˆ is the ladder operator, Ω is the drive amplitude, ρ a γ + 1 is the rate of one-particle gain, γ -1 is the rate of one-particle decay, γ 2 is the rate of two-particle decay, and D [ˆ]ˆ x ρ ≡ ˆˆˆ xρx † -{ ˆ ˆ x x, ρ † ˆ } / 2. In the Fock basis {| n , n 〉 = 0 1 2 , , , . . . } , Eq. (S5) can be written as

$$\cdots \quad \cdots \quad \cdots \quad \cdots \quad \cdots \\ \hat { \rho } _ { n, n ^ { \prime } } & = \underline { \in } \left ( \sqrt { n } \rho _ { n - 1, n ^ { \prime } } - \sqrt { n + 1 } \rho _ { n + 1, n ^ { \prime } } + \sqrt { n ^ { \prime } } \rho _ { n, n ^ { \prime } - 1 } - \sqrt { n ^ { \prime } + 1 } \rho _ { n, n ^ { \prime } + 1 } \right ) \\ & \quad + \gamma _ { 1 } ^ { + } \left [ \sqrt { n n ^ { \prime } } \rho _ { n - 1, n ^ { \prime } - 1 } - \left ( \frac { n + n ^ { \prime } } { 2 } + 1 \right ) \rho _ { n, n ^ { \prime } } \right ] \\ & \quad + \gamma _ { 1 } ^ { - } \left [ \sqrt { ( n + 1 ) ( n ^ { \prime } + 1 ) } \, \rho _ { n + 1, n ^ { \prime } + 1 } - \frac { n + n ^ { \prime } } { 2 } \, \rho _ { n, n ^ { \prime } } \right ] \\ & \quad + \gamma _ { 2 } \left [ \sqrt { ( n + 2 ) ( n + 1 ) ( n ^ { \prime } + 2 ) ( n ^ { \prime } + 1 ) } \, \rho _ { n + 2, n ^ { \prime } + 2 } - \frac { n ( n - 1 ) + n ^ { \prime } ( n ^ { \prime } - 1 ) } { 2 } \, \rho _ { n, n ^ { \prime } } \right ] \,.$$

Note that Ω couples neighboring elements of the density matrix, whereas the dissipative terms only couple elements in the same diagonal, i.e., with n -n ′ = constant. The response is given by 〈 ˆ a 〉 = Tr(ˆˆ) = aρ ∑ n √ nρ n,n -1 , where ρ is the steady-state density matrix. Since all coefficients in Eq. (S6) are real, ρ will be real and symmetric. Below we consider this steady-state response in various limits.

## A. Linear response for γ + 1 = γ -1 = 0

For γ ± 1 = 0 and weak drives Ω glyph[lessmuch] γ 2 , any two-particle excitation decays very rapidly. Hence, the dynamics are well approximated by retaining the lowest three oscillator levels, n = 0 1 2, in Eq. (S6), which yields , ,

$$\dot { \rho } _ { 0 0 } = 2 \gamma _ { 2 } \rho _ { 2 2 } - 2 \Omega \rho _ { 1 0 }$$

$$\rho _ { 0 0 } & = \epsilon _ { 1 2 } \rho _ { 2 2 } \overset { \sim } { \rho } \overset { \sim } { \rho } _ { 1 0 } \,, \\ \dot { \rho } _ { 1 1 } & = 2 \Omega ( \rho _ { 1 0 } - \sqrt { 2 } \rho _ { 2 1 } )$$

$$\dot { \rho } _ { 1 0 } = \Omega ( \rho _ { 0 0 } - \rho _ { 1 1 } )$$

√

$$\rho _ { 1 0 } & = \dots \rho _ { 0 0 } \quad \rho _ { 1 1 } \,, & \quad \langle \omega _ { 1 } \, \varkappa _ { 1 } \\ \dot { \rho } _ { 2 1 } & = - \gamma _ { 2 } \rho _ { 2 1 } + \Omega ( \sqrt { 2 } \rho _ { 1 1 } + \rho _ { 2 0 } - \sqrt { 2 } \rho _ { 2 2 } ) \, & \quad ( S 7 d )$$

$$\dot { \rho } _ { 2 0 } = - \gamma _ { 2 } \rho _ { 2 0 } + \Omega ( \sqrt { 2 } \rho _ { 1 0 } - \rho _ { 2 1 } )$$

where ρ 22 = 1 -ρ 00 -ρ 11 . Substituting Ω → 0 into these equations, we find the undriven steady state is given by ρ 22 = ρ 21 = ρ 20 = 0 and ρ 00 = ρ 11 = 1 2. / To linear order in Ω /γ 2 , Eqs. (S7b) and (S7d) yield ρ 10 = √ 2 ρ 21 = Ω /γ 2 . Thus, the linear response is given by 〈 ˆ a 〉 = ρ 10 + √ 2 ρ 21 = 2Ω /γ 2 , as evident in Fig. 1 of the main text.

## B. Nonmonotonic response for γ ± 1 glyph[lessmuch] γ 2

$$1. \ \gamma _ { 1 } ^ { + } = 0, \, 0 < \gamma _ { 1 } ^ { - } \ll \gamma _ { 2 }$$

Two-level approximation.As the system is fully damped, we can project the dynamics at weak drives onto the manifold spanned by n = 0 and 1. Then Eq. (S6) reduces to

˙

ρ

11

= 2Ω

ρ

10

-

γ

-

1

ρ

11

,

(S8a)

$$\dot { \rho } _ { 1 0 } = \Omega ( \rho _ { 0 0 } - \rho _ { 1 1 } ) - ( \gamma _ { 1 } ^ { - } / 2 ) \rho _ { 1 0 }$$

where ρ 00 = 1 -ρ 11 . In steady state, we find

$$\langle \hat { a } \rangle = \rho _ { 1 0 } = \frac { 2 \Omega \gamma _ { 1 } ^ { - } } { ( \gamma _ { 1 } ^ { - } ) ^ { 2 } + 8 \Omega ^ { 2 } } \.$$

Thus, the response grows linearly as 〈 ˆ a 〉 ≈ 2Ω /γ -1 for Ω glyph[lessorsimilar] γ -1 , but falls off as 〈 ˆ a 〉 ≈ γ -1 / (4Ω) for Ω glyph[greaterorsimilar] γ -1 . This is shown by the dashed black curve in Fig. S2. We see the linear susceptibility diverges as 2 /γ -1 for γ -1 → 0.

FIG. S2. Response 〈 ˆ a 〉 as a function of drive Ω for γ + 1 = 0 and γ -1 /γ 2 = 0 02. The perturbative result corresponds to Eq. (S17). .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-056035d0-4428-46c7-96cf-da978e31db90.jpg)

Three-level approximation.The two-level picture breaks down at large drives. As higher-energy modes come into play, the response grows again with Ω, as seen in Fig. S2. This revival can be captured by including the n = 2 level in the dynamics, yielding, with ρ 11 = 1 -ρ 00 -ρ 22 ,

$$\dot { \rho } _ { 0 0 } = 2 \gamma _ { 2 } \rho _ { 2 2 } - 2 \Omega \rho _ { 1 0 } + \gamma _ { 1 } ^ { - } \rho _ { 1 1 }$$

√

$$\rho _ { 0 0 } & = \lambda _ { 1 2 } \rho _ { 2 2 } \, \quad \lambda _ { 2 3 } \rho _ { 1 0 } \, + \, \lambda _ { 1 } \, \rho _ { 1 1 } \,, & \quad \langle \nu \nu \nu \nu \nu \\ \dot { \rho } _ { 1 0 } & = \Omega ( \rho _ { 0 0 } - \rho _ { 1 1 } ) + \gamma _ { 1 } ^ { - } ( \sqrt { 2 } \rho _ { 2 1 } - \rho _ { 1 0 } / 2 ) \,, & \quad \langle \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu$$

$$\dot { \rho } _ { 2 2 } = - 2 ( \gamma _ { 2 } + \gamma _ { 1 } ^ { - } ) \rho _ { 2 2 } + 2 \sqrt { 2 } \Omega \rho _ { 2 1 }$$

$$\Gamma _ { 2 1 } & = - \left \{ \sum _ { L } \alpha \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L l } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { R } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { \Lambda } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { l } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ \left \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \left \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L. } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ {L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ {. } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ \null \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L } \right \} _ { L. }$$

$$P ^ { \prime \prime } _ { 2 2 } & - \quad \langle \, \langle \, \rangle ^ { 2 } + \, \langle \, \rangle ^ { 1 } \, \rangle P ^ { \prime \prime } _ { 2 2 } \, \langle \, \rangle \, \langle \, \rangle ^ { \prime \prime } _ { 2 1 } \,, \\ \dot { \rho } _ { 2 1 } & = - ( \gamma _ { 2 } + 3 \gamma _ { 1 } ^ { - } / 2 ) \rho _ { 2 1 } + \Omega \left [ \sqrt { 2 } ( \rho _ { 1 1 } - \rho _ { 2 2 } ) + \rho _ { 2 0 } \right ] \,,$$

Let us solve for the steady state to lowest order in γ -1 /γ 2 and Ω /γ 2 . From the last three equations, we can write

$$\rho _ { 2 2 } \approx ( \sqrt { 2 } \Omega / \gamma _ { 2 } ) \rho _ { 2 1 } \,,$$

$$& \rho _ { 2 2 } \approx \langle \vee \alpha \lambda / \eta _ { 2 } \rangle \rho _ { 2 1 } \,, & \langle \vartheta _ { 1 1 } \alpha \rangle \\ & \rho _ { 2 1 } \approx ( \Omega / \gamma _ { 2 } ) \left [ \sqrt { 2 } ( \rho _ { 1 1 } - \rho _ { 2 2 } ) + \rho _ { 2 0 } \right ] \,, & \langle S 1 1 b ) \\ & \Lambda _ { \dots } \approx \Theta / \Lambda _ { \dots } \Lambda / \Lambda _ { \dots } - \Lambda _ { \dots } \, \Lambda & \Lambda _ { \dots } \, \Lambda$$

$$\rho _ { 2 0 } \approx ( \Omega / \gamma _ { 2 } ) ( \sqrt { 2 } \rho _ { 1 0 } - \rho _ { 2 1 } ) \,.$$

Substituting ρ 22 and ρ 20 into Eq. (S11b), we find

$$\mu \cdot \mu \cdot \mu \cdot \\ \rho _ { 2 1 } \approx ( \sqrt { 2 } \Omega / \gamma _ { 2 } ) [ \rho _ { 1 1 } + ( \Omega / \gamma _ { 2 } ) \rho _ { 1 0 } ] \,.$$

Using Eqs. (S11a) and (S12) in Eq. (S10a) yields

$$\rho _ { 1 0 } \approx \left ( \frac { 2 \Omega } { \gamma _ { 2 } } + \frac { \gamma _ { 1 } ^ { - } } { 2 \Omega } \right ) \rho _ { 1 1 } \,.$$

Substituting this expression back into Eqs. (S12) and (S11a) yields

$$\begin{array} {$$

Finally, substituting Eqs. (S13)-(S15) into Eq. (S10b) and using ρ 00 = 1 -ρ 11 -ρ 22 , we find

$$\rho _ { 1 1 } \approx \frac { 4 \Omega ^ { 2 }$$

which matches the two-level result. However, the response is modified as

$$\langle \hat { a } \rangle = \rho _ { 1 0 } + \$$

where we have used Eqs. (S13), (S14), and (S16). It reduces to the two-level expression in Eq. (S9) for Ω glyph[lessmuch] ( γ -1 γ 2 ) 1 / 2 . For Ω ∼ ( γ -1 γ 2 ) 1 / 2 glyph[greatermuch] γ 1 , we find 〈 ˆ a 〉 ≈ γ -1 / (4Ω) + 2Ω /γ 2 , which attains a minimum at Ω = ( γ -1 γ 2 ) 1 / 2 / (2 √ 2) and increases linearly with drive at large Ω. This nonmonotonic variation is evident in Fig. S2.

$$2. \ 0 \leq \gamma _ { 1 } ^ { + }, \gamma$$

In the presence of one-particle gain ( γ + 1 &gt; 0), the undriven steady state corresponds to a dynamic equilibrium where particles flow in and out of the oscillator. When γ 2 glyph[greatermuch] γ ± 1 , the dynamics are confined to the levels n = 0, 1, and 2 for weak drives. We consider this three-level system, governed by the equations of motion [from Eq. (S6)]

$$\dot { \rho } _ { 0 0 } = 2 \gamma _ { 2 } \$$

$$\rho _ { 1 0 } & = \ast \wp _ { 0 0 } \quad \rho _ { 1 1 / 2 } + \, \ r _ { 1 1 } \ \vee \, \rho _ { 2 1 } \quad \rho _ { 1 0 / 2 } \, \rho _ { 1 1 / 2 } \, \rho _ { 1 0 } \,, & \quad \, \langle \omega _ { 1 1 / 2 } \, \rangle, \\ \dot { \rho } _ { 2 2 } & = - 2 ( \gamma _ { 2 } + \gamma _ { 1 } ^ { - } ) \rho _ { 2 2 } + 2 \sqrt { 2 } \Omega \rho _ { 2 1 } + 2 \gamma _ { 1 } ^ { + } \rho _ { 1 1 } \,, & \quad \, \langle \omega _ { 2 2 } \, \rangle,$$

$$P _ { 0 0 } & = \Gamma _ { 1 2 } P Z Z \quad \Delta \Delta \Delta \Gamma _ { 1 0 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 0 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 0 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 0 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 0 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \\ \dot { \rho } _ { 1 0 } & = \Omega ( \rho _ { 0 0 } - \rho _ { 1 1 } ) + \gamma _ { 1 } ^ { - } ( \sqrt { 2 } \rho _ { 2 1 } - \rho _ { 1 0 } / 2 ) - ( 3 \gamma _ { 1 } ^ { + } / 2 ) \rho _ { 1 0 } \ \Gamma _ { 1 0 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 } \ \Gamma _ { 1 1 }$$

$$P ^ { \prime \prime } _ { \mathbb { Z } } & = - \varkappa / 2 + 1 / 1 / P ^ { \prime \prime } _ { \mathbb { Z } } + \varkappa \nu \varepsilon \varepsilon ^ { \prime \prime } _ { \mathbb { Z } } + \varkappa / 1 / P ^ { \prime } _ { \mathbb { Z } } 1 + \varkappa / 1 / P ^ { \prime } _ { \mathbb { Z } } 1 \,, & \quad \langle \cup \varkappa \nu \nu \\ \dot { \rho } _ { 2 1 } & = - ( \gamma _ { 2 } + 3 \gamma _ { 1 } ^ { - } / 2 ) \rho _ { 2 1 } + \Omega \left [ \sqrt { 2 } ( \rho _ { 1 1 } - \rho _ { 2 2 } ) + \rho _ { 2 0 } \right ] + \sqrt { 2 } \gamma _ { 1 } ^ { + } \rho _ { 1 0 } \,, & \quad ( S 1 8 d )$$

$$\Gamma _ { 2 1 } & = - \left \{ \Gamma _ { 2 } + \vartheta _ { 1 } \right \} \vartheta _ { 1 } \right \} \Gamma _ { 2 1 } + \vartheta _ { 1 } \left \{ \Gamma _ { 2 } + \vartheta _ { 1 } \right \} \Gamma _ { 2 2 } \right \} + \Gamma _ { 2 0 } \right \} + \vartheta _ { 1 } \left \{ \Gamma _ { 2 0 } \right \}, & \left \{ \chi \chi \chi \right \} \\ \dot { \rho } _ { 2 0 } & = - ( \gamma _ { 2 } + \gamma _ { 1 } ^ { - } ) \rho _ { 2 0 } + \Omega ( \sqrt { 2 } \rho _ { 1 0 } - \rho _ { 2 1 } ) \,, &$$

where ρ 11 = 1 -ρ 00 -ρ 22 . For Ω = 0, we find the steady state ρ 11 ≈ 1 -ρ 00 ≈ γ + 1 / (3 γ + 1 + γ -1 ). Thus, negative damping leads to a finite occupation of the n = 1 level. For Ω &gt; 0, the steady-state analysis can be carried out following the same steps as in the last section, which yields, to lowest order in γ + 1 /γ 2 , γ -1 /γ 2 , and Ω /γ 2 ,

$$\langle \hat { a } \rangle \approx \frac { 2 \Omega } { \gamma _ { 2 } } \left [ \frac { ( \gamma _ { 1 } ^ { + } + \gamma _ { 1 } ^ { - } ) \gamma _ { 2 } + 8 \Omega ^ { 2 } } { ( 3 \gamma _ { 1 } ^ { + } + \gamma _ { 1 } ^ { - } ) ^ { 2 } + 8 \Omega ^ { 2 } } \right ].$$

As before, we find a nonmonotonic response. For Ω glyph[lessmuch] 3 γ + 1 + γ -1 , the response is linear, 〈 ˆ a 〉 ≈ 2Ω( γ + 1 + γ -1 ) / (3 γ + 1 + γ -1 ) 2 , with a slope that diverges as γ ± 1 → 0. For 3 γ + 1 + γ -1 glyph[lessmuch] Ω glyph[lessmuch] [( γ + 1 + γ -1 ) γ 2 ] 1 / 2 , it falls off as 〈 ˆ a 〉 ≈ ( γ + 1 + γ -1 ) / (4Ω), exhibiting negative susceptibility. For Ω glyph[greatermuch] [( γ + 1 + γ -1 ) γ 2 ] 1 / 2 , it rises again as 〈 ˆ a 〉 ≈ 2Ω /γ 2 .

## C. Linear response for γ + 1 = γ -1 glyph[greatermuch] γ 2

Here, we derive the enhanced sensitivity for the critical condition γ + 1 = γ -1 ≡ Γ 1 glyph[greatermuch] γ 2 . We first find the undriven steady state, then consider perturbation at weak drives. At Ω = 0, the steady state has no coherence. The equation of motion for the populations p n ≡ ρ n,n is found by substituting n ′ = n in Eq. (S6), which yields

$$\dot { p } _ { n } = \Gamma _ { 1 } \left [ n p _ { n - 1 } + ( n + 1 ) p _ { n + 1 } - ( 2 n + 1 ) p _ { n } \right ] + \gamma _ { 2 } \left [ ( n + 1 ) ( n + 2 ) p _ { n + 2 } - n ( n - 1 ) p _ { n } \right ] \,.$$

At γ 2 = 0, the steady state is simply given by p n = constant, i.e., the system reaches an infinite-temperature state. For γ 2 &gt; 0, this distribution is curtailed at large n and p n is a slowly-varying function set by the scale ε ≡ √ γ / 2 Γ . 1 Thus, we can approximate p n by a continuous function p x ( ) with x = nε . Then Eq. (S20) can be written as

$$x p ( x - \varepsilon ) + ( x + \varepsilon ) p ( x + \varepsilon ) - ( 2 x + \varepsilon ) p ( x ) + \varepsilon \left [ ( x + \varepsilon ) ( x + 2 \varepsilon ) p ( x + 2 \varepsilon ) - x ( x - \varepsilon ) p ( x ) \right ] = 0 \. \quad \text{($S21$)}$$

The above equation is identically satisfied up to O ε ( ). Expanding to O ε ( 2 ), we find

$$x p ^ { \prime \prime } ( x ) + ( 2 x ^ { 2 } + 1 ) p ^ { \prime } ( x ) + 4 x p ( x ) = 0$$

where prime denotes d/d x . Requiring p x ( ) to be bounded at x = 0 leads to the unique solution p x ( ) = Ce -x 2 . The constant C can be determined from probability conservation ∑ n p n = 1, or (1 /ε ) ∫ ∞ 0 p x dx ( ) = 1, which gives

$$p ( x ) & = ( 2 \varepsilon / \sqrt { \pi } ) \, e ^ { - x ^ { 2 } }.$$

Thus, the undriven steady state has a Gaussian number distribution of width √ Γ 1 /γ 2 glyph[greatermuch] 1, as shown in Fig. S3(a). To find the linear response, we consider perturbation in the leading off-diagonal terms χ n ≡ ρ n,n -1 , which satisfy

$$\dot { \chi } _ { n } = \Gamma _ { 1 } [ \sqrt { n ( n - 1 ) } \chi _ { n - 1 } + \sqrt { n ( n + 1 ) } \chi _ { n + 1 } - 2 n \chi _ { n } ] + \gamma _ { 2 } [ \sqrt { n ( n + 2 ) } ( n + 1 ) \chi _ { n + 2 } - ( n - 1 ) ^ { 2 } \chi _ { n } ] + \Omega \sqrt { n } ( p _ { n - 1 } - p _ { n } ) \.$$

$$( S 2 4 )$$

This equation is simplified with the substitution q n ≡ √ nχ n . In steady state,

$$\Gamma _ { 1 } ( q _ { n + 1 } + q _ { n - 1 } - 2 q _ { n } ) + \gamma _ { 2 } \left [ ( n + 1 ) q _ { n + 2 } - ( n + 1 / n - 2 ) q _ { n } \right ] + \Omega ( p _ { n - 1 } - p _ { n } ) = 0 \.$$

Since q n will also be a slowly-varying function as p n , we rewrite Eq. (S25) in terms of the continuous variable x = nε ,

$$q ( x + \varepsilon ) + q ( x - \varepsilon ) - 2 q ( x ) + \varepsilon \left [ ( x + \varepsilon ) q ( x + 2 \varepsilon ) - ( x - 2 \varepsilon + \varepsilon ^ { 2 } / x ) q ( x ) \right ] + \eta \varepsilon \left [ \tilde { p } ( x - \varepsilon ) - \tilde { p } ( x ) \right ] = 0 \, \quad \ \ ( S 2 6 )$$

where η ≡ Ω Γ / 1 and ˜( p x ) ≡ p x /ε ( ) ∼ O (1). Expanding Eq. (S26) to O ε ( 2 ), we find

$$q ^ { \prime \prime } ( x ) + 2 x q ^ { \prime } ( x ) + 3 q ( x ) - \eta \tilde { p } ^ { \prime } ( x ) = 0 \,.$$

Substituting ˜ ( p ′ x ) = -(4 / √ π xe ) -x 2 from Eq. (S23) and using q x ( ) ≡ (4 η/ √ π e ) -x 2 h x ( ) yields

$$h ^ { \prime \prime } ( x ) - 2 x h ^ { \prime } ( x ) + h ( x ) + x = 0 \.$$

In addition, since q (0) = 0, we have the boundary condition h (0) = 0. Then the only bounded solution to Eq. (S28) is h x ( ) = x . Therefore,

$$q ( x ) = ( 4 \eta / \sqrt { \pi } ) x e ^ { - x ^ { 2 } }$$

which agrees well with numerics, as shown in Fig. S3(a). The linear response is given by

$$\langle \hat { a } \rangle = \sum _ { n = 0 } ^ { \infty } \sqrt { n } \chi _ { n } \approx \frac { 1 } { \varepsilon } \int _ { 0 } ^ { \infty } d x q ( x ) = \frac { 2 \eta } { \sqrt { \pi \varepsilon } } = \frac { 2 \Omega } { \sqrt { \pi \Gamma _ { 1 } \gamma _ { 2 } } } \,.$$

Hence, the linear susceptibility is χ ≈ 2 / √ π Γ 1 γ 2 . Figure S3(b) shows how the exact solution interpolates between this regime and that of Γ 1 glyph[lessmuch] γ 2 , where χ ≈ 1 / (4Γ ) [see Sec. SII B 2]. 1 In contrast, a fully damped 'passive' oscillator has a susceptibility χ p = 2 /γ -1 [see Sec. SII B 1]. Thus, operating the vdP oscillator in the regime γ + 1 = γ -1 glyph[greatermuch] γ 2 yields a sensitivity enhancement χ/χ p ≈ [ γ -1 / πγ ( 2 )] 1 / 2 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-70e4b8e2-9f94-4c2c-9194-a205e22e36ed.jpg)

FIG. S3. (a) Steady-state occupations p n ≡ ρ n,n and coherences q n ≡ √ nρ n,n -1 for Γ 1 /γ 2 = 1000 and Ω → 0. The coherences are rescaled to show on the same axes. Solid lines are obtained from exact numerics and dashed lines correspond to Eqs. (S23) and (S29). (b) Linear susceptibility χ as a function of Γ 1 /γ 2 . Solid curve shows exact result and dashed curves show aymptotes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-5d388647-7ed5-4e40-8877-fd95a004f285.jpg)

## D. Linear response for γ + 1 -γ -1 glyph[greatermuch] γ 2

As we discussed in the main text, the susceptibility is enhanced much further by operating in the limit-cycle phase. Here we derive an expression for the linear response in this regime. As in the last section, we start by writing down rate equations which govern the populations of the undriven system,

$$\dot { p } _ { n } = \gamma _ { 1 } ^ { + } \left [ n p _ { n - 1 } - ( n + 1 ) p _ { n } \right ] + \gamma _ { 1 } ^ { - } \left [ ( n + 1 ) p _ { n + 1 } - n p _ { n } \right ] + \gamma _ { 2 } \left [ ( n + 1 ) ( n + 2 ) p _ { n + 2 } - n ( n - 1 ) p _ { n } \right ]. \quad ( S 3 1 )$$

The steady state is set by the parameters ε ≡ ( γ /γ 2 + 1 2 1 ) / glyph[lessmuch] 1 and ζ ≡ γ -1 /γ + 1 &lt; 1. We expect, in steady state, p n will be peaked near the classical solution, n ∗ ≈ α 2 cl = (1 -ζ ) / (2 ε 2 ). To find this distribution, we approximate p n by a continuous function p x ( ) where x = nε . Then Eq. (S31) yields

$$x p ( x - \varepsilon ) - ( x + \varepsilon ) p ( x ) + \zeta \left [ ( x + \varepsilon ) p ( x + \varepsilon ) - x p ( x ) \right ] + \varepsilon \left [ ( x + \varepsilon ) ( x + 2 \varepsilon ) p ( x + 2 \varepsilon ) - x ( x - \varepsilon ) p ( x ) \right ] = 0 \,. \quad \text{($S3$)}$$

Since p x ( ) will be peaked near x ∗ = n ε ∗ ∼ O (1 /ε ), we shift the origin with the transformation x = β/ε + , where y β ∼ O (1). We also define u y ( ) ≡ p β/ε ( + ). Then expanding Eq. (S32) in powers of y ε , we find

$$\varepsilon \beta ( 2 \beta + \zeta - 1 ) u ^ { \prime } ( y ) + \varepsilon ^ { 2 } \left [ \beta ( 4 \beta + \zeta + 1 ) u ^ { \prime \prime } ( y ) / 2 + ( 4 \beta + \zeta - 1 ) \{ y u ^ { \prime } ( y ) + u ( y ) \} \right ] + O ( \varepsilon ^ { 3 } ) = 0 \. \quad \ \ ( S 3 3 )$$

The linear term implies β = (1 -ζ ) / 2, exactly matching the classical estimate. Substituting this value for β into the quadratic term leads to the differential equation

$$( 3 - \zeta ) u ^ { \prime \prime } ( y ) + 4 y u ^ { \prime } ( y ) + 4 u ( y ) = 0$$

which has a unique positive-definite solution, u y ( ) = Ae -2 y 2 / (3 -ζ ) . The integration constant A can be determined by requiring ∑ n p n = 1, or (1 /ε ) ∫ u y dy ( ) ≈ 1, which yields

$$u ( y ) \approx \varepsilon \sqrt { 2 / [ ( 3 - \zeta ) \pi ] } \, e ^ { - 2 y ^ { 2 } / ( 3 - \zeta ) } \,.$$

From the correspondence p n ≈ u nε ( -β/ε ) and the definitions ζ ≡ γ -1 /γ + 1 , ε ≡ ( γ /γ 2 + 1 2 1 ) / , we find the undriven steady state has a Gaussian number distribution with mean ( γ + 1 -γ -1 ) / (2 γ 2 ) and standard deviation [(3 γ + 1 -γ -1 ) / (4 γ 2 )] 1 / 2 . This result was also found in Ref. [S2] and agrees well with numerics, as shown in Fig. S4(a).

Next, we consider linear perturbation at weak drives. From Eq. (S6), the coherences χ n ≡ ρ n,n -1 satisfy

$$\dot { \chi } _ { n } & = \gamma _ { 1 } ^ { + } \left [ \sqrt { n ( n - 1 ) } \chi _ { n - 1 } - ( n + 1 / 2 ) \chi _ { n } \right ] + \gamma _ { 1 } ^ { - } \left [ \sqrt { n ( n + 1 ) } \chi _ { n + 1 } - ( n - 1 / 2 ) \chi _ { n } \right ] \\ & \quad + \gamma _ { 2 } \left [ \sqrt { n ( n + 2 ) } ( n + 1 ) \chi _ { n + 2 } - ( n - 1 ) ^ { 2 } \chi _ { n } \right ] + \Omega \sqrt { n } ( p _ { n - 1 } - p _ { n } )$$

where p n are the undriven occupations. To find the steady state, we again approximate χ n by a continuous function v y ( ) where χ n ≈ v nε ( -β/ε ). Following the same steps as before, we arrive at the solution

$$v ( y ) = 4 \eta ( 3 - \zeta ) ^ { - 3 / 2 } \sqrt { ( 1 - \zeta ) / \pi } \, e ^ { - 2 y ^ { 2 } / ( 3 - \zeta ) }$$

FIG. S4. (a) Steady-state populations p n ≡ ρ n,n and coherences χ n ≡ ρ n,n -1 for γ + 1 /γ 2 = 200, γ -1 /γ 2 = 20, and Ω → 0. The coherences have been rescaled to show on the same axes as p n . Solid lines correspond to exact numerics and dashed lines show perturbative results in Eqs. (S35) and (S37). (b) Linear susceptibility χ as a function of γ + 1 /γ -1 and γ -1 /γ 2 , rescaled by the susceptibility of a passive system, χ p = 2 /γ -1 . Solid lines show exact results and dotted lines are obtained from Eq. (S39).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-98e4b6a5-e6c0-4d6b-a149-40589044dbeb.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-abaff73d-0307-4e91-9b65-714cd3f1edcc.jpg)

where η ≡ Ω /γ + 1 . Thus, χ n has the same Gaussian profile as p n , as seen in Fig. S4(a). The linear response is given by

$$\langle \hat { a } \rangle = \sum _ { n } \sqrt { n } \chi _ { n } \approx \frac { 1 } { \varepsilon ^ { 2 } } \int _ { - \beta / \varepsilon } ^ { \infty } \sqrt { \beta + \varepsilon y } \, v ( y ) d y \approx \frac { 2 \eta } { 3 \varepsilon ^ { 2 } } \left [ 1 - \frac { 2 \zeta } { 3 } + O ( \zeta ^ { 2 } ) \right ] \,.$$

Expressing η , ε , and ζ in terms of the drive and dissipation rates, we find a susceptibility

$$\chi \approx \frac { 2 } { 3 \gamma _ { 2 } } \left ( 1 - \frac { 2 \gamma _ { 1 } ^ { - } } { 3 \gamma _ { 1 } ^ { + } } \right ) \,.$$

Thus, for γ -1 glyph[lessmuch] γ + 1 , χ is only limited by two-particle decay. Consequently, the sensitivity gain over a passive system, χ/χ p , scales as γ -1 /γ 2 [Recall: χ p = 2 /γ -1 ]. This scaling is illustrated in Fig. S4(b).

## SIII. MAPPING BETWEEN DENSITY MATRIX AND WIGNER FUNCTION

In this section, we review the Wigner representation of the density operator and derive an explicit mapping between the Wigner function and density matrix for an oscillator mode. The Wigner distribution was introduced in Ref. [S3] to enable the calculation of expectation values as integrals over a phase space, similar to classical ensemble averages. The Wigner function acts as a weight for such an integral representation of the density operator. For a single oscillator mode ˆ, it can be written as the expectation of a displaced parity operator [S4] a

$$\hat { \Pi } ( \alpha ) \equiv \hat { D } ( \alpha ) ( - 1 ) ^ { \hat { a } ^ { \dagger } \hat { a } } \hat { D } ^ { - 1 } ( \alpha ) \,,$$

where α is a point in phase space and ˆ ( D α ) ≡ e αa ˆ † -α a ∗ ˆ is a displacement operator. The Wigner function is given by

$$W ( \alpha ) = ( 2 / \pi ) \text{Tr} [ \hat { \rho } \hat { \Pi } ( \alpha ) ] \.$$

Conversely, the density operator ˆ can be obtained from the Wigner function as ρ

$$\hat { \rho } = 2 \, \int d ^ { 2 } \alpha \, W ( \alpha ) \, \hat { \Pi } ( \alpha ) \.$$

So, there is a one-to-one correspondence between ˆ and ρ W . Note the operators Π( ˆ α ) are Hermitian with eigenvalues ± 1. They can also be shown to constitute a basis for expanding any operator over the complex phase space [S5]. In particular, they satisfy the orthonormality and completeness relations

$$( 4 / \pi ) \, \text{Tr} \left [ \hat { \Pi } ( \alpha ), \hat { \Pi } ( \alpha ^ { \prime } ) \right ] = \delta ( \alpha - \alpha ^ { \prime } ) \hat { \mathbf i }$$

$$( 4 / \pi ) \int d ^ { 2 } \alpha \langle m | \hat { \Pi } ( \alpha ) | m ^ { \prime } \rangle \langle n | \hat { \Pi } ( \alpha ) | n ^ { \prime } \rangle ^ { * } = \delta _ { m, n } \delta _ { m ^ { \prime }, n ^ { \prime } }$$

where | n 〉 are Fock states. Since Π( ˆ α ) are parity operators, from Eq. (S41) we see that W is real valued and uniformly bounded, W α ( ) ∈ -[ 2 /π, 2 /π ]. Further, using the property Tr[ ˆ ( Π α )] = 1 / 2 in Eq. (S42), one finds

$$\int d ^ { 2 } \alpha W ( \alpha ) = 1 \.$$

Hence, W α ( ) can be interpreted as a quasiprobability distribution in phase space. In particular, one can compute the expectation of any symmetrically ordered product { (ˆ ) a † n ˆ a m } S as an ensemble average,

$$\langle \{ ( \hat { a } ^ { \dagger } ) ^ { n } \hat { a } ^ { m } \} _ { S } \rangle = \int d ^ { 2 } \alpha ( \alpha ^ { * } ) ^ { n } \alpha ^ { m } W ( \alpha ) \,,$$

which follows from using the relation ˆ D -1 ( α aD α )ˆ ˆ ( ) = ˆ + a α ˆ in Eq. (S42). 1 The linear response 〈 ˆ a 〉 corresponds to a special case of Eq. (S44) that amounts to measuring the center of mass of the Wigner distribution,

$$\langle \hat { a } \rangle = \int d ^ { 2 } \alpha \alpha W ( \alpha ) \,.$$

In principle, Eqs. (S41) and (S42) completely specify W α ( ) in terms of ˆ and ρ vice versa . However, one can also find an explicit mapping between the Wigner function and the density matrix in the Fock basis. To see this, note the parity operator ( -1) ˆ a a † ˆ reflects ˆ to a -ˆ, ( a -1) ˆ a a † ˆ ˆ( a -1) ˆ a a † ˆ = -ˆ. a It follows we have the identity [S4]

$$( - 1 ) ^ { \hat { a } ^ { \dagger } \hat { a } } \hat { D } ( \alpha ) ( - 1 ) ^ { \hat { a } ^ { \dagger } \hat { a } } = \hat { D } ( - \alpha ) = \hat { D } ^ { - 1 } ( \alpha ) \.$$

Using the above in Eq. (S40), we can write Π( ˆ α ) = ˆ (2 D α )( -1) ˆ a a † ˆ which, when used in Eq. (S41), yields

$$W ( \alpha ) = ( 2 / \pi ) \sum _ { n, n ^ { \prime } } \rho _ { n, n ^ { \prime } } \langle n ^ { \prime } | \hat { D } ( 2 \alpha ) | n \rangle ( - 1 ) ^ { n } \.$$

The matrix elements of ˆ D can be expressed in terms of associated Laguerre polynomials L ( p ) q ( x ) [S5]. Thus, we find

$$W ( \alpha ) = ( 2 / \pi ) e ^ { - 2 | \alpha | ^ { 2 } } \sum \nolimits _ { n, n ^ { \prime } } ( - 1 ) ^ { n } \sqrt { n! / n ^ { \prime }! } \, ( 2 \alpha ) ^ { n ^ { \prime } - n } \, L _ { n } ^ { ( n ^ { \prime } - n ) } ( 4 | \alpha | ^ { 2 } ) \, \rho _ { n, n ^ { \prime } } \,.$$

The mapping can be written in a more symmetric form in polar coordinates α = re iφ ,

$$W ( r, \phi ) = \frac { 2 } { \pi } e ^ { - 2 r ^ { 2 } } \left [ \sum _ { n = 0 } ^ { \infty } ( - 1 ) ^ { n } L _ { n } ^ { ( 0 ) } ( 4 r ^ { 2 } ) \, \rho _ { n, n } + \sum _ { j = 1 } ^ { \infty } ( 2 r ) ^ { j } \sum _ { n = 0 } ^ { \infty } ( - 1 ) ^ { n } \frac { L _ { n } ^ { ( j ) } ( 4 r ^ { 2 } ) } { \sqrt { ( n + 1 ) _ { j } } } ( e ^ { - i \phi } \rho _ { n + j, n } + e ^ { i j \phi } \rho _ { n, n + j } ) \right ],$$

where ( x ) j is the Pochhammer symbol, ( x ) j = x x ( +1) . . . ( x + j -1). We see the Wigner function has no angular variation for a purely diagonal density matrix, e.g., in the steady state of an undriven vdP oscillator. Further, elements in the j -th off diagonal contribute j units of angular momentum. Equation (S49) can be inverted using orthogonality properties of Laguerre polynomials and exponential functions to give ρ n,n ′ in terms of W r,φ ( ),

$$\rho _ { n, n ^ { \prime } } = ( - 1 ) ^ { n ^ { \prime } } \sqrt { \frac { n ^ { \prime }! } { n! } } \int _ { 0 } ^ { 2 \pi } d \phi e ^ { i ( n - n ^ { \prime } ) \phi } \int _ { 0 } ^ { \infty } d r \, e ^ { - 2 r ^ { 2 } } ( 2 r ) ^ { n - n ^ { \prime } + 1 } L _ { n ^ { \prime } } ^ { ( n - n ^ { \prime } ) } ( 4 r ^ { 2 } ) \, W ( r, \phi ) \.$$

## SIV. SIGNAL-TO-NOISE RATIO

We have focused on the response 〈 ˆ a 〉 of a vdP oscillator to an external drive. An experimental measurement of this response will be subject to quantum noise. Here we estimate the signal-to-noise ratio. As we showed in Sec. SIII, the 'signal' 〈 ˆ a 〉 as well as other physical quantities can be extracted from the Wigner function W α ( ) which represents a quasiprobability distribution in phase space. From Eq. (S45), 〈 ˆ a 〉 is given by the center of mass of W α ( ). The noise σ can be estimated from the spread about this center of mass,

$$\sigma ^ { 2 } = \, \int d ^ { 2 } \alpha \, | \alpha - \langle \hat { a } \rangle | ^ { 2 } W ( \alpha ) = - | \langle \hat { a } \rangle | ^ { 2 } + \, \int d ^ { 2 } \alpha \, | \alpha | ^ { 2 } W ( \alpha )$$

where we have used Eqs. (S43) and (S45). The last integral can be evaluated from Eq. (S44), yielding

$$\sigma = \sqrt { \langle \hat { a } ^ { \dagger } \hat { a } + \hat { a } \hat { a } ^ { \dagger } \rangle / 2 - | \langle \hat { a } \rangle | ^ { 2 } } = \sqrt { \langle \hat { n } \rangle + 1 / 2 - | \langle \hat { a } \rangle | ^ { 2 } }$$

where ˆ is the number operator, ˆ n n ≡ ˆ ˆ. a a † The signal-to-noise ratio is given by SNR = 〈 ˆ a /σ 〉 . We plot the response and the SNR as a function of drive and dissipation in Fig. S5, which shows one can have SNR glyph[greaterorsimilar] 1 while being in the quantum regime, 〈 ˆ a 〉 ∼ 1. In Fig. S5(a), we consider the case γ ± 1 glyph[lessmuch] γ 2 , where one finds a nonmonotonic response with increasing drive. Here SNR &lt; 1 as the oscillator is close to the vacuum state. Nonetheless, such states have been measured in experiments with close to 1% accuracy [S6-S12]. In Fig. S5(b), we find, for γ ± 1 glyph[greatermuch] γ 2 , one can have both SNR &gt; 1 as well as a susceptibility boost over a passive oscillator, χ &gt; χ p .

These noise estimates also show the advantage of measuring weak signals with a quantum vdP oscillator as opposed to simply using the classical limit of the critical oscillator. Although the susceptibility can be made infinitely large in this limit, the particle-number fluctuations also diverge as ∆ n ∼ ( γ ± 1 /γ 2 ) 1 / 2 [see Eq. (S23)]. Therefore, the SNR, in fact, vanishes for weak signals. Instead, by harnessing a quantum oscillator in the limit-cycle regime, one can enhance the 'minimum detectable signal per unit time,' which is an important figure-of-merit for sensors [S13].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-13fc6ba5-f425-4bba-b5e8-9b6f3c046069.jpg)

FIG. S5. Response 〈 ˆ a 〉 as a function of drive Ω and dissipation rates γ ± 1 . Dashed lines show contours of constant signal-to-noise ratio. (a) γ + 1 = 0. Solid curve encloses the region over which d ˆ 〈 a / 〉 dΩ &lt; 0. (b) γ -1 /γ 2 = 30. Below the solid curve, one finds a sensitivity enhancement over a passive system, i.e., d ˆ 〈 a / 〉 dΩ &gt; 2 /γ -1 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01002v3-figure-6a398e64-f1b2-4e10-a183-f9f9729075ae.jpg)

- [S2] V. V. Dodonov and S. S. Mizrahi, 'Exact stationary photon distributions due to competition between one-and two-photon absorption and emission,' J. Phys. A 30 , 5657 (1997).
- [S3] E. Wigner, 'On the quantum correction for thermodynamic equilibrium,' Phys. Rev. 40 , 749 (1932).
- [S4] K. E. Cahill and R. J. Glauber, 'Density operators and quasiprobability distributions,' Phys. Rev. 177 , 1882 (1969).
- [S5] K. E. Cahill and R. J. Glauber, 'Ordered expansions in boson amplitude operators,' Phys. Rev. 177 , 1857 (1969).
- [S6] K. Banaszek, C. Radzewicz, K. W´dkiewicz, o and J. S. Krasi´ski, 'Direct measurement of the Wigner function by photon n counting,' Phys. Rev. A 60 , 674 (1999).
- [S7] A. I. Lvovsky, H. Hansen, T. Aichele, O. Benson, J. Mlynek, and S. Schiller, 'Quantum state reconstruction of the single-photon Fock state,' Phys. Rev. Lett. 87 , 050402 (2001).
- [S8] P. Bertet, A. Auffeves, P. Maioli, S. Osnaghi, T. Meunier, M. Brune, J. M. Raimond, and S. Haroche, 'Direct measurement of the Wigner function of a one-photon Fock state in a cavity,' Phys. Rev. Lett. 89 , 200402 (2002).
- [S9] M. Hofheinz et al., 'Synthesizing arbitrary quantum states in a superconducting resonator,' Nature (London) 459 , 546 (2009).
- [S10] Y. Chu, P. Kharel, T. Yoon, L. Frunzio, P. T. Rakich, and R. J. Schoelkopf, 'Creation and control of multi-phonon Fock states in a bulk acoustic-wave resonator,' Nature (London) 563 , 666 (2018).
- [S11] D. T. Smithey, M. Beck, M. G. Raymer, and A. Faridani, 'Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomography: Application to squeezed states and the vacuum,' Phys. Rev. Lett. 70 , 1244 (1993).
- [S12] D. Leibfried, D. M. Meekhof, B. E. King, C. Monroe, W. M. Itano, and D. J. Wineland, 'Experimental determination of the motional quantum state of a trapped atom,' Phys. Rev. Lett. 77 , 4281 (1996).
- [S13] C. L. Degen, F. Reinhard, and P. Cappellaro, 'Quantum sensing,' Rev. Mod. Phys. 89 , 035002 (2017).
Shovan Dutta Nigel R. Cooper
Critical Response of a Quantum van der Pol Oscillator
Classical dynamical systems close to a critical point are known to act as
efficient sensors due to a strongly nonlinear response. We explore such systems
in the quantum regime by modeling a quantum version of a driven van der Pol
oscillator. We find the classical response survives down to one excitation
quantum. At very weak drives, genuine quantum features arise, including
diverging and negative susceptibilities. Further, the linear response is
greatly enhanced by using a strong incoherent pump. These results are largely
generic and can be probed in current experimental platforms suited for quantum
sensing.