A definition of the asymptotic phase for quantum nonlinear oscillators from the Koopman operator viewpoint

Consider a deterministic dynamical system

where $\bm{X}(t)\in\mathbb{R}^{N}$ is the system state at time $t$, ${\bm{A}}({\bm{X}})\in{\mathbb{R}}^{N}$ is a vector field representing the system dynamics, and $(\dot{})$ represents time derivative. We assume that this system has an exponentially stable limit-cycle solution ${\bm{X}}_{0}(t)$ with a natural period $T$ and frequency $\Omega_{c}=2\pi/T$, satisfying ${\bm{X}}_{0}(t+T)={\bm{X}}_{0}(t)$, and denote its basin of attraction as $B\subseteq{\mathbb{R}}^{N}$. The asymptotic phase $\Phi_{c}({\bm{X}}):{B}\to[0,2\pi)$ is defined such that ${\bm{A}}({\bm{X}})\cdot\nabla\Phi_{c}({\bm{X}})=\Omega_{c}$ is satisfied for $\forall{\bm{X}}\in B$, where $\nabla=\partial/\partial{\bm{X}}$ represents the gradient with respect to ${\bm{X}}$ winfree2001geometry ; kuramoto1984chemical ; pikovsky2001synchronization ; nakao2016phase ; ermentrout2010mathematical . The asymptotic phase $\phi(t)=\Phi_{c}({\bm{X}}(t))$ of the system state ${\bm{X}}(t)$ then obeys

i.e., $\phi$ always increases with a constant frequency $\Omega_{c}$ as ${\bm{X}}$ evolves in $B$. Thus, the asymptotic phase $\Phi_{c}$ gives a nonlinear transformation of the system state ${\bm{X}}$ to a phase value $\phi$ such that the dynamics of $\phi$ takes a simple linear form, $\phi(t)=\Omega_{c}t+const.$ The simplicity of the phase equation (2) has facilitated detailed studies of synchronization and collective dynamics in coupled-oscillator systems winfree2001geometry ; kuramoto1984chemical ; pikovsky2001synchronization ; nakao2016phase ; ermentrout2010mathematical . The level sets of $\Phi_{c}({\bm{X}})$ are called isochrons.

The linear operator ${A}={\bm{A}}({\bm{X}})\cdot\nabla$ in the definition of the asymptotic phase $\Phi_{c}({\bm{X}})$ is the infinitesimal generator of the Koopman operator describing the evolution of observables for the system described by Eq. (1) (see Refs. mauroy2013isostables ; mauroy2020koopman ; shirasaka2017phase ; kuramoto2019concept for details). The complex exponential $\Psi_{c}({\bm{X}})=e^{i\Phi_{c}({\bm{X}})}$ of $\Phi_{c}({\bm{X}})$ is an eigenfunction of ${A}$ with the eigenvalue $i\Omega_{c}$, namely, it satisfies the eigenvalue equation ${A}\Psi_{c}({\bm{X}})=i\Omega_{c}\Psi_{c}({\bm{X}})$. Therefore, the asymptotic phase $\Phi_{c}({\bm{X}})$ has a natural operator-theoretic interpretation as the argument of the Koopman eigenfunction $\Psi_{c}({\bm{X}})$ associated with the eigenvalue $i\Omega_{c}$, characterized by the fundamental frequency $\Omega_{c}$ of the oscillator mauroy2013isostables ; mauroy2020koopman ; shirasaka2017phase ; kuramoto2019concept . We can thus define the asymptotic phase by using the Koopman eigenfunction $\Psi_{c}({\bm{X}})$ of $A$ as

For stochastic oscillatory systems, we cannot use the deterministic limit-cycle solution in defining the asymptotic phase unless the noise is sufficiently weak. Thomas and Lindner thomas2014asymptotic defined the asymptotic phase for classical stochastic oscillators without relying on the deterministic limit cycle by using the eigenfunction with the slowest decay rate of the backward Fokker-Planck operator.

Consider a stochastic oscillator described by an Ito stochastic differential equation (SDE)

