A mathematical framework for amplitude and phase noise analysis of coupled oscillators

TL;DR

This paper develops a rigorous mathematical framework to analyze how noise affects the amplitude and phase synchronization of coupled oscillators, enabling phase-only modeling even under strong noise conditions.

Contribution

It introduces a novel set of equations for amplitude and phase deviations in noisy coupled oscillators, valid beyond the weak noise approximation, using Floquet theory for partial decoupling.

Findings

Decoupling of amplitude and phase dynamics via Floquet theory.

Reduced phase model accurately describes oscillator behavior under noise.

Insights into noise influence on oscillator frequency.

Abstract

Synchronization of coupled oscillators is a paradigm for complexity in many areas of science and engineering. Any realistic network model should include noise effects. We present a description in terms of phase and amplitude deviation for nonlinear oscillators coupled together through noisy interactions. In particular, the coupling is assumed to be modulated by white Gaussian noise. The equations derived for the amplitude deviation and the phase are rigorous, and their validity is not limited to the weak noise limit. We show that using Floquet theory, a partial decoupling between the amplitude and the phase is obtained. The decoupling can be exploited to describe the oscillator's dynamics solely by the phase variable. We discuss to what extent the reduced model is appropriate and some implications on the role of noise on the frequency of the oscillators.

Equations

$$\frac{d\mathbf{X}(t)}{dt}=\mathbf{a}(\mathbf{X}(t),\mathbf{\xi }(t))$$

$$\frac{d\mathbf{X}(t)}{dt}=\mathbf{a}(\mathbf{X}(t),0)+\frac{\partial \mathbf{a}(\mathbf{X}(t),0)}{\partial \mathbf{\xi }}\mathbf{\xi }(t)+\dots$$

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

$$\mathbf{X}(t)=\mathbf{X}(0)+{\int }_0^t\mathbf{a}(\mathbf{X}(s))ds+\epsilon {\int }_0^t\mathbf{B}(\mathbf{X}(s))d\mathbf{W}(s)$$

$$S_n^S=\sum _{i=1}^{n}B\left(\frac{X(t_i)+X(t_{i-1})}{2}\right)\left[W(t_i)-W(t_{i-1})\right]$$

$$S_n^I=\sum _{i=1}^{n}G(t_{i-1})\left[W(t_i)-W(t_{i-1})\right]$$

$$d\mathbf{X}={\mathbf{a}}^S(\mathbf{X})\text{leavevmode}dt+\epsilon {\mathbf{B}}^S(\mathbf{X})\circ d\mathbf{W}$$

$$d\mathbf{X}={\mathbf{a}}^I(\mathbf{X})\text{leavevmode}dt+\epsilon {\mathbf{B}}^I(\mathbf{X})d\mathbf{W}$$

$$a_i^I(\mathbf{X})=a_i^S(\mathbf{X})+\frac{{\epsilon }^2}{2}\sum _{j,k}\frac{\partial B_{ij}}{\partial x_k}{B}_{jk}$$

$${\mathbf{B}}^S(\mathbf{X})={\mathbf{B}}^I(\mathbf{X})$$

$$d{\mathbf{X}}_i=\left[{\mathbf{a}}_i({\mathbf{X}}_i)+\epsilon {\mathbf{c}}_i({\mathbf{X}}_1,\dots ,{\mathbf{X}}_N)\right]dt+\epsilon {\mathbf{B}}_i({\mathbf{X}}_1,\dots ,{\mathbf{X}}_N)d{\mathbf{W}}_{i}i=1,\dots ,N$$

$$\frac{d{\mathbf{x}}_i(t)}{dt}={\mathbf{a}}_i({\mathbf{x}}_i(t))$$

$${\mathbf{u}}_{1_i}(t)=\frac{{\mathbf{a}}_i({\mathbf{x}}_{S_i}(t))}{|{\mathbf{a}}_i({\mathbf{x}}_{S_i}(t))|}$$

$${\mathbf{u}}_{{\alpha }_i}^T(t){\mathbf{v}}_{{\beta }_i}(t)={\delta }_{\alpha \beta }$$

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

$$d{\theta }_i=$$

$$d{\mathbf{R}}_i=$$

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

$${\hat{a}}_{{\theta }_i}({\theta }_1\dots {\mathbf{R}}_N)$$

$$\displaystyle+\dfrac{1}{2}\bigg{(}\dfrac{\partial\boldsymbol{a}_{i}(\boldsymbol{x}_{S_{i}})}{\partial\theta_{i}}+\dfrac{\partial^{2}\boldsymbol{Y}_{i}}{\partial\theta_{i}^{2}}\,\boldsymbol{R}_{i}\bigg{)}\boldsymbol{B}_{\theta_{i}}\boldsymbol{B}_{\theta_{i}}^{T}\bigg{]}$$

$$c_{{\theta }_i}({\theta }_i,{\mathbf{R}}_i)=$$

$${\mathbf{B}}_{{\theta }_i}({\theta }_1\dots {\mathbf{R}}_N)$$

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

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

$${\hat{\mathbf{a}}}_{{\mathbf{R}}_i}({\theta }_1\dots {\mathbf{R}}_N)$$

$${\mathbf{c}}_{{\mathbf{R}}_i}({\theta }_1\dots {\mathbf{R}}_N)$$

$${\mathbf{B}}_{{\mathbf{R}}_i}({\theta }_1\dots {\mathbf{R}}_N)$$

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

$${\Phi }_i(t)={\mathbf{P}}_i(t)e^{\mathbf{D}t}{\mathbf{P}}_i^{-1}(0),$$

