Semiclassical Phase Reduction Theory for Quantum Synchronization

TL;DR

This paper introduces a semiclassical phase reduction framework for quantum limit-cycle oscillators, enabling analysis of quantum synchronization using classical methods and providing insights into quantum-classical relations.

Contribution

It develops a general semiclassical phase reduction theory for quantum oscillators, allowing phase dynamics analysis and synchronization control.

Findings

The framework accurately reconstructs quantum system properties from phase equations.

Applied to quantum van der Pol oscillator, it reveals synchronization behaviors under various conditions.

Provides a bridge between quantum and classical synchronization analysis.

Abstract

We develop a general theoretical framework of semiclassical phase reduction for analyzing synchronization of quantum limit-cycle oscillators. The dynamics of quantum dissipative systems exhibiting limit-cycle oscillations are reduced to a simple, one-dimensional classical stochastic differential equation approximately describing the phase dynamics of the system under the semiclassical approximation. The density matrix and power spectrum of the original quantum system can be approximately reconstructed from the reduced phase equation. The developed framework enables us to analyze synchronization dynamics of quantum limit-cycle oscillators using the standard methods for classical limit-cycle oscillators in a quantitative way. As an example, we analyze synchronization of a quantum van der Pol oscillator under harmonic driving and squeezing, including the case that the squeezing is strong and the oscillation is asymmetric. The developed framework provides insights into the relation between quantum and classical synchronization and will facilitate systematic analysis and control of quantum nonlinear oscillators.

Equations

$$\dot{\rho }=-i[H+ϵ\tilde{H}(t),\rho ]+\sum _{m=1}^n\mathcal{D}[L_m]\rho ,$$

$$\displaystyle\frac{\partial P(\bm{\alpha},t)}{\partial t}=\Big{[}-\sum_{j=1}^{2}\partial_{j}\{A_{j}(\bm{\alpha})+\epsilon\tilde{A}_{j}(\bm{\alpha},t)\}+\frac{1}{2}\sum_{j=1}^{2}\sum_{k=1}^{2}\partial_{j}\partial_{k}\{\epsilon D_{jk}(\bm{\alpha})\}\Big{]}P(\bm{\alpha},t),$$

$$d\mathbf{\alpha }$$

$$d\mathbf{X}=\{\mathbf{F}(\mathbf{X})+ϵ\mathbf{q}(\mathbf{X},t)\}dt+\sqrt{ϵ}\mathbf{G}(\mathbf{X})d\mathbf{W},$$

$$d\varphi =\left\{\omega +ϵ\mathbf{Z}(\varphi )\cdot \mathbf{q}(\varphi ,t)+ϵg(\varphi )\right\}dt+\sqrt{ϵ}\{\mathbf{G}(\varphi {)}^T\mathbf{Z}(\varphi )\}\cdot d\mathbf{W},$$

$$g(\varphi )=\frac{1}{2}\text{mbox}Tr\left\{\mathbf{G}(\varphi {)}^T\mathbf{Y}(\varphi )\mathbf{G}(\varphi )\right\},$$

$$\omega \frac{d}{d\varphi }\mathbf{Z}(\varphi )=-{\mathbf{J}}^T(\varphi )\mathbf{Z}(\varphi ),\mathbf{Z}(\varphi )\cdot \frac{d{\mathbf{X}}_0(\varphi )}{d\varphi }=1,$$

$$\tilde{\omega }=\omega +\frac{ϵ}{2\pi }{\int }_0^{2\pi }g({\psi }^{\prime })d{\psi }^{\prime },$$

$$d\psi =\{\tilde{\omega }+ϵ\mathbf{Z}(\psi )\cdot \mathbf{q}(\psi ,t)\}dt+\sqrt{ϵ}h(\psi )dW,$$

$$\frac{\partial }{\partial t}P(\varphi ,t)=-\frac{\partial }{\partial \varphi }\left\{\omega +ϵ\mathbf{Z}(\varphi )\cdot \mathbf{q}(\varphi ,t)+ϵg(\varphi )\right\}P(\varphi ,t)+\frac{ϵ}{2}\frac{{\partial }^2}{\partial {\varphi }^2}h(\varphi {)}^{2}P(\varphi ,t).$$

$$\rho \approx {\int }_0^{2\pi }d\varphi P(\varphi )|{\alpha }_0(\varphi )⟩⟨{\alpha }_0(\varphi )|,$$

$$\dot{\rho }=-i\left[-\Delta a^{†}a+iE(a-a^{†})+i\eta (a^2{e}^{-i\theta }-a^{†2}{e}^{i\theta }),\rho \right]+{\gamma }_1\mathcal{D}[a^{†}]\rho +{\gamma }_2\mathcal{D}[a^2]\rho ,$$

$$d\mathbf{X}=\left(\begin{array}{c}\frac{1}{2}{x}^{\prime }-{\Delta }^{\prime }{p}^{\prime }-{\gamma }_2^{\prime }{x}^{\prime }(x^{\prime 2}+p^{\prime 2})-ϵE^{\prime }-2\delta {\eta }^{\prime }(x^{\prime }\cos \theta +p^{\prime }\sin \theta )\\ \frac{1}{2}{p}^{\prime }+{\Delta }^{\prime }{x}^{\prime }-{\gamma }_2^{\prime }{p}^{\prime }(x^{\prime 2}+p^{\prime 2})+2\delta {\eta }^{\prime }(p^{\prime }\cos \theta -x^{\prime }\sin \theta )\end{array}\right)d{t}^{\prime }+\sqrt{ϵ}\mathbf{G}(\mathbf{X})d{\mathbf{W}}^{\prime },$$

$$\mathbf{G}(\mathbf{X})$$

$$H=-{\Delta }^{\prime }{a}^{\prime †}{a}^{\prime },ϵ\tilde{H}=ϵ\{i{E}^{\prime }(a^{\prime }-a^{\prime †})+i{\eta }^{\prime }(a^{\prime 2}{e}^{-i\theta }-a^{\prime †2}{e}^{i\theta })\}.$$

$${\mathbf{X}}_0(\varphi )=\sqrt{\frac{1}{2{\gamma }_2^{\prime }}}\left(\begin{array}{c}\cos \varphi \\ \sin \varphi \end{array}\right)$$

$$\mathbf{Z}(\varphi )=\sqrt{2{\gamma }_2^{\prime }}\left(\begin{array}{c}-\sin \varphi \\ \cos \varphi \end{array}\right),\mathbf{Y}(\varphi )=2{\gamma }_2^{\prime }\left(\begin{array}{cc}\sin 2\varphi & -\cos 2\varphi \\ -\cos 2\varphi & -\sin 2\varphi \end{array}\right).$$

$$d\varphi$$

$$\rho \approx {\int }_0^{2\pi }d\varphi P(\varphi )|\sqrt{\frac{{\gamma }_1}{2{\gamma }_2}}{e}^{i\varphi }⟩⟨\sqrt{\frac{{\gamma }_1}{2{\gamma }_2}}{e}^{i\varphi }|.$$

$$H=-{\Delta }^{\prime }{a}^{\prime †}{a}^{\prime }+i{\eta }^{\prime }(a^{\prime 2}{e}^{-i\theta }-a^{\prime †2}{e}^{i\theta }),ϵ\tilde{H}=ϵi{E}^{\prime }(a^{\prime }-a^{\prime †}).$$

$$S_{qm}(\omega )$$

$$R_{qm}(\tau )$$

$$S_{sc}(\omega )$$

$$R_{sc}(\tau )$$

$$\mathbf{D}(\mathbf{\alpha })=\mathbf{\beta }(\mathbf{\alpha })\mathbf{\beta }(\mathbf{\alpha }{)}^T=\left(\begin{array}{cc}{D}_{11}(\mathbf{\alpha })& D_{12}(\mathbf{\alpha })\\ D_{21}(\mathbf{\alpha })& D_{22}(\mathbf{\alpha })\end{array}\right)\in {\mathbb{C}}^{2\times 2},$$

$$\frac{\partial }{\partial t}P(\mathbf{X},t)=\left[-\frac{\partial }{\partial \mathbf{X}}\{\mathbf{F}(\mathbf{X})+ϵ\mathbf{q}(\mathbf{X},t)\}+\frac{1}{2}\frac{{\partial }^2}{\partial {\mathbf{X}}^2}\mathbf{D}(\mathbf{X})\right]P(\mathbf{X},t),$$

$$\frac{\partial }{\partial \alpha }=\frac{1}{2}\left(\frac{\partial }{\partial x}-i\frac{\partial }{\partial p}\right),\frac{\partial }{\partial {\alpha }^{\ast }}=\frac{1}{2}\left(\frac{\partial }{\partial x}+i\frac{\partial }{\partial p}\right).$$

$$\mathbf{D}(\mathbf{X})$$

$$=\frac{1}{2}\left(\begin{array}{cc}\text{mbox}Re{D}_{11}(\mathbf{\alpha })+D_{12}(\mathbf{\alpha })& \text{mbox}Im{D}_{11}(\mathbf{\alpha })\\ \text{mbox}Im{D}_{11}(\mathbf{\alpha })& -\text{mbox}Re{D}_{11}(\mathbf{\alpha })+D_{12}(\mathbf{\alpha })\end{array}\right)\in {\mathbb{R}}^{2\times 2}$$

$$\mathbf{D}(\mathbf{\alpha })$$

$${\lambda }_{\pm }(\mathbf{X})$$

$${\mathbf{v}}_{+}(\mathbf{X})$$

Full text

I Introduction

