Semiclassical Phase Reduction Theory for Quantum Synchronization

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 [1][2][3][4][5][6]. 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 [1][2][3][7][8][9]. By using this framework, we can systematically reduce multidimensional nonlinear dynamical equations describing weakly perturbed limit-cycle oscillators to a onedimensional phase equation that approximately describes the oscillator dynamics. The simple semilinear 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 [1][2][3][4][5][6]. 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 [1,2]. Recently, generalizations of the phase reduction theory to nonconventional physical systems, such as time-delayed oscillators [10,11], piecewise-smooth oscillators [12], collectively oscillating networks [13], and rhythmic spatiotemporal patterns [14,15], have also been discussed.
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 [28] and the precise measurement of the resistance standard with a superconducting device [42]; an application of the limit-cycle oscillation to analog memory in a quantum optical device [43] 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 FIG. 1. A schematic diagram of the semiclassical phase reduction for quantum synchronization. A quantum self-sustained oscillator, which has a stable limit-cycle solution in the classical limit, can be described by an approximate one-dimensional stochastic differential equation for a phase variable φ that characterizes the system state. The system state can be approximately reconstructed from the reduced phase equation.
distributions of quantum systems are commonly employed. In a pioneering study, Hamerly and Mabuchi [43] derived a phase equation from the stochastic differential equation (SDE) describing a truncated Wigner function of a quantum limitcycling 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 [43]. 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 [24,26,28,47]; 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, where the quantum dynamics can be approximately described by a SDE representing a system state in the phase space fluctuating along a deterministic classical trajectory due to small quantum noise. We derive a linearized multidimensional 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 [1][2][3][7][8][9]. 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.

A. 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 timescale 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 [51,52], where ρ is a density matrix representing the system state, H is a system Hamiltonian, H (t ) is a time-dependent Hamiltonian representing weak external perturbations applied to the system (0 < 1), n is the number of reservoirs, L m is the coupling operator between the system and the mth reservoir (m = 1, . . . , n), D[L]ρ = LρL † − (ρL † L + L † Lρ)/2 denotes the Lindblad form, and the Planck constant is set ash = 1. 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 multidimensional SDE by introducing a phase-space quasiprobability distribution, such as the P, Q, or Wigner representation [51,52]. 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 ρ = P(α)|α α|dα, where |α is a coherent state specified by a complex value α ∈ C, or equivalently by a two-dimensional complex vector α = (α, α * ) T ∈ C 2×1 , P(α) is a quasiprobability distribution of α, dα = dαdα * , 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 033012-2 where A j (α) andÃ j (α, t ) are the jth components of complex vectors A(α) = (A 1 (α), A * 1 (α)) T ∈ C 2×1 andÃ(α, t ) = (Ã 1 (α, t ),Ã * 1 (α, t )) T ∈ C 2×1 representing the system dynamics and perturbations, respectively, D jk (α) is the ( j, k)-th component of the symmetric diffusion matrix D(α) ∈ C 2×2 representing quantum fluctuations, and the complex partial derivatives are defined as ∂ 1 = ∂/∂α and . The drift term A(α) consists of terms arising from the system Hamiltonian H and the dissipation {L m }, Ã (α, t ) represents the small terms arising from the perturbation Hamiltonian H (t ), and the diffusion matrix D(α) represents the intensity of the small quantum noise, generally arising from all terms of H, 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, H (t ), and {L m } are given [51,52]. The external perturbation Ã (α, t ) and the diffusion matrix D(α) are assumed to be of the same order, O( ).
By introducing an appropriate complex matrix √ β(α) ∈ C 2×2 (see Appendix A for the explicit form), the diffusion matrix D(α) can be represented as D(α) = √ β(α)( √ β(α)) T and the Ito SDE corresponding to Eq. (2) for the phase-space variable α(t ) is given by where It should be noted that diffusion matrix of certain quantum systems in the P representation becomes negative definite for certain α [51,52]. 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. [40] 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.

