Semiclassical Phase Reduction Theory for Quantum Synchronization

# Semiclassical Phase Reduction Theory for Quantum Synchronization

Yuzuru Kato Department of Systems and Control Engineering, Tokyo Institute of Technology, Tokyo 152-8552, Japan    Naoki Yamamoto Department of Applied Physics and Physico-Informatics, Keio University, Kanagawa 223-8522, Japan    Hiroya Nakao Department of Systems and Control Engineering, Tokyo Institute of Technology, Tokyo 152-8552, Japan
July 23, 2019
###### 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.

## 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),

 ˙ρ=−i[H+ϵ~H(t),ρ]+n∑m=1D[Lm]ρ, (1)

where is a density matrix representing the system state, is a system Hamiltonian, is a time-dependent Hamiltonian representing weak external perturbations applied to the system (), is the number of reservoirs, is the coupling operator between the system and the th reservoir , and 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 is represented as , where is a coherent state specified by a complex value , or equivalently by a two-dimensional complex vector , is a quasiprobability distribution of , , the integral is taken over the entire space spanned by , and * indicates complex conjugate.

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

 ∂P(α,t)∂t=[−2∑j=1∂j{Aj(α)+ϵ~Aj(α,t)}+122∑j=12∑k=1∂j∂k{ϵDjk(α)}]P(α,t), (2)

where and are the th components of complex vectors , and representing the system dynamics and perturbations, respectively, is the -th component of the symmetric diffusion matrix representing quantum fluctuations, and the complex partial derivatives are defined as and (note that and ).

The drift term consists of terms arising from the system Hamiltonian and the dissipation , represents the small terms arising from the perturbation Hamiltonian , and the diffusion matrix represents the intensity of the small quantum noise, generally arising from all terms of , , and . These terms can be explicitly calculated from the master equation in Eq. (1) by using the standard calculus for phase-space representation when , , and are given Carmichael (2007); Gardiner and Haken (1991). The external perturbation and the diffusion matrix are assumed to be of the same order, .

By introducing an appropriate complex matrix (see Appendix A for the explicit form), the diffusion matrix can be represented as and the Ito SDE corresponding to Eq. (2) for the phase-space variable is given by

 dα ={A(α)+ϵ~A(α,t)}dt+√ϵβ(α)dW, (3)

where represents a vector of independent Wiener processes satisfying .

It should be noted that diffusion matrix of certain quantum systems in the P representation becomes negative definite for certain  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 representation. To this end, we define a real vector from the complex vector . The real-valued expression of Eq. (3) for is then given by an Ito SDE,

 dX={F(X)+ϵq(X,t)}dt+√ϵG(X)dW, (4)

where , , and are real-valued equivalent representations of the system dynamics , perturbation , and noise intensity in Eq. (3), respectively.

We assume that the system in the classical limit without perturbation and quantum noise, , has an exponentially stable limit-cycle solution with a natural period and frequency . 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 such that is satisfied for all system states in the basin of the limit cycle in the classical limit, where is the gradient of . Using this phase function, we define the phase of a system state as . It then follows that , i.e., always increases at a constant frequency with the evolution of . In the following formulation, we represent the system state on the limit cycle as as a function of the phase rather than the time . In this representation, is a -periodic function of , . Note that an identity is satisfied by the definition of .

When the noise and perturbations are sufficiently weak and the deviation of the state from the limit cycle is small, we can approximate by a state on the limit cycle as and derive a SDE for the phase in the lowest order approximation by using the Ito formula as (see Appendix B for details)

 dϕ={ω+ϵZ(ϕ)⋅q(ϕ,t)+ϵg(ϕ)}dt+√ϵ{G(ϕ)TZ(ϕ)}⋅dW, (5)

where the drift term is correct up to and the noise intensity is correct up to . Here, the inner product between two vectors and is defined as .

In the above phase equation, the gradient of at is approximately evaluated at on the limit cycle and is denoted as . We call this 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 on the limit cycle and they are denoted as and , respectively. The additional function in the drift term in Eq. (5) arises from the change of the variables and is given by

 g(ϕ)=12Tr{G(ϕ)TY(ϕ)G(ϕ)}, (6)

where is a Hessian matrix of the phase function also evaluated at on the limit cycle. All these functions are -periodic, as they are functions of .

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

 ωddϕZ(ϕ)=−JT(ϕ)Z(ϕ),Z(ϕ)⋅dX0(ϕ)dϕ=1, (7)

respectively, where is a Jacobian matrix of at on the limit cycle. It is also known that the Hessian matrix on the limit cycle can be calculated as a -periodic solution of an adjoint-type equation Suvak and Demir (2010); Takeshita and Feres (2010) with an appropriate constraint. These equations for are detailed in the Appendix B. In the numerical calculations, 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 can be obtained by a shooting method Suvak and Demir (2010).

Because of the additional term in Eq. (10), the effective frequency of the oscillator in the absence of the perturbation is given by

 ~ω=ω+ϵ2π∫2π0g(ψ′)dψ′, (8)

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

 dψ={~ω+ϵZ(ψ)⋅q(ψ,t)}dt+√ϵh(ψ)dW, (9)