where $\bm{X}(t)\in\mathbb{R}^{N}$ is the system state at time $t$, ${\bm{A}}({\bm{X}})\in{\mathbb{R}}^{N}$ is a drift vector representing the deterministic dynamics, ${\bm{B}}({\bm{X}})\in{\mathbb{R}}^{N\times N}$ is a matrix characterizing the effect of the noise, and ${\bm{W}}(t)\in{\mathbb{R}}^{N}$ is a $N$-dimensional Wiener process. This system is assumed to be oscillatory in the sense explained below. The transition probability density $p(\bm{X},t|\bm{Y},s)$ ($t\geq s$) of Eq. (4) obeys the forward and backward Fokker-Planck equations gardiner2009stochastic ,

and

respectively, where the (forward) Fokker-Planck operator is given by

and the backward Fokker-Planck operator is given by (in terms of ${\bm{X}}$)

Here, ${\bm{D}}({\bm{X}})={\bm{B}}({\bm{X}}){\bm{B}}({\bm{X}})^{\sf T}\in{\mathbb{R}}^{N\times N}$ is a matrix of diffusion coefficients (${\sf T}$ indicates matrix transposition). The forward and backward operators $L_{\bm{X}}$ and $L_{\bm{X}}^{*}$ are mutually adjoint, i.e., $\langle{{L}_{\bm{X}}G(\bm{X}),H(\bm{X})}\rangle_{\bm{X}}=\langle{G(\bm{X}),{L}^{*}_{\bm{X}}H(\bm{X})}\rangle_{\bm{X}}$, where the inner product is defined as $\langle{G(\bm{X}),H(\bm{X})}\rangle_{\bm{X}}=\int\overline{G(\bm{X})}H(\bm{X})d\bm{X}$ for two functions $G(\bm{X}),H(\bm{X}):{\mathbb{R}}^{N}\to{\mathbb{C}}$ (the overline indicates complex conjugate and the integration is taken over the whole range of ${\bm{X}}$).

As explained in Appendix A, the backward Fokker-Planck operator $L_{\bm{X}}^{*}$ in Eq. (8) is the infinitesimal generator of the Koopman operator for Eq. (4). We denote the probability density function of ${\bm{X}}$ at time $t$ as $p_{t}({\bm{X}})\in{\mathbb{R}}$, which also obeys the Fokker-Planck equation $\partial p_{t}({\bm{X}})/\partial t=L_{\bm{X}}p_{t}({\bm{X}})$, and an observable of the system at time $t$ as $g_{t}:{\mathbb{R}}^{N}\to{\mathbb{C}}$, which maps the system state ${\bm{X}}$ to a complex value. The evolution of the expectation $\langle g\rangle_{t}=\int p_{t}({\bm{X}})g_{t}({\bm{X}})d{\bm{X}}=\langle p_{t}({\bm{X}}),g_{t}({\bm{X}})\rangle_{\bm{X}}$of the observable $g$ at $t=t_{0}$ can be expressed as

where $g_{t}({\bm{X}})$ remains constant and $p_{t}({\bm{X}})$ evolves in the second expression, while $p_{t}({\bm{X}})$ remains constant and $g_{t}({\bm{X}})$ evolves in the last expression.

The linear differential operators ${L}_{\bm{X}}$ and ${L}_{\bm{X}}^{*}$ have a biorthogonal eigensystem $\{\lambda_{k},P_{k},Q_{k}\}_{k=0,1,2,...}$ of the eigenvalue $\lambda_{k}$ and eigenfunctions $P_{k}({\bm{X}})$ and $Q_{k}({\bm{X}})$ satisfying

where $k,l=0,1,2,\ldots$ and $\delta_{kl}$ represents the Kronecker delta gardiner2009stochastic . We assume that, among the eigenvalues, one eigenvalue $\lambda_{0}$ is zero, which is associated with the stationary state $p^{S}({\bm{X}})$ of the system satisfying $L_{\bm{X}}p^{S}({\bm{X}})=0$, and all other eigenvalues have negative real parts. It is assumed that the eigenvalues with the largest non-negative real part (or the smallest absolute real part, i.e., the slowest decay rate) are given by a complex-conjugate pair. We denote these eigenvalues as

where ${|\mu_{s}|}~{}(\mu_{s}<0)$ is the decay rate and $\Omega_{s}=\mbox{Im}\ \overline{\lambda_{1}}$ represents the fundamental oscillation frequency of the system. The oscillatory property of the system is embodied in this assumption. We also assume that this pair of principal eigenvalues are well separated from other branches of eigenvalues and this oscillatory mode is dominant in the system.

