Phase oscillator model for noisy oscillators

TL;DR

This paper derives a phase oscillator model for noisy oscillators, showing that phase noise behaves as a drift-diffusion process and that amplitude fluctuations influence the expected frequency, extending the Kuramoto model.

Contribution

It introduces a phase reduction for noisy oscillators that accounts for amplitude fluctuations and their impact on frequency shifts, providing a more accurate model than traditional approaches.

Findings

Phase noise is modeled as a drift-diffusion process.

Amplitude fluctuations affect the expected oscillator frequency.

Different approximation levels improve phase reduction accuracy.

Abstract

The Kuramoto model has become a paradigm to describe the dynamics of nonlinear oscillator under the influence of external perturbations, both deterministic and stochastic. It is based on the idea to describe the oscillator dynamics by a scalar differential equation, that defines the time evolution for the phase of the oscillator. Starting from a phase and amplitude description of noisy oscillators, we discuss the reduction to a phase oscillator model, analogous to the Kuramoto model. The model derived shows that the phase noise problem is a drift-diffusion process. Even in the case where the expected amplitude remains unchanged, the unavoidable amplitude fluctuations do change the expected frequency, and the frequency shift depends on the amplitude variance. We discuss different degrees of approximation, yielding increasingly accurate phase reduced descriptions of noisy oscillators.

Figures

Equations

$$d{\mathbf{X}}_t=\mathbf{a}({\mathbf{X}}_t)dt+\epsilon \mathbf{B}({\mathbf{X}}_t)d{\mathbf{W}}_t$$

$${\mathbf{u}}_1(t)=\frac{\mathbf{a}({\mathbf{x}}_s(t))}{|\mathbf{a}({\mathbf{x}}_s(t))|}$$

$$\mathbf{x}=\mathbf{h}(\theta ,\mathbf{R})={\mathbf{x}}_s(\theta (t))+\mathbf{Y}(\theta (t))\mathbf{R}(t)$$

$$d\theta =$$

$$d\mathbf{R}=$$

$$a_{\theta }(\theta ,\mathbf{R})=$$

$${\hat{a}}_{\theta }(\theta ,\mathbf{R})=$$

$${\mathbf{B}}_{\theta }(\theta ,\mathbf{R})=$$

$$\mathbf{L}(\theta )=$$

$${\mathbf{a}}_{\mathbf{R}}(\theta ,\mathbf{R})=$$

$${\hat{\mathbf{a}}}_{\mathbf{R}}(\theta ,\mathbf{R})=$$

$${\mathbf{B}}_{\mathbf{R}}(\theta ,\mathbf{R})=$$

$$I_{{\mathbf{x}}_s(0)}=\left\{\mathbf{x}(0)\in \text{mathds}{R}^n/{\mathbf{x}}_s:\underset{t\to +\infty }{\lim }||\mathbf{x}(t)-{\mathbf{x}}_s(t)||=0\right\}$$

$$\frac{d\mathbf{y}(t)}{dt}=\mathbf{A}(t)\mathbf{y}(t)$$

$$\big{[}r(t)\boldsymbol{u}_{1}(t),\boldsymbol{u}_{2}(t),\ldots,\boldsymbol{u}_{n}(t)\big{]}=\boldsymbol{P}(t)$$

$$d\theta =$$

$$d\mathbf{R}=$$

$$d\theta =dt+\epsilon {\mathbf{B}}_1(\theta ,0)d{\mathbf{W}}_t$$

$$d\theta=\big{(}1+\varepsilon^{2}\hat{a}_{1}(\theta,\mathbf{0})\big{)}dt+\varepsilon\boldsymbol{B}_{1}(\theta,\mathbf{0})\,d\boldsymbol{W}_{t}$$