$$[r_i{\mathbf{u}}_{1_i}(t),\dots ,{\mathbf{u}}_{n_i}(t)]={\mathbf{P}}_i(t),$$

$$d{\theta }_i=$$

$$d{\mathbf{R}}_i=$$

$$d\theta_{i}=\big{(}1+\varepsilon^{2}\hat{a}_{\theta_{i}}(\theta_{i})\big{)}dt+\varepsilon\boldsymbol{B}_{\theta_{i}}(\theta^{(1)},\ldots,\theta^{(N)})\,d\boldsymbol{W}_{i}$$

$$\frac{\partial u}{\partial t}=\mathcal{A}u,$$

$$\mathcal{A}f(\theta)=\big{(}1+\varepsilon^{2}\hat{a}_{\theta}(\theta)\big{)}\dfrac{\partial f(\theta)}{\partial\theta}+\dfrac{\varepsilon^{2}}{2}\,\big{(}\boldsymbol{B}_{\theta}(\theta)\boldsymbol{B}_{\theta}^{T}(\theta)\big{)}\,\dfrac{\partial^{2}f(\theta)}{\partial\theta^{2}}$$

$$E\left[{\int }_{t_0}^{t}f(\mathbf{X})d\mathbf{W}\right]=0$$

Full text

1 Introduction

Periodically driven oscillators and coupled oscillators are classical problems in nonlinear dynamics, with many relevant applications in physics, chemistry, biology and engineering [1, 2, 3]. To make the models more realistic, external inputs can be included, to represent the unavoidable random fluctuations that occur in real world systems, due to the physical properties of the oscillators or induced by the environment. Such disturbances can be modeled by stochastic forces applied to the oscillators, which are then described by stochastic differential equations [4, 5].

Corrupting noise can dramatically affect the performance of oscillators. This is of particular relevance, for instance, in the field of modern electronic devices. Phase noise in oscillators can produce distortion or complete loss of incoming information in traditional receivers, and high bit error rates in phase modulated applications. Traditionally, the action of noise on electronic oscillators has been described as purely diffusive process [7, 6]. It is commonly assumed that the effect of white noise on the spectrum of an oscillator is to produce a broadening of the oscillator’s spectrum without affecting the positions of the peaks. Recently, this assumption has been questioned by the analysis of some simple solvable models, and by the development of improved mathematical techniques [8, 9, 10]. These works have shown that the phase noise problem is best described as a convection–diffusion process, i.e. white noise may also be responsible for a shift in the oscillator’s angular frequency.

It may sound surprising that a random perturbation can produce some kind of coherent modification to the oscillator’s frequency. In fact one may expect that, as a result of their random nature, fluctuations have a null net effect and leave the oscillation frequency and amplitude unaffected. Random perturbations may produce a coherent modification to the oscillator’s frequency because of the peculiar characteristics of oscillators. First, the autonomous nature of oscillators implies that any time shifted version of a solution is still a solution. The consequence is that phase shifts are not absorbed, but rather they accumulate in time. Second, oscillators are nonlinear dynamical systems. Some directions are preferred to others, so that perturbations along some directions are amplified, while others are attenuated. The result is that coherent behavior can emerge from random excitations.

In the last few years the idea on the role of noise has also changed at a more fundamental level. For long time noise has been considered a nuisance to be reduced as much as possible. Only recently it was figured out that noise can play a constructive role in natural phenomena or engineering applications. Important examples are stochastic resonance, where a periodic signal is amplified by noise [11, 12, 13], and energy harvesting, where noise is used as a power source [14, 15]. It also been recognized that in particular situations, noise can favor the synchronization of oscillators. Synchronization is the result of two competing mechanisms. On the one hand, differences in the oscillators free running frequencies are destructive to synchronization. On the other hand, couplings between the oscillators favor the emergence of collective rhythms. Similarly, can at the same work both towards and against synchronization. The phase diffusion produced by noise obviously acts against the synchronicity, while the frequency shift produced by noise may decrease the frequency mismatch between oscillators thus enhancing the emergence of locked states [16, 17].

To study the balance between the phase diffusion and the frequency shift, one needs a mathematical model capable to capture the influence of noise on the phase of the oscillators. In this paper we present a mathematical framework to reduce a network of oscillators subject to white Gaussian noise described in terms of state variables to the equivalent amplitude and phase model. In section 2 we introduce some basic concepts about noisy oscillators and the theory of stochastic differential equations, that represents the ideal mathematical framework for the analysis of such problems. In section 3 we derive the main result of the paper, giving a rigorous description in terms of amplitude and phase variables for a network composed by nonlinear coupled oscillators of any order, subject to white Gaussian noise. We also show how the phase dynamics can be partially decoupled from the amplitude dynamics, thus suggesting the possibility to derive reduced order model analogous to the celebrated Kuramoto model [1]. In section 4 we discuss some implications of our model on the phase dynamics, with particular attention to the role of the noise on the expected frequencies of the oscillators. We discuss the physical origin of the frequency shift and why such an effect should not be neglected with respect to the phase diffusion process. In section 5 we present a simple example to show the application of the transformation to amplitude and phase variable. The example chosen admits an analytical solution to illustrate the influence of noise on a small network. Section 6 is devoted to conclusions.

2 Noisy oscillators and stochastic differential equations

Nonlinear oscillators can be conveniently described by the differential equation

[TABLE]

where $\boldsymbol{x}:\mathds{R}\mapsto\mathds{R}^{n}$ is the state of the oscillator, and $\boldsymbol{\xi}$ represents the unavoidable noise sources, both internal and external, that affect the oscillator. In most practical situation the noise level is expected to be small with respect to the oscillator state, so that we are legitimate to linearize equation (1) around the noiseless state

