Quantum asymptotic phase reveals signatures of quantum synchronization

Synchronization of quantum nonlinear oscillators has attracted much attention recently. To characterize the quantum oscillatory dynamics, we recently proposed a fully quantum-mechanical definition of the asymptotic phase, which is a key quantity in the synchronization analysis of classical nonlinear oscillators (Kato and Nakao 2022 Chaos 32 063133). In this work, we further extend this theory and introduce multiple asymptotic phases using the eigenoperators of the adjoint Liouville superoperator of the quantum nonlinear oscillator associated with different fundamental frequencies. We analyze a quantum van der Pol oscillator with Kerr effect in the strong quantum regime and show that the system has several different fundamental frequencies. By introducing order parameters and power spectra in terms of the associated quantum asymptotic phases, we reveal that phase locking of the system with a harmonic drive at several different frequencies, an explicit quantum signature observed only in the strong quantum regime, can be interpreted as synchronization on a torus rather than a simple limit cycle.

Synchronization of quantum nonlinear oscillators has attracted much attention recently. To characterize the quantum oscillatory dynamics, we recently proposed a fully quantum-mechanical definition of the asymptotic phase, which is a key quantity in the synchronization analysis of classical nonlinear oscillators [1]. In this work, we further extend this theory and introduce multiple asymptotic phases using the eigenoperators of the adjoint Liouville superoperator of the quantum nonlinear oscillator associated with different fundamental frequencies. We analyze a quantum van der Pol oscillator with Kerr effect in the strong quantum regime and show that the system has several different fundamental frequencies. By introducing order parameters and power spectra in terms of the associated quantum asymptotic phases, we reveal that phase locking of the system with a harmonic drive at several different frequencies, an explicit quantum signature observed only in the strong quantum regime, can be interpreted as synchronization on a torus rather than a simple limit cycle.

I. INTRODUCTION
Spontaneous rhythmic oscillations and synchronization are phenomena ubiquitously observed in various fields of science and technology [2][3][4][5][6][7]. Owing to the recent progress in nanotechnology, synchronization in micro-and nano-scale devices has been realized experimentally [8][9][10][11][12] and experimental demonstrations of quantum phase synchronization in spin-1 atoms [13], in nuclear spin systems [14], and on the IBM Q system [15] have been reported. A number of theoretical investigations have also been performed to reveal novel quantum signatures in synchronization [1, 62].
In the strong quantum regime where only a small number of energy states participate in the system dynamics, the discrete nature of the energy spectrum can give rise to explicit quantum signatures, such as multiple phase locking at several different frequencies [16] and synchronization blockade [17,18]. In characterizing quantum synchronization, several definitions of the system's oscillation phase have been proposed [16,[19][20][21][22][23][43][44][45][46].
In classical mechanics, nonlinear dissipative systems exhibiting spontaneous rhythmic dynamics can be modeled as limit-cycle oscillators. The asymptotic phase [2][3][4][5][6] is a fundamental quantity for the analysis of synchronization, which is defined by the oscillator's vector field and increasing with a constant frequency in the basin of the limit-cycle attractor. It provides the basis for phase reduction theory [2][3][4][5][6][7], a standard dimensionality-reduction method to derive phase equations describing weakly-coupled oscillators.
We recently introduced the asymptotic phase for quantum oscillatory systems [1] by extending the definition of the asymptotic phase for classical stochastic oscillatory systems [47] based on the Koopman operator theory [48]. This quantum asymptotic phase is defined fully quantum-mechanically in terms of the eigenoperator of the system's adjoint Liouville superoperator associated with the slowest decaying mode and hence applicable even in the strong quantum regimes where we cannot rely on the limit-cycle trajectory in the classical limit as in the semiclassical regime [27,28,49].
In this study, we further extend this definition and introduce multiple quantum asymptotic phases in terms of the eigenoperators with several different fundamental frequencies, i.e., with the eigenvalues possessing the smallest absolute imaginary part in the individual branches of the eigenvalues near the imaginary axis, which play dominant roles in synchronization dynamics. As an example, we analyze a quantum van der Pol oscillator with Kerr effect and show that it possesses several dominant eigenvalues with different fundamental frequencies in the strong quantum regime. By introducing the order parameters and power spectra in terms of the associated quantum asymptotic phases with respective fundamental frequencies, we reveal a torus-like structure of multiple-frequency phase locking of the system with a harmonic drive [16], which is observed only in the strong quantum regime.