$$\begin{array}[]{rl}\dfrac{\partial p}{\partial t}=&-\dfrac{\partial}{\partial\theta}\bigg{\{}\Big{[}1+a_{\theta}(\theta,\boldsymbol{R})+\varepsilon^{2}\hat{a}_{\theta}(\theta,\boldsymbol{R})\Big{]}p\bigg{\}}\\ &-\displaystyle{\sum_{i=1}^{n-1}}\dfrac{\partial}{\partial R_{i}}\bigg{\{}\Big{[}\sum_{j=1}^{n-1}L_{ij}(\theta)R_{j}+a_{R_{i}}(\theta,\boldsymbol{R})+\varepsilon^{2}\hat{a}_{R_{i}}(\theta,\boldsymbol{R})\Big{]}p\bigg{\}}+\dfrac{\varepsilon^{2}}{2}\dfrac{\partial^{2}}{\partial\theta^{2}}\Big{[}\boldsymbol{B}_{\theta}(\theta,\boldsymbol{R})\boldsymbol{B}_{\theta}^{T}(\theta,\boldsymbol{R})\,p\Big{]}\\[4.30554pt] &+\varepsilon^{2}\displaystyle\sum_{i=1}^{n-1}\dfrac{\partial^{2}}{\partial\theta\partial R_{i}}\Big{[}\boldsymbol{B}_{\theta}(\theta,\boldsymbol{R})\boldsymbol{B}_{R_{i}}^{T}(\theta,\boldsymbol{R})\,p\Big{]}+\dfrac{\varepsilon^{2}}{2}\displaystyle{\sum_{i,j=1}^{n-1}}\dfrac{\partial^{2}}{\partial R_{i}\partial R_{j}}\Big{[}\boldsymbol{B}_{R_{i}}(\theta,\boldsymbol{R})\boldsymbol{B}_{R_{j}}^{T}(\theta,\boldsymbol{R})\,p\Big{]}\end{array}$$

$$\frac{\partial p}{\partial t}=(L_0+\epsilon L_1+{\epsilon }^2{L}_2)p$$

$$L_{0}p$$

$$L_{1}p$$

$$+\sum _{i=1}^{n-1}\frac{{\partial }^2}{\partial \theta \partial {\tilde{R}}_i}\left({\mathbf{B}}_{\theta }{\mathbf{B}}_{R_i}^{T}p\right)$$

$$L_{2}p$$

$$d\theta$$

$$d\tilde{\mathbf{R}}$$

$$\frac{d\mathbf{\mu }}{d\theta }$$

$$\frac{d\mathbf{\sigma }}{d\theta }$$

$$p_{st}(\theta ,\tilde{\mathbf{R}})=\frac{1}{\sqrt{(2\pi {)}^n|\mathbf{\sigma }|}}\exp \left[-\frac{1}{|\mathbf{\sigma }|}\sum _{i,j=1}^{n-1}|{\mathbf{\sigma }}_{ij}|{\tilde{R}}_i{\tilde{R}}_j\right]$$

$$Pf(\theta ,\tilde{\mathbf{R}},t)=p_{st}(\theta ,\tilde{\mathbf{R}})\int d{R}_1\dots \int d{R}_{n-1}f(\theta ,\tilde{\mathbf{R}},t)$$

$$Pp(\theta ,\tilde{\mathbf{R}},t)=p_{st}(\theta ,\tilde{\mathbf{R}})\hat{p}(\theta ,t)$$

$$\frac{\partial Pp}{\partial t}=PLPp+PL{e}^{tQL}Qp+{\int }_0^{t}K(t-s)Ppds$$

$$\frac{\partial Pp(\theta ,\tilde{R},t)}{\partial t}\approx PLPp(\theta ,\tilde{R},t)$$

$$P{L}_{0}Pp$$

$$P{L}_{1}Pp$$

$$P{L}_{2}Pp$$

$$\frac{\partial \hat{p}}{\partial t}$$

Full text

1 Introduction

