Asymptotic phase and amplitude for classical and semiclassical stochastic oscillators via Koopman operator theory

The asymptotic phase is a fundamental quantity for the analysis of deterministic limit-cycle oscillators, and generalized definitions of the asymptotic phase for stochastic oscillators have also been proposed. In this article, we show that the asymptotic phase and also amplitude can be defined for classical and semiclassical stochastic oscillators in a natural and unified manner by using the eigenfunctions of the Koopman operator of the system. We show that the proposed definition gives appropriate values of the phase and amplitude for strongly stochastic limit-cycle oscillators, excitable systems undergoing noise-induced oscillations, and also for quantum limit-cycle oscillators in the semiclassical regime.

The asymptotic phase is a fundamental quantity for the analysis of deterministic limit-cycle oscillators, and generalized definitions of the asymptotic phase for stochastic oscillators have also been proposed. In this article, we show that the asymptotic phase and also amplitude can be defined for classical and semiclassical stochastic oscillators in a natural and unified manner by using the eigenfunctions of the Koopman operator of the system. We show that the proposed definition gives appropriate values of the phase and amplitude for strongly stochastic limit-cycle oscillators, excitable systems undergoing noise-induced oscillations, and also for quantum limit-cycle oscillators in the semiclassical regime.
Recently, the asymptotic phase and isochrons (level sets of the asymptotic phase), classical notions in the theory of nonlinear oscillations since Winfree [8] and Guckenheimer [9], have been studied from a viewpoint of the Koopman operator theory by Mauroy, Mezić, and Moehlis [10], and their relationship with the Koopman eigenfunction associated with the fundamental frequency of the oscillator has been clarified [10][11][12][13][14]. Moreover, they have shown that the (asymptotic) amplitude and isostables, which characterize deviation of the system state from the limit cycle and extend the Floquet coordinates [13,15,16] to the nonlinear regime, can be introduced naturally in terms of the Koopman eigenfunctions associated with the Floquet exponents with non-zero real parts [10][11][12][13][14]17]. By using the asymptotic phase and amplitude functions, we can obtain a reduced description of limit-cycle oscillators, which is useful for the analysis and control of synchronization dynamics of limit-cycle oscillators [18][19][20][21][22][23][24]. The theory can also be generalized to delay-differential systems [25] and spatially extended systems [26].
How to generalize the definition of the conventional asymptotic phase, which was essentially deterministic [8,9], to stochastic systems has been an intriguing problem [27][28][29][30][31][32][33][34][35][36]. When the stochasticity is sufficiently weak, the phase and also amplitude can be defined by using the drift term of the stochastic differential equation (SDE) describing the deterministic vector field of the oscillator. This approach can also be employed for quantum nonlinear oscillators in the semiclassical regime described by a quantum Fokker-Planck equation (FPE) [33,37]. However, this definition is no longer applicable to strongly stochastic oscillatory systems for which the deterministic vector field does not serve as a clear reference due to the strong effect of noise.
To cope with this problem, Schwabedal and Pikovsky [34] introduced a definition of the phase in terms of the mean first return time, and Thomas and Lindner [35] proposed a definition of the asymptotic phase in terms of the slowest decaying eigenfunction of the backward Fokker-Planck (Kolmogorov) operator describing the mean first passage time, both of which yield phase values that increase with a constant frequency on average for stochastic oscillations in a similar way to the ordinary asymptotic phase for deterministic oscillators. Recently, we pointed out that the definition of the stochastic asymptotic phase by Thomas and Lindner [35] can be seen as a natural extension of the deterministic definition from the viewpoint of the Koopman operator theory; namely, it is given by the argument of the Koopman eigenfunction associated with the fundamental frequency [38] (see also Reference [39]) and extended this idea to the definition of the asymptotic phase for quantum oscillatory systems.
In this article, based on the Koopman operator theory for stochastic systems, we propose a definition of the asymptotic phase and amplitude for strongly stochastic oscillators. They are introduced in terms of the eigenfunctions of the Koopman operator associated with the complex eigenvalues with the largest non-zero real part and with the largest non-zero real eigenvalue, respectively, which gives a natural extension of the definition in the deterministic case. The validity of the proposed definition is illustrated for stochastic limit-cycle oscillations and noise-induced oscillations of excitable systems using noisy Stuart-Landau [2,4] and FitzHugh-Nagumo [40,41] models as examples. Moreover, we apply the proposed definition of the stochastic phase and amplitude to a quantum limit-cycle oscillator in the semiclassical regime and show that they also yield appropriate results.

