General linearized theory of quantum fluctuations around arbitrary limit cycles

TL;DR

This paper develops a linearization method for analyzing quantum fluctuations around time-dependent limit cycles in nonlinear quantum systems, extending standard techniques to non-stationary states like spontaneous oscillations.

Contribution

It introduces a novel linearization scheme tailored for systems with time-dependent classical solutions, such as limit cycles, and demonstrates its effectiveness through the driven Van der Pol oscillator.

Findings

The new method accurately captures quantum fluctuations around limit cycles.

It maintains computational simplicity and scalability for large systems.

Comparison with numerical simulations validates the approach.

Abstract

The theory of Gaussian quantum fluctuations around classical steady states in nonlinear quantum-optical systems (also known as standard linearization) is a cornerstone for the analysis of such systems. Its simplicity, together with its accuracy far from critical points or situations where the nonlinearity reaches the strong coupling regime, has turned it into a widespread technique, which is the first method of choice in most works on the subject. However, such a technique finds strong practical and conceptual complications when one tries to apply it to situations in which the classical long-time solution is time dependent, a most prominent example being spontaneous limit-cycle formation. Here we introduce a linearization scheme adapted to such situations, using the driven Van der Pol oscillator as a testbed for the method, which allows us to compare it with full numerical simulations. On a conceptual level, the scheme relies on the connection between the emergence of limit cycles and the spontaneous breaking of the symmetry under temporal translations. On the practical side, the method keeps the simplicity and linear scaling with the size of the problem (number of modes) characteristic of standard linearization, making it applicable to large (many-body) systems.

Equations

$$\frac{d\hat{\rho }}{d\tau }=\left[\frac{F}{\sqrt{\gamma }}({\hat{a}}^{†}-\hat{a})+\mathrm{i}\Delta {\hat{a}}^{†}\hat{a},\hat{\rho }\right]+\frac{\gamma }{2}{\mathcal{D}}_{a^2}[\hat{\rho }]+{\mathcal{D}}_{a^{†}}[\hat{\rho }],$$

$$\dot{\beta }$$

$${\dot{\beta }}^{+}$$

$$\beta (\tau +\theta )=\bar{\beta }(\tau +\theta )+b(\tau +\theta ),$$

$${\beta }^{+}(\tau +\theta )={\bar{\beta }}^{\ast }(\tau +\theta )+b^{+}(\tau +\theta ).$$

$$\dot{\mathbf{b}}(\tau )+{\mathbf{p}}_0(\tau )\dot{\theta }(\tau )=\mathcal{L}(\tau )\mathbf{b}(\tau )+\sqrt{\gamma }\mathbf{n}(\tau )\text{,}$$

$$\mathcal{L}(\tau)=\left(\begin{array}[]{cc}1-2|\bar{\beta}(\tau)|^{2}+\mathrm{i}\Delta&-\bar{\beta}(\tau)^{2}\\ -\bar{\beta}^{*2}(\tau)&1-2|\bar{\beta}(\tau)|^{2}-\mathrm{i}\Delta\end{array}\right),$$

$$\mathcal{N}(\tau)=\left(\begin{array}[]{cc}-\bar{\beta}^{2}(\tau)&2\\ 2&-\bar{\beta}^{\ast 2}(\tau)\end{array}\right).$$

$${\dot{\mathbf{p}}}_j\left(\tau \right)$$

$${\dot{\mathbf{q}}}_j^{†}\left(\tau \right)$$

$$⟨[\theta (\tau )-\theta (0){]}^2⟩=\gamma {\int }_0^{\tau }d{\tau }^{\prime }{\mathbf{q}}_0^{†}({\tau }^{\prime })\mathcal{N}({\tau }^{\prime }){\mathbf{q}}_0^{\ast }({\tau }^{\prime }).$$

$$W(\mathbf{r},\tau +\theta )=\frac{e^{-\frac{1}{2}[\mathbf{r}-\bar{\mathbf{d}}(\tau +\theta ){]}^T{\bar{V}}^{-1}(\tau +\theta )[\mathbf{r}-\bar{\mathbf{d}}(\tau +\theta )]}}{2\pi \sqrt{\det \{\bar{V}(\tau +\theta )\}}},$$

