Quantum asymptotic phases reveal signatures of quantum synchronization

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 [51, 52],

$$\begin{align} \dot{\rho} = \mathcal{L}\rho = -i[H, \rho] + \sum_{j=1}^{n} \mathcal{D}[C_{j}]\rho, \end{align} \tag{ 1 }$$

where $\mathcal{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,\ldots,n)$, $[A,B] = AB-BA$ is the commutator, $\mathcal{D}[C]\rho = C \rho C^{\dagger} - (\rho C^{\dagger} C + C^{\dagger} C \rho)/2$ is the Lindblad form ($\dagger$ denotes Hermitian conjugate), and the reduced Planck's constant is set as $\hbar = 1$.

Introducing an inner product $\langle{X,Y}\,\rangle_{tr} = \textrm{Tr}\hspace{0.07cm}{( X^{\dagger}Y )}$ of linear operators X and Y, we define the adjoint superoperator $\mathcal{L}^{*}$ of $\mathcal{L}$ satisfying $\langle{\mathcal{L}^{*}X,Y}\,\rangle_{tr} = \langle{X,\mathcal{L}Y}\,\rangle_{tr}$:

$$\begin{align} \mathcal{L}^{*} X = i[H, X] + \sum_{j=1}^{n} \mathcal{D}^{*}[C_{j}]X, \end{align} \tag{ 2 }$$

where $\mathcal{D}^{*}[C]X = C^{\dagger} X C - (X C^{\dagger} C + C^{\dagger} C X)/2$ is the adjoint Lindblad form. This ${\cal L}^{*}$ describes the evolution of an observable F as:

$$\begin{align} \dot{F} = {\cal L}^{*} F. \end{align} \tag{ 3 }$$

In the Schrödinger picture, the density operator ρ evolves as in equation (1) while the observation operator F does not vary with time, and in the Heisenberg picture, F evolves as in equation (3) while ρ remains constant. The expectation value:

$$\begin{align} \langle F \rangle = \mbox{Tr}(\rho F) = \langle{\rho, F}\rangle_{tr} \end{align} \tag{ 4 }$$

of F with respect to ρ is kept the same in both pictures (note that ρ is self-adjoint).

We assume that the Liouville superoperator $\mathcal{L}$ has a set of eigensystem (an eigenvalue and right and left eigenoperators) $\{\lambda_{k}, U_{k}, V_{k} \}$ satisfying:

$$\begin{align} \mathcal{L} U_{k} = \lambda_{k} U_{k}, \quad \mathcal{L}^{*} V_{k} = \overline{\lambda_{k}} V_{k}, \quad \langle{V_{k}, U_{l}}\rangle_{tr} = \delta_{kl}, \end{align} \tag{ 5 }$$

for $k, l = 0, 1, 2, \ldots$, where the overline indicates complex conjugate [53]. We also assume that among $\{\lambda_k \}_{k \geqslant 0}$, one eigenvalue λ₀ is always 0, which corresponds to the stationary state ρ₀ of the system satisfying ${\cal L} \rho_0 = 0$, and all other eigenvalues have negative real parts. This assumption also means that the system has no decoherence-free subspace [54]. Considering the system's oscillatory dynamics, we also assume that the eigenvalues of ${\cal L}$ with the largest non-vanishing real part (i.e. with the slowest decay rate) are given by a complex-conjugate pair and denote them as Λ₁ and $\overline{\Lambda_1}$, where $\Omega_1 = \mbox{Im}\ \overline{\Lambda_1}$ gives the frequency of the slowest decay mode. One may also choose $\Omega_1 = \mbox{Im}\ \Lambda_1$, which reverses the direction of the asymptotic phase. The sign of Ω₁ can be chosen arbitrarily and will be fixed later.