A. Classical Definition of the Asymptotic Phase and Amplitude
In this section, we review the definition of the asymptotic phase and amplitude for deterministic limit-cycle oscillators and discuss their relationship with the Koopman eigenfunctions [10-14, 19, 26]. We consider a deterministic dynamical systemẊ where X(t) ∈ R N is the system state at time t, A : R N → R N is a sufficiently smooth vector field representing the system dynamics, and the dot (˙) represents the time derivative. We assume that the system has an exponentially stable limit-cycle solution X 0 (t) with a natural period T and frequency ω = 2π/T , satisfying X 0 (t + T ) = X 0 (t).
We denote this limit cycle as χ and its basin of attraction as B χ ⊆ R N . Instead of the time t, we can parameterize a point on the limit cycle χ using a phase φ ∈ [0, 2π) as χ(φ) = X 0 (ωt) (0 ≤ t < T ), where the phase value φ = ωt increases linearly with time t from 0 to 2π (2π is identified with 0), and the origin of the phase φ = 0 is assigned to the state X 0 (0) without loss of generality. The linear stability of χ is characterized by the Floquet exponents λ j ∈ C (j = 0, 1, ..., N − 1) [15,42], which we sort in decreasing order of their real parts, i.e., Re(λ 0 ) ≥ Re(λ 1 ) ≥ Re(λ 2 ) ≥ · · · ≥ Re(λ N −1 ). Here, the exponent λ 0 is zero and associated with the phase direction tangent to χ, and the other exponents λ 1 , . . . , λ N −1 possess negative real parts because χ is exponentially stable and is associated with the amplitude directions deviating from χ. We further assume that λ 1 is real and λ 1 Re(λ 2 ), namely, the relaxation of the slowest decaying mode is non-oscillatory and much slower than the other faster decaying modes. Such a situation often occurs in realistic models of limit-cycle oscillators. We can then focus only on the slowest decaying mode and introduce a single real amplitude associated with it.
The asymptotic phase function Φ 0 : B χ → [0, 2π) and amplitude function R 0 : B χ → R of the limit cycle χ are defined in the basin B χ of χ such that are satisfied for all X ∈ B χ [13]. Here, the inner product is defined as a · b = N j=1 a j b j (the overline denotes complex conjugate) and ∇ = ∂/∂X represents the gradient with respect to X. As stated above, we focus only on the phase associated with λ 0 = 0 and the slowest decaying amplitude associated with real λ 1 . In general, we can introduce N −1 amplitude variables associated with N − 1 exponents λ 1 , . . . , λ N −1 , which are in general complex, and obtain a closed set of equations for the phase and amplitudes [13]. The level sets of the phase function are called isochrons [8,9] and those of the amplitude function are called isostables [10].
By using the above definition, we can introduce the phase and amplitude variables for the oscillator state X(t) ∈ B χ at time t as φ(t) = Φ 0 (X(t)) and r(t) = R 0 (X(t)), which obeẏ that is, the phase φ always increases with a constant frequency ω and the amplitude r decays exponentially with the rate λ 1 as X evolves in B χ toward χ. Note that the phase function is determined only up to an arbitrary constant and the scale of the amplitude function R 0 (X) is also arbitrary, because Φ 0 (X)+c 1 and c 2 R 0 (X) with arbitrary constants c 1 , c 2 ∈ R also satisfy Equation (2). Suppose that the initial state is X 0 ∈ B χ at t = 0. If we assign the phase φ(0; X 0 ) and amplitude r(0; X 0 ) to the initial state X 0 , we obtain φ(t; X 0 ) = ωt + φ(0; X 0 ) and r(t; X 0 ) = r(0; X 0 ) exp(λ 1 t), whose dependence on X 0 is explicitly shown.
By focusing only on the asymptotic phase and amplitude, we can perform phase-amplitude reduction (or isochronisostable reduction) of a limit-cycle oscillator [12,14,18,19], in which we reduce the dimensionality of the system dynamics from N to 2 and approximately describe it by a simple set of two-dimensional phase and amplitude equations. The phase equation has been extensively used for the analysis of weakly coupled limit-cycle oscillators [1][2][3][4][5][6], and the amplitude equation has also been used recently for the analysis and control of limit-cycle oscillators [18][19][20][21]24].