$$\bar{\mathbf{d}}(\tau )$$

$$\bar{V}(\tau )$$

$$C(\tau )=\underset{\tau \to \infty }{\lim }{\int }_0^{\tau }d{\tau }^{\prime }{e}^{2{\mu }_1(\tau -{\tau }^{\prime })}{\mathbf{q}}_1^{†}({\tau }^{\prime })\mathcal{N}({\tau }^{\prime }){\mathbf{q}}_1^{\ast }({\tau }^{\prime }),$$

$$\bar{W}(\mathbf{r})={\int }_0^T\frac{d\theta }{T}W(\mathbf{r};\tau +\theta )={\int }_0^T\frac{d\theta }{T}W(\mathbf{r};\theta ).$$

$$\hat{\rho }(\tau )={\int }_{{\mathbb{C}}^2}{d}^4\alpha P(\alpha ;\tau )\underset{\hat{\Lambda }(\alpha )}{\underbrace{\frac{|\alpha ⟩⟨{\alpha }^{+\ast }|}{⟨{\alpha }^{+\ast }|\alpha ⟩}}},$$

$$⟨{\hat{a}}^{†m}{\hat{a}}^n⟩={\int }_{{\mathbb{C}}^2}{d}^4\alpha P(\alpha ){\alpha }^{+m}{\alpha }^n.$$

$$\hat{a}\hat{\Lambda }=\alpha \hat{\Lambda }\text{, }\hat{\Lambda }{\hat{a}}^{†}={\alpha }^{+}\hat{\Lambda }\text{, }\hat{\Lambda }\hat{a}=(\alpha +{\partial }_{{\alpha }^{+}})\hat{\Lambda }\text{, }{\hat{a}}^{†}\hat{\Lambda }=({\alpha }^{+}+{\partial }_{\alpha })\hat{\Lambda },$$

$${\int }_{{\mathbb{C}}^2}{d}^4\alpha \hat{\Lambda }(\alpha ){\partial }_{\tau }P(\alpha ;\tau )={\int }_{{\mathbb{C}}^2}{d}^4\alpha P(\alpha ;\tau )\left[\sum _{j=\alpha ,{\alpha }^{+}}{A}_j(\alpha ){\partial }_j+\frac{1}{2}\sum _{j,l=\alpha ,{\alpha }^{+}}{D}_{jl}(\alpha ){\partial }_j{\partial }_l\right]\hat{\Lambda }(\alpha ).$$

$$\mathbf{A}$$

$$\mathcal{D}$$

$${\partial }_{\tau }P(\alpha ;\tau )=\left[-\sum _{j=\alpha ,{\alpha }^{+}}{\partial }_j{A}_j(\alpha )+\frac{1}{2}\sum _{j,l=\alpha ,{\alpha }^{+}}{\partial }_j{\partial }_l{D}_{jl}(\alpha )\right]P(\alpha ;\tau ).$$

$$\dot{\mathbf{\alpha }}=\mathbf{A}+\mathcal{B}\mathbf{\eta }(\tau )\text{,}$$

$$⟨{\hat{a}}^{†m}{\hat{a}}^n⟩={\int }_{{\mathbb{C}}^2}{d}^4\alpha P(\alpha ;\tau ){\alpha }^{+m}{\alpha }^n=⟨{\alpha }^{+m}[\tau ;\mathbf{\eta }]{\alpha }^n[\tau ;\mathbf{\eta }]{⟩}_{\text{stochastic}},$$

$$\mathcal{B}=\left(\begin{array}[]{cccc}\mathrm{i}\sqrt{\gamma}\alpha&0&1&\mathrm{i}\\ 0&-\mathrm{i}\sqrt{\gamma}\alpha^{+}&1&-\mathrm{i}\end{array}\right),$$

$$\dot{\alpha }$$

$${\dot{\alpha }}^{+}$$

$$⟨\eta (\tau )\eta ({\tau }^{\prime })⟩=⟨{\eta }^{+}(\tau ){\eta }^{+}({\tau }^{\prime })⟩=⟨\xi (\tau ){\xi }^{\ast }({\tau }^{\prime })⟩=\delta (\tau -{\tau }^{\prime }),$$