Spontaneous rhythmic oscillations and synchronization arise in various science and technology fields, such as laser oscillations, electronic oscillators, and spiking neurons Winfree (2001); Kuramoto (1984); Ermentrout and Terman (2010); Pikovsky et al. (2001); Glass and Mackey (1988); Strogatz (1994). Various nonlinear dissipative systems exhibiting rhythmic dynamics can be modeled as limit-cycle oscillators. A standard theoretical framework for analyzing limit-cycle oscillators in classical dissipative systems is the phase reduction theory Winfree (2001); Kuramoto (1984); Nakao (2016); Ermentrout and Terman (2010); Ermentrout (1996); Brown et al. (2004). By using this framework, we can systematically reduce multi-dimensional nonlinear dynamical equations describing weakly-perturbed limit-cycle oscillators to a one-dimensional phase equation that approximately describes the oscillator dynamics. The simple semi-linear form of the phase equation, characterized only by the natural frequency and phase sensitivity function (PSF) of the oscillator, facilitates detailed theoretical analysis of the oscillator dynamics.

The phase reduction theory has been successfully used to analyze universal properties of limit-cycle oscillators in a systematic way, such as synchronization of oscillators with periodic forcing and mutual synchronization of coupled oscillators Winfree (2001); Kuramoto (1984); Ermentrout and Terman (2010); Pikovsky et al. (2001); Glass and Mackey (1988); Strogatz (1994). It has been essential in the understanding of synchronization phenomena in classical rhythmic systems, for example, the collective synchronization transition of a population of oscillators and oscillatory pattern dynamics in spatially extended chemical or biological systems Winfree (2001); Kuramoto (1984). Recently, generalizations of the phase reduction theory to non-conventional physical systems, such as time-delayed oscillators Kotani et al. (2012); Novičenko and Pyragas (2012), piecewise-smooth oscillators Shirasaka et al. (2017a), collectively oscillating networks Nakao et al. (2018), and rhythmic spatiotemporal patterns Kawamura and Nakao (2013); Nakao et al. (2014), have also been discussed.

Recent progress in experimental studies has revealed that synchronization can take place in coupled nonlinear oscillators with intrinsically quantum-mechanical origins, such as micro and nanomechanical oscillators Shim et al. (2007); Zhang et al. (2012); *zhang2015synchronization; Bagheri et al. (2013); Matheny et al. (2014, 2019), spin torque oscillators Kaka et al. (2005), and cooled atomic ensembles Weiner et al. (2017); Heimonen et al. (2018). Moreover, theoretical studies have been performed on the synchronization of nonlinear oscillators which explicitly show quantum signatures Ludwig and Marquardt (2013); Weiss et al. (2016); Amitai et al. (2017); Xu et al. (2014); Xu and Holland (2015); Lee and Sadeghpour (2013); Lee et al. (2014); Hush et al. (2015); Roulet and Bruder (2018a); *roulet2018quantum; *koppenhofer2019optimal; Nigg (2018); Lee and Sadeghpour (2013); Walter et al. (2014); Sonar et al. (2018); Lee et al. (2014); Walter et al. (2015); Lörch et al. (2016); Ishibashi and Kanamoto (2017); Amitai et al. (2018); Navarrete-Benlloch et al. (2017); Weiss et al. (2017); Hriscu and Nazarov (2013); Hamerly and Mabuchi (2015); Lee and Cross (2013); Mari et al. (2013); Ameri et al. (2015); Witthaut et al. (2017); Davis-Tilley et al. (2018); Lörch et al. (2017); de Mendoza et al. (2014), such as optomechanical oscillators Ludwig and Marquardt (2013); Weiss et al. (2016); Amitai et al. (2017), cooled atomic ensembles Xu et al. (2014); Xu and Holland (2015), trapped ions Lee and Sadeghpour (2013); Lee et al. (2014); Hush et al. (2015), spins Roulet and Bruder (2018a); *roulet2018quantum, and superconducting circuits Nigg (2018). In particular, a number of studies have analyzed the quantum van der Pol (vdP) oscillator Lee and Sadeghpour (2013), which is a typical model of quantum self-sustained oscillators, for example, synchronization of a quantum vdP oscillator by harmonic driving Walter et al. (2014); Amitai et al. (2017) or squeezing Sonar et al. (2018), mutual synchronization of coupled quantum vdP oscillators Lee et al. (2014); Walter et al. (2015), and quantum fluctuations around oscillating and locked states of a quantum vdP oscillator Navarrete-Benlloch et al. (2017); Weiss et al. (2017).

In addition to its fundamental importance as a novel physical phenomenon where nonlinear and quantum phenomena have combined effect, quantum synchronization may also be useful in developing metrological applications, such as the improvement of the measurement accuracy in the Ramsey spectroscopy for atomic clocks Xu and Holland (2015) and the precise measurement of the resistance standard with a superconducting device Hriscu and Nazarov (2013); an application of the limit-cycle oscillation to analog memory in a quantum optical device Hamerly and Mabuchi (2015) has also been considered.

Considering the importance of phase reduction for analyzing synchronization of classical nonlinear oscillators, we aim to develop a phase reduction theory also for quantum nonlinear oscillators. In the analysis of quantum synchronization, phase-space approaches using the quasiprobability distributions of quantum systems are commonly employed. In a pioneering study, Hamerly and Mabuchi Hamerly and Mabuchi (2015) derived a phase equation from the stochastic differential equation (SDE) describing a truncated Wigner function of a quantum limit-cycling system in a free-carrier cavity. However, it is not fully consistent with the classical phase reduction theory, because the notions of the asymptotic phase and PSF, which are essential in the classical theory, are not introduced. Consequently, the limit cycle needs to be approximately symmetric for the analysis of synchronization with periodic forcing Hamerly and Mabuchi (2015). Similar phenomenological phase equations, where the phase simply represents the geometric angle of a circular limit cycle, have also been used in several studies on quantum synchronization Ludwig and Marquardt (2013); Xu and Holland (2015); Witthaut et al. (2017); Amitai et al. (2017); however, a systematic phase reduction theory has not been established so far.

In this study, we formulate a general framework of the phase reduction theory for quantum synchronization under the semiclassical approximation. We derive a linearized multi-dimensional semiclassical SDE from a general master equation that describes weakly-perturbed quantum dissipative systems with a single degree of freedom exhibiting stable nonlinear oscillations, and subsequently reduce it to an approximate one-dimensional classical SDE for the phase variable of the system (see Fig. 1). The derived phase equation has a simple form, characterized by the natural frequency, PSF, and Hessian matrix of the limit cycle in the classical limit, and a noise term arising from quantum fluctuations around the limit cycle. The quantum-mechanical density matrix and power spectrum of the original system can be approximately reconstructed from the reduced phase equation.

On the basis of the reduced phase equation, synchronization dynamics of quantum nonlinear oscillators can be analyzed in detail by using standard techniques for classical nonlinear oscillators Winfree (2001); Kuramoto (1984); Nakao (2016); Ermentrout and Terman (2010); Ermentrout (1996); Brown et al. (2004). As an example, we analyze synchronization of a quantum vdP oscillator under harmonic driving and squeezing. In particular, we consider the case with strong squeezing, where the oscillation is asymmetric and the analytical solution is not available. It is shown that, even in such cases, we can numerically calculate the necessary quantities in the classical limit and use them to analyze the synchronization dynamics of the original quantum system, provided that the quantum noise and the perturbations given to the oscillator are sufficiently weak.

The rest of this paper is organized as follows; In Sec. II, the derivation of the approximate phase equation for a quantum limit-cycle oscillator subjected to weak perturbations is given. In Sec. III, we analyze a quantum vdP oscillator with harmonic driving and squeezing using the derived phase equation. Section  IV gives concluding remarks, and Appendices provide detailed derivations of the equations and discussions.

II Theory

II.1 Stochastic differential equation for phase-space variables

We consider quantum dissipative systems with a single degree of freedom interacting with linear and nonlinear reservoirs, which has a stable limit-cycle solution in the classical limit and is driven by weak perturbations. Under the assumption that correlation times of the reservoirs are significantly shorter than the time scale of the main system, a Markovian approximation of the reservoirs can be employed and the evolution of the system can be described by a quantum master equation Carmichael (2007); Gardiner and Haken (1991),

[TABLE]

where $\rho$ is a density matrix representing the system state, $H$ is a system Hamiltonian, $\epsilon{\tilde{H}}(t)$ is a time-dependent Hamiltonian representing weak external perturbations applied to the system ($0<\epsilon\ll 1$), $n$ is the number of reservoirs, $L_{m}$ is the coupling operator between the system and the $m$th reservoir $(m=1,\ldots,n)$, and $\mathcal{D}[L]\rho=L\rho L^{{\dagger}}-(\rho L^{{\dagger}}L-L^{{\dagger}}L\rho)/2$ denotes the Lindblad form. We consider a physical condition where the effects of the quantum noise and external perturbations are sufficiently weak and of the same order, and perturbatively analyze their effect on the semiclassical dynamics of the system.