Figure 1 shows typical examples of the eigenvalues of ${\cal L}$ of the quantum vdP oscillator (see the next section for details). The eigenvalues form several branches. In the semiclassical regime (figure 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 (figure 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 (figure 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.

Zoom In Zoom Out Reset image size

The density operator ρ can also be transformed into a quasiprobability distribution in the phase space [51, 52, 55]. We use the P-representation and describe ρ as:

$$\begin{align} \rho = \int p({\boldsymbol{\alpha}}) | \alpha \rangle \langle \alpha | d {\boldsymbol{\alpha}}, \end{align} \tag{ 6 }$$

where $| \alpha \rangle$ is a coherent state specified by a complex value α, or equivalently by a complex vector $\boldsymbol{\alpha} = (\alpha, \overline{\alpha})^{T}$, $p({\boldsymbol{\alpha}})$ is a quasiprobability distribution of α , $d{\boldsymbol{\alpha}} = d\alpha d\overline{\alpha}$, and the integral is taken over the entire complex plane. The observable F is also transformed into a function in the phase space as:

$$\begin{align} f({\boldsymbol{\alpha}}) = \langle \alpha | F | \alpha \rangle, \end{align} \tag{ 7 }$$

where the operator F is arranged in the normal order [51, 52, 55]. By introducing the L² inner product $\langle{g(\boldsymbol{\alpha}), h(\boldsymbol{\alpha})}\rangle_{\boldsymbol{\alpha}} = \int \overline{g(\boldsymbol{\alpha}) } h(\boldsymbol{\alpha}) d \boldsymbol{\alpha}$ of two functions $g(\boldsymbol{\alpha})$ and $h(\boldsymbol{\alpha})$, the expectation value of F with respect to ρ is expressed as:

$$\begin{align} \langle F \rangle = \mbox{Tr} (\rho F) = \int d{\boldsymbol{\alpha}} p({\boldsymbol{\alpha}}) f({\boldsymbol{\alpha}}) = \langle\,{p(\boldsymbol{\alpha}), f(\boldsymbol{\alpha})}\rangle_{\boldsymbol{\alpha}}. \end{align} \tag{ 8 }$$

The time evolution of $p({\boldsymbol{\alpha}})$ corresponding to equation (1) is described by a partial differential equation:

$$\begin{align} \partial_t{p}({\boldsymbol{\alpha}}) = L_{\boldsymbol{\alpha}} p({\boldsymbol{\alpha}}), \end{align} \tag{ 9 }$$

where the differential operator ${L}_{\boldsymbol{\alpha}}$ satisfies $\mathcal{L}\rho = \int {L}_{\boldsymbol{\alpha}} p({\boldsymbol{\alpha}}) | \alpha \rangle \langle \alpha | d{\boldsymbol{\alpha}}$. The explicit form of ${L}_{\boldsymbol{\alpha}}$ can be calculated from equation (1) by using the standard calculus for the phase-space representation [51, 52, 55]. The corresponding evolution of $f({\boldsymbol{\alpha}})$ in the Heisenberg picture is given by:

$$\begin{align} \partial_t f({\boldsymbol{\alpha}}) = {L}^{+}_{\boldsymbol{\alpha}} \,f({\boldsymbol{\alpha}}), \end{align} \tag{ 10 }$$

where the differential operator ${L}^{+}_{\boldsymbol{\alpha}}$ is the adjoint of ${L}_{\boldsymbol{\alpha}}$ with respect to the L² inner product, i.e. $\langle{{L}^{+}_{\boldsymbol{\alpha}} g(\boldsymbol{\alpha}), h(\boldsymbol{\alpha})}\rangle_{\boldsymbol{\alpha}} = \langle{g(\boldsymbol{\alpha}), {L}_{\boldsymbol{\alpha}} h(\boldsymbol{\alpha})}\rangle_{\boldsymbol{\alpha}}$, which satisfies ${L}^{+}_{\boldsymbol{\alpha}} f({\boldsymbol{\alpha}}) = \langle \alpha | \mathcal{L}^{*} F | \alpha \rangle$.

The differential operator ${L}_{\boldsymbol{\alpha}}$ also has a set of eigensystem (an eigenvalue and right and left eigenfunctions) $\{\lambda_{k}, {u}_{k}({\boldsymbol{\alpha}}), {v}_{k}({\boldsymbol{\alpha}}) \}$ which satisfies:

$$\begin{align} &{L}_{\boldsymbol{\alpha}} {u}_{k} = \lambda_{k} {u}_{k}, \quad {L}^{+}_{\boldsymbol{\alpha}} {v}_{k} = \overline{\lambda_{k}} {v}_{k}, \quad \langle{{v}_{k}, {u}_{l}}\rangle_{\boldsymbol{\alpha}} = \delta_{kl}. \end{align} \tag{ 11 }$$

This eigensystem has a one-to-one correspondence with equation (5), where the eigenvalues $\{\lambda_{k} \}_{k \geqslant 0}$ are the same as those of ${\cal L}$; the eigenfunctions ${u}_{k}$ and ${v}_{k}$ of ${L}_{\boldsymbol{\alpha}}$ and ${L}^{+}_{\boldsymbol{\alpha}}$ are related to the eigenoperators U_k and V_k of $\mathcal{L}$ and $\mathcal{L}^{*}$ as:

$$\begin{align} U_{k} = \int {u}_{k}({\boldsymbol{\alpha}}) | \alpha \rangle \langle \alpha | d{\boldsymbol{\alpha}}, \quad {v}_{k}({\boldsymbol{\alpha}}) = \langle \alpha | V_{k} | \alpha \rangle, \end{align} \tag{ 12 }$$

which follow from $\mathcal{L} U_{k} = \int {u}_{k}({\boldsymbol{\alpha}}) \left\{\mathcal{L} | \alpha \rangle \langle \alpha | \right\} d{\boldsymbol{\alpha}} = \int \left\{{L}_{\boldsymbol{\alpha}} {u}_{k}({\boldsymbol{\alpha}}) \right\} | \alpha \rangle \langle \alpha | d{\boldsymbol{\alpha}} = \lambda_{k} U_{k}$ and ${L}^{+}_{\boldsymbol{\alpha}} {v}_{k} = {L}^{+}_{\boldsymbol{\alpha}} \langle \alpha | V_{k} | \alpha \rangle = \langle \alpha | {\mathcal L}^* V_{k} | \alpha \rangle = \overline{\lambda_{k}} \langle \alpha | V_{k} | \alpha \rangle = \overline{\lambda_{k}} {v}_{k}$.

In [46], we defined the quantum asymptotic phase function Φ₁ of the system state ρ as the argument (polar angle) of the expectation of the eigenoperator V₁ associated with the eigenvalue Λ₁ satisfying ${\cal L}^* {V}_1 = \overline{\Lambda_1} {V}_1$, namely,

$$\begin{align} \Phi_1(\rho) = \arg \langle V_1 \rangle = \arg \langle \rho, V_1 \rangle_{tr} = \arg \langle\, p({\boldsymbol{\alpha}}), v_1({\boldsymbol{\alpha}}) \rangle_{\boldsymbol{\alpha}}. \end{align} \tag{ 13 }$$

We showed that this quantum asymptotic phase yields appropriate phase values, namely, it always increases with a constant frequency Ω₁ as $\dot{\Phi}_1(\rho) = \Omega_1$ with the evolution of ρ even in the strong quantum regime and reproduces the conventional asymptotic phase in the semiclassical regime [46]. We note that $\Phi_1(\rho)$ is not defined when $\langle V_1 \rangle = 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 [48] from the Koopman operator viewpoint [49].

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 figure 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 Λ₁ on the rightmost branch but also the eigenvalues with the slowest decay rates on the other branches play important roles; for comparison, see figure 1(b) for a typical example of the eigenvalue Λ₁ in the semiclassical regime, where only a single dominant branch of eigenvalues exist. We denote these eigenvalues by $\Lambda_1, \Lambda_2, \Lambda_3, \ldots$ and their imaginary parts, i.e. fundamental frequencies, by $\Omega_j = \mbox{Im}\ \overline{\Lambda_j }$ ($j \geqslant 1$), and call $\{\Lambda_{j} \}_{j \geqslant 1}$ the principal eigenvalues. Here, the first principal eigenvalue Λ₁ and the fundamental frequency Ω₁ are those introduced in the previous subsection. These principal eigenvalues on the individual branches are shown by red dots in figure 1(a).

In a similar manner to the phase function Φ₁, we introduce the jth quantum asymptotic phase function $\Phi_j$ ($j = 2, 3, \ldots$) of ρ as the argument of the P-representation of the eigenoperator ${V}_j$ associated with the principal eigenvalue $\Lambda_j$ satisfying ${\cal L}^* {V}_j = \overline{\Lambda_j} {V}_j$, namely,

$$\begin{align} \Phi_j(\rho) = \arg \langle V_j \rangle = \arg \langle \rho, V_j \rangle_{tr} = \arg \langle\, p({\boldsymbol{\alpha}}), v_j({\boldsymbol{\alpha}}) \rangle_{\boldsymbol{\alpha}}. \end{align} \tag{ 14 }$$

Note that $\Phi_j(\rho)$ is not defined when $\langle V_j \rangle = 0$. Since

$$\begin{align} \frac{d}{dt} \langle V_j \rangle = \langle \dot\rho, V_j \rangle_{tr} = \langle {\cal L} \rho, V_j \rangle_{tr} = \langle \rho, {\cal L}^* V_j \rangle_{tr} = \overline{\Lambda_j} \langle \rho, V_j \rangle_{tr} = \overline{\Lambda_j} \langle V_j \rangle, \end{align} \tag{ 15 }$$

we obtain $\langle V_j \rangle(t) = \langle V_j \rangle(0) e^{\overline{\Lambda_j} t}$ and therefore:

$$\begin{align} \frac{d}{dt} \Phi_j(\rho) = \frac{d}{dt} \arg \langle V_j \rangle = \mbox{Im}\ \overline{\Lambda_j} = \Omega_ j. \end{align} \tag{ 16 }$$

Thus, $\Phi_j$ always increases with a constant frequency $\Omega_j$ with the evolution of ρ and plays the role of the asymptotic phase for any $j = 1, 2, 3, \ldots.$

It should be noted that $\Phi_j(\rho)$ cannot be defined for the stationary state ρ₀, because ρ₀ is the eigenfunction of the Liouville superoperator $\mathcal{L}$ with the eigenvalue $\lambda_0 = 0$ and hence $\langle V_j \rangle = \langle \rho_0, V_j \rangle_{tr} = 0$ from the biorthogonality in equation (5). As will be shown later, for the stationary state of a periodically driven system, $\langle V_j \rangle$ 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, $\Phi_{j}$ ($j \geqslant 2$) is not independent from Φ₁ and does not provide additional information, because $\Omega_j$ (j ≠ 1) is equal to Ω₁ according to the Koopman operator theory [56–59]. In the semiclassical regime, $\Omega_j$ (j ≠ 1) differs only slightly from Ω₁ as shown in figure 1(b) due to the effect of weak quantum noise and the corresponding decay rate for $j \geqslant 2$ is much larger, so the corresponding mode does not play an important role. However, as we see in the next section, $\Phi_{j}$ ($j \geqslant 2$) is distinctly different from Φ₁ and the corresponding decay rate characterized by $\mbox{Re}\ \overline{\Lambda_j}$ ($j \geqslant 2$) is comparable to $\mbox{Re}\ \overline{\Lambda_1}$ in the strong quantum regime, hence the phase function $\Phi_{j}$ ($j \geqslant 2$) yields independent information from Φ₁.

As an example, we consider a quantum van der Pol oscillator with Kerr effect. The master equation is given by [15, 27, 46]:

$$\begin{align} \dot{\rho} = \mathcal{L}_0 \rho, \quad \mathcal{L}_0 \rho = - i \left[ H,\rho\right] + \gamma_{1} \mathcal{D}[a^{\dagger}]\rho + \gamma_{2}\mathcal{D}[a^{2}]\rho, \end{align} \tag{ 17 }$$

where $H = \omega_0 a^{\dagger}a + K a^{\dagger 2} a^2$, ω₀ is the natural frequency of the oscillator, K is the Kerr parameter, and γ₁ and γ₂ 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 γ₂ 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 $\gamma_1 = 0.1$ and $(\omega_0, \gamma_{2}, K)/\gamma_{1} = (300, 4, 100)$, which are the same as in our previous study [46]. In the numerical calculation, we approximately truncated the density operator as a large-dimensional N × N matrix and mapped it into a N²-dimensional vector in the double-ket notation [60]. We can then represent the Liouville operator by a $N^2 \times N^2$ matrix and obtain the asymptotic phase in equation(13) from the eigensystem of this matrix [46].

Figure 1(a) shows the eigenvalues of ${\cal 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 $\Lambda_1, \Lambda_2, \Lambda_3, \ldots$ whose fundamental frequencies $\Omega_1 = \mbox{Im}\ \overline{\Lambda_1}, \Omega_2 = \mbox{Im}\ \overline{\Lambda_2}, \Omega_3 = \mbox{Im}\ \overline{\Lambda_3}, \ldots$ are different from each other. Moreover, their decay rates, which are characterized by $\mbox{Re}\ \overline{\Lambda_1}, \mbox{Re}\ \overline{\Lambda_2}, \mbox{Re}\ \overline{\Lambda_3}, \ldots$ are comparable to each other. This indicates that not only the quantum asymptotic phase Φ₁ associated with the principle eigenvalue Λ₁ of the slowest decaying mode [46] but also those associated with the other principal eigenvalues $\Lambda_2, \Lambda_3, \ldots$ can play important roles in the dynamics. We choose a negative value for each $\Omega_{j}$ so that the corresponding $\Phi_j$ increases from 0 to 2π in the counterclockwise direction.

In [15], it is shown that, in the strong quantum regime, $|{m+1}\rangle\langle{m}|$ is an approximate eigenoperator of ${\cal L}_0$ with the eigenvalue:

$$\begin{align} \widetilde{\lambda}_m=i[-\omega_0-2mK]- \frac{1}{2} \left\{\gamma_{1}(2 m+3)+2 \gamma_{2} m^{2} \right\}, \end{align} \tag{ 18 }$$

where $m = 0, 1, 2, \ldots$, namely, $\mathcal{L}^{-1}_0 |{m+1}\rangle\langle{m}| \approx \widetilde{\lambda}^{-1}_m |{m+1}\rangle\langle{m}|$. As shown in figure 1(a), these eigenvalues correspond to the principal eigenvalues of ${\cal L}_0$, i.e. $\Lambda_{j} \approx \widetilde{\lambda}_{j-1}$, and thus $V_{j} \approx |{j}\rangle\langle{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 differences in the transition frequencies arises due to the unequal spacing of the energy levels characterized by the Kerr parameter in equation (18).

As a comparison, we next consider the semiclassical regime where γ₂ and K are sufficiently small and the semiclassical approximation can be taken [26]. We set the parameters as $\gamma_1 = 1$ and $(\omega_0, \gamma_{2}, K)/\gamma_{1} =$ $(0.1, 0.05, 0.025)$, which are the same as those used in [46]. In this regime, we can approximate equation (9) by a quantum Fokker–Planck equation for $p({\boldsymbol{\alpha}})$ and the system can be approximately described as a Stuart-Landau oscillator (Hopf normal form) subjected to small quantum noise [15, 27, 46]. (Therefore, the conventional quantum vdP oscillator is also called a quantum Stuart-Landau oscillator recently [36, 38] and a more appropriate model of the quantum van der Pol oscillator has also been proposed [38, 39].)

Figure 1(b) shows the eigenvalues of $\mathcal{L}_0$ near the imaginary axis, where the principal eigenvalues are shown by red dots on the individual branches; $\Lambda_1 = \mu_1 + i \Omega_1$ ($\mu_1 \lt 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 $\hat{\lambda}_n = i \Omega_1 n + \mu_1 n^2~ (n = 0, \pm 1, \pm 2, \ldots)$ passing through $\overline{\Lambda_1}$, is isolated from other branches of eigenvalues with faster decay rates; the relative decay rate $\mbox{Re}~\overline{\Lambda_2}/\mbox{Re}~\overline{\Lambda_1} \approx -0.996/(-0.0736) \approx 13.5$ is more than three times larger than $\mbox{Re}~\overline{\Lambda_2}/\mbox{Re}~\overline{\Lambda_1} \approx -0.65/(-0.15) \approx 4.33$ in the strong quantum regime in figure 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 Ω₁; the small differences in the fundamental frequencies arise from small quantum noise. Thus, it is sufficient to consider only Λ₁ and introduce a single phase function Φ₁ 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({\boldsymbol{\alpha}})$, which represents the Stuart-Landau oscillator and possesses a stable limit-cycle solution [15, 27, 46]. Figure 1(c) shows a schematic diagram of the eigenvalues of the differential operator ${L}^{+}_{\boldsymbol{\alpha}}$, which are equivalent to those of the backward Liouville operator $\mathcal{L}^{*}_0$, in the classical limit. They are given in the form $\lambda_c = m \kappa_1 + i n \omega_1~ (m = 0,1,2,\ldots, n = 0, \pm 1, \pm2)$, where κ₁ is the real part of the largest negative eigenvalue and ω₁ is the pure-imaginary eigenvalue with the smallest absolute imaginary part. Thus, $\Phi_{j}$ ($j \geqslant 2$) is identical to Φ₁ and does not provide additional information, because $\Omega_j$ (j ≠ 2) is identical to Ω₁.

The above results suggest that the phase function $\Phi_{j}$ yields independent information from Φ₁ 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 with multiple fundamental frequencies rather than a limit cycle with a single fundamental frequency, and that we need to consider the phase functions $\Phi_2, \Phi_3, \ldots$ associated with $\Lambda_2, \Lambda_3, \ldots$ in addition to Φ₁ and Λ₁. Figure 2 shows a schematic picture of the torus behavior of the system with two fundamental frequencies.

Zoom In Zoom Out Reset image size

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 $\Phi_j({\boldsymbol{\alpha}})$ for $j = 1, 2, 3,$ and 4 of the pure coherent state α on the complex plane $(x = \mbox{Re}\ \alpha, p = \mbox{Im}\ \alpha)$. The parameters are the same as in figure 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 $\Phi_j$ for $j \geqslant 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 [56–59], but the differences in their eigenfrequencies brought by the strong quantum effect take an important role in analyzing strong quantum signatures in synchronization.

Zoom In Zoom Out Reset image size

To demonstrate that these quantum asymptotic phase functions yield appropriate phase values, we consider free oscillatory relaxation of ρ from a pure coherent initial state $\rho = |{\alpha_0}\rangle\langle{\alpha_0}|$ with $\alpha_0 = 1$ at t = 0 and measured the evolution of the expectation values of V_j and their arguments, i.e. their asymptotic phases $\Phi_j(\rho) = \arg \langle V_j \rangle$ $(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 $\langle V_j \rangle$ exhibits exponentially damped harmonic oscillations as shown in figures 4(a), (c), (e) and (g), and correspondingly, each $\Phi_j(\rho)$ gives constantly varying phase values with frequency $\Omega_j$ as shown in figures 4(b), (d), (f) and (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, $\langle a \rangle$ shown in figure 4(i) exhibits more complex oscillatory dynamics and the simple argument $\arg \langle a \rangle$ does not vary constantly with time as shown in figure 4(j). Thus, $\arg \langle a \rangle$ 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 $\Omega_j$ and simple argument $\arg \langle a \rangle$ can be observed more clearly than in the semiclassical regime or in the classical limit [46].

Zoom In Zoom Out Reset image size

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 is:

$$\begin{align} \dot{\rho} = ( \mathcal{L}_0 + \mathcal{L}_1 ) \rho, \end{align} \tag{ 19 }$$

where ${\cal L}_0$ is given by equation (17) with $H = -\Delta a^{\dagger}a + K a^{\dagger 2} a^2,$ and the harmonic drive is represented by $\mathcal{L}_1 \rho = -i \left[ i E (a - a^{\dagger}),\rho\right],$ where $\Delta = \omega_{d} - \omega_{0}$ is the frequency detuning of the harmonic drive from the oscillator and E is the strength of the harmonic drive [15, 27, 46]. We use the same parameters as in figure 1(a) and vary the detuning parameter Δ by controlling ω_d while keeping the natural frequency ω₀ fixed.

In [15], Lörch et al showed that this system under strong Kerr effect exhibits multiple phase locking to the harmonic drive at several detuning frequencies $\Delta = 2mK$ ($m = 0, 1, 2, \ldots$), i.e. at $\omega_d = \omega_0 + 2m K$, 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 $\omega_d = (\,{p}/{q}) \omega_0$ with $p, q = 1, 2, 3, \ldots$ ($p/q \neq 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 [15]).

In [15], the following order parameter S_a and power spectrum P_a defined using the annihilation operator a are used to analyze the system:

$$\begin{align} S_a = {|S_a|} e^{i \theta_a} = \frac{\langle{a}\rangle }{\sqrt{\langle{a^{\dagger} a}\rangle}}, \end{align} \tag{ 20 }$$

$$\begin{align} P_{a}(\omega) = \int_{-\infty}^{\infty} d\tau e^{i\omega \tau} \left( \langle{a^{\dagger}( \tau ) a(0)}\rangle - \langle{a^{\dagger}( \tau )}\rangle \langle{a(0)}\rangle\right), \end{align} \tag{ 21 }$$

where the expectation is taken with respect to the steady-state density operator obtained from equation (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 ${\cal L}_0$ as:

$$\begin{align} S_j = {|S_j|} e^{i \theta_j} = \frac{\langle{V_j}\rangle }{\sqrt{\langle{V_j^{\dagger} V_j}\rangle}}, \end{align} \tag{ 22 }$$

$$\begin{align} P_j(\omega) = \int_{-\infty}^{\infty} d\tau e^{i\omega \tau}\left(\langle{V_j^{\dagger}( \tau ) V_j(0)}\rangle - \langle{V_j^{\dagger}( \tau )}\rangle \langle{V_j(0)}\rangle\right), \end{align} \tag{ 23 }$$

where $V_j^{\dagger}(\tau) = e^{{\cal L}_{0}^{*} \tau} V_j^{\dagger}(0) = e^{\Lambda_j \tau} V_j^{\dagger}(0)$ and $V_j^{\dagger}(0) = V_j^{\dagger}$ ($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 \omega_d a^{\dagger} a t} V_{j} e^{- i \omega_d a^{\dagger} a t} = {V}_{j} e^{i \Omega t}$ with a constant Ω (see appendix). The resulting power spectra in the rotating frame are simply shifted from the power spectra in the original frame by $-\Omega$. Using these two quantities in equations (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 figure 5(a) showing ${|S_a|}$, several Arnold tongues representing phase locking of the oscillator at different frequencies are observed [15]. Figures 5(b)–(e) show ${|S_j|}$ for $j = 1, 2, 3$ and 4, respectively. It is remarkable that the Arnold tongues in figure 5(a) are clearly decomposed into individual Arnold tongues around $\Delta = 2(j-1)K$ in figures 5(b)–(e). This indicates that the order parameters in terms of the asymptotic phases can capture the phase locking dynamics at $\omega_d = \omega_0 + 2(j-1)K$ for different j individually.

Zoom In Zoom Out Reset image size

Similarly, figure 6 shows the power spectra P_a and P_j for $\Delta = 0$. Multiple peaks of P_a in figure 6(a), which indicate multiple phase locking of the oscillator to the harmonic drive, are clearly decomposed into individual peaks around $\omega = \Delta - 2 (j-1) K$ in figures 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 $\mbox{Re}\ \Lambda_{j}$ is smaller. Though not shown, we can also detect even smaller tongues and peaks with $j \geqslant 5$.

Zoom In Zoom Out Reset image size

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 [61], 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 Λ₁ 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 $\Lambda_{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.