$$\dot{\beta }$$

$${\dot{\beta }}^{+}$$

$$\beta (\tau +\theta )=\bar{\beta }(\tau +\theta )+b(\tau +\theta ),{\beta }^{+}(\tau +\theta )={\bar{\beta }}^{\ast }(\tau +\theta )+b^{+}(\tau +\theta ),$$

$$\frac{d}{d\tau }\theta (\tau +\theta (\tau +\theta (...)))$$

$$\approx {\partial }_{\tau }\theta (\tau +\theta (\tau +\theta (...))),$$

Full text

Introduction. The advent of modern quantum technologies has triggered the discovery of a plethora of optical, atomic, and solid state devices working in the quantum regime Dowling03 (see also the starting paragraph of Benito16 ** **and the references therein). A first-principles approach leads to a description of such devices as open quantum systems evolving according to nonlinear Hamiltonians and incoherent processes like dissipation GardinerZoller ; BreuerPetruccione ; Carmichael1 ; GardinerZollerI . Mathematically, one has to face master equations for the state of the system or quantum Langevin equations for its operators, which are in general impossible to solve exactly.

On the other hand, quantum nonlinearities are very difficult to observe in the laboratory and therefore most experiments are well described by effective linear models. The most widespread method for obtaining such linear models starting from nonlinear ones is the so-called standard linearization Drummond81 ; Lugiato81 , which consists in a Gaussian-state ansatz centered at the solution of the system’s nonlinear equations in the classical limit CNB14 . The method combines incredible simplicity with pretty good accuracy in regions of the phase diagram where the system shows a finite number of well-spaced classical attraction points. However, it relies on two properties of the system’s state in the classical limit: It has to be stationary and stable along all directions of phase space. The first condition precludes its application to regions where the classical solutions are time dependent (such as limit cycles StrogatzBook ; GrimshawBook , ubiquitous to, e.g., lasing, second-harmonic generation, or optomechanical systems). The second condition excludes the possibility of applying it to systems which, being invariant under continuous transformations of some kind, have a classical solution which breaks that invariance via spontaneous symmetry breaking. This is because Goldstone’s theorem implies the existence of a zero eigenvalue of the linear stability matrix, and hence a direction of phase space which is not damped CNBthesis ; CNB08 ; CNB10 ; Perez06 ; Perez07 .

While standard linearization has been generalized to deal with spontaneous symmetry breaking of spatial, polarization, and phase symmetries CNBthesis ; CNB08 ; CNB10 ; Garcia09 ; Garcia10 ; Lane88 ; Perez06 ; Perez07 ; Reid88 ; Reid89 , an extension capable of dealing with limit cycles remains. In the case of spontaneous symmetry breaking the trick consists on using a phase-space representation of the state to keep track of the phase-space variable associated to the system’s invariance, which will carry the largest part of the fluctuations. Then, the theory can be linearized with respect to any other phase-space variable.

In this work we generalize standard linearization to regions where the classical long-time solution is time dependent, in particular describing a periodic orbit in phase space. Our idea relies on the connection between the emergence of such limit cycles, and the spontaneous breaking of a very particular symmetry: arbitrary translations in time.

For convenience, in this work we introduce the method for single-mode problems, using the driven quantum Van der Pol (VdP) oscillator Walter14 ; Weiss17 ; Lorch16 as an example. The simplicity of this model will allow for comparisons with full numerical simulations. The generalization to multi-mode problems is straightforward, and will be explored in the future for more practical and complex problems such as optomechanical cavities deep into the parametric instability regime Lorch14 ; Qian12 . Moreover, the complexity of the method scales only linearly with the number of modes, providing then an efficient route towards the analysis of many-body systems out of equilibrium such as optomechanical arrays Arrays1 ; Arrays2 ; Arrays3 ; Arrays4 ; ArraysWeiss in the self-sustained oscillations regime.

**Van der Pol model. **The quantum model for a driven VdP oscillator consists of a single bosonic mode with annihilation operator $\hat{a}$, whose state $\hat{\rho}$ evolves according to the master equation Walter14 ; Weiss17

[TABLE]