[TABLE]

Another common assumption is that the noise possesses some “nice” statistical properties. In particular we shall assume that any noise source can be modeled as a Gaussian white noise. This assumption is justified in a wide range of practical situations, e.g. molecular dynamics, thermal noise, shot noise and Johnson noise in electronics. In general, white noise is a good approximation to a colored noise process in the case where the typical time scales of the underlying deterministic dynamics are much smaller than the noise correlation time (quasi–white approximation). By the central limit theorem, it is reasonable to describe $\boldsymbol{\xi}(t)$ as Gaussian distributed. Equation (2) can be rewritten as the stochastic differential equation (SDE) [5, 4]

[TABLE]

where $\boldsymbol{X}:\mathds{R}\mapsto\mathds{R}^{n}$ is the state of the oscillator, $\boldsymbol{a}:\mathds{R}^{n}\mapsto\mathds{R}^{n}$ is a vector field that describes the oscillator dynamics, $\boldsymbol{B}:\mathds{R}^{n}\mapsto\mathds{R}^{n,m}$ is a modulating real valued matrix, and $\varepsilon$ is a parameter that measures the noise intensity. $\boldsymbol{W}:\mathds{R}\mapsto\mathds{R}^{m}$ is a vector of Wiener processes, also called Brownian motions, a continuous time stochastic process characterized by zero expectation value $E[\boldsymbol{W}]=0$, independent increments, and with a Gaussian distribution. The vector valued function $\boldsymbol{a}(\boldsymbol{X})$ is called the drift term, while the matrix $\boldsymbol{B}(\boldsymbol{X})$ is called diffusion the diffusion term. For matrices $\boldsymbol{B}$ with constant entries the noise is said additive, while for a state dependent matrices $\boldsymbol{B}(\boldsymbol{X})$ the noise is said multiplicative.

Stochastic processes are nowhere differentiable, consequently the SDE (3) should be interpreted as a shorthanded version of the integral equation

[TABLE]

The first integral on the right hand side is a Riemann integral, and it does not pose any particular problem. The second integral is a Riemann–Stieltjes type integral, but differently from common Riemann–Stieltjes integrals, the point at which the function is evaluated do matter. The two main interpretation schemes are Stratonovich and Itô. According to Stratonovich, the stochastic integral111Conforming to the standard notation, we use the symbol $\boldsymbol{B}(\boldsymbol{X})\circ d\boldsymbol{W}$ to denote Stratonovich integral, while we reserve notation $\boldsymbol{B}(\boldsymbol{X})d\boldsymbol{W}$ for Itô integral. $\int_{0}^{t}B(X(s))\circ dW(s)$ is defined as the mean square limit of the partial sum [5]

[TABLE]

By contrast, in Itô interpretation the stochastic integral $\int_{0}^{t}G(s)dW(s)$ is the mean square limit of [5]

[TABLE]

Both interpretations have their own pros and cons. The main advantage of Stratonovich interpretations is that traditional calculus rules apply. The drawback is that in each time interval, both the initial value $X(t_{i-1})$ and the final value $X(t_{i})$ of the stochastic process $X$ are required to solve the SDE. This feature is known as the “look in the future property” of Stratonovich integral. As a consequence the Stratonovich interpretation is not well suited for numerical integration schemes. Moreover, in Stratonovich interpretation the stochastic process and the noise increments are correlated, making the determination of stochastic expectations difficult. By contrast, Itô stochastic integral only requires the initial value of the stochastic process in each time interval. Therefore Itô interpretation is preferred in the implementation of numerical integration schemes. In Itô SDEs the stochastic process and the noise increments are uncorrelated, making the determination of stochastic expectations easier. The drawback of Itô view is that a new set of calculus rules, known as Itô calculus must be used.

The relevant consequence of the two different interpretations is that the same SDE has different solutions whether it is interpreted following Stratonovich or Itô. However, the two interpretations are linked by a transformation that converts any Stratonovich (respectively Itô) SDE into an equivalent Itô (respectively Stratonovich) SDE. By equivalent we mean a different SDE, interpreted with different rules, but that has the same solution [4, 5]. The equivalence opens the possibility to switch from one interpretation to the other to take advantage of the pros of both the definitions. The Stratonovich SDE (the apex $S$ and $I$ denote Stratonovich and Itô, respectively)

[TABLE]

is equivalent to the Itô SDE

[TABLE]

where ($a_{i}$ is the $i$-th component of $\boldsymbol{a}$, $B_{ij}$ the $(i,j)$ element of $\boldsymbol{B}$)

[TABLE]

and

[TABLE]

In the following we shall use Itô interpretation, omitting the apex $I$ for simplicity.

3 Amplitude and phase equations for a network of oscillators with noisy interactions

A network composed of $N$ weakly coupled nonlinear oscillators with noisy interactions can be described, similarly to (3), by the (Itô) SDEs

[TABLE]

where $\boldsymbol{X}_{i}:\mathds{R}\mapsto\mathds{R}^{n}$ is a stochastic process describing the state of the $i$–th oscillator, $\boldsymbol{a}_{i}:\mathds{R}^{n}\mapsto\mathds{R}^{n}$ is the $i$–th drift coefficient, $\boldsymbol{B}_{i}:\mathds{R}^{n\cdot N}\mapsto\mathds{R}^{n,m}$ is the $i$–th $n\times m$ diffusion matrix, and $\boldsymbol{W}_{i}:\mathds{R}\mapsto\mathds{R}^{m}$ is the $i$–th vector of Wiener processes. For the sake of simplicity, in equation (9) we assume that all oscillators are of the same order222Oscillators of different order do not pose any particular problem, they only make the notation more involved. The theorems 1 and 2 can be formulated, mutatis mutandis, for oscillator of different orders. ($\boldsymbol{X}_{i}\in\mathds{R}^{n}$, for all $i$), but we allow the interaction to vary for each oscillator both in the modulating matrix $\boldsymbol{B}_{i}$ and in the random fluctuation $\boldsymbol{W}_{i}$. Such a model may arise, for instance, if both $\boldsymbol{a}_{i}$ and $\boldsymbol{c}_{i}$ have some small stochastic components, and they are linearized around a noiseless state.