Thomas and Lindner thomas2014asymptotic proposed a definition of the stochastic asymptotic phase of the system described by Eq. (4) by using the argument of the complex conjugate of the eigenfunction ${Q_{1}({\bm{X}})}$ of $L^{*}_{\bm{X}}$ associated with $\overline{\lambda_{1}}$ as (in the notation used here)

This definition is natural from the Koopman operator viewpoint kato2021asymptotic , as $L_{\bm{X}}^{*}$ is the infinitesimal generator of the Koopman operator for Eq. (4) that formally goes back to the linear operator $A={\bm{A}}({\bm{X}})\cdot\nabla$, i.e., to the infinitesimal generator of the Koopman operator for the deterministic system Eq. (1) in the noiseless limit ${\bm{D}}({\bm{X}})\to 0$. The exponential average of $\Phi_{s}$ satisfies kato2021asymptotic

where $\mathbb{E}^{{\bm{X}}_{0}}$ represents the average over the stochastic trajectories of Eq. (4) starting from an initial point ${\bm{X}}_{0}\in{\mathbb{R}}^{N}$. Thus, the asymptotic phase defined by Eq. (13) increases with a constant frequency $\Omega_{s}$ on average with the stochastic evolution of the system and can be considered a natural generalization of the deterministic asymptotic phase in Eq. (3). It is noted that we may also choose the eigenvalue $\lambda_{1}$ and eigenfunction $\overline{Q_{1}(\bm{X})}$ to define the stochastic asymptotic phase, which reverses its direction of increase.

We consider quantum oscillatory systems with a single degree of freedom coupled to reservoirs. The system’s quantum state is represented by a Hermitian density operator $\rho$ ($=\rho^{\dagger}$) and the observable is described by an operator $F$, where ${\dagger}$ represents Hermitian conjugate. Introducing an inner product $\langle{X,Y}\rangle_{tr}={\rm Tr}\hskip 1.9919pt{(X^{{\dagger}}Y)}$ of two operators $X$ and $Y$, the expectation of the observable $F$ is expressed as

In the Schrödinger picture, the quantum state $\rho$ evolves with time while the observable $F$ remains constant. Assuming that the interactions of the system with the reservoirs are instantaneous and Markovian approximation can be employed, the evolution of the density operator $\rho_{t}$ at time $t$ obeys the quantum master equation carmichael2007statistical ; gardiner1991quantum ; breuer2002theory

where the Liouville operator ${\mathcal{L}}$ (sometimes called a superoperator because it acts on an operator) is given by

Here, $H(=H^{{\dagger}})$ is a system Hamiltonian, $C_{j}$ is a coupling operator between the system and $j$th reservoir $(j=1,\ldots,n)$, $[A,B]=AB-BA$ is the commutator, $\mathcal{D}[C]X=CXC^{{\dagger}}-(XC^{{\dagger}}C+C^{{\dagger}}CX)/2$ is the Lindblad form, and the reduced Planck’s constant is set as $\hbar=1$.

In the Heisenberg picture, the quantum state $\rho$ remains constant while the observable $F$ evolves with time as

where the time dependence of $F$ is explicitly shown. Here, ${\mathcal{L}}^{*}$ is the adjoint operator of ${\mathcal{L}}$ satisfying

and can be explicitly calculated as

where $\mathcal{D}^{+}[C]X=C^{{\dagger}}XC-(XC^{{\dagger}}C+C^{{\dagger}}CX)/2$ is the adjoint Lindblad form breuer2002theory . The evolution of the expectation $\langle F\rangle_{t}=\langle\rho_{t},F_{t}\rangle_{tr}$ of $F$ with respect to $\rho$ at $t=t_{0}$ can be expressed as

where $F$ remains constant and $\rho$ evolves in the second expression (Schrödinger picture), while $\rho$ remains constant and $F$ evolves in the last expression (Heisenberg picture). Equation (22) corresponds to Eq. (10) for the expectation of the observable for classical stochastic systems. Thus, the adjoint operator ${\mathcal{L}}^{*}$ is a counterpart of the backward Fokker-Planck operator $L_{\bm{X}}^{*}$ in Sec. II, namely, ${\mathcal{L}}^{*}$ corresponds to the infinitesimal generator of the Koopman operator.