B. Koopman Operator Viewpoint
The asymptotic phase and amplitude introduced in the previous subsection are closely related to the Koopman operator of the system [18][19][20][21]. The Koopman operator U τ , which describes the evolution of a general observable g of the system state X ∈ R N , is defined as where g : R N → C is the observable, i.e., an observation function that maps a system state to an observed value, and S τ : R N → R N is a flow of the system satisfying X(t + τ ) = S τ X(t) for τ ≥ 0. When the flow S τ is analytic, it can be expanded as for |τ | 1. Considering an analytic observable g, we can expand it as Therefore, the infinitesimal time evolution of g can be expressed as The operator which appeared in Equation (3), can thus be interpreted as an infinitesimal generator of the Koopman operator U τ . For the limit-cycle oscillator described by Equation (1), we can easily confirm that the complex exponential of the phase function Φ 0 (X), is an eigenfunction of the operator A with an eigenvalue is satisfied for X ∈ B χ . Therefore, from the viewpoint of the Koopman operator theory, the asymptotic phase can be introduced as the argument (polar angle) of the Koopman eigenfunction Ψ 0 (X) associated with the eigenvalue iω, which is determined by the natural frequency ω of the oscillator [18][19][20][21], as where Arg represents the principal argument of a complex number in the range [0, 2π). Moreover, the asymptotic amplitude function R 0 (X) is nothing but the eigenfunction of the linear operator A associated with the eigenvalue λ 1 for X ∈ B χ , i.e., Thus, the Koopman operator theory provides a natural and unified definition of the asymptotic phase and amplitude, and the simplified Equation (3) in the phase-amplitude coordinates can be interpreted as a global linearization of the nonlinear dynamics of the limit-cycle oscillator by using the Koopman eigenfunctions. In References [18,43], Mauroy and Mezić pointed out these facts and explicitly calculated the phase and amplitude functions for several models of limit-cycle oscillators.