For $\varepsilon=0$, the SDE (9) reduce to an ordinary differential equation (ODE) describing $N$ independent, noiseless oscillators. The $i$–th oscillator is described by the ODE

[TABLE]

We assume that the ODE (9) admits an asymptotically stable $T_{i}$–periodic solution, represented by a limit cycle $\boldsymbol{x}_{S_{i}}(t)$ in its state space. For each oscillator we define the vector

[TABLE]

$\boldsymbol{u}_{1_{i}}(t)$ is the unit vector that at each time instant is tangent to the limit cycle $\boldsymbol{x}_{S_{i}}(t)$. Together with $\boldsymbol{u}_{1_{i}}(t)$ we consider other $n-1$ vectors $\boldsymbol{u}_{2_{i}}(t),\ldots,\boldsymbol{u}_{n_{i}}(t)$, such that the set $\{\boldsymbol{u}_{1_{i}}(t),\ldots,\boldsymbol{u}_{n_{i}}(t)\}$ is a basis for $\mathds{R}^{n}$ for all $t$. Let $\boldsymbol{U}_{i}(t)=[\boldsymbol{u}_{1_{i}}(t),\ldots,\boldsymbol{u}_{n_{i}}(t)]$ be the matrix whose columns are $\boldsymbol{u}_{1_{i}}(t),\ldots,\boldsymbol{u}_{n_{i}}(t)$. Such a matrix is obviously invertible, and let $\boldsymbol{V}_{i}(t)=\boldsymbol{U}_{i}^{-1}(t)$ be the inverse. We define the reciprocal vectors $\boldsymbol{v}_{1_{i}}^{T}(t),\ldots,\boldsymbol{v}_{n_{i}}^{T}(t)$ to be the rows of $\boldsymbol{V}_{i}(t)$. By construction, the $\boldsymbol{u}_{i}$ and $\boldsymbol{v}_{i}$ vectors are bi–orthogonal, i.e.

[TABLE]

where $\delta_{\alpha\beta}$ is the Kronecker’s symbol. We shall also use the matrices $\boldsymbol{Y}_{i}(t)=[\boldsymbol{u}_{2_{i}}(t),\ldots,\boldsymbol{u}_{n_{i}}(t)]$, $\boldsymbol{Z}_{i}(t)=[\boldsymbol{v}_{2_{i}}(t),\ldots,\boldsymbol{v}_{n_{i}}(t)]$, and the modulus of the vector field $r_{i}(t)=|\boldsymbol{a}_{i}(\boldsymbol{x}_{S_{i}}(t))|$ evaluated over the limit cycle.

A crucial concept to be defined in the analysis of synchronization of oscillators is the phase concept. A phase function is intended to represent the projection of the oscillator’s state onto a reference trajectory, normally the unperturbed limit cycle. For each oscillator we introduce a phase function $\theta_{i}:\mathds{R}^{n}\mapsto[0,T_{i})$, interpreted as an elapsed time from an initial reference point. Consider a point $\boldsymbol{x}_{S_{i}}(0)$ on the limit cycle, and assign phase zero to this point, i.e. $\theta_{i}(\boldsymbol{x}_{S_{i}}(0))=0$. The phase of the point $\boldsymbol{x}_{S_{i}}(t)$ is $\theta_{i}(\boldsymbol{x}_{S_{i}}(t))=t,\mod T_{i}$. Thus, the phase represents a new parametrization of the limit cycle. Together with the phase function we shall consider an amplitude deviation function $\boldsymbol{R}_{i}:\mathds{R}^{n}\mapsto\mathds{R}^{n-1}$, with $\theta_{i},\boldsymbol{R}_{i}\in\mathcal{C}^{m}(R^{n})$, $m\geq 2$. The amplitude333We shall use the term “amplitude” instead of the more correct “amplitude deviation” for the sake of simplicity. function $\boldsymbol{R}_{i}$ is interpreted as an orbital deviation from the limit cycle $\boldsymbol{x}_{S_{i}}(t)$, see figure 1.

The following theorem establishes the amplitude and phase equation for the network

Theorem 1**.**

Consider the Itô diffusion (9), and the reciprocal bases $\{\boldsymbol{u}_{1_{i}}(t),\ldots,\boldsymbol{u}_{n_{i}}(t)\}$ and $\{\boldsymbol{v}_{1_{i}}(t),\ldots,\boldsymbol{v}_{n_{i}}(t)\}$, satisfying the bi–orthogonality condition $\boldsymbol{v}_{\alpha_{i}}^{T}\boldsymbol{u}_{\beta_{i}}=\delta_{\alpha\beta}$, for all $i=1,\ldots,N$. Consider the coordinate transformation

[TABLE]

Then in a neighborhood of the limit cycle $\boldsymbol{x}_{S_{i}}$ the phase $\theta_{i}(t)$ and the amplitude $\boldsymbol{R}_{i}(t)$ are Itô processes and satisfy

[TABLE]