We assume that the operators ${\mathcal{L}}$ and ${\mathcal{L}}^{*}$ have a biorthogonal eigensystem $\{\Lambda_{k},U_{k},V_{k}\}_{k=0,1,2,...}$ consisting of the eigenvalue $\Lambda_{k}$ and right and left eigenoperators $U_{k}$ and $V_{k}$, satisfying

for $k,l=0,1,2,\ldots$ li2014perturbative . Among the eigenvalues, one eigenvalue $\Lambda_{0}$ is always $0$, which is associated with the stationary state $\rho^{S}$ of the system satisfying ${\mathcal{L}}\rho^{S}=0$, and all other eigenvalues have negative real parts; This also indicates that the system has no decoherece free subspace lidar1998decoherence . We assume that, reflecting the system’s oscillatory dynamics, the eigenvalues with the largest non-vanishing real part (i.e., with the slowest decay rate) are given by a complex-conjugate pair. We denote these eigenvalues as

where $|\mu_{q}|$ $(\mu_{q}<0)$ is the decay rate and $\Omega_{q}=\mbox{Im}\ \overline{\Lambda_{1}}$ gives the fundamental frequency of the oscillation. As in the case of the stochastic oscillatory systems in Sec. II, the oscillatory property of the system is embodied in this assumption.

The density operator $\rho$ can also be represented by using quasiprobability distributions in the phase space such as the $P$, $Q$, and Wigner distributions  carmichael2007statistical ; gardiner1991quantum ; cahill1969density . We use the $P$ representation and express $\rho$ as

where $|\alpha\rangle$ is a coherent state specified by a complex value $\alpha\in\mathbb{C}$, or equivalently by a complex vector $\bm{\alpha}=(\alpha,\overline{\alpha})^{T}\in\mathbb{C}^{2}$, $p({\bm{\alpha}}):{\mathbb{C}}^{2}\to{\mathbb{R}}$ is a quasiprobability distribution of ${\bm{\alpha}}$, $d{\bm{\alpha}}=d\alpha d\overline{\alpha}$, and the integral is taken over $\mathbb{C}$. Defining the $P$ representation of an observable $F$ as

where $F$ is expressed in the normal order carmichael2007statistical ; gardiner1991quantum ; cahill1969density , the expectation of $F$ is expressed as

Here, we defined the $L^{2}$ inner product $\langle{g(\bm{\alpha}),h(\bm{\alpha})}\rangle_{\bm{\alpha}}=\int\overline{g(\bm{\alpha})}h(\bm{\alpha})d\bm{\alpha}$ of two functions $g(\bm{\alpha})$, $h(\bm{\alpha}):{\mathbb{C}}^{2}\to{\mathbb{C}}$, where the integral is taken over the complex plane.

In the Schrödinger picture, the time evolution of $p_{t}({\bm{\alpha}})$ (dependence on $t$ is explicitly denoted) corresponding to the master equation (16) is describied by a partial differential equation

where the differential operator ${L}_{\bm{\alpha}}$ is related to the Liouville operator ${\mathcal{L}}$ in Eq. (16) via

and can be explicitly calculated from Eq. (16) by using the standard calculus for the phase-space representation  carmichael2007statistical ; gardiner1991quantum ; cahill1969density .

The corresponding evolution of the $P$ representation $f_{t}({\bm{\alpha}})$ of the observable $F_{t}$ in the Heisenberg picture is given by

where the differential operator $L_{\bm{\alpha}}^{*}$ is the adjoint of $L_{\bm{\alpha}}$, i.e.,

and satisfies

Thus, $L_{\bm{\alpha}}^{*}$ is the generator of the Koopman operator in the $P$ representation describing the evolution of $f(\bm{\alpha})$, which corresponds to the adjoint Liouville operator ${\mathcal{L}}^{*}$ in Eq. (20).

Corresponding to the Liouville operators ${\mathcal{L}}$ and ${\mathcal{L}}^{*}$, the differential operators $L_{\bm{\alpha}}$ and $L^{*}_{\bm{\alpha}}$ also possess a biorthogonal eigensystem $\{\Lambda_{k},{u}_{k}({\bm{\alpha}}),{v}_{k}({\bm{\alpha}})\}_{k=0,1,2,...}$ of eigenvalue $\Lambda_{k}$ and eigenfunctions $u_{k}$ and $v_{k}$, satisfying