Nonlinear oscillators are key components of many modern electronic devices. They are used in communication systems for modulation and demodulation of data, and in digital systems to provide a time reference frame to synchronize operations. An ideal oscillator would exhibit a perfectly periodic behavior, represented by a limit cycle in its state space. However, the output of actual oscillators is always corrupted by different types of disturbances, such as internal noise sources, thermal noise, and interactions with the environment. Noise sources in electronic circuits cause diffusion (or jitter) of the oscillator phase, ultimately inducing phase noise and time jitter [1, 2]. In the frequency domain, the consequences are a broadening and a shift in the power spectrum peaks, that may eventually produce interference on nearby frequency bands in transmitters and receivers [3, 4].

Phase noise is a typical time domain phenomenon arising from the random walk process that the phase undergoes [5]. Therefore, time domain models are particularly well suited for the analysis of phase noise. The Kuramoto model has become a paradigm to describe the dynamics of nonlinear oscillator under the influence of external perturbations, both deterministic and stochastic [6, 7]. Kuramoto model is based on the idea to describe the state of a nonlinear oscillator by a single scalar variable, the phase, that represents the position of the system on a reference orbit. To derive the model, the solution of the perturbed system is written as a time shifted version of the unperturbed one, and an equation is derived for the time shift. The time shift is interpreted as a phase shift normalized by the free running frequency, and the equation for the phase shift is called phase oscillator. The amplitude dynamics is completely neglected, under the hypothesis that the limit cycle is strongly stable and that amplitude deviations are immediately absorbed. This assumption is reasonable for most electronic oscillators, but in other cases, e.g. biological oscillators, it is more taken for mathematical convenience than being physically plausible.

Phase oscillator models have become a paradigm to describe the synchronization of an oscillator with an external periodic signal [8, 9, 10], phase locking of coupled oscillators [11, 12, 13], and to investigate the phase noise in oscillators [2, 3, 14]. However, circuits and systems are usually described in terms of state variables and state equations, and the problem of deriving the phase oscillator model starting from the state equations arises. The exact phase oscillator model can be found only for few, trivial nonlinear oscillators. In all other cases, e.g. phase response curves, one must resort to numerical methods that are either approximate in nature, or unsuitable for oscillators of order higher than the second [15, 16]. In the seminal works [2, 3], the authors derived a phase oscillator model for nonlinear oscillators subject to weak white Gaussian noise. It was shown that, to the first order, the phase noise problem is a diffusion process. In [18, 19] it was shown that phase noise is best described as a drift–diffusion process, i.e. noise also produces a shift in the oscillation frequency. In [20, 21] a rigorous description of noisy oscillators is derived, in terms of phase and amplitude deviation. The resulting equations are exact, they are not restricted to the weak noise limit and they are valid for oscillators of any order. The phase and amplitude equations represent the ideal starting point for the derivation of phase oscillator models.

This paper describes a phase oscillator model for nonlinear oscillators subject to white Gaussian noise, analogous to the celebrated Kuramoto model. It is shown that in a neighborhood of the limit cycle, the phase variable introduced in the proposed model coincide with the phase function defined using the concept of isochrons. The novel phase oscillator model gives better approximations of the full dynamics than other reduced models previously proposed in literature. The relationships with these previous models are discussed. In particular, it is shown that previous models can be derived from the novel one, introducing increasing order of approximation. As an example, the theory is applied to Duffing–van der Pol system with modulated noise.

2 Phase and amplitude deviation equations for noisy oscillators

Different types of random disturbances in electronics and biology can be modeled as white Gaussian noise. As a consequence electronics and biological oscillators with noise can be conveniently described by stochastic differential equation (SDEs) [22]

[TABLE]