where $(\theta_{1}\ldots\boldsymbol{R}_{N})$ is a shorthanded notation for $(\theta_{1},\boldsymbol{R}_{1},\ldots,\theta_{N},\boldsymbol{R}_{N})$ and (we omit explicit dependence on $\theta_{i}$ and $t$ for simplicity)

[TABLE]

Proof: See appendix A.

The amplitude and phase equations (13a) and (13b) are exact, since no approximation is involved in their derivation, and they are valid for any value of the noise intensity $\varepsilon$ as long as the Jacobian matrices $D\boldsymbol{h}_{i}$ are regular. The amplitude and phase equations obtained crucially depend on the choice of the basis vectors $\boldsymbol{u}_{2_{i}},\ldots,\boldsymbol{u}_{n_{i}}$.

In general, the equations for the two Itô processes for the phase and for the amplitude are coupled together. It is possible to show that, making use of Floquet theory, a partial decoupling between the phase and the amplitude dynamics is obtained. Before introducing the theorem we recall the main results of the Floquet theory [6, 18]. Let $\boldsymbol{A}_{i}(t)=\frac{\partial\boldsymbol{a}_{i}(\boldsymbol{x}_{S_{i}})}{\partial\boldsymbol{x}_{i}}$ be the Jacobian matrix of the $i$–th oscillator evaluated on the limit cycle $\boldsymbol{x}_{S_{i}}(t)$, and let $\Phi_{i}(t)$ be the fundamental matrix of the variational equation

[TABLE]

Thus, from Floquet theory we get:

[TABLE]

where $\boldsymbol{P}_{i}(t)$ is a $T_{i}$–periodic matrix, and $\boldsymbol{D}_{i}=\text{diag}[\nu_{1_{i}},\dots,\nu_{n_{i}}]$ is a diagonal matrix whose diagonal entries are the Floquet characteristic exponents [6, 18].

Theorem 2**.**

If the vectors $\boldsymbol{u}_{2_{i}}(t),\ldots,\boldsymbol{u}_{n_{i}}(t)$ are chosen such that

[TABLE]

then the Itô processes (13a) and (13b) become

[TABLE]

where $\widetilde{\boldsymbol{D}}_{i}=\text{diag}[\nu_{2_{i}},\dots,\nu_{n_{i}}]$ and the Taylor series of $\tilde{a}_{\theta_{i}}(\theta_{i},\boldsymbol{R}_{i})$ and $\tilde{\boldsymbol{a}}_{\boldsymbol{R}_{i}}(\theta_{i},\boldsymbol{R}_{i})$ do not contain linear terms in $\boldsymbol{R}_{i}$.

Proof: See appendix B.

The asymptotic stability hypothesis of the limit cycle implies that the Floquet characteristic exponents $\nu_{2_{i}},\ldots,\nu_{n_{i}}$ have negative real parts, for all $i=1,\ldots,N$. As a consequence, in the limit $\varepsilon\ll 1$, the amplitude asymptotic dynamics is one order of magnitude slower than the phase dynamics. This observation suggests the idea to neglect the amplitude dynamics given in (41b), and to approximate the stochastic processes $\boldsymbol{R}_{i}$ in (41a) with some reasonable (possibly constant) estimate. This approach leads to the so called “phase reduced models”. For instance, assuming the unperturbed value $\boldsymbol{R}_{i}\approx 0$, and taking into account that $a_{\theta_{i}}(\theta_{i},0)=0$ we obtain the simple phase equation

[TABLE]

We remark that the assumptions leading to (15) is justified only if the amplitudes relax instantaneously to the unperturbed value, an assumption that is often made more for mathematical convenience than being physically plausible. In fact, this assumption relies on linear approximation of manifolds, and nonlinear effects will become stronger the further we move away from the limit cycle. Moreover, the presence of nearby invariant structures such as equilibrium points and invariant manifolds may result in trajectories spending long periods of time away from the limit cycle, thus nullifying the instantaneous relaxation hypothesis. As a consequence a better solution is to chose the approximation $\boldsymbol{R}_{i}\approx E[\boldsymbol{R}_{i}]$, provided that the expected amplitude can be computed [10].

4 Discussion

A full analysis of the amplitude and phase equations (13a), (13b) is a formidable problem and will not be treated here. We limit ourselves to some considerations and some general comments.

Equations (13a), (13b), or the phase reduced model (17) suggest that the drift effects due to $\varepsilon^{2}\hat{a}_{\theta_{i}}$ become negligible in the limit of vanishing small noise ($\varepsilon\rightarrow 0$). The phase model in this limit has been extensively studied both at the single oscillator and network level [7, 6, 19]. However, for small but finite values of $\varepsilon$, the drift effects may become significant if $\hat{a}_{\theta_{i}}$ become large enough. That the drift effects should not be neglected even for small values of $\varepsilon$ can be justified as follows. When dealing with stochastic processes, the single realization is not very much significant: it is often much more useful to look at expected quantities. To illustrate the point, it is sufficient to consider a single oscillator described by a phase reduced model, the full network model (13a), (13b) is conceptually analogous. Let $f(\theta)$ be an arbitrary scalar function of the phase, and let $u(t,\theta)=E[f(\theta)]$ be the expected value of this function, with initial value $u(0,\theta)=f(\theta)$. Then the time evolution of $u(t,\theta)$ is governed by the Kolmogorov Backward Equation [5]

[TABLE]

where $\mathcal{A}$ is the generator of the Itô diffusion

[TABLE]

Equation (19) shows that both the $\mathcal{O}(\varepsilon^{2})$ drift coefficient and the $\mathcal{O}(\varepsilon)$ diffusion coefficient in (17) contribute for $\mathcal{O}(\varepsilon^{2})$ terms to the evolution of expected quantities in (18), (19), and therefore we are not allowed to neglect one contribution with respect to the other.

