Critical Response of a Quantum van der Pol Oscillator

TL;DR

This paper investigates the quantum behavior of a driven van der Pol oscillator near critical points, revealing quantum-specific phenomena like negative susceptibilities and enhanced responses, with implications for quantum sensing.

Contribution

It introduces a quantum model of the van der Pol oscillator and demonstrates quantum effects in its critical response, extending classical nonlinear sensing concepts to the quantum domain.

Findings

Quantum response persists at single quantum level

Negative and diverging susceptibilities observed at weak drives

Strong incoherent pumping enhances linear response

Abstract

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.

Equations

$$\dot{\alpha }={\gamma }_1\alpha -{\gamma }_2|\alpha {|}^2\alpha +\Omega ,$$

$$\alpha =\{\begin{array}{ll}\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{array}$$

$$\dot{\hat{\rho }}=-\mathrm{i}[\hat{H},\hat{\rho }]+{\gamma }_1^{+}\mathcal{D}[{\hat{a}}^{†}]\hat{\rho }+{\gamma }_1^{-}\mathcal{D}[\hat{a}]\hat{\rho }+{\gamma }_2\mathcal{D}[{\hat{a}}^2]\hat{\rho },$$

$$⟨\dot{\hat{a}}⟩=[({\gamma }_1^{+}-{\gamma }_1^{-})/2]⟨\hat{a}⟩-{\gamma }_2⟨{\hat{a}}^{†}\hat{a}\hat{a}⟩+\Omega .$$

$$⟨\hat{a}⟩\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].$$

$$\chi\approx\textstyle\frac{2}{3}\hskip 0.85355pt\gamma_{2}^{-1}\big{[}1-2\gamma_{1}^{-}/(3\gamma_{1}^{+})\big{]}\;.$$

$$\dot{\alpha }={\gamma }_1\alpha -{\gamma }_2|\alpha {|}^2\alpha +\Omega ,$$

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

$$f(x)\equiv\bigg{(}\frac{1}{2}+\sqrt{\frac{1}{4}-\frac{x^{3}}{27}}\hskip 2.27626pt\bigg{)}^{1/3}+\hskip 1.42271pt\text{sgn}(x)^{2/3}\bigg{(}\frac{1}{2}-\sqrt{\frac{1}{4}-\frac{x^{3}}{27}}\hskip 2.27626pt\bigg{)}^{1/3}.$$

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

$$\dot{\hat{\rho }}=[\Omega ({\hat{a}}^{†}-\hat{a}),\hat{\rho }]+{\gamma }_1^{+}\mathcal{D}[{\hat{a}}^{†}]\hat{\rho }+{\gamma }_1^{-}\mathcal{D}[\hat{a}]\hat{\rho }+{\gamma }_2\mathcal{D}[{\hat{a}}^2]\hat{\rho },$$

$${\dot{\rho }}_{n,n^{\prime }}=$$

$$\displaystyle+\gamma_{1}^{+}\Big{[}\sqrt{nn^{\prime}}\hskip 0.85355pt\rho_{n-1,n^{\prime}-1}-\Big{(}\frac{n+n^{\prime}}{2}+1\Big{)}\rho_{n,n^{\prime}}\Big{]}$$

$$\displaystyle+\gamma_{1}^{-}\Big{[}\sqrt{(n+1)(n^{\prime}+1)}\hskip 1.42271pt\rho_{n+1,n^{\prime}+1}-\frac{n+n^{\prime}}{2}\hskip 0.85355pt\rho_{n,n^{\prime}}\Big{]}$$

$$\displaystyle+\gamma_{2}\hskip 1.13791pt\Big{[}\sqrt{(n+2)(n+1)(n^{\prime}+2)(n^{\prime}+1)}\hskip 1.42271pt\rho_{n+2,n^{\prime}+2}-\frac{n(n-1)+n^{\prime}(n^{\prime}-1)}{2}\hskip 0.85355pt\rho_{n,n^{\prime}}\Big{]}\;.$$

$${\dot{\rho }}_{00}$$

$${\dot{\rho }}_{11}$$

$${\dot{\rho }}_{10}$$

$${\dot{\rho }}_{21}$$

$${\dot{\rho }}_{20}$$

$${\dot{\rho }}_{11}$$

$${\dot{\rho }}_{10}$$

$$⟨\hat{a}⟩={\rho }_{10}=\frac{2\Omega {\gamma }_1^{-}}{({\gamma }_1^{-}{)}^2+8{\Omega }^2}.$$

$${\dot{\rho }}_{00}$$

$${\dot{\rho }}_{10}$$

$${\dot{\rho }}_{22}$$

$${\dot{\rho }}_{21}$$

$${\dot{\rho }}_{20}$$

$${\rho }_{22}$$

$${\rho }_{21}$$

$${\rho }_{20}$$

$${\rho }_{21}\approx (\sqrt{2}\Omega /{\gamma }_2)[{\rho }_{11}+(\Omega /{\gamma }_2){\rho }_{10}].$$

$${\rho }_{10}\approx \left(\frac{2\Omega }{{\gamma }_2}+\frac{{\gamma }_1^{-}}{2\Omega }\right){\rho }_{11}.$$

$${\rho }_{21}$$

$${\rho }_{22}$$

$${\rho }_{11}\approx \frac{4{\Omega }^2}{({\gamma }_1^{-}{)}^2+8{\Omega }^2},$$