First, we transform Eq. (1) into a multi-dimensional SDE by introducing a phase-space quasiprobability distribution, such as the P, Q, or Wigner representation Carmichael (2007); Gardiner and Haken (1991). In this paper, we use the P representation, because the density matrix and spectrum can be reconstructed using a simple and natural approximation. In the P representation, the density matrix $\rho$ is represented as $\rho=\int P({\bm{\alpha}})|\alpha\rangle\langle\alpha|d{\bm{\alpha}}$, where $|\alpha\rangle$ is a coherent state specified by a complex value $\alpha\in\mathbb{C}$, or equivalently by a two-dimensional complex vector $\bm{\alpha}=(\alpha,\alpha^{*})^{T}\in\mathbb{C}^{2\times 1}$, $P({\bm{\alpha}})$ is a quasiprobability distribution of ${\bm{\alpha}}$, $d{\bm{\alpha}}=d\alpha d\alpha^{*}$, the integral is taken over the entire space spanned by ${\bm{\alpha}}$, and * indicates complex conjugate.

The Fokker-Planck equation (FPE) equivalent to Eq. (1) can be written as

[TABLE]

where $A_{j}({\bm{\alpha}})$ and $\tilde{A}_{j}({\bm{\alpha}},t)$ are the $j$th components of complex vectors $\bm{A}(\bm{\alpha})=(A_{1}(\bm{\alpha})$, $A^{*}_{1}(\bm{\alpha}))^{T}\in{\mathbb{C}}^{2\times 1}$ and $\bm{\tilde{A}}(\bm{\alpha},t)=(\tilde{A}_{1}(\bm{\alpha},t),\tilde{A}^{*}_{1}(\bm{\alpha},t))^{T}\in{\mathbb{C}}^{2\times 1}$ representing the system dynamics and perturbations, respectively, $\epsilon D_{jk}({\bm{\alpha}})$ is the $(j,k)$-th component of the symmetric diffusion matrix $\epsilon{\bm{D}}({\bm{\alpha}})\in{\mathbb{C}}^{2\times 2}$ representing quantum fluctuations, and the complex partial derivatives are defined as $\partial_{1}=\partial/\partial\alpha$ and $\partial_{2}=\partial/\partial\alpha^{*}$ (note that $A_{2}(\bm{\alpha})=A^{*}_{1}(\bm{\alpha})$ and $\tilde{A}_{2}(\bm{\alpha},t)=\tilde{A}^{*}_{1}(\bm{\alpha},t)$).

The drift term $\bm{A}(\bm{\alpha})$ consists of terms arising from the system Hamiltonian $H$ and the dissipation $\{L_{m}\}$, $\epsilon\tilde{\bm{A}}(\bm{\alpha},t)$ represents the small terms arising from the perturbation Hamiltonian $\epsilon{\tilde{H}}(t)$, and the diffusion matrix $\epsilon{\bm{D}}({\bm{\alpha}})$ represents the intensity of the small quantum noise, generally arising from all terms of $H$, $\epsilon{\tilde{H}}(t)$, and $\{L_{m}\}$. These terms can be explicitly calculated from the master equation in Eq. (1) by using the standard calculus for phase-space representation when $H$, $\epsilon{\tilde{H}}(t)$, and $\{L_{m}\}$ are given Carmichael (2007); Gardiner and Haken (1991). The external perturbation $\epsilon\tilde{\bm{A}}({\bm{\alpha}},t)$ and the diffusion matrix $\epsilon{\bm{D}}({\bm{\alpha}})$ are assumed to be of the same order, $\mathcal{O}(\epsilon)$.

By introducing an appropriate complex matrix $\sqrt{\epsilon}\bm{\beta}(\bm{\alpha})\in\mathbb{C}^{2\times 2}$ (see Appendix A for the explicit form), the diffusion matrix $\epsilon{\bm{D}}({\bm{\alpha}})$ can be represented as $\epsilon\bm{D}(\bm{\alpha})=\sqrt{\epsilon}\bm{\beta}(\bm{\alpha})(\sqrt{\epsilon}\bm{\beta}(\bm{\alpha}))^{T}$ and the Ito SDE corresponding to Eq. (2) for the phase-space variable ${\bm{\alpha}}(t)$ is given by

[TABLE]

where $\bm{W}(t)=(W_{1}(t),W_{2}(t))^{T}\in\mathbb{R}^{2\times 1}$ represents a vector of independent Wiener processes $W_{i}(t)~{}(i=1,\ldots,2)$ satisfying $\langle{dW_{i}dW_{j}}\rangle=\delta_{ij}dt$.

It should be noted that diffusion matrix of certain quantum systems in the P representation becomes negative definite for certain ${\bm{\alpha}}$ Carmichael (2007); Gardiner and Haken (1991). For such systems, we need to employ, for example, the positive P representation with two additional nonclassical variables in place of the P representation, as used by Navarrete-Benlloch et al. Navarrete-Benlloch et al. (2017) in the Floquet analysis of quantum oscillations. In this study, to present the fundamental idea of the semiclassical phase reduction in its simplest form, we only consider the case for which the diffusion matrix is always positive semidefinite along the limit cycle and formulate the phase reduction theory in the two-dimensional phase space of classical variables.

II.2 Derivation of the phase equation

Our aim is to derive an approximate one-dimensional SDE for the phase variable of the system from the SDE in Eq. (3) in the $P$ representation. To this end, we define a real vector ${\bm{X}}=(x,p)^{T}=(\mbox{Re}\ \alpha,\mbox{Im}\ \alpha)^{T}\in{\mathbb{R}}^{2\times 1}$ from the complex vector ${\bm{\alpha}}$. The real-valued expression of Eq. (3) for ${\bm{X}}(t)$ is then given by an Ito SDE,

[TABLE]

where $\bm{F}({\bm{X}})\in{\mathbb{R}}^{2\times 1}$, ${\bm{q}}({\bm{X}},t)\in{\mathbb{R}}^{2\times 1}$, and ${\bm{G}}({\bm{X}})\in{\mathbb{R}}^{2\times 2}$ are real-valued equivalent representations of the system dynamics ${\bm{A}}({\bm{\alpha}})\in{\mathbb{C}}^{2\times 1}$, perturbation $\tilde{\bm{A}}({\bm{\alpha}},t)\in{\mathbb{C}}^{2\times 1}$, and noise intensity ${\bm{\beta}}({\bm{\alpha}})\in{\mathbb{C}}^{2\times 2}$ in Eq. (3), respectively.

We assume that the system in the classical limit without perturbation and quantum noise, $\dot{{\bm{X}}}={\bm{F}}({\bm{X}})$, has an exponentially stable limit-cycle solution ${\bm{X}}_{0}(t)=(x_{0}(t),p_{0}(t))^{T}={\bm{X}}_{0}(t+T)$ with a natural period $T$ and frequency $\omega=2\pi/T$. In the same way as the phase reduction for classical limit cycles Winfree (2001); Kuramoto (1984); Nakao (2016); Ermentrout and Terman (2010); Ermentrout (1996); Brown et al. (2004), we can introduce an asymptotic phase function $\Phi({\bm{X}}):{B}\subset{\mathbb{R}}^{2}\to[0,2\pi)$ such that $\nabla\Phi({\bm{X}})\cdot{\bm{F}}({\bm{X}})=\omega$ is satisfied for all system states ${\bm{X}}$ in the basin ${B}$ of the limit cycle in the classical limit, where $\nabla\Phi({\bm{X}})\in{\mathbb{R}}^{2\times 1}$ is the gradient of $\Phi({\bm{X}})$. Using this phase function, we define the phase of a system state ${\bm{X}}\in B$ as $\phi=\Phi({\bm{X}})$. It then follows that $\dot{\phi}=\dot{\Phi}({\bm{X}})={\bm{F}}({\bm{X}})\cdot\nabla\Phi({\bm{X}})=\omega$, i.e., $\phi$ always increases at a constant frequency $\omega$ with the evolution of ${\bm{X}}$. In the following formulation, we represent the system state ${\bm{X}}$ on the limit cycle as ${\bm{X}}_{0}(\phi)=(x_{0}(\phi),p_{0}(\phi))^{T}$ as a function of the phase $\phi$ rather than the time $t$. In this representation, ${\bm{X}}_{0}(\phi)$ is a $2\pi$-periodic function of $\phi$, ${\bm{X}}_{0}(\phi)={\bm{X}}_{0}(\phi+2\pi)$. Note that an identity $\Phi({\bm{X}}_{0}(\phi))=\phi$ is satisfied by the definition of $\Phi({\bm{X}})$.

When the noise and perturbations are sufficiently weak and the deviation of the state ${\bm{X}}$ from the limit cycle is small, we can approximate ${\bm{X}}(t)$ by a state ${\bm{X}}_{0}(\phi(t))$ on the limit cycle as ${\bm{X}}(t)\approx{\bm{X}}_{0}(\phi(t))$ and derive a SDE for the phase in the lowest order approximation by using the Ito formula as (see Appendix B for details)

[TABLE]

where the drift term is correct up to $\mathcal{O}(\epsilon)$ and the noise intensity is correct up to $\mathcal{O}(\sqrt{\epsilon})$. Here, the inner product between two vectors $\bm{a}=(a_{0},a_{2},\cdots,a_{N-1})^{T}\in{\mathbb{R}}^{N\times 1}$ and $\bm{b}=(b_{0},b_{2},\cdots,b_{N-1})^{T}\in{\mathbb{R}}^{N\times 1}$ is defined as $\bm{a}\cdot\bm{b}=\sum_{i=0}^{N-1}a_{i}b_{i}$.