Expected quantities can be determined using Itô calculus without solving the Kolmogorov Backward Equation (19). In fact, let $\boldsymbol{X}$ be the solution of an Itô SDE, and let $f$ be a non anticipating function (adapted process), then the zero expectation property of Itô stochastic integral holds

[TABLE]

Taking the stochastic expectation on both sides of equations (13a), (13b) and using the zero expectation property we can transform the SDEs for the amplitude and phase into a set of ODEs for the expectation values

[TABLE]

where the property $dE[\theta]/dt=E[d\theta/dt]$ has been used. The problem here is the nonlinear nature of the ODEs. In fact, to compute the expectation of the nonlinear functions one needs all the moments for the amplitudes and the phases. We illustrate the issue for the function $a_{\theta_{i}}$ for the simple case of a scalar amplitude $R_{i}$. Taking the Taylor series in the neighborhood of $\theta_{i}=0$, $R_{i}=0$ we have

[TABLE]

Therefore, to compute $E[\theta_{i}]$ and $E[R_{i}]$, one needs all the moments of $\theta_{i}$, $R_{i}$, i.e. the system is open. To close the system various approaches are available [20, 21]. Among others, moment closure techniques are procedures to approximate the exact (but open) moment dynamics with a closed (but approximate) system. A relatively simple closure technique amounts to assume that higher order moments can be expressed in terms of the lowest order ones, assuming that the stochastic processes satisfy certain distribution laws (e.g. Gaussian distribution). Equation (22) is also instrumental to show the limit of the phase reduction method. Although the amplitude is expected to remain small (E\big{[}R_{i}\big{]}\approx 0), higher order moments (e.g. the variance) may play a relevant role to modify the expected angular frequency. Obviously, the same consideration holds true for the amplitude variable.

We close this section with a final remark. One may argue that the drift terms $\hat{a}_{\theta_{i}}$ and $\hat{\boldsymbol{a}}_{\boldsymbol{R}_{i}}$ are artifacts due to Itô interpretation. However, it turns out that the frequency drift is also present if Stratonovich interpretation is used [8, 9, 10]. To clarify the point, consider the Stratonovich SDE describing the network of oscillators

[TABLE]

Taking into account that in Stratonovich interpretation traditional calculus rules apply, repeating the procedure used in the previous section the following reduced phase model is derived

[TABLE]

However, in Stratonovich interpretation is no longer true that the stochastic processes and the noise increments are uncorrelated. Because of the anticipating nature discussed in section 2, $E[\int_{t_{0}}^{t}f(\boldsymbol{X})\circ d\boldsymbol{W}]\neq 0$. To resolve the correlation, a Stratonovich SDE has to be transformed into its equivalent Itô SDE by the addition of the drift correction term [5, 4]. Here is where the drift coefficient, that arises naturally from the quadratic terms in Itô formula, comes into play [8, 9, 10].

5 Example

In this section we give an example to show the derivation of the phase and amplitude deviation equations starting from the network’s state equations. To keep everything reasonably simple, we consider a network composed by two oscillators written in polar coordinates

[TABLE]

for $i,j=1,2$, and $j\neq i$. The real parameters $\nu_{i}$ define the free running frequencies of the oscillators in absence of noise. Although most of the information concerning expectation values and phase locking can be directly obtained from equations (25a)–(25b), we first transform these equations into the equivalent amplitude and phase models using theorems 1 and 2 to show the application of the method.

For $\varepsilon=0$, equation (25a) admits two stationary solutions: $\rho_{i}=0$ that corresponds to an unstable equilibrium point, and $\rho_{i}=1$, that corresponds to an asymptotically stable limit cycle with angular frequency $d\phi_{i}/dt=\nu_{i}$. Without loss of generality, we shall assume $\nu_{1}>\nu_{2}$. To investigate the synchronization of the two oscillators we look at the phase difference $\psi=\phi_{1}-\phi_{2}$. The oscillators are phase locked if the phase difference remains constant in time. If the two oscillators have different free running frequencies $\nu_{1}\neq\nu_{2}$, the phase difference $\psi=(\nu_{1}-\nu_{2})t$ grows unboundedly large (in absolute value) as the time passes. Conversely, if the coupling effect is taken into account but the noise influence is ignored, the phase difference evolves according to

[TABLE]

Asymptotically the two oscillators become phase locked with phase difference

[TABLE]

Moreover, since $d\psi/dt>0$ for $\psi<\psi_{s}$ and $d\psi/dt<0$ for $\psi>\psi_{s}$, the phase locked state is asymptotically stable.

5.1 Amplitude and phase equations using an orthogonal frame

Each uncoupled oscillator of the network (25a)–(25b) admits the limit cycle

[TABLE]

It follows that the orthogonal basis is given by the tangent vector $\boldsymbol{u}_{1_{i}}(t)=[0,1]^{T}$ and the orthogonal unit vector $\boldsymbol{u}_{2_{i}}(t)=[1,0]^{T}$. Since the matrix $\boldsymbol{U}_{i}(t)=[\boldsymbol{u}_{1_{i}}(t),\boldsymbol{u}_{2_{i}}(t)]$ is orthogonal, we have $\boldsymbol{v}_{1_{i}}(t)=\boldsymbol{u}_{1_{i}}(t)$ and $\boldsymbol{v}_{2_{i}}(t)=\boldsymbol{u}_{2_{i}}(t)$. The change of coordinates (12) implies $\rho_{i}=1+R_{i}$ and $\phi_{i}=\nu_{i}\theta_{i}$. Since the Jacobian of the transformation is |D\boldsymbol{h}_{i}(\theta_{i},R_{i})\big{|}_{R_{i}=0}=\nu_{i}, the coordinate transformation holds for any value of the noise intensity. Using theorem 1 it is straightforward to derive the amplitude and phase equations