III. FOKKER-PLANCK EQUATION AND STOCHASTIC KOOPMAN OPERATOR A. Forward and Backward Fokker-Planck Equations
In the previous section, we considered deterministic limit-cycle oscillators and introduced the asymptotic phase and amplitude functions from the Koopman-operator viewpoint. Our aim in this study is to generalize the idea to stochastic oscillatory systems. In this section, we review some basic facts on the Fokker-Planck equations and stochastic Koopman operator for stochastic dynamical systems.
We consider a stochastic dynamical system described by a time-homogeneous SDE of Ito type [44][45][46], where X(t) ∈ R N is the system state at time t, A : R N → R N is a drift term representing the deterministic vector field of the oscillator, B : R N → R N ×N is a matrix characterizing the intensity of the noise, and W (t) is a Wiener process in R N representing the N -dimensional independent Gaussian-white noise. We assume that A and B satisfy the Lipschitz with some constant K for Equation (13) to possess a unique strong solution X(t) [44,46]. The FPE equivalent to the above SDE, describing the time evolution of the probability density function (PDF) p(X, t) : R N × R → R of X at time t is given by where is a (forward) Fokker-Planck operator. Here, the drift vector A : R N → R N is the same as in Equation (13) and D = BB T : R N → R N ×N is a symmetric diffusion matrix, where T indicates matrix transposition. We also assume that the functions A(X) and D(X) are smooth, satisfy the growth conditions with some constant M , and the uniform parabolicity (λ · D(X)λ) ≥ α|λ| 2 for all λ ∈ R N with a constant α > 0 in order that Equation (14) possesses a classical solution for t > 0 [44,[46][47][48].
The transition probability density p(X, t|Y , s), satisfying p(X, t) = p(X, t|Y , s) p(Y , s)dY for t > s and [44][45][46], obeys the forward FPE and also the corresponding backward FPE Here, the backward Fokker-Planck operator is the adjoint linear operator of L X with respect to the L 2 inner product of two functions G, H : R N → C, i.e., where the overline indicates complex conjugate and the integration is taken over the whole range of X here and hereafter.

B. Eigensystem of the Fokker-Planck Operators
The linear differential operators L X and L + X have the eigensystem {Λ k , P k , Q k } k≥0 consisting of the eigenvalue Λ k and eigenfunctions P k (X), Q k (X) satisfying and the biorthogonality conditions where k, l = 0, 1, 2, . . . and δ kl represents Kronecker's delta [35,45,49]. Here, Q k (X) is the complex conjugate of Q k (X), which is an eigenfunction of L + X associated with the eigenvalue Λ k , i.e., L + X Q k (X) = Λ k Q k (X). Because L X is a Fokker-Planck operator, the eigenvalue Λ 0 is zero and the associated eigenfunction P 0 (X) gives the stationary PDF of the FPE, i.e., L X P 0 (X) = 0, when appropriately normalized. All other eigenvalues have negative real parts and the associated eigenfunctions represent the relaxation eigenmodes of the FPE that eventually decay as t → ∞ [35,45,49].

C. Stochastic Koopman Operator
We here introduce the stochastic Koopman operator following Mezić [50] and discuss its relationship with the backward Fokker-Planck operator. Definition 1. For an observable g : R N → C and τ > 0, the stochastic Koopman operator U τ st is defined as represents the expectation over realizations of the stochastic flow S τ st of Equation (13) and p(Y , τ |X, 0) is the transition probability density satisfying Equation (16).
In the second expression of Equation (23), the expectation E[g(S τ st X)] is represented as an average over the transition probability density p(Y , τ |X, 0). The initial time can be taken as 0 without loss of generality because the process is time-homogeneous. We also introduce the infinitesimal generator of the stochastic Koopman operator.
Definition 2. For an observable g : R N → C, the infinitesimal generator A st of the stochastic Koopman operator U τ st is defined by From the above definitions, it can be shown that the infinitesimal generator of the stochastic Koopman operator is given by the backward Fokker-Planck operator. Thus, the infinitesimal generator of the stochastic Koopman operator is given by the backward Fokker-Planck operator, i.e., A st = L + X . Before proceeding to the definition of the asymptotic phase and amplitude, we show a result on the time evolution of the average of the eigenfunction Q k (k = 0, 1, 2, · · · ) of A st = L + X . Lemma 2. Let X(t) = S t st X 0 be a solution to Equation (13) with an initial condition X 0 ∈ R N , where S t st : R N → R N (t ≥ 0) is the stochastic flow of Equation (13). Then, the average for arbitrary X 0 , where E[·] represents the expectation over realizations of the stochastic flow S t st and p(X, t|X 0 , 0) is the transition probability density satisfying Equation (16).
We use the above result for discussing the evolution of the averaged phase and amplitude in the next section.

A. Stochastic Oscillatory Systems
The definitions of the phase and amplitude in Section II are based on the deterministic limit-cycle solution. These definitions are still applicable to noisy limit-cycle oscillators when the noise can be regarded as a weak perturbation [1][2][3][4][5]. However, they are no longer valid when the oscillator is subjected to stronger noise because we cannot rely on the deterministic limit-cycle solution in defining the phase and amplitude functions.
In Reference [35], Thomas and Lindner proposed a definition of the asymptotic phase for strongly stochastic oscillators without relying on the limit-cycle solution of the deterministic system, where they used the slowest decaying eigenfunction of the backward Fokker-Planck operator as the phase function based on the consideration of the mean first passage time. In this section, we show that their definition can be viewed as a natural extension of the deterministic definition in the sense that it is given by the argument of the Koopman eigenfunction associated with the fundamental frequency [38,39].

B. Assumptions on the Eigenvalues
Since we consider oscillatory stochastic systems, we introduce the following assumptions on the eigenvalues of the Fokker-Planck operator L X in Equation (15).
(i) We assume that the eigenvalues with the largest non-zero real part are given by a complex-conjugate pair, i.e., the slowest decaying eigenmode is oscillatory, and regard this eigenmode as the fundamental oscillation of the system. These eigenvalues are represented as where µ 1 < 0 and ω 1 > 0 characterize the decay rate and fundamental frequency of the oscillation, respectively.
(ii) We assume that the largest non-zero eigenvalue on the real axis, denoted as Λ 2 , is smaller than µ 1 , i.e., Λ 2 < Re Λ 1 , and consider that this eigenvalue characterizes the decay rate of the amplitude of the system, i.e., the deviation of the system state from the averaged oscillatory state.
When the system has a stable limit cycle as discussed in the previous section in the limit of vanishing noise intensity, these eigenvalues are expected to converge to iω and λ 1 of the deterministic limit cycle in the limit of vanishing noise, i.e., ω 1 → ω and Λ 2 → λ 1 .

C. Definition of the Asymptotic Phase Function
Thomas and Lindner [35] defined an asymptotic phase for the stochastic oscillatory system, Equation (13), by using the argument of the (complex conjugate of the) eigenfunction Q 1 (X) of the backward Fokker-Planck operator L + X associated with the eigenvalue Λ 1 characterized by the fundamental frequency ω 1 (in the notation of the present study), satisfying L + X Q 1 (X) = Λ 1 Q 1 (X), as and showed that this Φ(X) gives an appropriate phase value that increases with a constant frequency ω 1 with the evolution of X on average. They showed that, in the limit of vanishing noise where the system is described by the vector field A(X) possessing a stable limit-cycle solution, this definition of the asymptotic phase coincides with the deterministic definition in Section II [35]. In Reference [38], we pointed out that the above definition of the phase function by the backward Fokker-Planck operator can also be understood from the viewpoint of the Koopman operator theory. In what follows, we introduce the asymptotic phase and also the amplitude from the Koopmanoperator viewpoint.
Let us rephrase the above definition of the asymptotic phase Φ(X) for stochastic oscillators from the Koopmanoperator viewpoint.
Definition 3. We define the asymptotic phase Φ(X) of the oscillator state X ∈ R N described by Equation (13) Note that the above phase Φ(X) has a discontinuity at 2π, which causes difficulty in taking ensemble averages of Φ(X) over realizations of X. Rather, as in the standard convention in directional statistics [52], we consider a 'wrapped' distribution of the phase values and use the circular mean to calculate the average phase. This is accomplished by taking the ensemble average of Q 1 (X) over many realizations and then calculate its argument, rather than calculating the ensemble average of Arg Q 1 (X).
Definition 4. We define the averaged asymptotic phase of the stochastic oscillator described by Equation (13) at time t, started from an initial condition X 0 at time 0, as where E[·] represents the expectation over realizations of the stochastic flow S t st of Equation (13) and p(X, t|X 0 , 0) is the transition probability density satisfying Equation (16).
Let us confirm that the above definition of the phase function yields appropriate phase values on average. Lemma 3. The average asymptotic phase in Equation (31) increases with a constant frequency ω 1 , i.e., for arbitrary X 0 ∈ R N .
where Λ 1 = µ 1 + iω 1 . By integration, we obtain where we used that S 0 st X 0 = X 0 . The averaged asymptotic phase is thus given by which yields Equation (32) by differentiation by t.
Thus, the averaged asymptotic phase φ(t; X 0 ) of the oscillator satisfieṡ namely, φ(t; X 0 ) increases with a constant frequency ω 1 on average for any X 0 . This result indicates that the definition of the asymptotic phase in Equation (29) for the stochastic oscillators by Thomas and Lindner [35] is a natural extension of the definition in Equation (11) for the deterministic oscillators from the Koopman-operator viewpoint.

D. Definition of the Amplitude Function
We have seen that the definition of the stochastic asymptotic phase by using the backward Fokker-Planck operator can be naturally interpreted as a generalization of the deterministic definition from the viewpoint of the Koopman operator theory. Furthermore, as explained in Section II, the asymptotic amplitude can be naturally defined by using the Koopman eigenfunction associated with the largest non-zero real eigenvalue in deterministic systems. Therefore, to generalize the definition of the amplitude to stochastic oscillators, it appears natural to use the eigenfunction of the stochastic Koopman operator.
Definition 5. We define the amplitude of the oscillator state X ∈ R N described by Equation (13) by using the Koopman eigencfunction Q 2 (X) of L + X in Equation (18) associated with the largest non-zero real eigenvalue Λ 2 as Let us confirm that the above definition of the amplitude function yields appropriate amplitude values on average. Definition 6. We define the averaged amplitude of the stochastic oscillator described by Equation (13) at time t, started from the initial condition X(0) = X 0 at time 0, as where E[·] represents the expectation over realizations of the stochastic flow S t st of Equation (13) and p(X, t|X 0 , 0) is the transition probability density satisfying Equation (16).
Proof. From Lemma 2, Thus, the averaged amplitude r(t; X 0 ) obeysṙ and hence r(t; X 0 ) = e Λ2t r(0; X 0 ) as expected for any X 0 . Since Λ 2 is real and negative by assumption, r(t; X 0 ) decays exponentially with time. This result indicates that the definition of the amplitude in Equation (37) for stochastic oscillators is a natural extension of the definition in Equation (12) for the deterministic systems from the Koopman-operator viewpoint. As we illustrate in the next section with a few examples, the above definition yields an amplitude value that decays linearly with t on average and characterizes the deviation of the system state from the steady oscillation.

E. Limit of Vanishing Noise Intensity
Before proceeding to examples, we point out that the results for the stochastic oscillators formally reduce to the results for deterministic limit-cycle oscillators in the limit of vanishingly small noise.
If we assume that the noise does not exist, i.e., D(X) → 0 in the forward and backward Fokker-Planck Equations (14) and (17), we obtain the forward and backward Liouville equations [45,46,48,53], ∂ ∂t p(X, t|Y , s) = L X p(X, t|Y , s), where the forward Liouville operator is given by and the backward Liouville operator is given by Because the backward Liouville operator L + X coincides with the infinitesimal generator of the Koopman operator A in the deterministic case given in Equation (8), the Koopman eigenfunction Ψ 0 (X) of A in Equation (10) is an eigenfunction of L + X with an eigenvalue iω. Thus, the definition of the asymptotic phase for stochastic oscillators in Equation (30) can be considered a natural generalization of the definition of the asymptotic phase for deterministic oscillators in Equation (11). Similarly, the Koopman eigenfunction R 0 (X) of A in Equation (12) is an eigenfunction of L + X with an eigenvalue Λ 2 = λ 1 , so the definition of the amplitude for stochastic oscillators in Equation (37) also corresponds to that for deterministic oscillators in Equation (12).

V. EXAMPLES A. Numerical Methods
To demonstrate the validity of the phase and amplitude functions introduced in Section IV, we consider two classical examples of noisy limit-cycle oscillators, i.e., the Stuart-Landau model [2,4] and the FitzHugh-Nagumo model [40,41]. We numerically calculate the eigenvalues and eigenfunctions of the backward Fokker-Planck operator and evaluate the phase and amplitude functions. We also analyze a quantum limit-cycle oscillator in the semiclassical regime, which can be described by the same stochastic differential equations as those for the classical noisy limit-cycle oscillators.
In the numerical calculations, we truncated the state space and approximated it as a finite square domain −D ≤ x ≤ D, −D ≤ y ≤ D with a large enough value of D. In all models considered, the stationary PDF of the FPE rapidly decayed with the distance from the origin and took numerically negligible values at the edges of the domain. We discretized the domain into N × N grids and represented the PDF as a N 2 -dimensional vector. We then represented the operator L + X as a N 2 × N 2 matrix, calculated the eigenvalues and eigenvectors, and obtained the phase and amplitude functions.

B. Example 1: Noisy Stuart-Landau Model
As the first example, we consider the Stuart-Landau model (the normal form of the supercritical Hopf bifurcation [2,42]) under the effect of noise, described by where x and y are real variables, a, b, c, and d are real parameters, W x and W y are independent Wiener processes, and D x and D y represent the intensities of the noise acting on x and y, respectively. The noiseless system with D x = D y = 0 has a stable limit cycle with a natural frequency ω = b − ad/c and the largest non-zero Floquet exponent λ 1 = −2a when a > 0 and c > 0. For this system, we can explicitly calculate the limit cycle and the phase and amplitude functions as [4] (x 0 (φ), y 0 (φ)) T = a c (cos φ, sin φ) T , where C 0 is an arbitrary scalar constant. The basin B χ of this limit cycle χ is the whole complex plane except the origin. In the following numerical simulations, we set the parameter values as (a, b, c, d, D x , D y ) = (0.5, 1.5, 0.25, 0.25, 1, 1). The natural frequency and the largest non-zero Floquet exponent are ω = 1 and λ 1 = −1, respectively. It is noted that the fundamental frequency ω 1 = Im Λ 1 and decay rate Λ 2 under the effect of noise are generally different from these deterministic values. We used D = 3.6 and N = 151 for the numerical analysis. Figure 1a shows the eigenvalues of the Koopman operator L + X near the imaginary axis obtained numerically, where the eigenvalues Λ 1 = µ 1 + iω 1 and Λ 2 are shown by orange and red dots, respectively. The rightmost branch of the eigenvalues is approximately given by a parabolaλ n = iω 1 n − µ 1 n 2 (n = 0, ±1, ±2, . . .) passing through Λ 1 [35]. Figure 1b,c show the phase function Φ(x, y) and amplitude function |R(x, y)| associated with Λ 1 = µ 1 + iω 1 and Λ 2 , respectively. We can observe that a circular region representing the local minima of the amplitude exists along the limit-cycle solution in the deterministic case and the phase increases from 0 to 2π along this circle. In contrast to the deterministic case, Equation (46), the amplitude does not diverge at the unstable fixed point at the origin (x, y) = (0, 0), because the system state can escape from this point in a finite time due to the effect of noise.
To confirm that these functions yield appropriate values of the phase and amplitude on average, we obtained 10000 trajectories by direct numerical simulations of the Equation (45) from the initial point (x 0 , y 0 ) = (−1.5, −1.5) and calculated the averaged phase φ = Arg Q 1 (x, y) and amplitude |r| = |Q 2 (x, y)| , where · represents a sample average over all obtained trajectories. Figure 1d,e show that these values are in good agreement with the analytical solutions φ = ω 1 t + φ 0 and |r| = |r 0 | exp(Λ 2 t), where the fundamental frequency ω 1 = 0.728 and the decay rate Λ 2 = −1.680 are numerically evaluated from the eigenvalues plotted in Figure 1a. For comparison, we also show the analytical solutions for the deterministic case without noise (D x = D y = 0), namely, φ = ωt + φ 0 = t + φ 0 and |r| = |r 0 | exp(λ 1 t) = |r 0 | exp(−t). The averaged phase increases more slowly and the averaged amplitude decays more quickly than those in the deterministic case due to the effect of noise.

C. Example 2: Noisy FitzHugh-Nagumo Model
Next, we consider the FitzHugh-Nagumo model [40,41] subjected to noise, described by where x and y are real variables, a 1 , b 1 , and η 1 are real parameters, W x and W y are independent Wiener processes, and D x and D y represent the intensities of the noise, respectively. First, we consider parameter set (A): (a 1 , b 1 , η 1 , D x , D y ) = (1/3, 0.5, 0.5, 0.2, 0.2), at which the deterministic vector field possesses a stable limit-cycle solution with the natural frequency ω = 0.588 and the largest non-zero Floquet exponent λ 1 = −1.11. We used D = 4.2 and N = 151 for the numerical analysis. Figure 2a shows the eigenvalues of the Koopman operator L + X near the imaginary axis obtained numerically, where Λ 1 = µ 1 + iω 1 and Λ 2 are shown by orange and red dots, respectively. The rightmost branch of the eigenvalues is approximately a parabolaλ n = iω 1 n − µ 1 n 2 (n = 0, ±1, ±2, . . .) passing through Λ 1 , which is qualitatively similar to the one for the noisy Stuart-Landau model in Figure 1a. Figure 2b,c show the phase Φ(x, y) and amplitude |R(x, y)| associated with Λ 1 = µ 1 + iω 1 and Λ 2 . As in the case of the noisy Stuart-Landau model, a circular region corresponding to the local minima of the amplitude function exists around the deterministic limit-cycle solution and the phase increases along this region. The amplitude does not diverge at the unstable fixed point due to the effect of noise.
We calculated the time evolution of Φ(x, y) and |R(x, y)| by direct numerical simulations of Equation (47)  phase φ = Arg Q 1 (x, y) and amplitude |r| = |Q 2 (x, y)| . They are in good agreement with the analytical solutions φ = ω 1 t + φ 0 and |r| = |r 0 | exp(Λ 2 t), where the fundamental frequency ω 1 = 0.582 and the decay rate Λ 2 = −0.778 are numerically evaluated from the eigenvalues plotted in Figure 2a. For comparison, we also show the analytical solution for the deterministic case without noise (D x = D y = 0), namely, φ = ωt + φ 0 = 0.588t + φ 0 and |r| = |r 0 | exp(λ 1 t) = |r 0 | exp(−1.11t). The phase evolves more slowly and also the amplitude decays more slowly than those in the deterministic case due to the effect of noise.
Next, we consider parameter set (B): (a 1 , b 1 , η 1 , D x , D y ) = (1/3, 1.05, 0.25, 0.1, 0.1). In this case, the deterministic vector field does not have a stable limit-cycle, but the system is close to a supercritical Hopf bifurcation of a limit cycle. Thus, relatively regular noise-induced oscillations occur even though the system does not have a deterministic limit cycle, a phenomenon known as the coherence resonance [54][55][56]. We used D = 3.9 and N = 151 for the numerical analysis. Figure 3a shows the eigenvalues of the Koopman operator L + X near the imaginary axis obtained numerically, where Λ 1 = µ 1 + iω 1 and Λ 2 are shown by orange and red dots, respectively. Figure 3b,c show the phase Φ(x, y) and amplitude |R(x, y)| associated with Λ 1 = µ 1 + iω 1 and Λ 2 , respectively. Interestingly, although the deterministic system does not have a limit-cycle solution, we can still observe in Figure 3b,c a circular region representing the local minima of the amplitude function. This region corresponds to the noise-induced oscillations and the phase increases along this circular region. The averaged phase φ = Arg Q 1 (x, y) and amplitude r = |Q 2 (x, y)| show good agreement with the analytical solutions φ = ω 1 t + φ 0 and |r| = |r 0 | exp(Λ 2 t), where the fundamental frequency ω 1 = 0.287 and the decay rate Λ 2 = −0.815 are numerically evaluated from the eigenvalues in Figure 3a.
Thus, we can introduce the phase and amplitude functions also in this case without a deterministic limit cycle by using the present definition using the Koopman eigenfunctions.

D. Example 3: Semiclassical Stuart-Landau Model
Finally, we apply the proposed definition of the stochastic phase and amplitude functions to a quantum limit-cycle oscillator in the semiclassical regime. As an example, we use the quantum Stuart-Landau model [57,58] (also known as the quantum van der Pol model [59]) with a Kerr effect [33,60] in quantum optics.
Employing the phase space approach [61,62], the system state can be represented by a Wigner function W (x, y) (see Reference [33,60] for details). In the semiclassical regime, the quantum noise is sufficiently weak and W (x, y) approximately obeys a quantum FPE, which has the same form as the ordinary FPE for classical systems. Thus, we can derive the corresponding Ito SDE from the quantum FPE as where β(x, y) = γ1 2 +2γ 2 x 2 + y 2 − 1 2 , ω 0 is a frequency parameter of the oscillator, K represents the Kerr parameter, and γ 1 and γ 2 are the decay rates for the negative damping and nonlinear damping, respectively, and W x and W y are independent Wiener processes. The semiclassical approximation is valid when γ 2 and K are sufficiently small [33,60]. The deterministic part of Equation (48) is equivalent to the Stuart-Landau model used in Example 1 and only the coefficient of the noise term differs. Therefore, the deterministic limit cycle and the phase and amplitude functions can be obtained from the results for Equation (46), where the parameters are given by a = γ1+2γ2 2 , b = ω 0 +2K, c = γ 2 , d = 2K. We set the parameter values as (γ 1 , γ 2 , ω 0 , K) = (1, 0.05, 1, 0.025). With these values, the deterministic vector field of Equation (48) possesses a stable limit-cycle solution with the natural frequency ω = ω 0 −Kγ 1 /γ 2 = 0.5 and the largest non-zero Floquet exponent λ 1 = γ 1 + 2γ 2 = −1.1. We used D = 6.7 and N = 151 for the numerical analysis. Figure 4a shows the eigenvalues of the Koopman operator L + X , where Λ 1 = µ 1 + iω 1 and Λ 2 are indicated by orange and red dots, respectively. Figure 4b,c show the phase function Φ(x, y) and amplitude function R(x, y) associated with Λ 1 = µ 1 + iω 1 and Λ 2 . For comparison, we also show in Figure 4d the deterministic phase function Φ 0 (x, y) in Equation (46). The stochastic phase function Φ(x, y) in Figure 4b is slightly different from the deterministic phase function Φ 0 (x, y) in Figure 4d, in particular near the origin, because the stochastic phase function includes the effect of weak quantum noise. Figure 4c shows that the amplitude function takes the minima around the deterministic limit cycle. We note that the amplitude function does not diverge at the origin in contrast to the deterministic case. In (e,f ), averaged results over 10,000 trajectories (orange and red thin lines) and analytical solutions for the semiclassical case (blue-dotted lines) and results for the deterministic case (green-dotted lines) are shown.
Here, we used a color map with the maximum value of 0.005 to enhance the local minima of the amplitude function, and the region near the origin where the amplitude is larger than 0.005 is shown in the same color (the true maximum amplitude is 0.0876 at the origin). Figure 4e,f show the time evolution of the values of the phase and amplitude averaged over 10,000 trajectories obtained by direct numerical simulations of the semiclassical Equation (45) from the initial point (x 0 , y 0 ) = (2.5, 0). They show good agreement with the analytical solutions φ = ω 1 t + φ 0 and |r| = |r 0 | exp(Λ 2 t), where the values of the fundamental frequency ω 1 = 0.496 and the decay rate Λ 2 = −0.798 are numerically evaluated from the eigenvalues shown in Figure 4a.
Thus, the present definition of the phase and amplitude functions is also applicable to a quantum Stuart-Landau model in the semiclassical regime and yields reasonable values.

VI. DISCUSSION
We proposed a definition of the asymptotic phase and amplitude functions for stochastic oscillatory systems by generalizing the definition for deterministic limit-cycle oscillators on the basis of the Koopman operator theory, motivated by the definition of the stochastic asymptotic phase introduced by Thomas and Lindner [35].
The proposed asymptotic phase and amplitude for strongly stochastic oscillatory systems may be used for systematic and quantitative analysis of synchronization phenomena in stochastic oscillators. We may also be able to develop a phase-amplitude reduction theory for strongly stochastic oscillators by using the phase and amplitude functions, which allows us to reduce the system dynamics subjected to weak external inputs to a simple two-dimensional set of equations. Such theories will facilitate detailed analysis, control, and optimization of noise-induced oscillatory phenomena, including the coherence-resonance and self-induced-stochastic-resonance oscillations [63,64].
It will also be interesting to introduce amplitude functions for strongly quantum oscillatory systems on the basis of the Koopman operator theory [38]. We have recently defined the quantum asymptotic phase of strongly quantum oscillators by using the eigenoperator of the adjoint Liouville superoperator and found that it can largely differ from the asymptotic phase in the classical limit [38]. The amplitude function, which can be defined similarly, may also be very different from the classical counterpart and characterize the quantum signatures of synchronization observed in the strong quantum regime.

VII. CONCLUSIONS
We proposed a definition of the asymptotic phase and amplitude functions for stochastic oscillatory systems. The proposed phase and amplitude functions are introduced in terms of the backward Fokker-Planck operator, which can be interpreted as the Koopman operator for the stochastic system. The validity of the phase and amplitude functions was numerically demonstrated for noisy Stuart-Landau and FitzHugh-Nagumo models and also for a quantum Stuart-Landau model in the semiclassical regime.
Note added -During the preparation of this article , we noticed a new, closely related study by Pérez-Cervera, Lindner, and Thomas [65], which introduced the isostables (level sets of the asymptotic amplitude function) for stochastic oscillators on the basis of the backward Kolmogorov equation and analyzed the examples of a spiral sink, a noisy Stuart-Landau oscillator, and a noisy heteroclinic oscillator. Our present results differ from Reference [65] in that we explicitly discussed the relationship with the Koopman operator theory and analyzed a noisy excitable system and quantum limit-cycle oscillator in the semiclassical regime, in addition to noisy limit-cycle oscillators. We thus believe our results provide different insights and are complementary to Reference [65].