where $\mathcal{D}_{J}[\hat{\rho}]=2\hat{J}\hat{\rho}\hat{J}^{\dagger}-\hat{J}^{\dagger}\hat{J}\hat{\rho}-\hat{\rho}\hat{J}^{\dagger}\hat{J}$ and the bosonic operators satisfy canonical commutation relations $[\hat{a},\hat{a}]=0$ and $[\hat{a},\hat{a}^{\dagger}]=1$. The Hamiltonian includes a coherent drive at rate $F/\sqrt{\gamma}>0$ detuned by $\Delta$ with respect to the natural oscillation frequency of the oscillator (note that we work in a picture rotating at the driving frequency). The model contains two incoherent terms as well, the first one corresponding to pairs of excitations lost irreversibly at rate $\gamma$ (nonlinear losses), and the second one to linear pumping. The rate of the latter is used to normalize the rest of rates and frequencies, while its inverse normalizes time, so that $\tau$, $\gamma$, $F$, and $\Delta$ are dimensionless. We show later that with these choices the classical phase diagram of the system is determined uniquely by $F$ and $\Delta$, while $\gamma$ determines the strength of the quantum fluctuations.

The method is best introduced by mapping the master equation to a set of stochastic equations. This can be done with the help of phase-space quasiprobability distributions GardinerZoller ; GardinerZollerI ; Carmichael1 ; SchleichBook such as standard Wigner, Husimi, or Glauber-Sudharsan representations. Here we choose the positive P representation Drummond80 ; CarmichaelBook2 ; GardinerZoller ; GardinerZollerI because, unlike the previous representations, it always leads to stochastic equations equivalent to the master equation without any approximation. This representation associates two independent stochastic variables that we denote by $\beta/\sqrt{\gamma}$ and $\beta^{+}/\sqrt{\gamma}$ with the annihilation and creation operators $\hat{a}$ and $\hat{a}^{\dagger}$, respectively, in such a way that normally-ordered quantum expectation values and stochastic averages are related by $\langle\hat{a}^{\dagger m}\hat{a}^{n}\rangle=\langle\beta^{+m}\beta^{n}\rangle/\gamma^{(m+n)/2}$, with $m,n\in\mathbb{N}$. Using standard techniques Drummond80 ; CarmichaelBook2 ; GardinerZoller ; GardinerZollerI ; GardinerBook , we show in SupMat that the stochastic amplitudes evolve according to

[TABLE]

where $\eta(\tau)$, $\eta^{+}(\tau)$, and $\xi(\tau)$ are independent white Gaussian noises (real the first two, and complex the last one).

**Limit cycles in the classical limit. **Coming from a normally ordered representation (where vacuum noise is already taken into account in the ordering), the equations above predict a large-amplitude coherent state for $\gamma\rightarrow 0$. We talk then about the classical limit. The remaining deterministic equation $\dot{\beta}=F+(1+\mathrm{i}\Delta-|\beta|^{2})\beta$ is a paradigm for synchronization phenomena Weiss17 , and its phase diagram is well known (we provide an overview of it in SupMat ). In general terms, its stationary solutions, corresponding to solutions oscillating at the driving frequency, are stable only provided a strong enough drive is fed; otherwise, the oscillations are not synchronized to the drive, so that for long times the system ends up in a nontrivial stable periodic solution $\bar{\beta}(\tau)=\bar{\beta}(\tau+T)$ which we call *limit cycle *or periodic orbit StrogatzBook ; GrimshawBook . In Fig. 1 we show an example of such solution, where it can be appreciated that it describes a closed curve in phase space (a), with an absolute value and a phase that oscillate periodically (b). Note that analytical solutions for these limit cycles exist only in limited cases, and therefore one needs to find them numerically in general.

Linearization around limit cycles. We are now able to introduce the linearization technique for quantum fluctuations around limit cycles. We start by expanding the stochastic amplitudes as

[TABLE]