[TABLE]

As it was expected from theorem 2, in the phase equation the drift coefficient contains a linear term in $R_{i}$.

5.2 Amplitude and phase equations using a Floquet frame

The Jacobian matrix of (25a), (25b) for $\varepsilon=0$ evaluated over the limit cycle for $\varepsilon=0$ is

[TABLE]

with eigenvalues $\lambda_{1_{i}}=0$, $\lambda_{2_{i}}=-1$, for all $i$. The corresponding eigenvectors are the Floquet vectors $\boldsymbol{u}_{1_{i}}(t)=[0,1]^{T}$ and $\boldsymbol{u}_{2_{i}}(t)=[1,-\nu_{i}]^{T}$. Inverting the matrix $\boldsymbol{U}_{i}(t)=[\boldsymbol{u}_{1_{i}}(t),\boldsymbol{u}_{2_{i}}(t)]$ we find the Floquet co–vectors $\boldsymbol{v}_{1_{i}}(t)=[\nu_{i},1]^{T}$ and $\boldsymbol{v}_{2_{i}}(t)=[1,0]^{T}$. The relation between the old and the new coordinates is $\rho_{i}=1+R_{i}$ and $\phi_{i}=\nu_{i}(\theta_{i}-R_{i})$. As before the Jacobian matrix $D\boldsymbol{h}_{i}$ is regular on the whole plane $\theta_{i},R_{i}$. Consequently the phase and amplitude equations in the new basis are

[TABLE]

Comparing (28a) to (30a) we observe that according to theorem 2, the latter has a drift coefficient that starts with a quadratic term in $R_{i}$. The use of Floquet basis also emphasizes the role played by higher order moments, e.g. the variance, on the angular frequencies of the oscillators.

In this particular example $E[R_{i}]$ and $E[R_{i}^{2}]$ can be determined analytically, because equation (30b) depends on $R_{i}$ only. Thus we can write a one dimensional Fokker–Planck equation (see [4]) for the probability density function (PDF) of the amplitude variable

[TABLE]

In the limit $t\rightarrow+\infty$ it admits the stationary solution

[TABLE]

where $\mathcal{N}$ is normalization constant that can be determined through the requirement $\int_{-\infty}^{+\infty}\overline{p}_{i}(R_{i})\,dR_{i}=1$. Using the stationary PDF we can compute the expectation value for an arbitrary function of the amplitude through

[TABLE]

Figure 2 shows the stationary PDF $\overline{p}_{i}(R_{i})$ for different values of the parameter $\varepsilon$. For small noise intensity the PDF is well approximated by a Gaussian distribution. Increasing the noise intensity, we observe an increase in the variance (due to diffusion) and a shift in the mode (the most probable value of the amplitude). We also observe the PDF becomes asymmetric with respect to the maximum value, thus indicating that higher order moments become more and more relevant. Figure 3 shows the first two moments for the amplitude deviation. It is seen the quadratic dependence on the noise intensity as predicted by theorem 1.

Taking stochastic expectations of (30a) and using the zero expectation property of Itô integral we have

[TABLE]

Multiplying444We recall that $\theta_{i}$ is a normalized phase variable, multiplication for $\nu_{i}$ is necessary to retrieve the non normalized phase variable. by $\nu_{i}$ it is straightforward to obtain the equation for the expected phase difference

[TABLE]

Since the amplitude equation is the same for all the oscillators555The case where the noise intensity $\varepsilon$ is not equal for all the oscillators can be treated similarly. Obviously $E[R_{i}]\neq E[R_{j}]$ and $E[R_{i}^{2}]\neq E[R_{j}^{2}]$ and the solution of (35) would be more complicated., we have $E[R_{1}]=E[R_{2}]$ and $E[R_{1}^{2}]=E[R_{2}^{2}]$, then

[TABLE]

Asymptotically the oscillators converge to the phase locked state

[TABLE]

which is different from the phase locked state in absence of noise (26). The phase difference in presence of noise is compared with that obtained without noise in figure 4. On the left we can see the asymptotic expected phase difference versus the noise intensity, while on the right it is shown the phase difference versus time for a specific realization of the noise process. It can be seen how noise operates to actively reduce the phase difference between the oscillators.

6 Conclusions

We have considered networks of coupled nonlinear oscillators subject to white Gaussian noise. We have shown that the network can be conveniently described by stochastic differential equations. The advantages and disadvantages of the two most popular interpretations, i.e. Itô and Stratonovich, have been briefly outlined.

Using projection techniques and Itô calculus, we have derived a rigorous mathematical description for the network dynamics, in terms of the phases and amplitude deviations of the oscillators. We have shown that using Floquet theory a partial decoupling between the phase and the amplitude dynamics can be obtained. This idea leads to the development of phase reduced models analogous to the celebrated Kuramoto model.

The amplitude and phase description highlights the influence of noise on the phases of the oscillators. It represents a good starting point for the analysis of the role of noise on synchronization. It is shown that noise can prevent phase locking of oscillators, by producing phase diffusion. It also shown that noise can favor synchronization. In fact the oscillators may adjust their frequency in response to noise intensity, and as a consequence noise can actively contribute to the synchronization by decreasing the frequency mismatch. A simple example has been used to illustrate the derivation of the amplitude and phase equations and their analysis.

Acknowledgement

This work was supported by the Università Italo Francese in the framework Galileo 2014–2015 (project number G14-145).