which has one-to-one correspondence with Eq. (23). The eigenvalues $\{\Lambda_{k}\}_{k=0,1,2,...}$ are the same as those of ${\mathcal{L}}$, and the eigenfunctions ${u}_{k}$ and ${v}_{k}$ of $L_{\bm{\alpha}}$ are related to the eigenoperators $U_{k}$ and $V_{k}$ of ${\mathcal{L}}$ via

which follow from

and

Generalizing the definition for classical stochastic oscillatory systems in Sec. II, we here propose a definition of the quantum asymptotic phase. We note that, in quantum systems, the system state is given by the density operator $\rho$ and individual trajectories as in the classical stochastic systems cannot be considered.

First, we define the quantum asymptotic phase $\Phi_{q}({\bm{\alpha}})$ of the coherent state $\bm{\bm{\alpha}}$ in the $P$ representation. Considering the definition Eq. (13) of the asymptotic phase in terms of the eigenfunction $Q_{1}({\bm{X}})$ of the backward Fokker-Planck operator $L^{*}_{\bm{X}}$ in the classical stochastic case, we define the quantum asymptotic phase $\Phi_{q}(\bm{\alpha})$ as the argument of the complex conjugate of the eigenfunction $v_{1}({\bm{\alpha}})$ in the $P$ representation associated with the principal eigenvalue $\overline{\Lambda_{1}}$ as

Next, considering that the general quantum state $\rho$ is represented as a superposition of coherent states with the weight $p({\bm{\alpha}})$, Eq. (25), we define the asymptotic phase of $\rho$ as

It can be shown that the asymptotic phase $\Phi_{q}(\rho)$ evolves with a constant frequency $\Omega_{q}$ as the quantum state $\rho$ evolves according to the master equation (16). Defining

where the dependence on time $t$ is explicitly denoted, we have

or, equivalently,

Integrating by time, we obtain

where $\rho_{0}$ is the initial state at $t=0$, and hence the asymptotic phase is given by

Differentiating by $t$, we obtain

namely, the asymptotic phase $\Phi_{q}(\rho_{t})$ increases with a constant frequency $\Omega_{q}$ with the evolution of the quantum state $\rho_{t}$.

Thus, by using the eigenfunction $v_{1}({\bm{\alpha}})$ of the adjoint linear operator $L^{*}_{\bm{\alpha}}$ or equivalently the eigenoperator $V_{1}$ of the adjoint operator ${\mathcal{L}}^{*}$ associated with the eigenvalue $\overline{\Lambda_{1}}$, we can define the asymptotic phase $\Phi_{q}(\rho)$ of the quantum state $\rho$. The quantum master equation (16), adjoint Liouville operator ${\mathcal{L}}^{*}$ (or adjoint differential operator $L^{*}_{\bm{\alpha}}$ in the $P$ representation), and eigenoperator $V_{1}$ (or the eigenfunction $v_{1}({\bm{\alpha}})$ in the $P$ representation) with the eigenvalue $\overline{\Lambda_{1}}$ correspond to the forward Fokker-Planck equation (5), backward Fokker-Planck operator $L_{\bm{X}}^{*}$, and eigenfunction ${Q_{1}(\bm{X})}$ with the eigenvalue $\overline{\lambda_{1}}$ in the classical stochastic system discussed in Sec. II, respectively.

We stress, however, that the system state is generally represented by the density operator $\rho$ or quasiprobability distribution $p({\bm{\alpha}})$ and individual trajectories cannot be considered in the quantum case. As in the classical stochastic case, we may also choose the eigenvalue $\Lambda_{1}$, eigenfunction $\overline{v_{1}(\bm{\alpha})}$, and eigenoperator $V_{1}^{{\dagger}}$ to define the quantum asymptotic phase, which reverses its direction of increase.

As an example of quantum oscillatory systems, we consider the quantum van der Pol model with quantum Kerr effect. The system’s density operator $\rho$ obeys the master equation