In the above phase equation, the gradient $\nabla\Phi$ of $\Phi({\bm{X}})$ at ${\bm{X}}$ is approximately evaluated at ${\bm{X}}(\phi)$ on the limit cycle and is denoted as $\bm{Z}(\phi)=\nabla\Phi|_{{\bm{X}}={\bm{X}}_{0}(\phi)}\in{\mathbb{R}}^{2\times 1}$. We call this ${\bm{Z}}(\phi)$ the phase sensitivity function (PSF) of the limit cycle, which characterizes the linear response property of the oscillator phase to given perturbations Kuramoto (1984); Nakao (2016). Similarly, the perturbation and noise intensity can also be evaluated approximately at ${\bm{X}}={\bm{X}}_{0}(\phi)$ on the limit cycle and they are denoted as ${\bm{q}}(\phi,t)={\bm{q}}({\bm{X}}_{0}(\phi),t)$ and $\bm{G}(\phi)=\bm{G}({\bm{X}}_{0}(\phi))$, respectively. The additional function $g(\phi)$ in the drift term in Eq. (5) arises from the change of the variables and is given by

[TABLE]

where ${\bm{Y}}(\phi)=\nabla^{T}\nabla\Phi|_{{\bm{X}}={{\bm{X}}_{0}(\phi)}}$ is a Hessian matrix of the phase function $\Phi({\bm{X}})$ also evaluated at ${\bm{X}}={\bm{X}}_{0}(\phi)$ on the limit cycle. All these functions are $2\pi$-periodic, as they are functions of ${\bm{X}}_{0}(\phi)$.

It is well known in the classical phase reduction theory that the PSF can be obtained as a $2\pi$-periodic solution to the following adjoint equation and an additional normalization condition Ermentrout (1996); Brown et al. (2004); Nakao (2016):

[TABLE]

respectively, where ${\bm{J}}(\phi)={\bm{J}}({\bm{X}}_{0}(\phi))\in{\mathbb{R}}^{2\times 2}$ is a Jacobian matrix of ${\bm{F}}({\bm{X}})$ at ${\bm{X}}={\bm{X}}_{0}(\phi)$ on the limit cycle. It is also known that the Hessian matrix ${\bm{Y}}(\phi)$ on the limit cycle can be calculated as a $2\pi$-periodic solution of an adjoint-type equation Suvak and Demir (2010); Takeshita and Feres (2010) with an appropriate constraint. These equations for ${\bm{Y}}(\phi)$ are detailed in the Appendix B. In the numerical calculations, ${\bm{Z}}(\phi)$ can easily be obtained by the backward integration of the adjoint equation with occasional normalization as proposed by Ermentrout Ermentrout and Terman (2010), and then the Hessian ${\bm{Y}}(\phi)$ can be obtained by a shooting method Suvak and Demir (2010).

Because of the additional term $g(\phi)$ in Eq. (10), the effective frequency $\tilde{\omega}=\langle d\phi\rangle/dt$ of the oscillator in the absence of the perturbation ${\bm{q}}(\phi,t)$ is given by

[TABLE]

which is slightly different from the natural frequency of the oscillator $\omega$ in the classical limit. Though not used in the present study, we can further introduce a new phase variable $\psi$ that is only slightly different from $\phi$ by a near-identity transform as $\phi=\psi+\epsilon n(\psi)$, where $n(\psi)$ is a $2\pi$-periodic function with $n(0)=0$, and eliminate the additional function $g(\phi)$ in Eq. (5) by renormalizing it into the frequency term. The new phase $\psi$ then obeys a simpler SDE of the form

[TABLE]

where $h(\psi)=\sqrt{\sum_{i=1}^{2}\left\{{\bm{G}}(\psi)^{T}{\bm{Z}}(\psi)\right\}_{i}^{2}}$ and $W(t)$ is a one-dimensional Wiener process. As before, the drift term is correct up to $\mathcal{O}(\epsilon)$ and the noise intensity is correct up to $\mathcal{O}(\sqrt{\epsilon})$. See Appendix C for the details. In this study, we use the original phase equation in Eq. (5) for numerical simulations and verify its validity. We also note here that the phase equation derived in Ref. Hamerly and Mabuchi (2015) does not contain a term with the Hessian matrix, because the order of the noise intensity is implicitly assumed to be $\mathcal{O}(\epsilon)$ in Hamerly and Mabuchi (2015).

From the reduced SDE in Eq. (5), we can derive a corresponding FPE describing the probability density function $P(\phi,t)$ of the phase variable $\phi$ as

[TABLE]

Using this FPE, we can obtain the stationary distribution and transition probability of the phase variable $\phi$ and use them to reconstruct the density matrix and power spectrum.

II.3 Reconstruction of the density matrix

From the reduced phase equation, we can approximately reconstruct the quantum state as follows. Using the phase variable $\phi$, the oscillator state in the classical limit can be approximated as ${\bm{X}}\approx{\bm{X}}_{0}(\phi)=(x_{0}(\phi),p_{0}(\phi))^{T}$, or ${\bm{\alpha}}\approx{\bm{\alpha}}_{0}(\phi)=(\alpha_{0}(\phi),\alpha_{0}(\phi)^{*})^{T}=(x_{0}(\phi)+ip(\phi),x_{0}(\phi)-ip(\phi))^{T}$ in the original complex representation. Therefore, the quantum state at phase $\phi$ is approximately described as $\left|\alpha_{0}(\phi)\right\rangle$ and the density matrix $\rho$ is approximately represented by using the probability density function $P(\phi)$ of the phase variable $\phi$, obtained from the SDE in Eq. (5) or FPE in Eq. (10), as

[TABLE]

which is simply a mixture of coherent states weighted by the distribution of the phase on the classical limit cycle. Thus, we can approximately reconstruct the density matrix of the original quantum oscillator from the classical SDE for the phase variable $\phi$, which is characterized by the natural frequency $\omega$, PSF ${\bm{Z}}(\phi)$, Hessian matrix ${\bm{Y}}(\phi)$, and noise intensity ${\bm{G}}(\phi)$ that represents quantum fluctuations around the limit cycle.

The derivation of the phase equation in Eq. (5) from the original quantum-mechanical master equation in Eq. (1) and reconstruction of the quantum-mechanical density matrix from the approximate phase equation, Eq. (11), are the main result of the present work. A schematic diagram of the proposed method is illustrated in Fig. 1. The reduced phase equation is essentially the same as that for the classical limit-cycle oscillator driven by noise, and synchronization dynamics of the weakly perturbed quantum nonlinear oscillator in the semiclassical regime can be analyzed on the basis of the reduced phase equation by using the standard methods for the classical limit-cycle oscillator.

III Examples

III.1 Quantum van der Pol oscillator with harmonic driving and squeezing

As an example, we consider a quantum vdP oscillator subjected to harmonic driving and squeezing. We assume that the harmonic driving is sufficiently weak and treat it as a perturbation. As for the squeezing, we consider two cases; (i) the squeezing is sufficiently weak and can also be treated as a perturbation, and (ii) the squeezing is relatively strong and cannot be treated as a perturbation.

We denote by $\omega_{0}$, $\omega_{d}$, and $\omega_{sq}$ the frequencies of the oscillator, harmonic driving, and pump beam of squeezing, respectively. We consider the case where the squeezing is generated by a degenerate parametric amplifier and assume $\omega_{sq}=2\omega_{d}$ Gardiner and Haken (1991). In the rotating coordinate frame of frequency $\omega_{d}$, the master equation is given by Walter et al. (2014); Sonar et al. (2018)

[TABLE]

where $\Delta=\omega_{d}-\omega_{0}$ is the frequency detuning of the harmonic driving from the oscillator, $E$ is the intensity of the harmonic driving, $\eta e^{i\theta}$ is the squeezing parameter, $\gamma_{1}$ and $\gamma_{2}$ are the decay rates for negative damping and nonlinear damping, respectively, and the Planck constant is set as $\hbar=1$. The harmonic driving is represented by a constant $E$, because a coordinate frame rotating with the driving frequency $\omega_{d}$ is used.

We assume that $\gamma_{2}$ is sufficiently small and of $\mathcal{O}(\epsilon)$, for which the semiclassical approximation is valid, and represent $\gamma_{2}$ as $\gamma_{2}=\epsilon\gamma_{1}\gamma_{2}{{}^{\prime}}$ using a dimensionless parameter $\gamma_{2}{{}^{\prime}}$ of $\mathcal{O}(1)$. In this setting, the size of the stable limit-cycle solution in Eq. (12) in the classical limit is $\mathcal{O}(1/\sqrt{\epsilon})$, while we have implicitly assumed it to be $\mathcal{O}(1)$ in the derivation of Eq. (5). Therefore, we introduce a rescaled annihilation operator $a^{\prime}$ and the corresponding classical variable $\alpha^{\prime}$ (${\bm{\alpha}}^{\prime}=(\alpha^{\prime},\alpha^{\prime*})$ in the vector representation) as $a^{\prime}|\alpha^{\prime}\rangle=\sqrt{\epsilon}a|\sqrt{\epsilon}\alpha\rangle$, and represent the parameters as $\Delta=\gamma_{1}\Delta{{}^{\prime}},E=\sqrt{\epsilon}\gamma_{1}E{{}^{\prime}},\eta=\delta\gamma_{1}\eta{{}^{\prime}}$, where $\Delta^{\prime},E^{\prime}$, and $\eta^{\prime}$ are dimensionless parameters of $\mathcal{O}(1)$. By this rescaling, the size of the limit cycle becomes $\mathcal{O}(1)$ and the parameter $\delta$ determines the relative intensity of the squeezing.

The real-valued representation $\bm{X}=(x^{\prime},p^{\prime})^{T}=(\mbox{Re}~{}\alpha^{\prime},\mbox{Im}~{}\alpha^{\prime})^{T}$ of Eq. (4) after rescaling is then obtained as