B. 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 X = (x, p) T = (Re α, Im α) T ∈ R 2×1 from the complex vector α. The real-valued expression of Eq. (3) for X (t ) is then given by an Ito SDE, where F(X ) ∈ R 2×1 , q(X , t ) ∈ R 2×1 , and G(X ) ∈ R 2×2 are real-valued equivalent representations of the system dynamics A(α) ∈ C 2×1 , perturbationÃ(α, t ) ∈ C 2×1 , and noise intensity β(α) ∈ C 2×2 in Eq. (3), respectively. We assume that the system in the classical limit without perturbation and quantum noise,Ẋ = F(X ), has an exponentially stable limit-cycle solution X 0 (t ) = (x 0 (t ), p 0 (t )) T = X 0 (t + T ) with a natural period T and frequency ω = 2π/T .
In the same way as the phase reduction for classical limit cycles [1][2][3][7][8][9], we can introduce an asymptotic phase function (X ) : B ⊂ R 2×1 → [0, 2π ) such that ∇ (X ) · F(X ) = ω is satisfied for all system states X in the basin B of the limit cycle in the classical limit, where ∇ (X ) ∈ R 2×1 is the gradient of (X ). Using this phase function, we define the phase of a system state X ∈ B as φ = (X ). It then follows thaṫ φ =˙ (X ) = F(X ) · ∇ (X ) = ω, i.e., φ always increases at a constant frequency ω with the evolution of X . Here, the inner product between two vectors a = (a 0 , a 2 , · · · , a N−1 ) T ∈ In the following formulation, we represent the system state X on the limit cycle as X 0 (φ) = (x 0 (φ), p 0 (φ)) T as a function of the phase φ rather than the time t. In this representation, X 0 (φ) is a 2π -periodic function of φ, X 0 (φ) = X 0 (φ + 2π ). Note that an identity (X 0 (φ)) = φ is satisfied by the definition of (X ).
When the noise and perturbations are sufficiently weak and the deviation of the state X from the limit cycle is small, we can approximate X (t ) by a state X 0 (φ(t )) on the limit cycle as X (t ) ≈ X 0 (φ(t )) and derive a SDE for the phase in the lowest order approximation by using the Ito formula as (see Appendix B for details) where the drift term is correct up to O( ) and the noise intensity is correct up to O( √ ). In the above phase equation, the gradient ∇ of (X ) at X is approximately evaluated at X (φ) on the limit cycle and is denoted as Z(φ) = ∇ | X =X 0 (φ) ∈ R 2×1 . We call this Z(φ) the phase sensitivity function (PSF) of the limit cycle, which characterizes the linear response property of the oscillator phase to given perturbations [2,7]. Similarly, the perturbation and noise intensity can also be evaluated approximately at X = X 0 (φ) on the limit cycle and they are denoted as q(φ, t ) = q(X 0 (φ), t ) and G(φ) = G(X 0 (φ)), respectively. The additional function g(φ) in the drift term in Eq. (5) arises from the change of the variables and is given by where is a Hessian matrix of the phase function (X ) also evaluated at X = X 0 (φ) on the limit cycle. All these functions are 2π -periodic, as they are functions of X 0 (φ). It is well known in the classical phase reduction theory that the PSF can be obtained as a 2π -periodic solution to the following adjoint equation and an additional normalization condition [7][8][9]: respectively, where J(φ) = J(X 0 (φ)) ∈ R 2×2 is a Jacobian matrix of F(X ) at X = X 0 (φ) on the limit cycle. It is also known that the Hessian matrix Y (φ) on the limit cycle can be calculated as a 2π -periodic solution to an adjoint-type equation [53,54] with an appropriate constraint. These equations for Y (φ) are detailed in the Appendix B. In the numerical calculations, Z(φ) can easily be obtained by the backward integration of the adjoint equation with occasional normalization as proposed by Ermentrout [3], and then the Hessian Y (φ) can be obtained by a shooting method [53]. Because of the additional term g(φ) in Eq. (10), the effective frequencyω = dφ /dt of the oscillator in the absence of the perturbation q(φ, t ) is given bỹ 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 φ = ψ + n(ψ ), where n(ψ ) is a 2π -periodic function with n(0) = 0, and eliminate the additional function g(φ) in Eq. (5) by renormalizing it into the frequency term.
The new phase ψ then obeys a simpler SDE of the form where is a onedimensional Wiener process. As before, the drift term is correct up to O( ) and the noise intensity is correct up to O( √ ). 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. [43] does not contain a term with the Hessian matrix, because the order of the noise intensity is implicitly assumed to be O( ) in Ref. [43].
From the reduced SDE in Eq. (5), we can derive a corresponding FPE describing the probability density function P(φ, t ) of the phase variable φ as 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.