Appendix A: Proof of theorem 1

First we show that a neighborhood of the limit cycle exists, where $\theta_{i}$ and $\boldsymbol{R}_{i}$ are Itô processes. The Jacobian matrix of the coordinate transformation (12) evaluated on the limit cycle is

[TABLE]

Since $\{\boldsymbol{u}_{1_{i}}(t),\ldots,\boldsymbol{u}_{n_{i}}(t)\}$ is a basis for $\mathds{R}^{n}$, it follows that \det D\boldsymbol{h}_{i}(\theta_{i},\boldsymbol{R}_{i})\big{|}_{\boldsymbol{R}_{i}=0}\neq 0. Then by the inverse function theorem a neighborhood of $\boldsymbol{R}_{i}=0$ exists, where $\boldsymbol{h}_{i}$ is invertible. Moreover, if $\boldsymbol{h}_{i}$ is of class $\mathcal{C}^{k}$ then its inverse is also of class $\mathcal{C}^{k}$. Taking the inverse of $\boldsymbol{h}_{i}$ we can write $\theta_{i}=\theta_{i}(\boldsymbol{x}_{i})$ and $\boldsymbol{R}_{i}=\boldsymbol{R}(\boldsymbol{x}_{i})$, and if the basis vectors are smooth enough it follows from Itô formula that $\theta_{i}$ and $\boldsymbol{R}_{i}$ are Itô processes.

Now we prove that $\theta_{i}$ and $\boldsymbol{R}_{i}$ satisfy equations (13a) and (13b). Using Itô formula and eq. (9), $\boldsymbol{x}_{i}=\boldsymbol{h}_{i}(\theta_{i},\boldsymbol{R}_{i})$ implies

[TABLE]

where $\partial\boldsymbol{h}_{i}/\partial\boldsymbol{R}_{i}$ and $\partial^{2}\boldsymbol{h}_{i}/\partial\boldsymbol{R}_{i}^{2}$ are the Jacobian and the Hessian matrices, respectively. Introducing (12) yields

[TABLE]

Multiplying (38) to the left by $\boldsymbol{v}^{T}_{1_{i}}$ and using the bi–orthogonality condition we get

[TABLE]

Multiplying (38) to the left by $\boldsymbol{Z}_{i}^{T}$ gives

[TABLE]

Since $\theta_{i}$ and $\boldsymbol{R}_{i}$ are Itô processes they satisfy relations of type $d\theta_{i}=\alpha_{i}\,dt+\boldsymbol{\beta_{i}}\,d\boldsymbol{W}_{i}$, and $d\boldsymbol{R}_{i}=\boldsymbol{\gamma}_{i}\,dt+\boldsymbol{\sigma}_{i}\,d\boldsymbol{W}_{i}$, respectively. By Itô lemma $(d\theta_{i})^{2}=\boldsymbol{\beta}_{i}\boldsymbol{\beta}_{i}^{T}\,dt$, and $d\theta_{i}\,d\boldsymbol{R}_{i}=\boldsymbol{\sigma}_{i}\boldsymbol{\beta}_{i}^{T}\,dt$. Introducing these results in (39), (40) and equating terms in $d\boldsymbol{W}_{i}$ we obtain

[TABLE]

Finally, using (41a), (41b) together with $(d\theta_{i})^{2}=\boldsymbol{\beta}_{i}\boldsymbol{\beta}_{i}^{T}\,dt$, $d\theta_{i}\,d\boldsymbol{R}_{i}=\boldsymbol{\sigma}_{i}\boldsymbol{\beta}_{i}^{T}\,dt$ in (39), $\eqref{eq17}$, and rearranging the terms we get the thesis. $\square$

Appendix B: Proof of theorem 2

First of all, we observe that by hypothesis the columns of matrix $\boldsymbol{P}_{i}(t)$ are linearly independent for any $t$, and therefore they can be chosen as a basis for $\mathds{R}^{n}$. Moreover, from equation (11), we have $r_{i}\boldsymbol{u}_{1_{i}}(t)=d\boldsymbol{x}_{S_{i}}/dt$. Therefore $\nu_{1_{i}}=0$, since $d\boldsymbol{x}_{S_{i}}/dt$ is the solution of the variational equation associated to the structural Floquet exponent. Furthermore, equation (15) implies $\boldsymbol{P}_{i}(t)=\Phi_{i}(t)\boldsymbol{P}_{i}(0)e^{-\boldsymbol{D}_{i}t}.$ Taking the derivatives of (15) yields

[TABLE]

and taking into account that $\Phi_{i}(t)$ is a fundamental matrix of the variational equation we obtain: $d\boldsymbol{P}_{i}/dt=\boldsymbol{A}_{i}\boldsymbol{P}_{i}-\boldsymbol{P}_{i}\boldsymbol{D}_{i}$ and consequently $d\boldsymbol{Y}_{i}/dt=\boldsymbol{A}_{i}\boldsymbol{Y}_{i}-\boldsymbol{Y}_{i}\widetilde{\boldsymbol{D}}_{i}.$ Substituting the expression of $d\boldsymbol{Y}_{i}/dt$ in (14a), (14e) and (14f), taking the Taylor series $\boldsymbol{a}_{i}(\boldsymbol{x}_{S_{i}}+\boldsymbol{Y}_{i}\boldsymbol{R}_{i})=\boldsymbol{a}_{i}(\boldsymbol{x}_{S_{i}})+\boldsymbol{A}_{i}\boldsymbol{Y}_{i}\boldsymbol{R}_{i}+\dots$, and using the bi-orthogonality condition the thesis follows. $\square$