where and is a one-dimensional Wiener process. As before, the drift term is correct up to and the noise intensity is correct up to . 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 in Hamerly and Mabuchi (2015).

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

 ∂∂tP(ϕ,t)=−∂∂ϕ{ω+ϵZ(ϕ)⋅q(ϕ,t)+ϵg(ϕ)}P(ϕ,t)+ϵ2∂2∂ϕ2h(ϕ)2P(ϕ,t). (10)

Using this FPE, we can obtain the stationary distribution and transition probability of the phase variable 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 , the oscillator state in the classical limit can be approximated as , or in the original complex representation. Therefore, the quantum state at phase is approximately described as and the density matrix is approximately represented by using the probability density function of the phase variable , obtained from the SDE in Eq. (5) or FPE in Eq. (10), as

 ρ≈∫2π0dϕP(ϕ)|α0(ϕ)⟩⟨α0(ϕ)|, (11)

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 , which is characterized by the natural frequency , PSF , Hessian matrix , and noise intensity 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 , , and 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 Gardiner and Haken (1991). In the rotating coordinate frame of frequency , the master equation is given by Walter et al. (2014); Sonar et al. (2018)

 ˙ρ=−i[−Δa†a+iE(a−a†)+iη(a2e−iθ−a†2eiθ),ρ]+γ1D[a†]ρ+γ2D[a2]ρ, (12)

where is the frequency detuning of the harmonic driving from the oscillator, is the intensity of the harmonic driving, is the squeezing parameter, and are the decay rates for negative damping and nonlinear damping, respectively, and the Planck constant is set as . The harmonic driving is represented by a constant , because a coordinate frame rotating with the driving frequency is used.

We assume that is sufficiently small and of , for which the semiclassical approximation is valid, and represent as using a dimensionless parameter of . In this setting, the size of the stable limit-cycle solution in Eq. (12) in the classical limit is , while we have implicitly assumed it to be in the derivation of Eq. (5). Therefore, we introduce a rescaled annihilation operator and the corresponding classical variable ( in the vector representation) as , and represent the parameters as , where , and are dimensionless parameters of . By this rescaling, the size of the limit cycle becomes and the parameter determines the relative intensity of the squeezing.

The real-valued representation of Eq. (4) after rescaling is then obtained as

 dX=(12x′−Δ′p′−γ′2x′(x′2+p′2)−ϵE′−2δη′(x′cosθ+p′sinθ)12p′+Δ′x′−γ′2p′(x′2+p′2)+2δη′(p′cosθ−x′sinθ))dt′+√ϵG(X)dW′, (13)

where and . The noise intensity matrix is explicitly given by

 G(X) =⎛⎜ ⎜ ⎜ ⎜⎝√(1+R′1)2cosχ′12√(1−R′1)2sinχ′12√(1+R′1)2sinχ′12−√(1−R′1)2cosχ′12⎞⎟ ⎟ ⎟ ⎟⎠, (14)

with . 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 . The rescaled system and perturbation Hamiltonians are given by

 H=−Δ′a′†a′,ϵ~H=ϵ{iE′(a′−a′†)+iη′(a′2e−iθ−a′†2eiθ)}. (15)

For this system, we obtain . The perturbation is represented by . Note that the vector field in this case is simply a normal form of the supercritical Hopf bifurcation. A classical nonlinear oscillator described by this 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

 (16)

as a function of phase , where the frequency is given by . The basin of this limit cycle is the whole -plane except . The phase function of this limit cycle can be expressed as  Nakao (2016), which gives . The PSF and Hessian matrix can be obtained by calculating the gradients of the phase function at on the limit cycle as

 Z(ϕ)=√2γ′2(−sinϕcosϕ),Y(ϕ)=2γ′2(sin2ϕ−cos2ϕ−cos2ϕ−sin2ϕ). (17)

In this case, the additional term in Eq. (5) and therefore the frequency shift in Eq. (8) vanishes, i.e., . The terms in the noise intensity 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

 dϕ ={Δ′+√2ϵ√γ′2E′sinϕ+2ϵη′sin(2ϕ−θ)}dt′+√ϵ√3γ′22dW′ (18)

in the lowest-order approximation, where . Using the probability density function of the phase described by the FPE (10) corresponding to Eq. (18), the approximate density matrix, Eq. (11), is explicitly given by

 ρ≈∫2π0dϕP(ϕ)∣∣∣√γ12γ2eiϕ⟩⟨√γ12γ2eiϕ∣∣∣. (19)

### iii.3 Strong squeezing

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

 H=−Δ′a′†a′+iη′(a′2e−iθ−a′†2eiθ),ϵ~H=ϵiE′(a′−a′†). (20)

We obtain with extra terms due to squeezing, characterized by the parameter . When (i.e., ), this vector field possesses a stable limit-cycle solution 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 , this limit cycle disappears via a saddle-node bifurcation on invariant circle. The perturbation is given by .