C. 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 approx- T in the original complex representation. Therefore the quantum state at phase φ is approximately described as |α 0 (φ) and the density matrix ρ is approximately represented by using the probability density function P(φ) of the phase variable φ, obtained from the SDE in Eq. (5) or FPE in Eq. (10), as 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 Z(φ), Hessian matrix Y (φ), and noise intensity G(φ) 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.

A. 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 ω 0 , ω d , and ω 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 ω sq = 2ω d [52]. In the rotating coordinate frame of frequency ω d , the master equation is given by [34,35] where = ω d − ω 0 is the frequency detuning of the harmonic driving from the oscillator, E is the intensity of the harmonic driving, ηe iθ is the squeezing parameter and γ 1 and γ 2 are the decay rates for negative damping and nonlinear damping, respectively. The harmonic driving is represented by a constant E , because a coordinate frame rotating with the driving frequency ω d is used. Note that the Lindblad term with the quadratic annihilation operator, D[a 2 ], is essentially important in giving rise to the limit-cycle oscillations.
We assume that γ 2 is sufficiently small and of O( ), for which the semiclassical approximation is valid, and represent γ 2 as γ 2 = γ 1 γ 2 using a dimensionless parameter γ 2 of O (1). In this setting, the size of the stable limit-cycle solution in Eq. (12) in the classical limit is O(1/ √ ), while we have implicitly assumed it to be O(1) in the derivation of Eq. (5). Therefore we introduce a rescaled annihilation operator a and the corresponding classical variable α (α = (α , α * ) in the vector representation) as a |α = √ a| √ α , and represent the parameters as = γ 1 , E = √ γ 1 E , and η = δγ 1 η , where , E , and η are dimensionless parameters of O (1). By this rescaling, the size of the limit cycle becomes O(1) and the parameter δ determines the relative intensity of the squeezing.

033012-4
The real-valued representation X = (x , p ) T = (Re α , Im α ) T of Eq. (4) after rescaling is then obtained as where dt = γ 1 dt and dW = √ γ 1 dW . The noise intensity matrix is explicitly given by with . Further details of the derivation can be found in the Appendix D.