where $\boldsymbol{X}_{t}:\mathds{R}\mapsto\mathds{R}^{n}$ is a stochastic process describing the state of the oscillator, $\boldsymbol{a}:\mathds{R}^{n}\mapsto\mathds{R}^{n}$ is a vector valued function that defines the oscillator dynamics, $\boldsymbol{B}:\mathds{R}^{n}\mapsto\mathds{R}^{n,m}$ is a real valued $n\times m$ matrix, $\varepsilon$ is a parameter that measures the noise intensity and $\boldsymbol{W}_{t}:\mathds{R}\mapsto\mathds{R}^{m}$ is a vector of Brownian motion components (the integral of a white noise). We shall assume $\boldsymbol{a}\in\mathcal{C}^{k}(\Omega\subseteq\mathds{R}^{n})$, with $k\geq 1$ and that $\boldsymbol{B}$ satisfy a Lipschitz condition, to guarantee the existence and uniqueness of the solution [22]. Following [2], we shall interpret (1) as an Itô SDE.

In absence of noise (for $\varepsilon=0$) eq. (1) reduces to a system of ordinary differential equations (ODEs). We assume that these ODEs admit an asymptotically stable $T$–periodic limit cycle denoted by $\boldsymbol{x}_{s}(t)$. We define the unit vector tangent to the limit cycle

[TABLE]

Together with $\boldsymbol{u}_{1}(t)$ we consider other $n-1$ linear independent vectors $\boldsymbol{u}_{2}(t),\ldots,\boldsymbol{u}_{n}(t)$, such that the set $\{\boldsymbol{u}_{1}(t),\ldots,\boldsymbol{u}_{n}(t)\}$ is a basis for $\mathds{R}^{n}$, for all $t$. Given the matrix $\boldsymbol{U}(t)=[\boldsymbol{u}_{1}(t),\ldots,\boldsymbol{u}_{n}(t)]$, we define the reciprocal vectors $\boldsymbol{v}_{1}^{T}(t),\ldots,\boldsymbol{v}_{n}^{T}(t)$ to be the rows of the inverse matrix $\boldsymbol{V}(t)=\boldsymbol{U}^{-1}(t)$. Thus $\{\boldsymbol{v}_{1}(t),\ldots,\boldsymbol{v}_{n}(t)\}$ also span $\mathds{R}^{n}$ and the bi–orthogonality condition $\boldsymbol{v}_{i}^{T}\boldsymbol{u}_{j}=\boldsymbol{u}_{i}^{T}\boldsymbol{v}_{j}=\delta_{ij}$ for all $t$, holds. We shall also use the matrices $\boldsymbol{Y}(t)=[\boldsymbol{u}_{2}(t),\ldots,\boldsymbol{u}_{n}(t)]$, $\boldsymbol{Z}(t)=[\boldsymbol{v}_{2}(t),\ldots,\boldsymbol{v}_{n}(t)]$, and the modulus of the vector field evaluated on the limit cycle, $r(t)=|\boldsymbol{a}(\boldsymbol{x}_{s}(t))|$.

The phase of a nonlinear oscillator can be interpreted as an elapsed time from an initial reference point. Thus, the phase represents a new parametrization of the trajectories. Together with the phase function $\theta:\mathds{R}\mapsto\mathds{R}$ we shall consider an amplitude deviation function $\boldsymbol{R}:\mathds{R}\mapsto\mathds{R}^{n-1}$, with $\theta,\boldsymbol{R}\in\mathcal{C}^{m}(\mathds{R})$, $m\geq 2$. The amplitude deviation $\boldsymbol{R}(\boldsymbol{x}(t))$ is interpreted as an orbital deviation from the limit cycle. In [20, 21] the following theorem was established, that gives the phase and amplitude deviation equations corresponding to the noisy oscillator (1)

Theorem 1

Consider the Itô SDEs (1) such that the ODEs obtained setting $\varepsilon=0$ admit a $T$–periodic limit cycle $\boldsymbol{x}_{s}(t)$. Let $\{\boldsymbol{u}_{1}(t),\ldots,\boldsymbol{u}_{n}(t)\}$ and $\{\boldsymbol{v}_{1}(t),\ldots,\boldsymbol{v}_{n}(t)\}$ be two reciprocal bases such that $\boldsymbol{u}_{1}(t)$ is given by (2) and such that the bi–orthogonality condition $\boldsymbol{v}_{i}^{T}\boldsymbol{u}_{j}=\boldsymbol{u}_{i}^{T}\boldsymbol{v}_{j}=\delta_{ij}$ holds. Consider the coordinate transformation