where $a$ and $a^{\dagger}$ are annihilation and creation operators, $H=\omega_{0}a^{{\dagger}}a+Ka^{{\dagger}2}a^{2}$ is the Hamiltonian, $\omega_{0}$ is the frequency parameter of the oscillator, $K$ is the Kerr parameter, and $\gamma_{1}$ and $\gamma_{2}$ are the decay rates for negative damping and nonlinear damping, respectively kato2020semiclassical ; lorch2016genuine . By using the standard rule of calculus carmichael2007statistical ; gardiner1991quantum ; cahill1969density , we can derive the differential operator $L_{\bm{\alpha}}$ in Eq. (28) describing the evolution of the $P$ representation $p({\bm{\alpha}})$ of $\rho$ from the Liouville operator ${\mathcal{L}}$, which consists of third- and higher-order derivative with respect to $\alpha$ lorch2016genuine .

To evaluate the fundamental frequency $\Omega_{q}$ and the asymptotic phase $\Phi_{q}(\rho)$ of the system, we numerically calculate the eigenvalues and eigenfunctions of the Liouville and adjoint Liouville operators ${\mathcal{L}}$ and ${\mathcal{L}}^{*}$. To this end, we approximately truncate the number representation of the density operator as a large $N\times N$ matrix and map it to a $N^{2}$-dimensional vector of the double-ket notation albert2018lindbladians . We can then approximately represent the Liouville operators ${\mathcal{L}}$ and ${\mathcal{L}}^{*}$ by $N^{2}\times N^{2}$ matrices and calculate their eigensystem to obtain the asymptotic phase in Eq. (36).

We first consider the semiclassical regime where $\gamma_{2}$ and $K$ are sufficiently small. In this case, as explained in Appendix B, we can approximate Eq. (28) by a Fokker-Planck equation for $p({\bm{\alpha}})$, namely, the system state is approximately equivalent to a classical stochastic system. Furthermore, in the classical limit where the quantum noise vanishes, the system is described by a single complex variable $\alpha\in{\mathbb{C}}$ obeying a deterministic ordinary differential equation

This equation represents the Stuart-Landau oscillator (normal form of the supercritical Hopf bifurcation) kuramoto1984chemical and possesses a stable limit-cycle solution ${\alpha}_{0}(\phi)=Re^{i\phi}$, which is represented as a function of the phase $\phi=\Omega_{c}t+const.$ with a natural frequency $\Omega_{c}=-\omega_{0}-K\gamma_{1}/\gamma_{2}$ and radius $R=\sqrt{{\gamma_{1}}/{2\gamma_{2}}}$. The basin $B$ of this limit cycle is the whole complex plane except the origin. The classical asymptotic phase $\Phi_{c}$ of this system is given by kato2019semiclassical

and satisfies $\dot{\Phi}_{c}({\bm{\alpha}})=\Omega_{c}$ as ${\alpha}$ evolves in $B$ under Eq. (47) nakao2016phase . As stated in Sec. II, this $\Phi_{c}({\bm{\alpha}})$ is the argument of the Koopman eigenfunction $\Psi_{c}({\bm{\alpha}})$ of the system associated with the eigenvalue $i\Omega_{c}$. In Ref. kato2019semiclassical , we used this $\Phi_{c}$ for the phase-reduction analysis of quantum synchronization in the semiclassical regime with weak quantum noise. It is expected that the quantum asymptotic phase $\Phi_{q}(\bm{\alpha})$ of the coherent state $\bm{\alpha}$ is close to the classical asymptotic phase $\Phi_{c}(\bm{\alpha})$ when the quantum noise is sufficiently small.

Figure 1(a) shows the eigenvalues of ${\mathcal{L}}$ near the imaginary axis obtained numerically, where the principal eigenvalue $\overline{\Lambda_{1}}=\mu_{q}+i\Omega_{q}$ is shown by a red dot ($\mu_{q}<0$), and Figs. 1(b) and 1(c) compare the quantum-mechanical phase $\Phi_{q}(\bm{\alpha})$ with the corresponding classical phase $\Phi_{c}({\bm{\alpha}})$. Here, we adopt a negative value for $\Omega_{q}$ so that the resulting phase $\Phi_{q}$ increases in the counterclockwise direction from $0$ to $2\pi$ on the complex plane, i.e., $\Phi_{q}$ satisfies $\oint_{c}\nabla\Phi_{q}({\bm{x}})\cdot d{\bm{x}}=2\pi$ where ${\bm{x}}=(x,p)=(\mbox{Re}~{}\alpha,\mbox{Im}~{}\alpha)$ and $C$ is a circle around $0$.