B. Weak squeezing
First, we consider the case of weak squeezing with δ = . The rescaled system and perturbation Hamiltonians are given by Note that the vector field F(X ) in this case is simply a normal form of the supercritical Hopf bifurcation. A classical nonlinear oscillator described by this F(X ) is known as the Stuart-Landau (SL) oscillator [2] (which is different from the classical vdP oscillator) and it is analytically solvable.
The stable limit cycle of the SL oscillator is given by as a function of phase φ = ωt, where the frequency is given by ω = . The basin B of this limit cycle is the whole (x , p )-plane except (0, 0). The phase function (X ) of this limit cycle can be expressed as ( Hessian matrix Y (φ) can be obtained by calculating the gradients of the phase function (X ) at X = X 0 (φ) on the limit cycle as In this case, the additional term g(φ) in Eq. (5)  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 in the lowest-order approximation, where dW = √ γ 1 dW .
Using the probability density function P(φ) of the phase φ described by the FPE (10) corresponding to Eq. (18), the approximate density matrix, Eq. (11), is explicitly given by

C. Strong squeezing
Next, we consider the case of strong squeezing with δ = 1 and incorporate it into the system Hamiltonian. The rescaled system and perturbation Hamiltonians are given by We obtain F(X ) = (x /2 − p − γ 2 x (x 2 + p 2 ) − 2η (x cos θ + p sin θ ), p /2 + x − γ 2 p (x 2 + p 2 ) + 2η (p cos θ − x sin θ )) T with extra terms due to squeezing, characterized by the parameter η . When > 2η (i.e., > 2η), this vector field F(X ) possesses a stable limit-cycle solution 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 nonuniform. At = 2η , this limit cycle disappears via a saddle-node bifurcation on invariant circle. The perturbation is given by q(X , t ) = (−E , 0).
In this case, the system is not analytically solvable, but we can numerically obtain the limit-cycle solution X 0 (φ) = (x 0 (φ), p 0 (φ)) T , natural frequency ω, PSF Z(φ), and Hessian matrix Y (φ), and use them in the phase equation in Eq. (5). The density matrix can be approximately reconstructed from Eq. (11), where α 0 (φ) = x 0 (φ) + ip 0 (φ). 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 Figs. 3(d) and 3(e). The effective frequency is evaluated asω = 0.7743 at the parameter values given in Fig. 3, which is slightly different from the natural frequency ω = 0.7746 of the system in the classical limit without noise. From the phase equation, we can obtain the stationary phase distribution P(φ) by solving the corresponding FPE and reconstruct the density matrix as a mixture of the coherent states on the limit cycle.

D. Reconstruction of density matrices
To test the validity of the reduced phase equation, we compare the density matrix ρ sc , which is reconstructed from Eq. (11) by using P(φ) obtained from the FPE in Eq. (10) associated with the reduced phase equation in Eq. (5), with the true density matrix ρ qm , which is obtained by direct numerical simulation of master equation in Eq. (12), in the steady state of the system. We use the fidelity F = Tr [ √ ρ sc ρ qm √ ρ sc ] [55] to quantify the similarity between ρ sc and ρ qm . Numerical simulations of the master equation have been performed by using QuTiP [56] numerical toolbox. Figures 2(a)-2(d) show the steady-state Wigner distributions corresponding to ρ sc and ρ 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 Figs. 2(a) and 2(b), only the harmonic driving is given as the perturbation (η = 0), while in Figs. 2(c) and 2(d), only the squeezing is given as the perturbation (E = 0). It can be seen that the true density matrix ρ qm is accurately approximated by the density matrix ρ 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 Figs. 2(a) and 2(b), which indicates that there is a 1:1 phase locking [4] 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 Figs. 2(c) and 2(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 ω d and the locked phase rotates with frequency ω d in the original coordinate. Figures 3(a) and 3(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 ω d , and constant driving is applied on 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.  Fig. 3(i). The origin of the phase φ = 0 is chosen as the intersection of the limit cycle and the x axis with x > 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 ρ sc can be reconstructed from the phase distribution P(φ) obtained numerically. As shown in Figs. 3(a) and 3(b), the true density matrix ρ qm is well approximated by ρ sc with fidelity F = 0.976. In Figs. 3(a) and 3(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. Note that the reconstructed density matrix ρ sc is slightly more concentrated than the original density matrix ρ qm in Figs. 2 and 3. This is because ρ sc is approximated as a weighted mixture of coherent states with minimum uncertainty along the limit cycle.
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 ω qm,sc on the parameter , where determines the natural frequency of the limit cycle in the classical limit, is shown in Figs. 4(d) and 4(e). It can be confirmed that ω qm is accurately approximated by ω 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 ω 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 > 2η. We cannot analyze synchronization with the harmonic driving as a stationary problem by using a rotating coordinate frame of frequency ω d , because the limit cycle is asymmetric and the variation in does not correspond directly to the variation in ω d . We thus explicitly apply harmonic driving with periodic amplitude modulation E cos ω e t of frequency ω e and measure ω qm and ω sc as functions of In this case, we obtain a periodic (cyclostationary) solution of period T e = 2π/ω e instead of a stationary solution. As shown in Fig. 5(a), the quantum-mechanical averages x and p of the position and momentum operators x = (a + a † )/2 and p = −i(a − a † )/2 exhibit steady periodic dynamics after the initial transient. Here, the initial condition is a coherent state |α 0 (φ = 0) , where α 0 (φ = 0) is a point on the limit cycle with φ = 0. Figures 5(b)-5(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 [57]). 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 t e < T e ) of the steady state oscillation as R t e qm,sc (τ ), where R t e qm (τ ) is calculated by using a density matrix ρ qm (t e ) at time t e and R t e sc (τ ) is calculated by using a phase distribution P sc (φ, t e ) at time t e , respectively, in the steadily oscillating state. Then we use the averaged power spectraS qm,sc (ω) = ∞ −∞ dτ e iωτ T e 0 dt e R t e qm,sc (τ )/T e to evaluate the observed frequencies relative to the frequency of the amplitude modulation asω qm,sc = arg max ωS qm,sc (ω) − ω e . Figures 4(c) and 4(f) compare the averaged spectraS qm,sc (ω) and observed frequenciesω 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 [24,29,30,47,48].
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 033012-8 the formulation. Detailed comparisons of the results between different representations, including the positive-P representation which is necessary to treat negative-definite diffusion matrices [51], 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 Wigner quasiprobability distributions [37,41] as observed in the steady state of a quantum vdP oscillator with a strong Kerr drive and external drive [37], 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 [58][59][60][61][62][63]. In classical dissipative systems, the phase reduction theory has already been used in technical applications of synchronization such as in the ring laser gyroscope [64,65], phase-locked loop [4,66], and Josephson voltage standard [4,67,68]. The quantum version of these applications, as well as the recent demonstrations [28,43], 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. 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 where D 22 (α) = D * 11 (α) and D 12 (α) = D 21 (α). The nondiagonal element D 12 (α) = D 21 (α) is real and positive, because it is a constant of cross diffusion described by ∂ 2 P(α, t )/∂α∂α * 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(X , t ) with X = (x, p) T = (Re α, Im α) T , as where The real-valued diffusion matrix D(X ) in the above FPE and the complex-valued diffusion matrix D(α) are related as and By denoting the matrix components of D(α) in the polar representation as D 11 (α) = R 11 (α)e iχ (α) and D 12 (α) = R 12 (α), where R 11 (α), R 22 (α) 0 and χ (α) ∈ [0, 2π ), the eigenvalues λ ± (X ) and eigenvectors v ± (X ) of D(X ) can be expressed as and D(X ) can be decomposed as Thus G(X ) is given by and β(α) is obtained from G(X ) as The assumption in the main text that the diffusion matrix is always positive semidefinite along the limit cycle is equivalent to the assumption that λ − (X 0 (φ)) 0, that is, R 12 (α 0 (φ)) R 11 (α 0 (φ)) is satisfied for all φ, because λ + (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 033012-9 ρ = |α 0 (φ 0 ) α 0 (φ 0 )| at a given phase point φ 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 (X ) : in the basin B of the limit cycle, where ∇ ∈ R 2×1 indicates the gradient of with respect to X . Using this (X ), we define the phase φ of the oscillator state X as φ = (X ). As long as X evolves in B,φ =˙ (X ) =Ẋ · ∇ (X ) = F(X ) · ∇ (X ) = ω holds. Recently, it has been shown that this phase function is closely related to an eigenfunction of the Koopman operator of the systemẊ = F(X ) associated with the eigenvalue iω [69].
When X obeys the Ito SDE in Eq. (4), we obtain an Ito SDEs for the phase φ as dφ = [(∇ (X )) · (F(X ) + q(X , t )) where the third term in the drift part arises from the change of the variables by the Ito formula and ∇ T ∇ ∈ R 2×2 represents the Hessian matrix of (X ) with respect to X . This equation is still not closed in the phase variable φ, because each term on the right-hand side depends on X . When the perturbation and quantum noise are weak, the deviation of the system state X from the limit cycle is small and of the order of O( √ ) because the limit cycle is exponentially stable and the system state is subjected to Gaussianwhite noise. Thus, in the lowest-order approximation, we can approximate the state X by a state X 0 (φ) on the limit cycle as X (t ) = X 0 (φ(t )) + O( √ ). We then obtain an Ito SDE for the phase variable φ, which is correct up to O( ) in the drift term and up to O( √ ) in the noise intensity, where represents the effect of the perturbation on φ, represents the effect of the quantum noise on φ, and g(φ) = 1 2 Tr{G(X 0 (φ)) T (∇ T ∇ | X =X 0 (φ) )G(X 0 (φ))} (B6) represents a term arising from the change of the variables, respectively.
We denote the gradient vector (PSF) and Hessian matrix of the phase function (X ) evaluated at X = X 0 (φ) on the limit cycle as Z(φ) = ∇ | X =X 0 (φ) and Y (φ) = ∇ T ∇ | X =X 0 (φ) , respectively. The components of the PSF and Hessian matrix ∇ T ∇ | X =X 0 (φ) are explicitly given by for i, j = 1, 2, respectively. It is well known in the classical phase reduction theory [3,[7][8][9] that Z(φ) is given by a 2π -periodic solution to the following adjoint equation and normalization condition: It is also known [53,54] that the Hessian matrix Y (φ) of the phase function, evaluated at X = X 0 (φ) on the limit cycle, is given by a 2π -periodic solution to a differential equation which satisfies a constraint In the above equations, J(φ) ∈ R 2×2 is a Jacobian matrix of F(X ) at X = X 0 (φ) and K(φ) ∈ R 2×2×2 is a third order tensor, respectively, whose components are given by and the matrix components of the product Z(φ) • K(φ) ∈ R 2×2 are given by 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 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 [70]. As Eq. (B13) is an Ito SDE, using the property of the Wiener process, the noise term can be rewritten as where h(φ) = 2 i=1 (h(φ)) 2 i 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 O( 2 ) in the drift term and O( ) in the noise intensity, respectively. Therefore, the error in the mean of φ from the true value grows with time as O( 2 t ), and the error in the variance of φ grows as O( 2 t ). Thus, these errors in the phase dynamics remain O( ) up to t = O(1/ ).

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 φ in the lowest-order approximation, it has an additional function g(φ) in the drift term, which adds tiny periodic fluctuations to the deterministic part. By further introducing a new phase ψ that is only slightly different from φ, we can eliminate this term and obtain a simpler SDE, where f (ψ, t ) = Z(ψ ) · q(ψ, t ), W (t ) is a one-dimensional Wiener process, and h(ψ ) is a 2π -periodic function of ψ.
Here, the new phase ψ is defined from φ by a near-identity transform as φ = ψ + n(ψ ), where n(ψ ) is a 2π -periodic function with n(0) = 0. Using this transformation, the additional term g(φ) in Eq. (5) can be renormalized into the frequency term as whereω is the effective frequency of the system. As is assumed to be sufficiently small, the transformation between the two variables φ and ψ 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 O( ) 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. (C1) from the original phase equation, Eq. (5) or (B13), by using the near-identity transform [71]. In Eq. (B13), the function g(φ) contains the Hessian matrix Y (φ) of (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 φ to a new phase ψ, φ = ψ + n(ψ ), where the transformation function n(ψ ) is a smooth 2πperiodic function of ψ satisfying n(0) = 0, and assume that ψ obeys an Ito SDE of the form in the lowest-order approximation, which does not contain a term corresponding to g(φ) but has a small shift in the frequency. From this SDE, we obtain an Ito SDE for φ by using the Ito formula as where we omitted the tiny terms of O( 2 ) in the drift term and O( 3/2 ) in the noise intensity. The replacement of φ by ψ in the functions f and h also results in errors of O( 2 ) and O( 3/2 ) in the drift term and noise intensity, respectively, which can also be neglected. The above equation coincides with the original Eq. (B13) if n(ψ ) satisfies + ωn (ψ ) = g(φ). (C6) As g(φ) = g(ψ ) + O( ), the equation for n(ψ ) is obtained at the lowest order as which gives where n(0) = 0 is used. Moreover, as n(ψ ) is 2π -periodic, n(2π ) = n(0) = 0 should hold, which determines the frequency shift as = 2π 2π 0 dψ g(ψ ).
Thus, by introducing the near-identity transform, we obtain an averaged phase equation whereω = ω + is a renormalized, effective frequency. This corresponds to Eq. (C1). The orders of errors caused by the above near-identity transformation are O( 2 ) in the drift term and O( 3/2 ) in the noise intensity. Therefore, the phase equations in Eq. (B13) and (C10) are equally correct in the lowest-order approximation and valid up to t = O(1/ ).
The frequency shift can be evaluated by numerically calculating the Hessian matrix of (X ) in g(ψ ) and integrating Eq. (C9), or alternatively by measuringω by numerically evolving the SDE in Eq. (3) or Eq. (4) without perturbations.
In the examples used in the main text, the frequency shift 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ω may not be required (only the frequency difference matters). In such cases, one may simply assumeω ≈ ω and avoid the calculation of .

033012-11
In the case of strong squeezing with δ = 1, the rescaled system Hamiltonian and the perturbation Hamiltonian are given by H = − a † a +iη (a 2 e −iθ −a †2 e iθ ), respectively. The deterministic part F(X ) gives an asymmetric limit cycle when η > 0, which is difficult to solve analytically. However, we can still obtain the limit cycle X 0 (φ) numerically and use it to evaluate the PSF Z(φ), Hessian matrix Y (φ), and the noise intensity G(φ), and use these quantities in the phase equation. The PSF Z(φ) can be numerically calculated by the adjoint method, and the Hessian matrix Y (φ) 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.