Here, $\theta$ determines at which point of the cycle the solution $\bar{\beta}(\tau+\theta)$ starts for $\tau=0$, and it is precisely the parameter which is not fixed by the classical equations of motion: $\bar{\beta}(\tau+\theta)$ is a solution of the equations for any choice of $\theta$. Owed to this symmetry, quantum fluctuations cannot be considered small in arbitrary points and directions of phase space, as nothing prevents them from acting on $\theta$ without resistance. Hence, in order for any linearized theory of quantum fluctuations to work, $\theta$ has to be taken as a variable itself (making it time dependent in the expansion above) and only then the fluctuations $b$ and $b^{+}$ can be taken as small quantities. In addition, $\dot{\theta}$ can be taken as small quantity as well, since variations of $\theta$ are induced by quantum noise, which is weak in the region of interest. Introducing (3) in (2), to first order in the small variables (including noise) we then get SupMat

[TABLE]

where $\mathbf{b}=(b,b^{+})^{T}$, $\mathbf{p}_{0}=(\partial_{\tau}\bar{\beta},\partial_{\tau}\bar{\beta}^{\ast})^{T}$, $\mathbf{n}(\tau)=[\sqrt{2}\xi(\tau)+\mathrm{i}\bar{\beta}(\tau)\eta(\tau),\sqrt{2}\xi(\tau)-\mathrm{i}\bar{\beta}^{\ast}(\tau)\eta^{+}(\tau)]^{T}$, and

[TABLE]

is the linear stability matrix. Note that the noise correlations can be written in the compact form $\langle n_{j}(\tau)n_{l}(\tau^{\prime})\rangle=N_{jl}(\tau)\delta(\tau-\tau^{\prime})$, where $N_{jl}$ are the elements of the diffusion matrix

[TABLE]

As we will see, the introduction of $\theta(\tau)$ as an explicit variable will allow us to describe properly spontaneous temporal symmetry breaking and its associated undamped phase-space direction.

Floquet method and eigenvectors. The main difference of Eq. (4) with respect to the linearized Langevin equations found in previous linearization methods is the time periodicity of $\textbf{p}_{0}(\tau)$ and $\mathcal{L}(\tau)$. We deal with this by applying Floquet theory CodingtonBook ; GrimshawBook as we explain next.

Let us define the fundamental matrix $\mathcal{R}(\tau)$, which satisfies the initial value problem $\dot{\mathcal{R}}(\tau)=\mathcal{L}(\tau)\mathcal{R}(\tau)$ with $\mathcal{R}(0)=\mathcal{I}$, the latter being the identity matrix. From it, we further define the matrix $\mathcal{M}$ through $\exp(\mathcal{M}T)=\mathcal{R}(T)$, and the $T$*-*periodic matrix $\mathcal{P}(\tau)=\mathcal{R}(\tau)\exp(-\mathcal{M}\tau)$. Given the eigensystem $\{\mathbf{v}_{j},\mathbf{w}_{j};\mu_{j}\}_{j=0,1}$ of $\mathcal{M}$, composed of right and left orthogonal ($\mathbf{w}_{j}^{\dagger}\mathbf{v}_{l}=\delta_{jl}$) eigenvectors satisfying $\mathcal{M}\mathbf{v}_{j}=\mu_{j}\mathbf{v}_{j}$ and $\mathbf{w}_{j}^{\dagger}\mathcal{M}=\mu_{j}\mathbf{w}_{j}^{\dagger}$, we introduce the Floquet eigenvectors $\mathbf{p}_{j}(\tau)=\mathcal{P}(\tau)\mathbf{v}_{j}$ and $\mathbf{q}_{j}^{\dagger}(\tau)=\mathbf{w}_{j}^{\dagger}\mathcal{P}^{-1}(\tau)$. As we show along the next sections, knowledge of these vectors is enough to derive the linearized quantum properties of the system. To this aim, it is also convenient to point out that they satisfy the initial value problems

[TABLE]

and the orthogonality conditions $\mathbf{q}_{j}^{\dagger}(\tau)\mathbf{p}_{l}(\tau)=\delta_{jl}\;\forall\tau$, as easily proven from their definition.