$$⟨\hat{a}⟩={\rho }_{10}+\sqrt{2}{\rho }_{21}\approx \left(\frac{4\Omega }{{\gamma }_2}+\frac{{\gamma }_1^{-}}{2\Omega }\right){\rho }_{11}\approx \frac{2\Omega }{{\gamma }_2}\left[\frac{{\gamma }_1^{-}{\gamma }_2+8{\Omega }^2}{({\gamma }_1^{-}{)}^2+8{\Omega }^2}\right],$$

$${\dot{\rho }}_{00}$$

Full text

A key insight from the theory of phase transitions is that systems close to a critical point are highly sensitive to perturbations Kardar (2007). 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 Strogatz (1994). A particularly important prototype is an oscillator with nonlinear damping, which can transition from a dormant to a limit-cycle phase across a Hopf bifurcation Marsden and McCracken (1976). Generically, the response of such a van der Pol (vdP) oscillator van der Pol (1926) to a resonant drive $\Omega$ grows as $\Omega^{1/3}$ at the critical point Duke and Jülicher (2008). This nonlinearity is what enables the ear and other biological sensors to detect very weak stimuli and maximally process environmental inputs Duke and Jülicher (2008); Eguíluz et al. (2000); Mora and Bialek (2011); Munõz (2018); Roli et al. (2018).

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 Degen et al. (2017); Braun et al. (2018). 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 Sieberer et al. (2016); Rotter and Bird (2015), led by advances in synthetic experimental platforms with unprecedented control and tunability Müller et al. (2012); Houck et al. (2012); Aspelmeyer et al. (2014); Noh and Angelakis (2017); Georgescu et al. (2014).

At the semiclassical level, the vdP oscillator is realized in the physics of lasers Vahala et al. (2009); Sargent III et al. (1974). Recent theoretical studies have also examined a quantum vdP oscillator, largely in the context of synchronization Lee and Sadeghpour (2013); Walter et al. (2014a); Lörch et al. (2016); Lee et al. (2014); Walter et al. (2014b); Ishibashi and Kanamoto (2017); Scarlatella et al. (2019). 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 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.

The classical vdP model describes a harmonic oscillator with nonlinear damping, $\ddot{x}=-\omega^{2}x+(2\gamma_{1}-8\gamma_{2}x^{2})\dot{x}$. Here, $\omega$ is the natural frequency, $\gamma_{1}$ is a negative damping which may arise from an energy source, and $\gamma_{2}$ is the leading nonlinear term which damps the system at large $x$. For $|\gamma_{1}|\ll\omega$, the small-amplitude oscillations are sinusoidal with a slowly-varying amplitude, $x=\text{Re}[\alpha(t)e^{{\rm i}\omega t}]$, where $\dot{\alpha}=\gamma_{1}\alpha-\gamma_{2}|\alpha|^{2}\alpha$ van der Pol (1926). Thus, if $\gamma_{1}<0$, the system is fully damped and all oscillations die out, whereas for $\gamma_{1}>0$, they reach a stable amplitude with $|\alpha|=\sqrt{\gamma_{1}/\gamma_{2}}$. When subjected to a resonant drive $F\cos(\omega t+\phi)$, the equation of motion gains a forcing term,

[TABLE]

where $\Omega=-{\rm i}e^{{\rm i}\phi}F/(2\omega)$. We choose $\phi=\pi/2$ so that $\Omega$ is real and nonnegative. From Eq. (1), the steady-state response at the critical point $\gamma_{1}=0$ is purely nonlinear, $\alpha=(\Omega/\gamma_{2})^{1/3}$. Hence, the susceptibility $\chi_{\text{cl}}\equiv\text{d}\alpha/\text{d}\Omega$ is divergent at zero drive, $\chi_{\text{cl}}\propto\Omega^{-2/3}$. Away from criticality, the response is of the form $\alpha=\smash{\tilde{\Omega}^{1/3}f(\tilde{\gamma}_{1}/\tilde{\Omega}^{2/3})}$ where $\tilde{\gamma}_{1}\equiv\gamma_{1}/\gamma_{2}$ and $\tilde{\Omega}\equiv\Omega/\gamma_{2}$ (details in Supplement 111See 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.). In particular, the limiting behavior is such that

[TABLE]

which are consistent with a diverging susceptibility, $\chi_{\text{cl}}\propto 1/|\gamma_{1}|$, at the critical point $\gamma_{1}\to 0$.