II. ASYMPTOTIC PHASES FOR QUANTUM OSCILLATORY SYSTEMS
In Ref.
[1], we proposed a definition of the asymptotic phase for quantum oscillatory systems in terms of the eigenoperator of the adjoint Liouville superoperator associated with the slowest decaying mode, inspired by the definition of the asymptotic phase for classical stochastic oscillators [47,48]. In this section, we briefly review this definition and extend it to introduce multiple quantum asymptotic phases by using the eigenoperators associated with several eigenvalues with different fundamental frequencies.

A. Quantum master equation
We consider a quantum nonlinear oscillatory system with a single degree of freedom coupled to reservoirs, and assume that the interactions of the system with the reservoirs are instantaneous and Markovian approximation can be employed. The time evolution of the system's density operator ρ is described by a quantum master equation [50,51], where L is a Liouville superoperator representing the evolution of ρ, H is a system Hamiltonian, C j is a coupling operator between the system and jth reservoir (j = 1, . . . , n), [A, B] = AB − BA is the commutator, D[C]ρ = CρC † − (ρC † C + C † Cρ)/2 is the Lindblad form ( † denotes Hermitian conjugate), and the reduced Planck's constant is set as = 1.
Introducing an inner product X, Y tr = Tr (X † Y ) of linear operators X and Y , we define the adjoint superoperator L * of L satisfying L * X, Y tr = X, LY tr where D * [C]X = C † XC − (XC † C + C † CX)/2 is the adjoint Lindblad form. This L * describes the evolution of an observable F asḞ In the Schrödinger picture, the density operator ρ evolves as in Eq.
(1) while the observation operator F does not vary with time, and in the Heisenberg picture, F evolves as in Eq. (3) while ρ remains constant. The expectation value of F with respect to ρ is kept the same in both pictures (note that ρ is self-adjoint). We assume that the Liouville superoperator L has a set of eigensystem (an eigenvalue and right and left eigenop- for k, l = 0, 1, 2, . . ., where the overline indicates complex conjugate [52]. We assume that among {λ k } k≥0 , one eigenvalue λ 0 is always 0, which corresponds to the stationary state ρ 0 of the system satisfying Lρ 0 = 0, and all other eigenvalues have negative real parts. This assumption also means that the system has no decoherence-free subspace [53]. Considering the system's oscillatory dynamics, we also assume that the eigenvalues of L with the largest nonvanishing real part (i.e., with the slowest decay rate) are given by a complex-conjugate pair and denote them as Λ 1 and Λ 1 , where Ω 1 = Im Λ 1 gives the frequency of the slowest decay mode. One may also choose Ω 1 = Im Λ 1 , which reverses the direction of the asymptotic phase. The sign of Ω 1 can be chosen arbitrarily and will be fixed later. Figure 1 shows typical examples of the eigenvalues of L of the quantum vdP oscillator (see the next section for details). The eigenvalues form several branches. In the semiclassical regime ( Fig. 1(b)), the rightmost branch is far apart from the other branches and thus it is the only dominant branch, while in the strong quantum regime ( Fig. 1(a)), the branches are closer to each other and the system possess several comparably important branches. It can also be seen in the semiclassical regime ( Fig. 1(b)) that the imaginary part of the eigenvalues are approximately integer multiples of the fundamental frequency (the smallest absolute imaginary part with the slowest decay rate) in each individual branch.

B. Phase-space representation
The density operator ρ can also be transformed into a quasiprobability distribution in the phase space [50,51,54]. We use the P -representation and describe ρ as where |α is a coherent state specified by a complex value α, or equivalently by a complex vector α = (α, α) T , p(α) is a quasiprobability distribution of α, dα = dαdα, and the integral is taken over the entire complex plane. The observable F is also transformed into a function in the phase space as where the operator F is arranged in the normal order [50,51,54]. By introducing the L 2 inner product g(α), h(α) α = g(α)h(α)dα of two functions g(α) and h(α), the expectation value of F with respect to ρ is expressed as The time evolution of p(α) corresponding to Eq. (1) is described by a partial differential equation where the differential operator L α satisfies Lρ = L α p(α)|α α|dα. The explicit form of L α can be calculated from Eq. (1) by using the standard calculus for the phase-space representation [50,51,54]. The corresponding evolution of f (α) in the Heisenberg picture is given by where the differential operator L + α is the adjoint of L α with respect to the L 2 inner product, i.e., L + α g(α), h(α) α = g(α), L α h(α) α , which satisfies L + α f (α) = α|L * F |α . The differential operator L α also has a set of eigensystem (an eigenvalue and right and left eigenfunctions) This eigensystem has a one-to-one correspondence with Eq. (5), where the eigenvalues {λ k } k≥0 are the same as those of L; the eigenfunctions u k and v k of L α and L + α are related to the eigenoperators U k and V k of L and L * as C. Quantum asymptotic phase associated with the slowest decaying mode In Ref.
[1], we defined the quantum asymptotic phase function Φ 1 of the system state ρ as the argument (polar angle) of the expectation of the eigenoperator V 1 associated with the eigenvalue Λ 1 satisfying We showed that this quantum asymptotic phase yields appropriate phase values, namely, it always increases with a constant frequency Ω 1 asΦ 1 (ρ) = Ω 1 with the evolution of ρ even in the strong quantum regime and reproduces the conventional asymptotic phase in the semiclassical regime [1]. We note that Φ 1 (ρ) is not defined when V 1 = 0. This quantum asymptotic phase is a natural extension of the asymptotic phase for stochastic limit-cycle oscillators defined in terms of the slowest decaying eigenfunction of the backward Kolmogorov (Fokker-Planck) operator [47] from the Koopman operator viewpoint [48].

D. Multiple quantum asymptotic phases associated with different fundamental frequencies
In this study, we further extend the definition in the previous subsection and introduce multiple quantum asymptotic phases associated with several different fundamental frequencies by using the principal eigenvalues on different eigenvalue branches near the imaginary axis, and use them to characterize quantum signatures of synchronization in the strong quantum regime.
As shown in Fig. 1(a), in the strong quantum regime, multiple branches of the eigenvalues with different fundamental frequencies can exist near the imaginary axis, suggesting that not only the eigenvalue Λ 1 on the rightmost branch but also the eigenvalues with the slowest decay rates on the other branches play important roles; for comparison, see Fig. 1(b) for a typical example of the eigenvalue Λ 1 in the semiclassical regime, where only a single dominant branch of eigenvalues exist. We denote these eigenvalues by Λ 1 , Λ 2 , Λ 3 , . . . and their imaginary parts, i.e., fundamental frequencies, by Ω j = Im Λ j (j ≥ 1), and call {Λ j } j≥1 the principal eigenvalues. Here, the first principal eigenvalue Λ 1 and the fundamental frequency Ω 1 are those introduced in the previous subsection. These principal eigenvalues on the individual branches are shown by red dots in Fig. 1(a).
In a similar manner to the phase function Φ 1 , we introduce the j-th quantum asymptotic phase function Φ j (j = 2, 3, . . .) of ρ as the argument of the P -representation of the eigenoperator V j associated with the principal eigenvalue Λ j satisfying L * V j = Λ j V j , namely, Note that Φ j (ρ) is not defined when V j = 0. Since we obtain V j (t) = V j (0)e Λj t and therefore Thus, Φ j always increases with a constant frequency Ω j with the evolution of ρ and plays the role of the asymptotic phase for any j = 1, 2, 3, . . .. It should be noted that Φ j (ρ) cannot be defined for the stationary state ρ 0 , because ρ 0 is the eigenfunction of the Liouville superoperator L with the eigenvalue λ 0 = 0 and hence V j = ρ 0 , V j tr = 0 from the biorthogonality in Eq. (5). As will be shown later, for the stationary state of a periodically driven system, V j takes a non-zero value and the above definition of the phase functions can be used for the analysis of phase locking.
We stress that in the classical deterministic limit, Φ j (j ≥ 2) is not independent from Φ 1 and does not provide additional information, because Ω j (j = 1) is equal to Ω 1 according to the Koopman operator theory [55][56][57][58]. In the semiclassical regime, Ω j (j = 1) differs only slightly from Ω 1 as shown in Fig. 1(b) due to the effect of weak quantum noise and the corresponding decay rate for j ≥ 2 is much larger, so the corresponding mode does not play an important role. However, as we see in the next section, Φ j (j ≥ 2) is distinctly different from Φ 1 and the corresponding decay rate Re Λ j (j ≥ 2) is comparable to Re Λ 1 in the strong quantum regime, hence the phase function Φ j (j ≥ 2) yields independent information from Φ 1 .

III. QUANTUM VAN DER POL OSCILLATOR WITH KERR EFFECT
A. Eigenvalues of the Liouville superoperator As an example, we consider a quantum van der Pol oscillator with Kerr effect. The master equation is given by [1, 16,28] where H = ω 0 a † a + Ka †2 a 2 , ω 0 is the natural frequency of the oscillator, K is the Kerr parameter, and γ 1 and γ 2 are the decay rates for negative damping and nonlinear damping, respectively. We added a subscript 0 to the Liouville operator as we will introduce an additional external drive later. We first consider a strong quantum regime with large γ 2 and K, where only a small number of energy states participate in the system dynamics and the discrete nature of the energy spectrum plays important roles in the dynamics. We set the parameters as γ 1 = 0.1 and (ω 0 , γ 2 , K)/γ 1 = (300, 4, 100), which are the same as in our previous study [1]. In the numerical calculation, we approximately truncated the density operator as a large-dimensional N ×N matrix and mapped it into a N 2 -dimensional vector in the double-ket notation [59]. We can then represent the Liouville operator by a N 2 × N 2 matrix and obtain the asymptotic phase in Eq.(13) from the eigensystem of this matrix [1]. Figure 1(a) shows the eigenvalues of L 0 near the imaginary axis obtained numerically. We can identify several branches of the eigenvalues near the imaginary axis characterized by the principal eigenvalues Λ 1 , Λ 2 , Λ 3 , . . . whose fundamental frequencies Ω 1 = Im Λ 1 , Ω 2 = Im Λ 2 , Ω 3 = Im Λ 3 , . . . are different from each other. Moreover, their decay rates, which are characterized by Re Λ 1 , Re Λ 2 , Re Λ 3 , . . . are comparable to each other. This indicates that not only the quantum asymptotic phase Φ 1 associated with the principle eigenvalue Λ 1 of the slowest decaying mode [1] but also those associated with the other principal eigenvalues Λ 2 , Λ 3 , . . . can play important roles in the dynamics. We choose a negative value for each Ω j so that the corresponding Φ j increases from 0 to 2π in the counterclockwise direction.
In Ref. [16], it is shown that, in the strong quantum regime, |m + 1 m| is an approximate eigenoperator of L 0 with the eigenvalue where m = 0, 1, 2, . . ., namely, L −1 0 |m + 1 m| ≈ λ −1 m |m + 1 m|. As shown in Fig. 1(a), these eigenvalues correspond to the principal eigenvalues of L 0 , i.e., Λ j ≈ λ j−1 , and thus V j ≈ |j j − 1| (j = 1, 2, 3, 4). This indicates that that periodic transitions between the adjacent discrete energy states play important roles in the oscillatory behavior in the strong quantum regime, where the difference in the transition frequencies arises due to the unequal spacing of the energy levels characterized by the Kerr parameter in Eq. (18).
As a comparison, we next consider the semiclassical regime where γ 2 and K are sufficiently small and the semiclassical approximation can be taken [27]. We set the parameters as γ 1 = 1 and(ω 0 , γ 2 , K)/γ 1 = (0.1, 0.05, 0.025), which are the same as those used in Ref. [1]. In this regime, we can approximate Eq. (9) by a quantum Fokker-Planck equation for p(α) and the system can be approximately described as a Stuart-Landau oscillator (Hopf normal form) subjected to small quantum noise [1, 16,28]. (Therefore, the conventional quantum vdP oscillator is also called a quantum Stuart-Landau oscillator recently [37,39] and a more appropriate model of the quantum van der Pol oscillator has also been proposed [39,40].) Figure 1(b) shows the eigenvalues of L 0 near the imaginary axis, where the principal eigenvalues are shown by red dots on the individual branches; Λ 1 = µ 1 + iΩ 1 (µ 1 < 0) is on the rightmost light-blue branch. In contrast to the strong quantum regime, the rightmost branch of the eigenvalues, approximately given by a parabolaλ n = iΩ 1 n + µ 1 n 2 (n = 0, ±1, ±2, . . .) passing through Λ 1 , is isolated from other branches of eigenvalues with faster decay rates; the relative decay rate Re Λ 2 /Re Λ 1 ≈ −0.996/(−0.0736) ≈ 13.5 is more than three times larger than Re Λ 2 /Re Λ 1 ≈ −0.65/(−0.15) ≈ 4.33 in the strong quantum regime in Fig. 1(a), indicating that only the rightmost branch is dominant in the semiclassical regime. Also, the fundamental frequencies of the other branches, defined as the smallest absolute imaginary part of the eigenvalues, are approximately equal to Ω 1 ; the small differences in the fundamental frequencies arise from small quantum noise. Thus, it is sufficient to consider only Λ 1 and introduce a single phase function Φ 1 in this regime.
The system in the classical limit, i.e., in the limit of vanishing quantum noise, is described by the drift term of the approximate quantum Fokker-Planck equation for p(α), which represents the Stuart-Landau oscillator and possesses a stable limit-cycle solution [1, 16,28]. Figure 1(c) shows a schematic diagram of the eigenvalues of the differential operator L + α , which are equivalent to those of the backward Liouville operator L * 0 , in the classical limit. They are given in the form λ c = mκ 1 + inω 1 (m = 0, 1, 2, . . . , n = 0, ±1, ±2), where κ 1 is the real part of the largest negative eigenvalue and ω 1 is the pure-imaginary eigenvalue with the smallest absolute imaginary part. Thus, Φ j (j ≥ 2) is identical to Φ 1 and does not provide additional information, because Ω j (j = 2) is identical to Ω 1 .
The above results suggest that the phase function Φ j yields independent information from Φ 1 only in the strong quantum regime. The existence of several dominant fundamental frequencies in the strong quantum regime suggests that the system behaves like a torus rather than a limit cycle with a single fundamental frequency, and that we need to consider the phase functions Φ 2 , Φ 3 , . . . associated with Λ 2 , Λ 3 , . . . in addition to Φ 1 and Λ 1 . Figure 2 shows a schematic picture of the torus behavior of the system with two fundamental frequencies.

B. Asymptotic phase functions in the strong quantum regime
In this subsection, we examine the validity of the quantum asymptotic phase functions in the strong quantum regime. Figure 3 shows the asymptotic phase functions Φ j (α) for j = 1, 2, 3, and 4 of the pure coherent state α on the complex plane (x = Re α, p = Im α). The parameters are the same as in Fig. 1(a). These asymptotic phase functions look similar, but they are associated with different fundamental frequencies and slightly different from each other near the origin as shown in the enlarged figures. Thus, they capture different oscillatory dynamics of the system. Though not shown, we may also draw similar asymptotic phase functions Φ j for j ≥ 5, which are less dominant. Here, the asymptotic phase functions for different values of j look similar to each other since they are identical in the classical limit [55][56][57][58], but the differences in their eigenfrequencies brought by the strong quantum effect take an important role in analyzing strong quantum signatures in synchronization.
To demonstrate that these quantum asymptotic phase functions yield appropriate phase values, we consider free oscillatory relaxation of ρ from a pure coherent initial state ρ = |α 0 α 0 | with α 0 = 1 at t = 0 and measured the evolution of the expectation values of V j and their arguments, i.e., their asymptotic phases Φ j (ρ) = arg V j (j = 1, 2, 3, 4), as well as those of the annihilation operator a for comparison. It is noted that, though the initial condition is a pure coherent state, the system state quickly becomes mixed due to the interaction with the reservoirs.
We can confirm that each V j exhibits exponentially damped harmonic oscillations as shown in Fig. 4(a, c, e, g), and correspondingly, each Φ j (ρ) gives constantly varying phase values with frequency Ω j as shown in Fig. 4(b, d, f, h), verifying the validity of the definition of the quantum asymptotic phases. It is noted that different asymptotic phases independently capture different oscillation modes in the evolution of ρ. In contrast, a shown in Fig. 4(i) exhibits more complex oscillatory dynamics and the simple argument arg a does not vary constantly with time as shown in Fig. 4(j). Thus, arg a cannot be considered the asymptotic phase, although it is often used to define power spectra in the analysis of quantum synchronization. It is noted that the asymptotic phase is quantitatively different from the geometric angle also in the classical limit when the Kerr effect is added. In the strong quantum regime, due to the strong Kerr effect, the difference between the asymptotic phase Ω j and simple argument arg a can be observed more clearly than in the semiclassical regime or in the classical limit  C. Revealing multiple phase-locking structure We now analyze quantum synchronization of the vdP oscillator with a harmonic drive in the strong quantum regime using the proposed asymptotic phases.
The master equation in the rotating frame of the frequency ω d of the harmonic drive iṡ where L 0 is given by Eq. (17) with H = −∆a † a + Ka †2 a 2 and the harmonic drive is represented by is the frequency detuning of the harmonic drive from the oscillator and E is the strength of the harmonic drive [16,28]. We use the same parameters as in Fig. 1(a) and vary the detuning parameter ∆ by controlling ω d while keeping the natural frequency ω 0 fixed. In Ref. [16], Lörch et al. showed that this system under strong Kerr effect exhibits multiple phase locking to the harmonic drive at several detuning frequencies ∆ = 2mK (m = 0, 1, 2, . . .), i.e., at ω d = ω 0 + 2mK, observed as multiple sharp Arnold tongues, while the corresponding classical system exhibits only a single broad Arnold tongue. It is stressed that K is the Kerr parameter, hence this is not the ordinary higher-harmonic phase locking at the frequencies ω d = (p/q)ω 0 with p, q = 1, 2, 3., ... (p/q = 1); the multiple Arnold tongues are related to the periodic transitions between the adjacent discrete energy states. It is also noted that such multiple phase locking is robust to thermal noise under the finite temperature environment (See Supplemental Material in Ref. [16]).
In Ref. [16], the following order parameter S a and power spectrum P a defined using the annihilation operator a are used to analyze the system: where the expectation is taken with respect to the steady-state density operator obtained from Eq. (19).
Here, in addition to these quantities, we introduce the order parameters and the power spectra in terms of the quantum asymptotic phases (or, more precisely, the corresponding principal Koopman eigenoperators) defined in this study. They are defined using the left eigenoperators V j of L 0 as where V † j (τ ) = e L * 0 τ V † j (0) = e Λj τ V † j (0) and V † j (0) = V † j (j = 1, 2, 3 and 4). Here, |S a | and |S j | quantify the phase coherence of the system, while θ a and θ j characterize the averaged phases of the system relative to the harmonic drive. We note here that we can use V j , which is defined in the original frame, in the rotating frame with the external drive since the eigenoperators V j changes only its phase factor by the coordinate rotation, i.e., e iω d a † at V j e −iω d a † at = V j e iΩt with a constant Ω (see Appendix A). The resulting power spectra in the rotating frame are simply shifted from the power spectra in the original frame by −Ω. Using these two quantities in Eqs. (22) and (23), we can analyze both the phase coherence and frequency characteristics of the oscillator in quantum synchronization. Figure 5 shows the dependence of the order parameters |S a | and |S j | on the detuning ∆ and strength E of the harmonic drive. In Fig. 5(a) showing |S a |, several Arnold tongues representing phase locking of the oscillator at different frequencies are observed [16]. Figures 5(b)-(e) show |S j | for j = 1, 2, 3 and 4, respectively. It is remarkable that the Arnold tongues in Fig. 5(a) are clearly decomposed into individual Arnold tongues around ∆ = 2(j − 1)K in Figs. 5(b)-5(e). This indicates that the order parameters in terms of the asymptotic phases can capture the phase locking dynamics at ω d = ω 0 + 2(j − 1)K for different j individually.
Similarly, Fig. 6 shows the power spectra P a and P j for ∆ = 0. Multiple peaks of P a in Fig. 6(a), which indicate multiple phase locking of the oscillator to the harmonic drive, are clearly decomposed into individual peaks around ω = ∆ − 2(j − 1)K in Figs. 6(b)-(e) for P j (j = 1, 2, 3, and 4). The Arnold tongue and power spectrum are sharper when the decay rate characterized by Re Λ j is smaller. Though not shown, we can also detect even smaller tongues and peaks with j ≥ 5.
The above results reveal that, in the strong quantum regime, the system behaves like a torus with several fundamental frequencies and each of the associated oscillating mode individually exhibits phase locking to the harmonic drive at the respective frequency [60], resulting in the multiple Arnold tongues and spectral peaks. Such a torus-like behavior cannot be observed in the semiclassical regime where the system behaves like a noisy limit-cycle oscillator with a single fundamental frequency, because only the principal eigenvalue Λ 1 is dominant. In contrast, in the strong quantum regime, due to the effect of quantum noise, visible differences between the imaginary part of the principal eigenvalues Λ j (j = 1, 2, 3, and 4) arise and also the different branches of eigenvalues become closer to each other, resulting in the torus-like behavior.

IV. CONCLUSION
In this study, we defined multiple quantum asymptotic phases of quantum nonlinear oscillators in terms of the eigenoperators of the adjoint Liouville superoperator associated with several fundamental frequencies from the Koopman operator viewpoint, which extends our previous definition of the quantum asymptotic phase associated with the slowest decaying mode of the system. We introduced the order parameters and power spectra in terms of the proposed quantum asymptotic phases and applied them to the analysis of a quantum van der Pol oscillator with Kerr effect. We successfully revealed that the multiple phase locking of the system with a harmonic drive at several different frequencies [16], which is an explicit quantum signature observed only in the strong quantum regime, can be interpreted as synchronization on a torus rather than on a simple limit cycle.
Though not discussed in this paper, it will also be interesting to introduce the amplitude functions in addition to the quantum asymptotic phases by extending the definition for classical stochastic oscillators [48,61], which can be defined using the eigenfunctions associated with the eigenvalues on the real axis. The phase-amplitude description that uses both the phases and amplitudes may be applied to the analysis of quantum complete synchronization [19]. Also, investigating the mutual synchronization between two quantum nonlinear oscillators [17,29,62] by using quantum asymptotic phase functions is also a future work. We may be able to introduce a new measure for quantum synchronization between two oscillators, e.g., the normalized correlator [20], by using quantum asymptotic phase functions. It may be also interesting to extend our definition of the asymptotic phase for quantum synchronization in non-Markov open quantum systems [63].
We expect that the proposed definition of the quantum asymptotic phases will serve as a fundamental tool for analyzing strong quantum effects in synchronization [16,17] and will be useful for future applications of quantum synchronization in the growing fields of quantum technologies.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available upon reasonable request from the authors.
Appendix A: Proof of e iω d a † at Vje −iω d a † at = Vje iΩt The time evolution ofṼ j (t) = e iω d a † at V j e −iω d a † at is given by Because the adjoint superoperator L * 0 couples only the elements of the form |n + k m + k| to the elements of the form |n + k m + k | for an arbitrary pair of integers {n, m} and integers k, k [66,67], the eigenoperators V j can be explicitly written as v j (n)|n n + k|. (A2)