[TABLE]

where $dt^{\prime}=\gamma_{1}dt$ and $d\bm{W}^{\prime}=\sqrt{\gamma_{1}}d\bm{W}$. The noise intensity matrix is explicitly given by

[TABLE]

with $R^{\prime}_{1}e^{i\chi^{\prime}_{1}}=-(\gamma^{\prime}_{2}(x^{\prime}+ip^{\prime})^{2}+2\delta\eta^{\prime}e^{i\theta})$. Further details of the derivation can be found in the Appendix D.

III.2 Weak squeezing

First, we consider the case of weak squeezing with $\delta=\epsilon$. The rescaled system and perturbation Hamiltonians are given by

[TABLE]

For this system, we obtain ${\bm{F}}({\bm{X}})=({x}^{\prime}/2-\Delta^{\prime}{p}^{\prime}-\gamma^{\prime}_{2}{x}^{\prime}({x}^{\prime 2}+{p}^{\prime 2}),\ {p}^{\prime}/2+\Delta^{\prime}{x}^{\prime}-\gamma^{\prime}_{2}{p}^{\prime}({x}^{\prime 2}+{p}^{\prime 2}))^{T}$. The perturbation is represented by ${\bm{q}}({\bm{X}},t)=(-E^{\prime}-2\eta^{\prime}(x^{\prime}\cos\theta+p^{\prime}\sin\theta),\ 2\eta^{\prime}(p^{\prime}\cos\theta-x^{\prime}\sin\theta))^{T}$. Note that the vector field ${\bm{F}}({\bm{X}})$ in this case is simply a normal form of the supercritical Hopf bifurcation. A classical nonlinear oscillator described by this ${\bm{F}}({\bm{X}})$ is known as the Stuart-Landau (SL) oscillator Kuramoto (1984) (which is different from the classical vdP oscillator) and it is analytically solvable.

The stable limit cycle of the SL oscillator is given by

[TABLE]

as a function of phase $\phi=\omega t$, where the frequency is given by $\omega=\Delta^{\prime}$. The basin $B$ of this limit cycle is the whole $(x^{\prime},p^{\prime})$-plane except $(0,0)$. The phase function $\Phi(\bm{X})$ of this limit cycle can be expressed as $\Phi(x^{\prime},p^{\prime})=\tan^{-1}(p^{\prime}/x^{\prime})$ Nakao (2016), which gives $\dot{\phi}=\dot{\Phi}(x^{\prime},p^{\prime})=\omega$. The PSF ${\bm{Z}}(\phi)$ and Hessian matrix ${\bm{Y}}(\phi)$ can be obtained by calculating the gradients of the phase function $\Phi({\bm{X}})$ at ${\bm{X}}={\bm{X}}_{0}(\phi)$ on the limit cycle as

[TABLE]

In this case, the additional term $g(\phi)$ in Eq. (5) and therefore the frequency shift in Eq. (8) vanishes, i.e., $\tilde{\omega}=\omega$. The $\mathcal{O}(\epsilon\sqrt{\epsilon})$ terms in the noise intensity ${\bm{G}}(\phi)$ given by Eq. (14) are neglected.

From these results, the phase equation in Eq. (5) for the quantum vdP oscillator driven by weak harmonic driving and squeezing is explicitly given by

[TABLE]

in the lowest-order approximation, where $dW^{\prime}=\sqrt{\gamma_{1}}dW$. Using the probability density function $P(\phi)$ of the phase $\phi$ described by the FPE (10) corresponding to Eq. (18), the approximate density matrix, Eq. (11), is explicitly given by

[TABLE]

III.3 Strong squeezing

Next, we consider the case of strong squeezing with $\delta=1$ and incorporate it into the system Hamiltonian. The rescaled system and perturbation Hamiltonians are given by

[TABLE]

We obtain ${\bm{F}}({\bm{X}})=({x}^{\prime}/2-\Delta^{\prime}{p}^{\prime}-\gamma^{\prime}_{2}{x}^{\prime}({x}^{\prime 2}+{p}^{\prime 2})-2\eta^{\prime}({x}^{\prime}\cos\theta+{p}^{\prime}\sin\theta),\;{p}^{\prime}/2+\Delta^{\prime}{x}^{\prime}-\gamma^{\prime}_{2}{p}^{\prime}({x}^{\prime 2}+{p}^{\prime 2})+2\eta^{\prime}({p}^{\prime}\cos\theta-{x}^{\prime}\sin\theta))^{T}$ with extra terms due to squeezing, characterized by the parameter $\eta^{\prime}$. When $\Delta^{\prime}>2\eta^{\prime}$ (i.e., $\Delta>\ 2\eta$), this vector field ${\bm{F}}({\bm{X}})$ possesses a stable limit-cycle solution ${\bm{X}}_{0}(t)$ in the classical limit. Due to the strong squeezing, this limit cycle is asymmetric and the angular velocity of the oscillator state is non-uniform. At $\Delta^{\prime}=2\eta^{\prime}$, this limit cycle disappears via a saddle-node bifurcation on invariant circle. The perturbation is given by ${\bm{q}}({\bm{X}},t)=(-E^{\prime},0)$.

In this case, the system is not analytically solvable, but we can numerically obtain the limit-cycle solution ${\bm{X}}_{0}(\phi)=(x_{0}(\phi),p_{0}(\phi))^{T}$, natural frequency $\omega$, PSF $\bm{Z}(\phi)$, and Hessian matrix ${\bm{Y}}(\phi)$, and use them in the phase equation in Eq. (5). The density matrix can be approximately reconstructed from Eq. (11), where $\alpha_{0}(\phi)=x_{0}(\phi)+ip_{0}(\phi)$. In this case, the frequency shift does not vanish generally and the effective frequency $\tilde{\omega}$ is slightly different from $\omega$ in the classical limit without noise.

An example of the limit cycle in the classical limit is shown in Fig. 3(c), and the PSF is shown in Fig. 3(d) and (e). The effective frequency is evaluated as $\tilde{\omega}=0.7743$ at the parameter values given in Fig. 3, which is slightly different from the natural frequency $\omega=0.7746$ of the system in the classical limit without noise. From the phase equation, we can obtain the stationary phase distribution $P(\phi)$ by solving the corresponding FPE and reconstruct the density matrix as a mixture of the coherent states on the limit cycle.

III.4 Reconstruction of density matrices

To test the validity of the reduced phase equation, we compare the density matrix $\rho_{sc}$, which is reconstructed from Eq. (11) by using $P(\phi)$ obtained from the FPE in Eq. (10) associated with the reduced phase equation in Eq. (5), with the true density matrix $\rho_{qm}$, which is obtained by direct numerical simulation of the original quantum master equation in Eq. (12), in the steady state of the system. We use the fidelity $F={\rm Tr}\hskip 1.9919pt[\sqrt{\sqrt{\rho_{sc}}\rho_{qm}\sqrt{\rho_{sc}}}]$ Nielsen and Chuang (2000) to quantify the similarity between $\rho_{sc}$ and $\rho_{qm}$. Numerical simulations of the master equation have been performed by using QuTiP Johansson et al. (2012); *johansson2013qutip numerical toolbox.

Figure 2(a)-(d) show the steady-state Wigner distributions corresponding to $\rho_{sc}$ and $\rho_{qm}$ under the weak harmonic driving or the squeezing. In both cases, the distribution is localized along the limit cycle in the classical limit, where the width of the distribution is determined by the intensity of the quantum noise. In Fig. 2(a) and (b), only the harmonic driving is given as the perturbation $(\eta=0)$, while in Fig. 2(c) and (d), only the squeezing is given as the perturbation ($E=0$). It can be seen that the true density matrix $\rho_{qm}$ is accurately approximated by the density matrix $\rho_{sc}$ reconstructed from the phase equation in both cases. The fidelity is $F=0.963$ in the former case and $F=0.982$ in the latter case.

It is notable that the Wigner distribution is localized around one phase point on the limit cycle in Fig. 2(a) and (b), which indicates that there is a 1:1 phase locking Pikovsky et al. (2001) between the oscillator and the harmonic driving; In the classical limit, the phase is locked to the point where the deterministic part of Eq. (18) vanishes, and thus the Wigner distribution takes large values around such a point. Similarly, the Wigner distribution is localized around two phase points on the cycle shown in Fig. 2(c) and (d), because the frequency of the squeezing is twice that of the harmonic driving and 1:2 phase locking occurs, as can be expected from the third term in the deterministic part of Eq. (18) representing the effect of the squeezing. Note that Fig. 2 is depicted in the rotating coordinate frame of frequency $\omega_{d}$ and the locked phase rotates with frequency $\omega_{d}$ in the original coordinate.

Figure 3(a) and (b) show the Wigner distributions in the case of strong squeezing and weak harmonic driving, where all quantities are calculated numerically. In this case, the system exhibits a stable limit cycle in the rotating coordinate frame of frequency $\omega_{d}$, and constant driving is applied on the the system as in Eq. (13). The limit cycle in the classical limit is shown in Fig. 3(c), the $x$ and $p$ components of the PSF obtained from Eq. (7) are shown in Figs. 3(d) and (e), the $xx$, $pp$, $xp$ components of the Hessian matrix are shown in Fig. 3(f),(g), and (h) (the $px$ component is equal to the $xp$ component), and the additional term $g(\phi)$ is shown in Fig. 3(i). The origin of the phase $\phi=0$ is chosen as the intersection of the limit cycle and the ${x}^{\prime}$ axis with ${x}^{\prime}>0$.