In this case, the system is not analytically solvable, but we can numerically obtain the limit-cycle solution , natural frequency , PSF , and Hessian matrix , and use them in the phase equation in Eq. (5). The density matrix can be approximately reconstructed from Eq. (11), where . In this case, the frequency shift does not vanish generally and the effective frequency is slightly different from 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 at the parameter values given in Fig. 3, which is slightly different from the natural frequency of the system in the classical limit without noise. From the phase equation, we can obtain the stationary phase distribution 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 , which is reconstructed from Eq. (11) by using obtained from the FPE in Eq. (10) associated with the reduced phase equation in Eq. (5), with the true density matrix , 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  Nielsen and Chuang (2000) to quantify the similarity between and . 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 and 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 , while in Fig. 2(c) and (d), only the squeezing is given as the perturbation (). It can be seen that the true density matrix is accurately approximated by the density matrix reconstructed from the phase equation in both cases. The fidelity is in the former case and 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 and the locked phase rotates with frequency 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 , 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 and components of the PSF obtained from Eq. (7) are shown in Figs. 3(d) and (e), the , , components of the Hessian matrix are shown in Fig. 3(f),(g), and (h) (the component is equal to the component), and the additional term is shown in Fig. 3(i). The origin of the phase is chosen as the intersection of the limit cycle and the axis with .

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 can be reconstructed from the phase distribution obtained numerically. As shown in Fig. 3(a) and (b), the true density matrix is well approximated by with fidelity . 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 of the original quantum system in the steady state is defined as

 Sqm(ω) =∫∞−∞dτeiωτRqm(τ), (21) Rqm(τ) =⟨a†(τ)a(0)⟩qm−⟨a†(τ)⟩qm⟨a(0)⟩qm, (22)

where is the autocovariance and represents the expectation value of an operator with respect to the steady state density matrix 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

 Ssc(ω) =∫∞−∞dτeiωτRsc(τ), (23) Rsc(τ) =⟨α∗0(ϕ2(τ))α0(ϕ1(0))⟩sc−⟨α∗0(ϕ2(τ))⟩sc⟨α0(ϕ1(0))⟩sc. (24)

Here, is the autocovariance reconstructed from the phase equation, the mean of a -periodic function is given by , and the autocorrelation is given by , where is a steady phase distribution and is a transition probability. Both of these probability distributions can be calculated from Eq. (10). The observed frequency of the original system and its approximation by the phase reduction can be evaluated from the maxima of the spectra as , respectively.

First, we consider the cases with weak squeezing. Figure 4(a) shows the two power spectra and 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 can be accurately approximated by the reconstructed spectrum . The dependence of the observed frequencies on the parameter , where 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 is accurately approximated by 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 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 . We cannot analyze synchronization with the harmonic driving as a stationary problem by using a rotating coordinate frame of frequency , because the limit cycle is asymmetric and the variation in does not correspond directly to the variation in . We thus explicitly apply harmonic driving with periodic amplitude modulation of frequency and measure and as functions of for , where .

In this case, we obtain a periodic (cyclo-stationary) solution of period instead of a stationary solution. As shown in Fig. 5(a), the quantum-mechanical averages and of the position and momentum operators and exhibit steady periodic dynamics after the initial transient. Here, the initial condition is a coherent state , where is a point on the limit cycle with . Figure 5(b)-(e) show snapshots of the Wigner distributions in the periodic state, where the system evolves as (b) (c) (d) (e) (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 () of the steady state oscillation as , where is calculated by using a density matrix at time and is calculated by using a phase distribution at time , respectively, in the steadily oscillating state. Then we use the averaged power spectra to evaluate the observed frequencies relative to the frequency of the amplitude modulation as . Figure 4(c) and (f) compare the averaged spectra and observed frequencies 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 β(α)

In this section, we derive an explicit expression of in Eq. (3). The diffusion matrix of the FPE in Eq. (2) in the complex representation is given by

 D(α)=β(α)β(α)T=(D11(α)D12(α)D21(α)D22(α))∈C2×2, (25)

where and . The non-diagonal element is real and positive, because it is a constant of cross diffusion described by 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 with , as

 ∂∂tP(X,t)=[−∂∂X{F(X)+ϵq(X,t)}+12∂2∂X2D(X)]P(X,t), (26)

where

 ∂∂α=12(∂∂x−i∂∂p),∂∂α∗=12(∂∂x+i∂∂p). (27)

The real-valued diffusion matrix in the above FPE and the complex-valued diffusion matrix are related as

 D(X) =14(11−ii)D(α)(1−i1i) (28) =12(Re D11(α)+D12(α)Im D11(α)Im D11(α)−Re D11(α)+D12(α))∈R2×2 (29)

and

 D(α) =(1i1−i)D(X)(11i−i). (30)

By denoting the matrix components of in the polar representation as and , where and , the eigenvalues and eigenvectors of can be expressed as

 λ±(X) =12(R12(α)±R11(α)), (31) v+(X) =⎛⎜⎝cosχ(α)2sinχ(α)2⎞⎟⎠,v−(X)=⎛⎜⎝sinχ(α)2−cosχ(α)2⎞⎟⎠, (32)

and can be decomposed as

 D(X)=(