[TABLE]

Then a neighborhood of the limit cycle $\boldsymbol{x}_{s}(t)$ exists, where the phase $\theta(t)$ and the amplitude $\boldsymbol{R}(t)$ are Itô processes and satisfy

[TABLE]

where (explicit dependence on $\theta$ and $t$ is omitted for simplicity)

[TABLE]

Proof: See [20, 21].

The phase and amplitude equations given in theorem 1 are exact, because no approximation is used in their derivation, and their validity is not limited to weak noise. However, the resulting phase and amplitude depend on the choice of the basis $\{\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{n}\}$. The phase of a nonlinear oscillator can be unambiguously defined using the concept of isochrons [14, 23, 19]. Isochrons are $(n-1)$–dimensional manifolds transverse to the limit cycle. Given a reference point $\boldsymbol{x}_{s}(0)$ on the limit cycle, the isochron transverse to the cycle at $\boldsymbol{x}_{s}(0)$ is

[TABLE]

That is, $I_{\boldsymbol{x}_{s}(0)}$ is the set of initial conditions $\boldsymbol{x}(0)$ such that the trajectories leaving from $\boldsymbol{x}(0)$ meet asymptotically on the limit cycle at $\boldsymbol{x}_{s}(t)$. The phase of a nonlinear oscillator is introduced by assigning the same phase to all points on the same isochron, i.e. isochrons are the level sets of the phase function. The following theorem, given in [20, 21], shows that if the Floquet’s basis is used to derive eqs. (4), the resulting phase locally coincides in a neighborhood of the limit cycle, with the phase defined through isochrons. It was also shown that the resulting phase dynamics can be partially decoupled from the amplitude deviation dynamics.

Theorem 2

Consider the variational equation for the noiseless ($\varepsilon=0$) oscillator (1)

[TABLE]

where $\boldsymbol{A}(t)=\frac{\partial\boldsymbol{a}(\boldsymbol{x}_{s}(t))}{\partial\boldsymbol{x}}$ is the Jacobian matrix evaluated on $\boldsymbol{x}_{s}(t)$. Let $\boldsymbol{\Phi}(t)=\boldsymbol{P}(t)e^{\boldsymbol{D}\,t}\boldsymbol{S}_{0}$ be the fundamental matrix solution of (7), where $\boldsymbol{P}(t)$ is a $T$–periodic matrix, $\boldsymbol{S}_{0}=\boldsymbol{P}^{-1}(0)$, and $\boldsymbol{D}=\textrm{diag}[\nu_{1},\ldots,\nu_{n}]$ is a diagonal matrix whose diagonal entries are the Floquet’s characteristic exponents. If the basis vectors $\{\boldsymbol{u}_{1}(t),\ldots,\boldsymbol{u}_{n}(t)\}$ are chosen such that

[TABLE]

then, the Itô processes for the phase and amplitude reduce to

[TABLE]

where $\widetilde{\boldsymbol{D}}=\textrm{diag}[\nu_{2},\ldots,\nu_{n}]$, and the Taylor series of $\widetilde{a}_{\theta}(\theta,\boldsymbol{R})$, $\widetilde{\boldsymbol{a}}_{\boldsymbol{R}}(\theta,\boldsymbol{R})$, do not contain linear terms in $\boldsymbol{R}$.

Proof: See [20].

Theorem 2 shows that, if Floquet basis is used, then up to the first order, the phase equation (8a) is independent on the amplitude. In this sense Floquet basis represents a privileged basis. The relationship between the phase $\theta$ and the phase defined through isochrons can be inferred from (3). It is well known that at any time instant the Floquet vectors $\boldsymbol{u}_{2}(t),\ldots,\boldsymbol{u}_{n}(t)$ span the hyperplane tangent to the isochron [6, 14, 19]. For small deviations, in the neighborhood of the limit cycle, isochrons can be approximated by the hyperplane spanned by $\boldsymbol{u}_{2}(t),\ldots,\boldsymbol{u}_{n}(t)$. Therefore, any point lying on the isochron can be approximated by (3).