It can be seen that the limit cycle in the classical limit is asymmetric due to the effect of the strong squeezing. The density matrix $\rho_{sc}$ can be reconstructed from the phase distribution $P(\phi)$ obtained numerically. As shown in Fig. 3(a) and (b), the true density matrix $\rho_{qm}$ is well approximated by $\rho_{sc}$ with fidelity $F=0.976$. In Fig. 3(a) and (b), the Wigner distribution is concentrated around the stable phase point where the deterministic part of the phase equation vanishes. Thus, the reduced phase equation well reproduces the density matrix of the original quantum system also in this case.

III.5 Reconstruction of spectra and observed frequencies

The power spectrum $S_{qm}$ of the original quantum system in the steady state is defined as

[TABLE]

where $R_{qm}$ is the autocovariance and $\langle{A}\rangle_{qm}={\rm Tr}\hskip 1.9919pt{[A\rho_{qm}]}$ represents the expectation value of an operator $A$ with respect to the steady state density matrix $\rho_{qm}$ obtained from the master equation in Eq. (12). From the reduced phase equation, using the correspondence between the operators and c-numbers in the P representation, the power spectrum in Eq. (21) under the semiclassical approximation can be reconstructed as

[TABLE]

Here, $R_{sc}$ is the autocovariance reconstructed from the phase equation, the mean of a $2\pi$-periodic function $B(\phi)$ is given by $\langle{B(\phi)}\rangle_{sc}=\int_{0}^{2\pi}d\phi B(\phi)P_{sc}(\phi)$, and the autocorrelation is given by $\langle{B(\phi_{2}(\tau))B(\phi_{1}(0))}\rangle_{sc}$ $=\int_{0}^{2\pi}d\phi_{1}\int_{0}^{2\pi}d\phi_{2}(B(\phi_{2}(\tau))B(\phi_{1}(0)))$ $P(\phi_{2},\tau|\phi_{1},0)P_{sc}(\phi_{1})$, where $P_{sc}(\phi)$ is a steady phase distribution and $P(\phi_{2},t_{2}|\phi_{1},t_{1})$ is a transition probability. Both of these probability distributions can be calculated from Eq. (10). The observed frequency $\omega_{qm}$ of the original system and its approximation $\omega_{sc}$ by the phase reduction can be evaluated from the maxima of the spectra as $\omega_{qm,sc}=\mathop{\rm arg~{}max}\limits_{\omega}S_{qm,sc}(\omega)$, respectively.

First, we consider the cases with weak squeezing. Figure 4(a) shows the two power spectra $S_{qm}$ and $S_{sc}$ for the case where only the harmonic driving is given, and Fig. 4(b) shows the spectra for the case with squeezing only. In both cases, the true spectrum $S_{qm}$ can be accurately approximated by the reconstructed spectrum $S_{sc}$. The dependence of the observed frequencies $\omega_{qm,sc}$ on the parameter $\Delta$, where $\Delta$ determines the natural frequency of the limit cycle in the classical limit, is shown in Fig. 4(d) and (e). It can be confirmed that $\omega_{qm}$ is accurately approximated by $\omega_{sc}$ in both cases. The oscillator strictly synchronizes to the external driving when the frequency of the oscillator vanishes in the classical limit, because the harmonic driving acts as a constant force in the rotating frame. Here, strict synchronization is prevented by the quantum noise and the observed frequencies $\omega_{qm,sc}$ do not vanish completely; however, the tendency toward synchronization can be clearly seen from the decrease in the observed frequency compared to that of the unperturbed case.

Next, we consider the case with strong squeezing, where the system exhibits asymmetric limit cycle in the classical limit when $\Delta>2\eta$. We cannot analyze synchronization with the harmonic driving as a stationary problem by using a rotating coordinate frame of frequency $\omega_{d}$, because the limit cycle is asymmetric and the variation in $\Delta$ does not correspond directly to the variation in $\omega_{d}$. We thus explicitly apply harmonic driving with periodic amplitude modulation $E\cos\omega_{e}t$ of frequency $\omega_{e}$ and measure $\omega_{qm}$ and $\omega_{sc}$ as functions of $\Delta_{e}=\omega-\omega_{e}$ for $0\leq\Delta_{e}\leq 0.1~{}(\omega-0.1\leq\omega_{e}\leq\omega)$, where $\omega=0.7746$.

In this case, we obtain a periodic (cyclo-stationary) solution of period $T_{e}=2\pi/\omega_{e}$ instead of a stationary solution. As shown in Fig. 5(a), the quantum-mechanical averages $\langle{x}\rangle$ and $\langle{p}\rangle$ of the position and momentum operators ${x}=(a+a^{{\dagger}})/2$ and ${p}=-i(a-a^{{\dagger}})/2$ exhibit steady periodic dynamics after the initial transient. Here, the initial condition is a coherent state $|\alpha_{0}(\phi=0)\rangle$, where ${\bm{\alpha}}_{0}(\phi=0)$ is a point on the limit cycle with $\phi=0$. Figure 5(b)-(e) show snapshots of the Wigner distributions in the periodic state, where the system evolves as (b) $\to$ (c) $\to$ (d) $\to$ (e) $\to$ (b) (see Supplemental Video for the continuous evolution sup ). The tendency toward synchronization can be clearly observed from the existence of the dense region co-rotating with the external forcing.

We denote the quantum and approximated autocovariance functions at a given time $t_{e}$ ($0\leq t_{e}<T_{e}$) of the steady state oscillation as $R_{qm,sc}^{t_{e}}(\tau)$, where $R^{t_{e}}_{qm}(\tau)$ is calculated by using a density matrix $\rho_{qm}(t_{e})$ at time $t_{e}$ and $R^{t_{e}}_{sc}(\tau)$ is calculated by using a phase distribution $P_{sc}(\phi,t_{e})$ at time $t_{e}$, respectively, in the steadily oscillating state. Then we use the averaged power spectra $\bar{S}_{qm,sc}(\omega)=\int_{-\infty}^{\infty}d\tau e^{i\omega\tau}\int_{0}^{T_{e}}dt_{e}R_{qm,sc}^{t_{e}}(\tau)/T_{e}$ to evaluate the observed frequencies relative to the frequency of the amplitude modulation as $\bar{\omega}_{qm,sc}=\mathop{\rm arg~{}max}\limits_{\omega}\bar{S}_{qm,sc}(\omega)-\omega_{e}$. Figure 4(c) and (f) compare the averaged spectra $\bar{S}_{qm,sc}(\omega)$ and observed frequencies $\bar{\omega}_{qm,sc}$ obtained by direct numerical simulation of the original master equation and by the approximate phase equation, respectively. It can be seen that the spectrum and observed frequency obtained from the original master equation are accurately reproduced by those obtained from the approximate phase equation. Thus, by using the reduced phase equation, we can approximately reconstruct the spectrum and observed frequency of the original system also in this case.

IV Concluding remarks

We have developed a general framework of the phase reduction theory for quantum limit-cycle oscillators under the semiclassical approximation and confirmed its validity by analyzing synchronization dynamics of the quantum vdP model. The proposed framework can approximately characterize the dynamics of a quantum nonlinear oscillator by using a simple classical phase equation, which would serve as a starting point for analyzing synchronization of quantum nonlinear oscillators under the semiclassical approximation. Although we have only analyzed a single-oscillator problem with a single degree of freedom in this study, the developed framework can be directly extended to two or more quantum oscillators with weak coupling by using standard methods from the classical phase reduction theory. Analysis of large many-body systems and the study of their collective dynamics are of particular interest Lee et al. (2014); Witthaut et al. (2017); Ludwig and Marquardt (2013); Lee and Sadeghpour (2013); Davis-Tilley et al. (2018).

In this study, we have employed the P-representation for formulating the semiclassical phase reduction theory; however, other quasiprobability distributions can also be used for the formulation. Detailed comparisons of the results between different representations, including the positive-P representation which is necessary to treat negative-definite diffusion matrices Carmichael (2007), will be discussed in our forthcoming studies. Also, analysis on the genuine quantum signature of a quantum limit-cycle oscillator, which, for instance, can be measured by the negativity of a Wigner quasiprobability distribution Weiss et al. (2017); Lörch et al. (2016), could be performed via an extended version of the developed phase reduction theory.

Recently, the phase reduction theory has been applied to control and optimization of synchronization dynamics in classical nonlinear oscillators  Harada et al. (2010); Zlotnik and Li (2012); Zlotnik et al. (2013); Pikovsky (2015); Watanabe et al. (2019); Monga et al. (2018). In classical dissipative systems, the phase reduction theory has already been used in technical applications of synchronization such as in the ring laser gyroscope Macek and Davis Jr (1963); Cresser et al. (1982a); *cresser1982quantum2; *cresser1982quantum3, phase-locked loop Best (1984); Pikovsky et al. (2001), and Josephson voltage standard Josephson (1962); Shapiro (1963); Pikovsky et al. (2001). The quantum version of these applications, as well as the recent demonstrations Xu and Holland (2015); Hamerly and Mabuchi (2015), could be systematically investigated via the semiclassical phase reduction theory developed in the present study. These subjects will also be discussed in our forthcoming studies.

Acknowledgements.

This research was supported by the JSPS KAKENHI Grant Numbers JP16K13847, JP17H03279, 18K03471, and JP18H03287.

Appendix A Explicit form of $\bm{\beta}(\bm{\alpha})$

In this section, we derive an explicit expression of $\bm{\beta}(\bm{\alpha})$ in Eq. (3). The diffusion matrix of the FPE in Eq. (2) in the complex representation is given by