Let us now comment on the general properties of this eigensystem, which we prove in detail in SupMat . There always exists a null eigenvalue, say $\mu_{0}=0$, with related (right) Floquet eigenvector $\mathbf{p}_{0}(\tau)$. This property is a byproduct of the spontaneous temporal symmetry breaking generated by the limit cycle (Goldstone theorem). In the single-mode case, there is only one other eigenvalue, which is given by $\mu_{1}=\int_{0}^{T}\frac{d\tau}{T}\text{tr}\{\mathcal{L}(\tau)\}$, and has associated (left) Floquet eigenvector $\mathbf{q}_{1}(\tau)=(-\mathrm{i}\partial_{\tau}\bar{\beta},\mathrm{i}\partial_{\tau}\bar{\beta}^{*})^{T}\exp\left\{\int_{0}^{\tau}d\tau^{\prime}\text{tr}\{\mathcal{L}(\tau^{\prime})\}-\mu_{1}\tau\right\}$. This vector is the temporal counterpart of the linear or angular momentum found in previous works which deal with spatial symmetries CNBthesis .

Note that $\mathbf{p}_{0}(\tau)$ and $\mathbf{q}_{1}(\tau)$ are, respectively, the tangent and normal vectors of the limit cycle’s trajectory, see Fig. 1(a). We haven’t found explicit expressions of the other Floquet eigenvectors in terms of the $\bar{\beta}(\tau)$, but they can always be found numerically in an efficient fashion, as we do for Fig. 1(a).

**Diffusion of the temporal pattern. **As a first physical consequence of the properties above, we now show that $\theta$ is diffusing due to quantum noise, and hence quantum fluctuations smear off the classical periodic orbit.

In order to show this, we just need to apply $\mathbf{q}_{0}^{\dagger}(\tau)$ on (4), obtaining $\frac{d}{d\tau}(\mathbf{q}_{0}^{\dagger}\mathbf{b}+\theta)=\sqrt{\gamma}\mathbf{q}_{0}^{\dagger}(\tau)\mathbf{n}(\tau)$. Note that by taking $\theta$ as a variable in (3) we introduced a redundancy in the number of variables, which is now consistently removed by setting $\mathbf{q}_{0}^{\dagger}\mathbf{b}=0$ (in other words, introducing $\theta$ simply allowed us to track and give physical meaning to this part of the quantum fluctuations). The previous equation turns then into a diffusion equation for $\theta$, leading to a variance

[TABLE]

Note that the kernel is periodic, and therefore, the coarse-grained dynamics of $\theta$ corresponds to a diffusion process, with a variance increasing linearly with time, making $\theta$ fully undetermined asymptotically as shown in Fig. 1(c).

Steady state as a mixture of Gaussians. The above considerations imply that the steady state is formed by a balanced mixture of Gaussian states, one for each value of $\theta$. As we prove below, the Wigner functions of these Gaussian states CNBbook are given by

[TABLE]

where $\mathbf{r}=(x,p)^{T}$ is the coordinate vector in phase space, and the mean vector and covariance matrix are given by

[TABLE]

\mathcal{U}=\tiny{\left(\begin{array}[]{cc}1&1\\ -\mathrm{i}&\mathrm{i}\end{array}\right)} is the matrix that connects the complex representation of the bosonic mode to its real representation in phase space, and

[TABLE]

is a $T$-periodic function.

Let us now prove the expressions above. First, we introduce the quadrature vector $\mathbf{R}=\mathcal{U}(\beta,\beta^{+})^{T}/\sqrt{\gamma}$. Within the positive P representation the elements of the long-time mean vector $\bar{\mathbf{d}}$ and and covariance matrix $\bar{V}$** **are found as $\bar{d}_{m}(\tau)=\lim_{\tau\rightarrow\infty}\langle R_{m}(\tau)\rangle$ and $\bar{V}_{mn}(\tau)=\delta_{mn}+\lim_{\tau\rightarrow\infty}\langle\delta R_{m}(\tau)\delta R_{n}(\tau)\rangle$, where $\delta R_{m}=R_{m}-\langle R_{m}\rangle$ CNBbook . Next, note that the condition $\mathbf{q}_{0}^{\dagger}(\tau)\mathbf{b}(\tau)=0$ allows us to write the quantum fluctuations as $\mathbf{b}(\tau)=c_{1}(\tau)\mathbf{p}_{1}(\tau)$, where we define the projection $c_{1}(\tau)=\mathbf{q}_{1}^{\dagger}(\tau)\mathbf{b}(\tau)$. Using the expansion (3), we can then write the quadrature vector as $\sqrt{\gamma}\mathbf{R}(\tau)=\mathcal{U}[\bar{\beta}(\tau),\bar{\beta}^{+}(\tau)]^{T}+c_{1}(\tau)\mathcal{U}\mathbf{p}_{1}(\tau)$, whose stochastic properties are all then concentrated on $c_{1}(\tau)$. On the other hand, applying $\mathbf{q}_{1}^{\dagger}(\tau)$ on (4) we find $\dot{c}_{1}=\mu_{1}c_{1}+\sqrt{\gamma}\mathbf{q}_{1}^{\dagger}(\tau)\mathbf{n}(\tau)$, whose solution leads to the moments $\lim_{\tau\rightarrow\infty}\langle c_{1}(\tau)\rangle=0$ and $\lim_{\tau\rightarrow\infty}\langle c_{1}^{2}(\tau)\rangle=\gamma C(\tau)$, which provide the mean vector and covariance matrix in (10).