3 Phase oscillator model of noisy oscillators

A simple phase oscillator model can be derived from (4), by substituting to the amplitude deviation the noiseless value $\boldsymbol{R}=\mathbf{0}$. Such an assumption is made under the hypothesis that amplitude fluctuations instantaneously relax to the stable limit cycle. Taking into account that $a_{\theta}(\theta,\mathbf{0})=0$, and if $\mathcal{O}(\varepsilon^{2})$ terms are neglected, the following phase oscillator model, equivalent to those presented in [2, 3]111Please note that in the present work the angular frequency is normalized by the oscillator’s free running frequency, is obtained

[TABLE]

The expect angular frequency is readily found from eq. (9). Taking the stochastic expectation on both sides of (9) and using the zero expectation property of Itô integral, it follows that $E[d\theta/dt]=1$. However, the prediction that the noise does not modify the expected angular frequency of the oscillator turns out to be incorrect. Retaining $\mathcal{O}(\varepsilon^{2})$ terms yields the more accurate model

[TABLE]

that correspond to those given in [18, 4, 19]. The resulting expected angular frequency is $E[d\theta/dt]=1+\varepsilon^{2}E[\hat{a}_{1}(\theta,\mathbf{0})]$. Eq. (10) gives better predictions than eq. (9) for moderate values of $\varepsilon$, but the effect of amplitude deviations is still completely neglected. We shall now derive an improved phase oscillator model that takes into account the effect of noise induced amplitude fluctuations on the phase, yet remaining reasonably simple.

Let $p(\theta,\boldsymbol{R},t)$ be the probability density function (PDF) for the solution of the phase–amplitude model (4). The time evolution of the PDF is described by the Fokker–Planck Equation (FPE) given by [26]

[TABLE]

where $R_{i}$ denotes the $i^{\textrm{th}}$ component of $\boldsymbol{R}$, $L_{ij}$ is the $i,j$ element of the matrix $\boldsymbol{L}$, $a_{R_{i}}$ is the $i^{\textrm{th}}$ component of $\boldsymbol{a}_{\boldsymbol{R}}$, and $\boldsymbol{B}_{R_{i}}$ is the $i^{\textrm{th}}$ row of matrix $\boldsymbol{B}_{\boldsymbol{R}}$.

We write the amplitude deviation as a small deviation from the noiseless value in the form $\boldsymbol{R}=\varepsilon\widetilde{\boldsymbol{R}}$, and we expand the coefficients of (11) in Taylor series. Taking into account that $a_{\theta}(\theta,\mathbf{0})=0$, $\boldsymbol{a}_{\boldsymbol{R}}(\theta,\mathbf{0})=\mathbf{0}$ and that $\frac{\partial}{\partial R_{i}}=\frac{\partial}{\varepsilon\partial\widetilde{R}_{i}}$, the FPE (11) can be rewritten in compact form

[TABLE]

with operators

[TABLE]

In eqs. (13), the functions $a_{\theta}$, $\boldsymbol{a}_{R}$, $\hat{a}_{\theta}$, $\hat{\boldsymbol{a}}_{\boldsymbol{R}}$, $\boldsymbol{B}_{\theta}$, $\boldsymbol{B}_{\boldsymbol{R}}$ and their derivatives are evaluated at $(\theta,\mathbf{0})$. The zeroth order equation $\frac{\partial p}{\partial t}=L_{0}p$ is the FPE corresponding to the SDEs

[TABLE]

Eq. (14a) implies $\theta=t$, while (14b) is a time dependent Ornstein–Uhlenbeck process. It stems that $\widetilde{\boldsymbol{R}}(\theta)$ is a vector with Gaussian distributed components, and the moments are obtained solving