We study a quantum version of the vdP model, where a quantum harmonic oscillator is subjected to several dissipative processes, as in Refs. Lee and Sadeghpour (2013); Walter et al. (2014a); Lörch et al. (2016); Lee et al. (2014); Walter et al. (2014b); Ishibashi and Kanamoto (2017); Scarlatella et al. (2019). These are (i) one-particle loss with rate $\smash{\gamma_{1}^{-}}$, (ii) one-particle gain with rate $\smash{\gamma_{1}^{+}}$, and (iii) two-particle loss with rate $\gamma_{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 $\Omega$, given by the Hamiltonian $\hat{H}=\omega\hat{a}^{\dagger}\hat{a}+\Omega({\rm i}e^{{\rm i}\omega t}\hat{a}^{\dagger}+\text{H.c.})$. Here, $\hat{a}$ annihilates a particle in the oscillator mode and $\hbar=1$. It is convenient to work in the rotating frame where the Hamiltonian reads $\hat{H}={\rm i}\Omega(\hat{a}^{\dagger}-\hat{a})$. In general, the full dynamics of such a driven-dissipative system are governed by a master equation for the density matrix $\hat{\rho}$, obtained by tracing out the environment degrees of freedom Puri and Lawande (1977). However, for typical atomic/photonic setups, the environment relaxes to equilibrium on optical timescales, several orders of magnitude faster than the system dynamics Daley (2014). Under such a routine Markov approximation Rivas et al. (2010, 2014), the master equation reduces to the Lindblad form

[TABLE]

where $\mathcal{D}[\hat{x}]\hat{\rho}\equiv\hat{x}\hat{\rho}\hat{x}^{\dagger}-\{\hat{x}^{\dagger}\hat{x},\hat{\rho}\}/2$ Breuer and Petruccione (2002); Gardiner and Zoller (2004). Using Eq. (3), or operator equations of motion Sachdev (1984), the expectation value $\langle\hat{a}\rangle$ can be shown to obey

[TABLE]

Replacing $\hat{a}$ with a complex number $\alpha$ in Eq. (4) leads to the classical limit in Eq. (1) with $\gamma_{1}=(\gamma_{1}^{+}-\gamma_{1}^{-})/2$. Hence, the critical case corresponds to $\gamma_{1}^{+}=\gamma_{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. Lee and Sadeghpour (2013); Walter et al. (2014a); Lörch et al. (2016); Lee et al. (2014); Walter et al. (2014b); Ishibashi and Kanamoto (2017); Scarlatella et al. (2019).

We find the response by numerically solving Eq. (3) for the steady-state density matrix, then computing $\langle\hat{a}\rangle=\text{Tr}(\hat{a}\hat{\rho})$. It is most straightforward to use the Fock basis, $\hat{\rho}=\sum_{nn^{\prime}}\rho_{n,n^{\prime}}|n\rangle\langle n^{\prime}|$, where $n=0,1,2,\dots$ denotes the number of particles in the oscillator mode. The response is given by $\smash{\langle\hat{a}\rangle=\sum_{n}\sqrt{n}\rho_{n,n-1}}$. Past studies have found closed-form expressions for $\rho_{n,n^{\prime}}$ in special cases Simaan and Loudon (1975); Dodonov and Mizrahi (1997); Kheruntsyan (1999).

First we consider the simplest critical case, $\gamma_{1}^{+}=\gamma_{1}^{-}=0$. Figure 1 shows the susceptibility $\chi\equiv\text{d}\langle\hat{a}\rangle/\text{d}\Omega$ as a function of drive. Note that it coincides with the classical result for $\langle\hat{a}\rangle\gtrsim 1$. However, at weaker drives, the classical divergence is cut off and $\chi$ saturates, producing a linear response. We can understand this low-energy cutoff as follows. At $\Omega=0$, both $|0\rangle$ and $|1\rangle$ are steady states as there is only two-particle decay. A nonzero drive couples these neighboring Fock states, yielding $\rho_{10},\rho_{21}\sim\Omega$ and $\rho_{22}\sim\Omega^{2}$. Hence, to linear order, one can restrict the dynamics to the lowest three levels, which gives $\langle\hat{a}\rangle=2\Omega/\gamma_{2}$ (see Supplement Note (1)). Stronger drives inject more particles, inducing a crossover to the classical limit.

We now consider departures from $\gamma_{1}^{\pm}=0$ by allowing either $\gamma_{1}^{+}$ or $\gamma_{1}^{-}$ to be nonzero, spanning the transition from quiescent ($\gamma_{1}^{+}=0$, $\gamma_{1}^{-}>0$) to the limit-cycle regime ($\gamma_{1}^{+}>0$, $\gamma_{1}^{-}=0$). In Fig. 2 we plot the zero-drive susceptibility across the transition, which, remarkably, shows a divergence as $1/\gamma_{1}^{\pm}$ at small damping or antidamping. This divergence can be understood from a competition of energy scales. For $\gamma_{1}^{+}=0$ and $\gamma_{1}^{-}>0$, the system is fully damped and the undriven steady state is $|0\rangle$. For weak drives $\Omega<\gamma_{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 $\langle\hat{a}\rangle=2\Omega/\gamma_{1}^{-}$, matching the classical result. The response for $\gamma_{1}^{+}>0$ and $\gamma_{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 Li et al. (2014); Albert et al. (2016); Konopik and Lutz (2019). However, for $\gamma_{1}^{+}\ll\gamma_{2}$, the dynamics are confined to the lowest three levels, as shown in Fig. 2 inset, which yields $\langle\hat{a}\rangle=2\Omega/(9\gamma_{1}^{+})$ (see Supplement Note (1)).

At first sight, the results in Figs. 1 and 2 might appear to be contradictory: at $\gamma_{1}^{+}=\gamma_{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 $\chi(\Omega\to 0,\gamma_{1}^{\pm}\to 0)$ depends on the order of the two limits. This is better understood by plotting the full response as a function of $\Omega$ for small $\gamma_{1}^{\pm}$, as shown in Fig. 3. The response is strongly nonlinear for $\Omega\gtrsim\gamma_{1}^{\pm}$. Thus, although the linear susceptibility diverges, the linear region itself shrinks to zero as $\gamma_{1}^{\pm}\to~{}0$.

The full response exhibits four different regimes: (i) a linear response for $\smash{\Omega\lesssim\gamma_{1}^{\pm}}$, (ii) negative susceptibility for $\gamma_{1}^{\pm}\lesssim\Omega\lesssim(\gamma_{1}^{\pm}\gamma_{2})^{1/2}$, (iii) an extended quantum response for $(\gamma_{1}^{\pm}\gamma_{2})^{1/2}\lesssim\Omega\lesssim\gamma_{2}$, and (iv) classical response for $\smash{\Omega\gtrsim\gamma_{2}}$.

The first three regions are well reproduced by a three-level approximation, which yields, to lowest order in $\gamma_{1}^{\pm}/\gamma_{2}$ and $\Omega/\gamma_{2}$ (see Supplement Note (1) for a derivation),

[TABLE]

Hence, in the linear region, $\chi\approx 2(\gamma_{1}^{+}+\gamma_{1}^{-})/(3\gamma_{1}^{+}+\gamma_{1}^{-})^{2}$, which diverges for $\gamma_{\pm}\to 0$, whereas for $\Omega\sim(\gamma_{1}^{\pm}\gamma_{2})^{1/2}$, $\langle\hat{a}\rangle\approx 2\Omega/\gamma_{2}+(\gamma_{1}^{+}+\gamma_{1}^{-})/(4\Omega)$, 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 Gardiner and Zoller (2004). The Wigner function is widely used as an integral representation of the density matrix Note (1).

Physically, the separate regions originate from an interplay between drive and dissipation within a few-level manifold. Consider the case $\smash{\gamma_{1}^{+}=0}$ and $\smash{0<\gamma_{1}^{-}\ll\gamma_{2}}$. A weak drive couples the steady state $|0\rangle$ to $|1\rangle$, producing a coherence $\rho_{10}\sim\Omega/\gamma_{1}^{-}$. However, $\rho_{11}$ also grows with $\Omega$ and saturates at $1/2$. This saturation acts as a negative feedback for $\rho_{10}$, as seen from the equation of motion $\dot{\rho}_{11}=2\hskip 0.85355pt\Omega\hskip 0.85355pt\rho_{10}-\gamma_{1}^{-}\rho_{11}$. For $\rho_{11}\approx 1/2$, $\rho_{10}$ falls off as $\gamma_{1}^{-}/(4\Omega)$, 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 $\gamma_{1}^{\pm}=0$, with a purely quantum regime for $\langle\hat{a}\rangle\lesssim 1$ and a classical regime for $\langle\hat{a}\rangle\gtrsim 1$ (cf. Fig. 1). A similar variation with drive is found for $\smash{0<\gamma_{1}^{+}\ll\gamma_{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 driven-dissipative 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 $\smash{\dot{\hat{a}}=\Omega-(\gamma_{1}^{-}/2)\hat{a}}$, so the steady state is $\langle\hat{a}\rangle=2\Omega/\gamma_{1}^{-}$ and the passive susceptibility is $\chi_{p}=2/\gamma_{1}^{-}$. To discuss whether pumping enhances sensitivity, we define the sensitivity gain at weak drives, $G_{0}\equiv\chi/\chi_{p}|_{\Omega\to 0}$.

One regime where there is a gain, $G_{0}>1$, is the classical limit. Then the response is only a function of $\gamma_{1}^{+}-\gamma_{1}^{-}$, and one can simply set $\gamma_{1}^{+}=\gamma_{1}^{-}$ to make the system infinitely sensitive to weak signals for any damping. Quantum mechanically, however, the response also depends on the average $\Gamma_{1}\equiv(\gamma_{1}^{+}+\gamma_{1}^{-})/2$. Figure 4 shows how the classical limit emerges with increasing $\Gamma_{1}$ for the critical case, $\gamma_{1}^{+}=\gamma_{1}^{-}$. We see that $G_{0}$ scales as $\smash{\sqrt{\Gamma_{1}/\gamma_{2}}}$ and one can have $\chi>\chi_{p}$ if $\gamma_{1}^{-}\gtrsim 3\gamma_{2}$, which holds for typical experimental systems et al. (2015). The relatively slow growth of $G_{0}$ with $\Gamma_{1}$ stems from large number fluctuations at criticality (which also limits the sensing efficiency, see below). Since $\gamma_{1}^{+}=\gamma_{1}^{-}$, the undriven system would reach an infinite-temperature state in the limit $\gamma_{2}\to 0$. A nonzero $\gamma_{2}$ results in a Gaussian number distribution of width $\smash{\sqrt{\Gamma_{1}/\gamma_{2}}}$ Dodonov and Mizrahi (1997), leading to the linear susceptibility $\chi\approx 2/\sqrt{\pi\Gamma_{1}\gamma_{2}}$ (derivation in Supplement Note (1)).

The sensitivity can be enhanced much further by operating deep in the limit-cycle phase, $\gamma_{1}^{+}-\gamma_{1}^{-}\gg\gamma_{2}$, as in Fig. 5. Here, a vanishingly small drive breaks the phase-rotation symmetry, yielding a nonzero classical response. For the quantum oscillator, the coherence $\langle\hat{a}\rangle$ builds up at a finite rate (see Supplement for a derivation Note (1))

[TABLE]

Hence, the susceptibility is only limited by $\gamma_{2}$ and one can obtain a large enhancement. In particular, for $\gamma_{1}^{+}\gg\gamma_{1}^{-}$, $G_{0}\approx\gamma_{1}^{-}/(3\gamma_{2})\gg 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 $\Omega\gtrsim(\gamma_{1}^{+}\gamma_{2})^{1/2}$, as shown in Fig. 5. Such a behavior is reminiscent of sensing protocols based on switching in optically bistable systems Aldana et al. (2014).

Apart from the “signal” $\langle\hat{a}\rangle$, the efficiency of sensors is limited by the signal-to-noise ratio (SNR) Degen et al. (2017); Lau and Clerk (2018). Estimating the noise from the spread of the Wigner function, we find one can have $G_{0}>1$ and $\text{SNR}\gtrsim 1$ in the quantum domain, $\langle\hat{a}\rangle\sim 1$ (see Supplement Note (1)). 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 Lee and Sadeghpour (2013); Walter et al. (2014a), 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 $\Delta E^{-1}$, where $\Delta E$ is the level spacing between $|g\rangle$ and $|e\rangle$, which dominates all other energy scales Daley (2014). Hence, the Markov approximation holds. Additionally, one must have resolved sidebands and suppress off-resonant excitations. As shown in Ref. Lee and Sadeghpour (2013), 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 Lvovsky and Raymer (2009) and more direct mapping of the Wigner function Banaszek and Wódkiewicz (1996); Banaszek et al. (1999); Leibfried et al. (1996); Lutterbach and Davidovich (1997); Bertet et al. (2002).

Similar techniques can be used for phonon modes of an optomechanical membrane Walter et al. (2014a) or photon modes in microwave resonators Johansson et al. (2014). One can also realize an incoherent pump through spontaneous emission from a set of inverted qubits Lebreuilly et al. (2016); Ma et al. (2017). Strong two-body loss can result from Josephson junctions in superconducting circuits et al. (2015) or polariton blockade in an optical cavity Jia et al. (2018).

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. Lörch et al. (2017); Nigg (2018); Viennot et al. (2018); Rips et al. (2012). 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 Chandrasekar et al. (2007); Bender et al. (2016); Ruby et al. (2012).

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 Roli et al. (2018) and quantum-to-classical crossover Davidovich (2016); Schlosshauer (2019).

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.

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

[TABLE]

where $\alpha$ is the complex amplitude, $\gamma_{1}$ is the negative damping, $\gamma_{2}$ is the nonlinear damping, and $\Omega$ is the drive. The phase of the drive is chosen such that $\Omega$ is real and positive. Hence, $\alpha$ is real in steady state. To solve for the response, we introduce the rescaled quantities $\tilde{\gamma}_{1}\equiv\gamma_{1}/\gamma_{2}$, $\tilde{\Omega}\equiv\Omega/\gamma_{2}$, and $\tilde{\alpha}\equiv\alpha/\tilde{\Omega}^{1/3}$. In steady state,

[TABLE]

The only stable solution to Eq. (S2) is given by $\tilde{\alpha}=f(\tilde{\gamma}_{1}/\tilde{\Omega}^{2/3})$, where

[TABLE]

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., $\alpha=(\Omega/\gamma_{2})^{1/3}$ for $\gamma_{1}=0$. This is the characteristic nonlinear response utilized in biological sensors Eguíluz et al. (2000). Conversely, for large $x$, $f(x)$ can be approximated as

[TABLE]

The former limit is realized for $\gamma_{1}>0$ and $\Omega/\gamma_{1}\ll\sqrt{\gamma_{1}/\gamma_{2}}$, where the steady state corresponds to a perturbed limit cycle, \alpha=\sqrt{\gamma_{1}/\gamma_{2}}+\Omega/(2\gamma_{1})+O\big{[}(\Omega/\gamma_{1})^{2}\sqrt{\gamma_{2}/\gamma_{1}}\hskip 0.85355pt\big{]}. The latter limit is realized for $\gamma_{1}<0$ and $\Omega/|\gamma_{1}|\ll\sqrt{|\gamma_{1}|/\gamma_{2}}$. Here the system is fully damped and the response grows linearly with drive, $\alpha=\Omega/|\gamma_{1}|+O[(\Omega/\gamma_{1})^{3}(\gamma_{2}/\gamma_{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),

[TABLE]

where $\hat{\rho}$ is the density matrix, $\hat{a}$ is the ladder operator, $\Omega$ is the drive amplitude, $\gamma_{1}^{+}$ is the rate of one-particle gain, $\gamma_{1}^{-}$ is the rate of one-particle decay, $\gamma_{2}$ is the rate of two-particle decay, and $\mathcal{D}[\hat{x}]\hat{\rho}\equiv\hat{x}\hat{\rho}\hat{x}^{\dagger}-\{\hat{x}^{\dagger}\hat{x},\hat{\rho}\}/2$. In the Fock basis $\{|n\rangle,n=0,1,2,\dots\}$, Eq. (S5) can be written as

[TABLE]

Note that $\Omega$ couples neighboring elements of the density matrix, whereas the dissipative terms only couple elements in the same diagonal, i.e., with $n-n^{\prime}=$ constant. The response is given by $\langle\hat{a}\rangle=\text{Tr}(\hat{a}\hat{\rho})=\sum_{n}\sqrt{n}\rho_{n,n-1}$, where $\rho$ is the steady-state density matrix. Since all coefficients in Eq. (S6) are real, $\rho$ will be real and symmetric. Below we consider this steady-state response in various limits.

SII.1 Linear response for $\gamma_{1}^{+}=\gamma_{1}^{-}=0$

For $\gamma_{1}^{\pm}=0$ and weak drives $\Omega\ll\gamma_{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

[TABLE]

where $\rho_{22}=1-\rho_{00}-\rho_{11}$. Substituting $\Omega\to 0$ into these equations, we find the undriven steady state is given by $\rho_{22}=\rho_{21}=\rho_{20}=0$ and $\rho_{00}=\rho_{11}=1/2$. To linear order in $\Omega/\gamma_{2}$, Eqs. (S7b) and (S7d) yield $\rho_{10}=\sqrt{2}\rho_{21}=\Omega/\gamma_{2}$. Thus, the linear response is given by $\langle\hat{a}\rangle=\rho_{10}+\sqrt{2}\rho_{21}=2\Omega/\gamma_{2}$, as evident in Fig. 1 of the main text.

SII.2 Nonmonotonic response for $\gamma_{1}^{\pm}\ll\gamma_{2}$

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

[TABLE]

where $\rho_{00}=1-\rho_{11}$. In steady state, we find

[TABLE]

Thus, the response grows linearly as $\langle\hat{a}\rangle\approx 2\Omega/\gamma_{1}^{-}$ for $\Omega\lesssim\gamma_{1}^{-}$, but falls off as $\langle\hat{a}\rangle\approx\gamma_{1}^{-}/(4\Omega)$ for $\Omega\gtrsim\gamma_{1}^{-}$. This is shown by the dashed black curve in Fig. S2. We see the linear susceptibility diverges as $2/\gamma_{1}^{-}$ for $\gamma_{1}^{-}\to 0$.

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

[TABLE]

Let us solve for the steady state to lowest order in $\gamma_{1}^{-}/\gamma_{2}$ and $\Omega/\gamma_{2}$. From the last three equations, we can write

[TABLE]

Substituting $\rho_{22}$ and $\rho_{20}$ into Eq. (S11b), we find

[TABLE]

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

[TABLE]

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

[TABLE]

Finally, substituting Eqs. (S13)–(S15) into Eq. (S10b) and using $\rho_{00}=1-\rho_{11}-\rho_{22}$, we find

[TABLE]

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

[TABLE]

where we have used Eqs. (S13), (S14), and (S16). It reduces to the two-level expression in Eq. (S9) for $\Omega\ll(\gamma_{1}^{-}\gamma_{2})^{1/2}$. For $\Omega\sim(\gamma_{1}^{-}\gamma_{2})^{1/2}\gg\gamma_{1}$, we find $\langle\hat{a}\rangle\approx\gamma_{1}^{-}/(4\Omega)+2\Omega/\gamma_{2}$, which attains a minimum at $\Omega=(\gamma_{1}^{-}\gamma_{2})^{1/2}/(2\sqrt{2})$ and increases linearly with drive at large $\Omega$. This nonmonotonic variation is evident in Fig. S2.

SII.2.2 $0\leq\gamma_{1}^{+},\gamma_{1}^{-}\ll\gamma_{2}$

In the presence of one-particle gain ($\gamma_{1}^{+}>0$), the undriven steady state corresponds to a dynamic equilibrium where particles flow in and out of the oscillator. When $\gamma_{2}\gg\gamma_{1}^{\pm}$, 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)]

[TABLE]

where $\rho_{11}=1-\rho_{00}-\rho_{22}$. For $\Omega=0$, we find the steady state $\rho_{11}\approx 1-\rho_{00}\approx\gamma_{1}^{+}/(3\gamma_{1}^{+}+\gamma_{1}^{-})$. Thus, negative damping leads to a finite occupation of the $n=1$ level. For $\Omega>0$, the steady-state analysis can be carried out following the same steps as in the last section, which yields, to lowest order in $\gamma_{1}^{+}/\gamma_{2}$, $\gamma_{1}^{-}/\gamma_{2}$, and $\Omega/\gamma_{2}$,

[TABLE]

As before, we find a nonmonotonic response. For $\Omega\ll 3\gamma_{1}^{+}+\gamma_{1}^{-}$, the response is linear, $\langle\hat{a}\rangle\approx 2\Omega(\gamma_{1}^{+}+\gamma_{1}^{-})/(3\gamma_{1}^{+}+\gamma_{1}^{-})^{2}$, with a slope that diverges as $\gamma_{1}^{\pm}\to 0$. For $3\gamma_{1}^{+}+\gamma_{1}^{-}\ll\Omega\ll[(\gamma_{1}^{+}+\gamma_{1}^{-})\gamma_{2}]^{1/2}$, it falls off as $\langle\hat{a}\rangle\approx(\gamma_{1}^{+}+\gamma_{1}^{-})/(4\Omega)$, exhibiting negative susceptibility. For $\Omega\gg[(\gamma_{1}^{+}+\gamma_{1}^{-})\gamma_{2}]^{1/2}$, it rises again as $\langle\hat{a}\rangle\approx 2\Omega/\gamma_{2}$.

SII.3 Linear response for $\gamma_{1}^{+}=\gamma_{1}^{-}\gg\gamma_{2}$

Here, we derive the enhanced sensitivity for the critical condition $\gamma_{1}^{+}=\gamma_{1}^{-}\equiv\Gamma_{1}\gg\gamma_{2}$. We first find the undriven steady state, then consider perturbation at weak drives. At $\Omega=0$, the steady state has no coherence. The equation of motion for the populations $p_{n}\equiv\rho_{n,n}$ is found by substituting $n^{\prime}=n$ in Eq. (S6), which yields

[TABLE]

At $\gamma_{2}=0$, the steady state is simply given by $p_{n}=$ constant, i.e., the system reaches an infinite-temperature state. For $\gamma_{2}>0$, this distribution is curtailed at large $n$ and $p_{n}$ is a slowly-varying function set by the scale $\varepsilon\equiv\sqrt{\gamma_{2}/\Gamma_{1}}$. Thus, we can approximate $p_{n}$ by a continuous function $p(x)$ with $x=n\varepsilon$. Then Eq. (S20) can be written as

[TABLE]

The above equation is identically satisfied up to $O(\varepsilon)$. Expanding to $O(\varepsilon^{2})$, we find

[TABLE]

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 $\sum_{n}p_{n}=1$, or $(1/\varepsilon)\int_{0}^{\infty}p(x)dx=1$, which gives

[TABLE]

Thus, the undriven steady state has a Gaussian number distribution of width $\sqrt{\Gamma_{1}/\gamma_{2}}\gg 1$, as shown in Fig. S3(a).

To find the linear response, we consider perturbation in the leading off-diagonal terms $\chi_{n}\equiv\rho_{n,n-1}$, which satisfy

[TABLE]

This equation is simplified with the substitution $q_{n}\equiv\sqrt{n}\chi_{n}$. In steady state,

[TABLE]

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\varepsilon$,

[TABLE]

where $\eta\equiv\Omega/\Gamma_{1}$ and $\tilde{p}(x)\equiv p(x)/\varepsilon\sim O(1)$. Expanding Eq. (S26) to $O(\varepsilon^{2})$, we find

[TABLE]

Substituting $\tilde{p}\hskip 0.56917pt^{\prime}(x)=-(4/\sqrt{\pi})xe^{-x^{2}}$ from Eq. (S23) and using $q(x)\equiv(4\eta/\sqrt{\pi})e^{-x^{2}}h(x)$ yields

[TABLE]

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,

[TABLE]

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

[TABLE]

Hence, the linear susceptibility is $\chi\approx 2/\sqrt{\pi\Gamma_{1}\gamma_{2}}$. Figure S3(b) shows how the exact solution interpolates between this regime and that of $\Gamma_{1}\ll\gamma_{2}$, where $\chi\approx 1/(4\Gamma_{1})$ [see Sec. SII.2.2]. In contrast, a fully damped “passive” oscillator has a susceptibility $\chi_{p}=2/\gamma_{1}^{-}$ [see Sec. SII.2.1]. Thus, operating the vdP oscillator in the regime $\gamma_{1}^{+}=\gamma_{1}^{-}\gg\gamma_{2}$ yields a sensitivity enhancement $\chi/\chi_{p}\approx[\gamma_{1}^{-}/(\pi\gamma_{2})]^{1/2}$.

SII.4 Linear response for $\gamma_{1}^{+}-\gamma_{1}^{-}\gg\gamma_{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,

[TABLE]

The steady state is set by the parameters $\varepsilon\equiv(\gamma_{2}/\gamma_{1}^{+})^{1/2}\ll 1$ and $\zeta\equiv\gamma_{1}^{-}/\gamma_{1}^{+}<1$. We expect, in steady state, $p_{n}$ will be peaked near the classical solution, $n^{*}\approx\alpha_{\text{cl}}^{2}=(1-\zeta)/(2\varepsilon^{2})$. To find this distribution, we approximate $p_{n}$ by a continuous function $p(x)$ where $x=n\varepsilon$. Then Eq. (S31) yields

[TABLE]

Since $p(x)$ will be peaked near $x^{*}=n^{*}\varepsilon\sim O(1/\varepsilon)$, we shift the origin with the transformation $x=\beta/\varepsilon+y$, where $\beta\sim O(1)$. We also define $u(y)\equiv p(\beta/\varepsilon+y)$. Then expanding Eq. (S32) in powers of $\varepsilon$, we find

[TABLE]

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

[TABLE]

which has a unique positive-definite solution, $u(y)=Ae^{-2y^{2}/(3-\zeta)}$. The integration constant $A$ can be determined by requiring $\sum_{n}p_{n}=1$, or $(1/\varepsilon)\int u(y)dy\approx 1$, which yields

[TABLE]

From the correspondence $p_{n}\approx u(n\varepsilon-\beta/\varepsilon)$ and the definitions $\zeta\equiv\gamma_{1}^{-}/\gamma_{1}^{+}$, $\varepsilon\equiv(\gamma_{2}/\gamma_{1}^{+})^{1/2}$, we find the undriven steady state has a Gaussian number distribution with mean $(\gamma_{1}^{+}-\gamma_{1}^{-})/(2\gamma_{2})$ and standard deviation $[(3\gamma_{1}^{+}-\gamma_{1}^{-})/(4\gamma_{2})]^{1/2}$. This result was also found in Ref. Dodonov and Mizrahi (1997) and agrees well with numerics, as shown in Fig. S4(a).

Next, we consider linear perturbation at weak drives. From Eq. (S6), the coherences $\chi_{n}\equiv\rho_{n,n-1}$ satisfy

[TABLE]

where $p_{n}$ are the undriven occupations. To find the steady state, we again approximate $\chi_{n}$ by a continuous function $v(y)$ where $\chi_{n}\approx v(n\varepsilon-\beta/\varepsilon)$. Following the same steps as before, we arrive at the solution

[TABLE]

where $\eta\equiv\Omega/\gamma_{1}^{+}$. Thus, $\chi_{n}$ has the same Gaussian profile as $p_{n}$, as seen in Fig. S4(a). The linear response is given by

[TABLE]

Expressing $\eta$, $\varepsilon$, and $\zeta$ in terms of the drive and dissipation rates, we find a susceptibility

[TABLE]

Thus, for $\gamma_{1}^{-}\ll\gamma_{1}^{+}$, $\chi$ is only limited by two-particle decay. Consequently, the sensitivity gain over a passive system, $\chi/\chi_{p}$, scales as $\gamma_{1}^{-}/\gamma_{2}$ [Recall: $\chi_{p}=2/\gamma_{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. Wigner (1932) 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 $\hat{a}$, it can be written as the expectation of a displaced parity operator Cahill and Glauber (1969a)

[TABLE]

where $\alpha$ is a point in phase space and $\hat{D}(\alpha)\equiv e^{\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a}}$ is a displacement operator. The Wigner function is given by

[TABLE]

Conversely, the density operator $\hat{\rho}$ can be obtained from the Wigner function as

[TABLE]

So, there is a one-to-one correspondence between $\hat{\rho}$ and $W$. Note the operators $\hat{\Pi}(\alpha)$ are Hermitian with eigenvalues $\pm 1$. They can also be shown to constitute a basis for expanding any operator over the complex phase space Cahill and Glauber (1969b). In particular, they satisfy the orthonormality and completeness relations

[TABLE]

where $|n\rangle$ are Fock states. Since $\hat{\Pi}(\alpha)$ are parity operators, from Eq. (S41) we see that $W$ is real valued and uniformly bounded, $W(\alpha)\in[-2/\pi,2/\pi]$. Further, using the property $\text{Tr}[\hat{\Pi}(\alpha)]=1/2$ in Eq. (S42), one finds

[TABLE]

Hence, $W(\alpha)$ can be interpreted as a quasiprobability distribution in phase space. In particular, one can compute the expectation of any symmetrically ordered product $\{(\hat{a}^{\dagger})^{n}\hat{a}^{m}\}_{S}$ as an ensemble average,

[TABLE]

which follows from using the relation $\hat{D}^{-1}(\alpha)\hat{a}\hat{D}(\alpha)=\hat{a}+\alpha\hat{\mathds{1}}$ in Eq. (S42). The linear response $\langle\hat{a\rangle}$ corresponds to a special case of Eq. (S44) that amounts to measuring the center of mass of the Wigner distribution,

[TABLE]

In principle, Eqs. (S41) and (S42) completely specify $W(\alpha)$ in terms of $\hat{\rho}$ 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)^{\hat{a}^{\dagger}\hat{a}}$ reflects $\hat{a}$ to $-\hat{a}$, $(-1)^{\hat{a}^{\dagger}\hat{a}}\hat{a}(-1)^{\hat{a}^{\dagger}\hat{a}}=-\hat{a}$. It follows we have the identity Cahill and Glauber (1969a)

[TABLE]

Using the above in Eq. (S40), we can write $\hat{\Pi}(\alpha)=\hat{D}(2\alpha)(-1)^{\hat{a}^{\dagger}\hat{a}}$ which, when used in Eq. (S41), yields

[TABLE]

The matrix elements of $\hat{D}$ can be expressed in terms of associated Laguerre polynomials $L^{(p)}_{q}(x)$ Cahill and Glauber (1969b). Thus, we find

[TABLE]

The mapping can be written in a more symmetric form in polar coordinates $\alpha=re^{i\phi}$,

[TABLE]

where $(x)_{j}$ is the Pochhammer symbol, $(x)_{j}=x(x+1)\dots(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 $\rho_{n,n^{\prime}}$ in terms of $W(r,\phi)$,

[TABLE]

SIV Signal-to-noise ratio

We have focused on the response $\langle\hat{a}\rangle$ 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” $\langle\hat{a}\rangle$ as well as other physical quantities can be extracted from the Wigner function $W(\alpha)$ which represents a quasiprobability distribution in phase space. From Eq. (S45), $\langle\hat{a}\rangle$ is given by the center of mass of $W(\alpha)$. The noise $\sigma$ can be estimated from the spread about this center of mass,

[TABLE]

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

[TABLE]

where $\hat{n}$ is the number operator, $\hat{n}\equiv\hat{a}^{\dagger}\hat{a}$. The signal-to-noise ratio is given by $\text{SNR}=\langle\hat{a}\rangle/\sigma$. We plot the response and the SNR as a function of drive and dissipation in Fig. S5, which shows one can have $\text{SNR}\gtrsim 1$ while being in the quantum regime, $\langle\hat{a}\rangle\sim 1$. In Fig. S5(a), we consider the case $\gamma_{1}^{\pm}\ll\gamma_{2}$, where one finds a nonmonotonic response with increasing drive. Here $\text{SNR}<1$ as the oscillator is close to the vacuum state. Nonetheless, such states have been measured in experiments with close to 1% accuracy Banaszek et al. (1999); Lvovsky et al. (2001); Bertet et al. (2002); et al. (2009); Chu et al. (2018); Smithey et al. (1993); Leibfried et al. (1996). In Fig. S5(b), we find, for $\gamma_{1}^{\pm}\gg\gamma_{2}$, one can have both $\text{SNR}>1$ as well as a susceptibility boost over a passive oscillator, $\chi>\chi_{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 $\Delta n\sim(\gamma_{1}^{\pm}/\gamma_{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 Degen et al. (2017).