[TABLE]

where $D_{22}(\bm{\alpha})=D_{11}^{*}(\bm{\alpha})$ and $D_{12}(\bm{\alpha})=D_{21}(\bm{\alpha})$. The non-diagonal element $D_{12}(\bm{\alpha})=D_{21}(\bm{\alpha})$ is real and positive, because it is a constant of cross diffusion described by $\partial^{2}P({\bm{\alpha}},t)/\partial\alpha\partial\alpha^{*}$ and it can be obtained as an absolute value of a complex variable.

We rewrite the FPE in Eq. (2) corresponding to the SDE in Eq. (4) in the real-valued representation, i.e., for the quasiprobability distribution $P({\bm{X}},t)$ with $\bm{X}=(x,p)^{T}=(\mbox{Re}~{}\alpha,\mbox{Im}~{}\alpha)^{T}$, as

[TABLE]

where

[TABLE]

The real-valued diffusion matrix ${\bm{D}}({\bm{X}})$ in the above FPE and the complex-valued diffusion matrix ${\bm{D}}({\bm{\alpha}})$ are related as

[TABLE]

and

[TABLE]

By denoting the matrix components of ${\bm{D}}({\bm{\alpha}})$ in the polar representation as $D_{11}(\bm{\alpha})=R_{11}(\bm{\alpha})e^{i\chi({\bm{\alpha}})}$ and $D_{12}(\bm{\alpha})=R_{12}(\bm{\alpha})$, where $R_{11}({\bm{\alpha}}),R_{22}({\bm{\alpha}})\geq 0$ and $\chi({\bm{\alpha}})\in[0,2\pi)$, the eigenvalues $\lambda_{\pm}(\bm{X})$ and eigenvectors $\bm{v}_{\pm}(\bm{X})$ of $\bm{D}(\bm{X})$ can be expressed as

[TABLE]

and ${\bm{D}}({\bm{X}})$ can be decomposed as

[TABLE]

Thus, ${\bm{G}}({\bm{X}})$ is given by

[TABLE]

and $\bm{\beta}(\bm{\alpha})$ is obtained from ${\bm{G}}({\bm{X}})$ as

[TABLE]

The assumption in the main text that the diffusion matrix is always positive semidefinite along the limit cycle is equivalent to the assumption that $\lambda_{-}(\bm{X}_{0}(\phi))\geq 0$, that is, $R_{12}(\bm{\alpha}_{0}(\phi))\geq R_{11}(\bm{\alpha}_{0}(\phi))$ is satisfied for all $\phi$, because $\lambda_{+}(\bm{X})$ is always positive. With this assumption, if the initial state is given in the form of Eq. (11), for instance, by a pure coherent state $\rho=\left|\alpha_{0}(\phi_{0})\right\rangle\left\langle\alpha_{0}(\phi_{0})\right|$ at a given phase point $\phi_{0}$ on the limit cycle, the state always remains in the two-dimensional phase space of the classical variables.

Appendix B Derivation of the phase equation

In this section, we give a detailed derivation of the phase equation in Eq. (5). The asymptotic phase function $\Phi({\bm{X}}):B\subset{\mathbb{R}}^{2\times 1}\to[0,2\pi)$ introduced in the main text satisfies

[TABLE]

in the basin $B$ of the limit cycle, where $\nabla\Phi\in{\mathbb{R}}^{2\times 1}$ indicates the gradient of $\Phi$ with respect to ${\bm{X}}$. Using this $\Phi({\bm{X}})$, we define the phase $\phi$ of the oscillator state ${\bm{X}}$ as $\phi={\Phi}({\bm{X}})$. As long as ${\bm{X}}$ evolves in $B$, $\dot{\phi}=\dot{\Phi}({\bm{X}})=\dot{\bm{X}}\cdot\nabla\Phi({\bm{X}})=\bm{F}({\bm{X}})\cdot\nabla\Phi({\bm{X}})=\omega$ holds. Recently, it has been shown that this phase function is closely related to an eigenfunction of the Koopman operator of the system $\dot{\bm{X}}={\bm{F}}({\bm{X}})$ associated with the eigenvalue $i\omega$ Mauroy et al. (2013).

When ${\bm{X}}$ obeys the Ito SDE in Eq. (4), we obtain an Ito SDEs for the phase $\phi$ as

[TABLE]

where the third term in the drift part arises from the change of the variables by the Ito formula and $\nabla^{T}\nabla\Phi\in{\mathbb{R}}^{2\times 2}$ represents the Hessian matrix of $\Phi({\bm{X}})$ with respect to ${\bm{X}}$. This equation is still not closed in the phase variable $\phi$, because each term on the right-hand side depends on ${\bm{X}}$.

When the perturbation and quantum noise are weak, the deviation of the system state ${\bm{X}}$ from the limit cycle is small and of the order of $\mathcal{O}(\sqrt{\epsilon})$ because the limit cycle is exponentially stable and the system state is subjected to Gaussian-white noise. Thus, in the lowest-order approximation, we can approximate the state ${\bm{X}}$ by a state ${\bm{X}}_{0}(\phi)$ on the limit cycle as ${\bm{X}}(t)={\bm{X}}_{0}(\phi(t))+\mathcal{O}(\sqrt{\epsilon})$. We then obtain an Ito SDE for the phase variable $\phi$,

[TABLE]

which is correct up to $\mathcal{O}(\epsilon)$ in the drift term and up to $\mathcal{O}(\sqrt{\epsilon})$ in the noise intensity, where

[TABLE]

represents the effect of the perturbation on $\phi$,

[TABLE]

represents the effect of the quantum noise on $\phi$, and

[TABLE]

represents a term arising from the change of the variables, respectively.

We denote the gradient vector (PSF) and Hessian matrix of the phase function $\Phi({\bm{X}})$ evaluated at ${\bm{X}}={\bm{X}}_{0}(\phi)$ on the limit cycle as $\bm{Z}(\phi)=\nabla\Phi|_{{\bm{X}}={\bm{X}}_{0}(\phi)}$ and $\bm{Y}(\phi)=\nabla^{T}\nabla\Phi|_{{\bm{X}}={\bm{X}}_{0}(\phi)}$, respectively. The components of the PSF and Hessian matrix $\nabla^{T}\nabla\Phi|_{{\bm{X}}={{\bm{X}}_{0}(\phi)}}$ are explicitly given by

[TABLE]

for $i,j=1,2$, respectively.

It is well known in the classical phase reduction theory Nakao (2016); Ermentrout and Terman (2010); Ermentrout (1996); Brown et al. (2004) that ${\bm{Z}}(\phi)$ is given by a $2\pi$-periodic solution to the following adjoint equation and normalization condition:

[TABLE]

It is also known Suvak and Demir (2010); Takeshita and Feres (2010) that the Hessian matrix ${\bm{Y}}(\phi)$ of the phase function, evaluated at ${\bm{X}}={\bm{X}}_{0}(\phi)$ on the limit cycle, is given by a $2\pi$-periodic solution to a differential equation

[TABLE]

which satisfies a constraint

[TABLE]

In the above equations, ${\bm{J}}(\phi)\in{\mathbb{R}}^{2\times 2}$ is a Jacobian matrix of ${\bm{F}}({\bm{X}})$ at ${\bm{X}}={\bm{X}}_{0}(\phi)$ and ${\bm{K}}(\phi)\in{\mathbb{R}}^{2\times 2\times 2}$ is a third order tensor, respectively, whose components are given by

[TABLE]

and the matrix components of the product ${\bm{Z}}(\phi)\circ{\bm{K}}(\phi)\in{\mathbb{R}}^{2\times 2}$ are given by

[TABLE]

for $i,j,k=1,2$.

Thus, when the noise and perturbations are sufficiently weak, we obtain an approximate Ito SDE for the phase variable as

[TABLE]

at the lowest order, which corresponds to Eq. (5) in the main text. It can be shown that the amplitude effect does not enter the phase dynamics at the lowest order Shirasaka et al. (2017b). As Eq. (55) is an Ito SDE, using the property of the Wiener process, the noise term can be rewritten as

[TABLE]

where $h(\phi)=\sqrt{\sum_{i=1}^{2}({\bm{h}}(\phi))_{i}^{2}}$ and $W(t)$ is a one-dimensional Wiener process.

The errors in the evolution of the phase variable resulting from the lowest-order approximation above are $\mathcal{O}(\epsilon^{2})$ in the drift term and $\mathcal{O}(\epsilon)$ in the noise intensity, respectively. Therefore, the error in the mean of $\phi$ from the true value grows with time as $\mathcal{O}(\epsilon^{2}t)$, and the error in the variance of $\phi$ grows as $\mathcal{O}(\epsilon^{2}t)$. Thus, these errors in the phase dynamics remain $\mathcal{O}(\epsilon)$ up to $t=\mathcal{O}(1/\epsilon)$.

Appendix C Averaged phase equation

In this section, we derive the averaged phase equation, Eq. (9), by using the near-identity transform. Although Eq. (5) is a correct phase equation for the phase $\phi$ in the lowest-order approximation, it has an additional function $g(\phi)$ in the drift term, which adds tiny periodic fluctuations to the deterministic part. By further introducing a new phase $\psi$ that is only slightly different from $\phi$, we can eliminate this term and obtain a simpler SDE,

[TABLE]