[TABLE]

Without loss of generality we can assume null initial conditions. Therefore (15a) has the simple solution $\boldsymbol{\mu}(\theta)=\boldsymbol{0}$ for all $\theta$. The stationary distribution is

[TABLE]

where $|\boldsymbol{\sigma}|$ is the determinant of the covariance matrix and $|\boldsymbol{\sigma}_{ij}|$ is the cofactor. The stationary distribution is used to define a projection operator $P$, defined through its action over a generic function $f(\theta,\widetilde{\boldsymbol{R}},t)$

[TABLE]

Although not indicated explicitly, the integrals in (17) extend over the whole set of admissible values of $R_{i}$. The application of the projector to the PDF $p(\theta,\widetilde{\boldsymbol{R}},t)$ gives the marginal PDF for the phase variable

[TABLE]

Together with the projector $P$ we introduce the complementary operator $Q$, such that $P+Q=\boldsymbol{1}$, where $\boldsymbol{1}$ is the identity operator. Applying the standard Mori–Zwanzig projection procedure (see [4]) the reduced FPE for the marginal probability density function is obtained

[TABLE]

where $L=L_{0}+\varepsilon L_{1}+\varepsilon^{2}L_{2}$ and $K(t-s)=PLe^{(t-s)QL}QL$ is the memory kernel. The first term on the right hand side $PLPp(\theta,\widetilde{R},t)$ is called the Markovian term, since its contribution is completely determined by the initial value of the marginal PDF $\hat{p}(\theta,t)$. The second term is the noise term and summarizes the contribution of the initial distribution of the amplitude variables $Qp(\theta,\widetilde{\boldsymbol{R}},0)$. If the initial value $p(\theta,\widetilde{\boldsymbol{R}},0)$ is random, the function $Qp(\theta,\widetilde{\boldsymbol{R}},0)$ is a stochastic process. If the initial distribution is deterministic the noise term can be eliminated by a proper choice of the initial condition, because $Qp(\theta,\widetilde{\boldsymbol{R}},0)$ in the null space of $P$. The last term is the memory term, because it collects the dependence on the whole past history of the marginal PDF through the time integral. The non Markovian behavior of the reduced equation is a consequence of the attempt to describe the system evolution with a reduced number of variables [4].

The analysis of (19) is a very challenging problem, but it can be simplified with suitable approximations. In this work we adopt a short time approximation, assuming that the memory kernel decays so fast that we can drop the integral term. The analysis of the influence of the memory term is left to future works. Thus we consider the simplified FPE

[TABLE]

Using equations (13) and (16) it is straightforward to find

[TABLE]

Introducing (18) and (21) into (20) yields the reduced FPE for the marginal PDF $\hat{p}(\theta,t)$

[TABLE]

that corresponds to the reduced SDE for the phase variable

[TABLE]

With respect to the reduced phase model (10), equation (22) includes an additional term of order $\mathcal{O}(\varepsilon^{2})$. The additional frequency shift depends on the covariance of the amplitude deviation, that can be found solving (15b). Taking the stochastic expectation on both sides of (22) and using the zero expectation property of Itô SDE, the expected angular frequency of the reduced phase model is readily obtained

[TABLE]

4 Example

As an example of application, we consider a Duffing–van der Pol system with multiplicative noise, described by the state equations

[TABLE]

In the noiseless limit ($\varepsilon=0$), the analytic solution for this system cannot be expressed in terms of simple elementary functions [24]. Thus we use numerical simulations to determine the limit cycle. The results is then used to compute the vectors $\boldsymbol{u}_{1}(t)$, $\boldsymbol{u}_{2}(t)$, $\boldsymbol{v}_{1}(t)$, and $\boldsymbol{v}_{2}(t)$, for both the orthogonal and the Floquet basis. Floquet vectors are computed using the formulas given in [25].