The steady state associated to the expansion (3) of the stochastic variables is then given by the balanced mixture

[TABLE]

In Fig. 2 we compare the Wigner function (12) with the one obtained by exact simulation CNBnumericsNotes of the master equation (1). We find very good agreement even for relatively large $\gamma$, where quantum fluctuations are still quite relevant, as can be appreciated.

This Wigner function has a very suggestive interpretation, see Fig. 2. First, (10a) tells us that the Gaussian states are centered along the points of the limit cycle’s trajectory, as expected. As for quantum fluctuations, note that the eigenvalues of the covariance matrix $\bar{V}(\theta)$ are 1 and $\mathrm{det}\{\bar{V}(\theta)\}$, which inform us about the variance along the principal axes of the uncertainty ellipse. It is easy to check that the directions of these principal axes follow the vectors $\mathcal{U}\mathbf{q}_{0}(\theta)$ and $\mathcal{U}\mathbf{p}_{1}(\theta)$ for the 1 and $\mathrm{det}\{\bar{V}(\theta)\}$ eigenvalues, respectively (see Fig. 2). Hence, the quadrature of the Gaussian state which goes in the direction of $\mathbf{q}_{0}(\theta)$ (Goldstone mode) carries vacuum fluctuations, which one can trace back to the condition $\mathbf{q}_{0}^{\dagger}(\theta)\mathbf{b}(\theta)=0$ that the method naturally demands. On the other hand, since in principle all physical covariance matrices satisfy $\mathrm{det}\{\bar{V}\}\geq 1$ (uncertainty principle) CNBbook , this seems to suggest that the quadrature going in the direction of $\mathbf{p}_{1}(\theta)$ carries fluctuations above the shot noise limit. While this is indeed the case for the VdP oscillator studied here, our experience with other nonlinear systems CNBthesis tells us that we could find $\mathrm{det}\{\bar{V}(\theta)\}<1$ (squeezing below shot noise) without violating the uncertainty principle. This is because the two quadratures of each Gaussian state are not conjugate variables, but they are both conjugate to the diffusing variable $\theta$ CNBthesis , which is completely undetermined in the steady state.

Conclusions. In this Letter we have introduced a linearization method capable of dealing with quantum nonlinear systems in the regime where they show spontaneous limit-cycle formation. The technique keeps the simplicity of standard linearization around stationary solutions. It requires finding the fundamental matrix of the Floquet method over a period of the cycle by solving a linear initial value problem with time-periodic coefficients. Only two equations are added with each mode that is introduced in the problem, giving the method a linear scaling with the size of the system that makes it suitable for complex driven-dissipative many-body problems such as optomechanical arrays Arrays1 ; Arrays2 ; Arrays3 ; Arrays4 ; ArraysWeiss . Moreover, the linearity of the equations should give efficient access also to dynamical objects such as multi-time correlation functions, which are of crucial relevance for experiments GardinerZoller ; Carmichael1 ; CarmichaelBook2 ; GardinerZollerI and the emergent field of quantum synchronization Lorch16 ; Lee13 ; Walter14 ; Weiss17 ; Walter15 ; Mari13 .