where $f(\psi,t)=\bm{Z}(\psi)\cdot{\bm{q}}(\psi,t)$, $W(t)$ is a one-dimensional Wiener process, and $h(\psi)$ is a $2\pi$-periodic function of $\psi$. Here, the new phase $\psi$ is defined from $\phi$ by a near-identity transform as $\phi=\psi+\epsilon n(\psi)$, where $n(\psi)$ is a $2\pi$-periodic function with $n(0)=0$. Using this transformation, the additional term $g(\phi)$ in Eq. (5) can be renormalized into the frequency term as

[TABLE]

where $\tilde{\omega}$ is the effective frequency of the system. As $\epsilon$ is assumed to be sufficiently small, the transformation between the two variables $\phi$ and $\psi$ is invertible. Thus, the qualitative properties of the dynamics predicted by the two-phase equations, such as whether synchronization occurs or not, are invariant. In the classical phase-reduction theory, the $\mathcal{O}(\epsilon)$ difference between the phase variables due to the near-identity transformation or averaging is often neglected and both phases are considered to be the same. Below, we derive the simplified phase equation in Eq. (57) from the original phase equation, Eq. (5) or (55), by using the near-identity transform Sanders and Verhulst (1985).

In Eq. (55), the function $g(\phi)$ contains the Hessian matrix ${\bm{Y}}(\phi)$ of $\Phi({\bm{X}})$ on the limit cycle, which is typically not included in the phase equation for classical limit-cycle oscillators and gives a tiny but complex periodic contribution to the phase dynamics. To eliminate this term, we renormalize it into the frequency term. For this purpose, we consider a near-identity transform from the original phase $\phi$ to a new phase $\psi$,

[TABLE]

where the transformation function $n(\psi)$ is a smooth $2\pi$-periodic function of $\psi$ satisfying $n(0)=0$, and assume that $\psi$ obeys an Ito SDE of the form

[TABLE]

in the lowest-order approximation, which does not contain a term corresponding to $g(\phi)$ but has a small shift $\epsilon\Omega$ in the frequency. From this SDE, we obtain an Ito SDE for $\phi$ by using the Ito formula as

[TABLE]

where we omitted the tiny terms of $\mathcal{O}(\epsilon^{2})$ in the drift term and $\mathcal{O}(\epsilon^{3/2})$ in the noise intensity. The replacement of $\phi$ by $\psi$ in the functions $f$ and $h$ also results in errors of $\mathcal{O}(\epsilon^{2})$ and $\mathcal{O}(\epsilon^{3/2})$ in the drift term and noise intensity, respectively, which can also be neglected.

The above equation coincides with the original Eq. (55) if $n(\psi)$ satisfies

[TABLE]

As $g(\phi)=g(\psi)+\mathcal{O}(\epsilon)$, the equation for $n(\psi)$ is obtained at the lowest order as

[TABLE]

which gives

[TABLE]

where $n(0)=0$ is used. Moreover, as $n(\psi)$ is $2\pi$-periodic, $n(2\pi)=n(0)=0$ should hold, which determines the frequency shift $\Omega$ as

[TABLE]

Thus, by introducing the near-identity transform, we obtain an averaged phase equation

[TABLE]

where $\tilde{\omega}=\omega+\epsilon\Omega$ is a renormalized, effective frequency. This corresponds to Eq. (57). The orders of errors caused by the above near-identity transformation are $\mathcal{O}(\epsilon^{2})$ in the drift term and $\mathcal{O}(\epsilon^{3/2})$ in the noise intensity. Therefore, the phase equations in Eq. (55) and  (68) are equally correct in the lowest-order approximation and valid up to $t=\mathcal{O}(1/\epsilon)$.

The frequency shift $\epsilon\Omega$ can be evaluated by numerically calculating the Hessian matrix of $\Phi({\bm{X}})$ in $g(\psi)$ and integrating Eq. (67), or alternatively by measuring $\tilde{\omega}$ by numerically evolving the SDE in Eq. (3) or Eq. (4) without perturbations. In the examples used in the main text, the frequency shift $\epsilon\Omega$ is zero in the case of Eq. (18) with the symmetric limit cycle with weak squeezing, and takes a tiny value in the case with strong squeezing. In other applications, for example, in the analysis of coupled identical limit-cycle oscillators without external forcing, the precise value of $\tilde{\omega}$ may not be required (only the frequency difference matters). In such cases, one may simply assume $\tilde{\omega}\approx\omega$ and avoid the calculation of $\epsilon\Omega$.

Appendix D Phase-space representation of a quantum vdP oscillator with harmonic driving and squeezing

D.1 Weak squeezing

Here, we derive a phase equation for a quantum vdP oscillator with harmonic driving and squeezing. In the case of weak squeezing with $\delta=\epsilon$, the rescaled system Hamiltonian and the perturbation Hamiltonian are given by

[TABLE]

respectively, where the squeezing term is included in the perturbation. The functions $\bm{A}(\bm{\alpha^{\prime}})$, $\epsilon\bm{A}^{\prime}(\bm{\alpha^{\prime}})$, and $\epsilon\bm{D}(\bm{\alpha^{\prime}})$ in the quantum FPE are calculated as

[TABLE]

and

[TABLE]

where the tiny terms of $\mathcal{O}(\epsilon^{2})$ in $\epsilon{\bm{D}}({\bm{\alpha}}^{\prime})$ are dropped. The explicit form of $\bm{\beta}(\bm{\alpha}^{\prime})$ given by Eq. (39) can be obtained from Eq. (71) as

[TABLE]

where the modulus $R^{\prime}$ and argument $\delta^{\prime}$ of $\alpha^{\prime}$ is introduced as $\alpha^{\prime}=R^{\prime}e^{i\delta^{\prime}}$. In the real-valued representation with ${\bm{X}}=(x^{\prime},p^{\prime})^{T}=(\mbox{Re}\ \alpha,\mbox{Im}\ \alpha)^{T}$, the functions ${\bm{F}}({\bm{X}})$, $\epsilon{\bm{q}}({\bm{X}})$, and $\sqrt{\epsilon}{\bm{G}}({\bm{X}})$ are given by

[TABLE]

and

[TABLE]

respectively.

As discussed in the main text, the deterministic part of this equation, $\dot{\bm{X}}={\bm{F}}({\bm{X}})$, is a normal form of the supercritical Hopf bifurcation, also known as the Stuart-Landau oscillator, and it is analytically solvable. The limit cycle of this system in the classical limit can be obtained as ${\bm{X}}_{0}(\phi)=\sqrt{\frac{1}{2\gamma^{\prime}_{2}}}(\cos\phi,\sin\phi)^{T}$ with $\phi=\omega t$, or ${\bm{\alpha}}_{0}^{\prime}(\phi)=\sqrt{\frac{1}{2\gamma^{\prime}_{2}}}(e^{i\phi},e^{-i\phi})^{T}$ in the complex-valued representation, where the natural frequency is given by $\omega=\Delta^{\prime}$, and the frequency shift $\epsilon\Omega$ vanishes. From Eq. (71), the eigenvalues of matrix ${\bm{D}}({\bm{\alpha}})$ can be calculated as

[TABLE]

By plugging the limit-cycle solution ${\bm{X}}_{0}(\phi)$ into this equation, it can be seen that $\lambda_{-}(\bm{X}_{0}(\phi))=\frac{1}{4}>0$ is satisfied for any $\phi$ on the limit cycle and the diffusion matrix is always positive semidefinite along the limit cycle, because the magnitudes of the squeezing and nonlinear damping, which can cause negative diffusion, are assumed to be sufficiently small.

D.2 Strong squeezing

In the case of strong squeezing with $\delta=1$, the rescaled system Hamiltonian and the perturbation Hamiltonian are given by

[TABLE]

respectively, where the squeezing term is included in the system Hamiltonian. The functions $\bm{A}(\bm{\alpha^{\prime}})$, $\bm{A}^{\prime}(\bm{\alpha^{\prime}})$, and $\bm{D}(\bm{\alpha^{\prime}})$ in the phase-space representation are given by

[TABLE]

and

[TABLE]

The explicit form of $\bm{\beta}(\bm{\alpha}^{\prime})$ in this case is given by

[TABLE]

where $R^{\prime}_{2}e^{i\chi^{\prime}_{2}}=-(\gamma^{\prime}_{2}\alpha^{\prime 2}+2\eta^{\prime}e^{i\theta})$. In the real-valued representation with ${\bm{X}}=(x^{\prime},p^{\prime})^{T}=(\mbox{Re}\ \alpha,\mbox{Im}\ \alpha)^{T}$, the functions ${\bm{F}}({\bm{X}})$, $\epsilon{\bm{q}}({\bm{X}})$, and $\sqrt{\epsilon}{\bm{G}}({\bm{X}})$ are given by

[TABLE]

and

[TABLE]

respectively.

The deterministic part ${\bm{F}}({\bm{X}})$ gives an asymmetric limit cycle when $\eta>0$, which is difficult to solve analytically. However, we can still obtain the limit cycle ${\bm{X}}_{0}(\phi)$ numerically and use it to evaluate the PSF ${\bm{Z}}(\phi)$, Hessian matrix ${\bm{Y}}(\phi)$, and the noise intensity ${\bm{G}}(\phi)$, and use these quantities in the phase equation. The PSF ${\bm{Z}}(\phi)$ can be numerically calculated by the adjoint method, and the Hessian matrix ${\bm{Y}}(\phi)$ can be calculated by using a shooting-type numerical algorithm.

When the squeezing is too strong, the diffusion matrix can generally be negative definite on the limit cycle. We choose parameter settings where the diffusion matrix is always positive semidefinite along the limit cycle in the main text.