Figures 1 and 2 show the limit cycle, and the vectors $\boldsymbol{u}_{1}(t)$, $\boldsymbol{u}_{2}(t)$ obtained with the orthogonal and the Floquet basis, respectively. Some isochrons (computed using the algorithm described in [15], pag. 490) are also shown for reference. As expected, when Floquet basis is used the vector $\boldsymbol{u}_{2}(t)$ is locally tangent to the isochron on the limit cycle. It is worth noting that for this example isochrons cannot be expressed in terms of elementary functions, therefore, an analysis based on isochrons would require their numerical computation (together with their first and second spatial derivatives, [4]) for all points of the state space, a very demanding task.

The limit cycle $\boldsymbol{x}_{s}$ and the vectors $\boldsymbol{u}_{1}$, $\boldsymbol{u}_{2}$ are used to determine the coordinate transformation (3). Figure 3 shows the value of the Jacobian for the coordinate transformation of the Duffing–van der Pol oscillator. The thick black lines represent where the determinant is null, and thus the coordinate transformation to phase–amplitude variables is singular. It is evident that large values of the amplitude deviation (and thus of the noise intensity) can be reached before the amplitude–phase model is no longer valid.

The vectors $\boldsymbol{u}_{1}$, $\boldsymbol{u}_{2}$, together with the covectors $\boldsymbol{v}_{1}$, $\boldsymbol{v}_{2}$ are used to compute the phase and amplitude equations (4). For computational purposes, in particular for the derivatives that appear in (5), it is convenient to consider a truncated Fourier series for $\boldsymbol{u}_{2}$

[TABLE]

where $\omega_{0}$ is the fundamental frequency and $N$ is a sufficient high number of harmonics. The phase and amplitude equations (4) corresponding to the Duffing–van der Pol system (24) have been integrated numerically to determine the expected angular frequency, using both Euler–Maruyama and Milstein numerical integration schemes. The expected normalized angular frequency has been determined from the solution of the phase equation using

[TABLE]

for $t_{2}\gg t_{1}$. This is justified under the hypothesis that the system is ergodic, so that ensemble averages coincide with time averages. The expressions for $\boldsymbol{x}_{s}$, $\boldsymbol{u}_{1}$, $\boldsymbol{u}_{2}$, $\boldsymbol{v}_{1}$ and $\boldsymbol{v}_{2}$ also allow to derive and solve (15b). Note that since the Duffing–van der Pol is a second order system, the variance $\sigma(\theta)$ is just a scalar. The variance is then introduced into (23) to find the expected angular frequency given by the reduced phase oscillator model. Figure 4 shows the normalized expected angular frequency versus the noise intensity for the example under investigation. For simplicity, only Floquet basis is considered. The result obtained using numerical simulations (solid red line) is compared with the phase oscillator model given by (10) (dashed line) and with the improved model (23) (solid lines). It is clear that the proposed model (23) provides a much more accurate estimate of the expected angular frequency.

5 Conclusions

In this paper, a novel phase reduced model for nonlinear oscillators subject to white Gaussian noise has been presented. Starting from a recently developed description of the noisy oscillator dynamics in terms of amplitude and phase variables, a single scalar stochastic differential equation for the phase variable is given. The resulting phase oscillator model is analogous to the celebrated Kuramoto model.

The reduction procedure leading to the phase oscillator model is based on the use of projection operator techniques and stochastic calculus. It is shown that, with the help of Floquet theory, the phase variable obtained coincides in the vicinity of a limit cycle with the asymptotic phase defined using the concept of isochrons.

The phase equation shows a quadratic dependence of the phase dynamics on the noise intensity. It is shown that the expected angular frequency depends on the covariance matrix of amplitude variable. In particular, even in the case where the expected oscillation amplitude remains unchanged, the unavoidable amplitude fluctuations do modify the expected oscillation frequency.

As an example, a Duffing–van der Pol oscillator is analyzed both theoretically and numerically. The validity of the proposed phase oscillator model is confirmed by the numerical analysis. It is also shown that the proposed model provides more accurate prediction that other models previously proposed in literature.