As the quantum noise is small, the quantum frequency $\Omega_{q}$ and the asymptotic phase $\Phi_{q}(\bm{\alpha})$ obtained numerically from ${\mathcal{L}}$ are close to the classical frequency $\Omega_{c}$ and asymptotic phase $\Phi_{c}({\bm{\alpha}})$, respectively. The difference between $\Omega_{q}$ and $\Omega_{c}$ arises from the small quantum noise; in the limit of vanishing quantum noise, the eigenfunction ${v}_{1}({\bm{\alpha}})$ of $L_{\bm{\alpha}}$ in the $P$ representation coincides with the Koopman eigenfunction of the deterministic system Eq. (47) with the eigenvalue $i\Omega_{c}$ and therefore $\Phi_{q}$ reproduces the classical phase $\Phi_{c}$ (see Appendix C). Here, we note that the principal eigenvalues $i\Omega_{q}$ and $i\Omega_{c}$ are well separated from other branches of eigenvalues with faster decay rates and the corresponding oscillatory modes become quickly dominant.

To demonstrate that the present definition of the quantum asymptotic phase yields appropriate values, we consider free oscillatory relaxation of $\rho_{t}$ from a coherent initial state $\rho^{\alpha_{0}}=|{\alpha_{0}}\rangle\langle{\alpha_{0}}|$ with $\alpha_{0}=1$ at $t=0$ and measure $\Phi_{q}(\rho_{t})$. For comparison, we also measure the argument $\mbox{arg}\ \langle a\rangle_{t}$ of the expectation $\langle a\rangle_{t}=\langle\rho_{t},a\rangle_{tr}$ of the annihilation operator $a$, which simply gives the polar angle of $\langle a\rangle_{t}$ on the complex plane. We note that, although the system state $\rho$ starts from a pure coherent state, it soon becomes a mixed state due to the coupling with the reservoirs and eventually relaxes to the stationary state $\rho^{S}$.

Figures 1(d) and (e) plot the evolution of the expectation $\Psi_{q}(\rho_{t})=\langle\rho_{t},V_{1}\rangle_{tr}$ of the eigenoperator $V_{1}$ and the quantum asymptotic phase $\Phi_{q}(\rho_{t})=\arg\Psi_{q}(\rho_{t})$, respectively, and Figs. 1(f) and (g) show the evolution of the expectation $\langle a\rangle$ and its polar angle $\arg\langle a\rangle$, respectively. The asymptotic phase $\Phi_{q}(\rho_{t})$ increases with a constant frequency $\Omega_{q}$ and appropriately yields isochronous phase values. In contrast, the polar angle $\arg\langle a\rangle$ does not increase constantly with time, in particular in the transient process before $t=10$, as shown in Fig. 1(g); as the limit cycle in the classical limit is rotationally symmetric in this model, the polar angle also yields almost constantly increasing phase values after relaxation.

Thus, the quantum asymptotic phase $\Phi_{q}$ increases with a constant frequency $\Omega_{q}$ in the semiclassical regime of a quantum van der Pol model with the quantum Kerr effect.

Next, we consider a strong quantum regime with relatively large $\gamma_{2}$ and $K$, where only a small number of energy states participates in the system dynamics and the semiclassical description is not valid. The eigenvalues of ${\mathcal{L}}$ obtained numerically are shown in Fig. 2(a). Figures 2(b) and 2(c) show the quantum-mechanical phase $\Phi_{q}$ and the corresponding classical phase $\Phi_{c}$. Because the system is in the strong quantum regime, $\Phi_{c}$ is distinctly different from $\Phi_{q}$ and the classical frequency $\Omega_{c}$ also differs largely from the true quantum frequency $\Omega_{q}$. Here, we again note that the principal eigenvalues have much smaller (less than half of the second largest) decay rate than the eigenvalues in the other branches.

We consider free oscillatory relaxation of $\rho$ from a coherent initial state $\rho=|{\alpha_{0}}\rangle\langle{\alpha_{0}}|$ with $\alpha_{0}=1$ at $t=0$ and measure the evolution of the asymptotic phase $\Phi_{q}(\rho_{t})=\langle V_{1}\rangle_{t}$ of the system state $\rho_{t}$. For comparison, we also measure the polar angle $\mbox{arg}\langle a\rangle_{t}$ of $\langle a\rangle_{t}$.

Figures 2(d) and (e) plot the evolution of the expectation $\Psi_{q}(\rho_{t})=\langle\rho_{t},V_{1}\rangle_{tr}$ of the eigenoperator $V_{1}$ and the quantum asymptotic phase $\Phi_{q}(\rho_{t})=\arg\Psi_{q}(\rho_{t})$, respectively, and Figs. 2(f) and (g) show the evolution of the expectation $\langle a\rangle$ and its polar angle $\arg\langle a\rangle$, respectively. As expected, the asymptotic phase $\Phi_{q}(\rho_{t})$ appropriately gives constantly varying phase values with the frequency $\Omega_{q}$. In contrast, the polar angle $\arg\langle a\rangle_{t}$ does not vary constantly with time and is not isochronous. This is because the transition between relatively a small number of energy levels takes part in the system dynamics and the discreteness of the energy spectra can play important roles in this strong quantum regime.

Thus, the quantum asymptotic phase $\Phi_{q}$ gives appropriate phase values that increase with a constant frequency $\Omega_{q}$ even in the strong quantum regime of a quantum van der Pol model with quantum Kerr effect, even though the strong quantum effect strongly alters the dynamics of the system from the classical limit.

Although we introduced the definition of the quantum asymptotic phase $\Phi_{q}$ by analogy with the classical deterministic phase $\Phi_{c}$ and classical stochastic asymptotic phase $\Phi_{s}$, a fundamental difference exists between the quantum and classical cases. Specifically, in the quantum case, the system state is described by the density operator $\rho$ and the asymptotic phase $\Phi_{q}$ assigns a phase value on each $\rho$, while in the classical case, the phase function $\Phi_{c}$ or $\Phi_{s}$ assigns a phase value on each individual state ${\bm{X}}$, although the probability density function $p({\bm{X}})$ is used in defining the stochastic asymptotic phase $\Phi_{s}$. This difference arises because the system state can be described by a SDE representing a single stochastic trajectory of the noisy oscillator in the classical stochastic case, whereas the system state can only be characterized by the density operator $\rho$ representing the statistical state of the system in the quantum case.

Though we have taken the viewpoint that the asymptotic phase is defined for a single stochastic trajectory in the classical stochastic case in this study, we may also consider that the asymptotic phase is defined for the probability density $p({\bm{X}})$ rather than for the state ${\bm{X}}$ in the classical stochastic case. Then, the quantum asymptotic phase $\Phi_{q}$ as a function of the density matrix $\rho$ corresponds to the stochastic asymptotic phase $\Phi_{s}(p({\bm{X}}))=\arg\int p({\bm{X}})Q_{1}({\bm{X}})d{\bm{X}}$ as a function of the probability density $p({\bm{X}})$ in the classical stochastic case in Sec. II. We regarded this quantity as the exponential average of the asymptotic phase values defined for individual states in Eq. (14).

In the present approach, it is always possible to formally introduce a ’phase function’ for any oscillatory system as long as it has the decaying mode with a non-zero imaginary part by using the associated eigenoperator. However, this phase function does not necessarily play the role of a quantum asymptotic phase unless the system exhibits limit-cycle oscillations and synchronization. As an example, in Appendix D, we introduce a phase function of a damped quantum harmonic oscillator, which cannot be considered the quantum asymptotic phase since the system does not exhibit synchronization. Only when the system is nonlinear and exhibits synchronization, the asymptotic phase captures the synchronization dynamics of the system.

The asymptotic phase is the basis for developing phase reduction theory for classical nonlinear oscillators. For quantum oscillatory systems, however, even if we can introduce the asymptotic phase as described in this study, it does not necessarily mean that we can develop the phase reduction theory. This is because the quantum state $\rho$ may not be appropriately localized in the phase-space representation and hence it may not be reconstructed from the phase value even approximately. Still, as we showed for the case of the quantum vdP oscillator in the semiclassical regime kato2019semiclassical , we may approximately describe the dynamics of the quantum state by using the reduced phase equation and perform a detailed synchronization analysis in appropriate physical situations. It should be noted that we also encounter a similar problem in developing phase reduction theory for classical oscillators under strong noise. We also note that, even if we cannot derive a reduced phase description for quantum nonlinear oscillators, the asymptotic phase can be used to characterize the peculiar properties of quantum synchronization kato2021koopman2 .

The data that supports the findings of this study are available within the article.