Phase descriptions of a multidimensional Ornstein-Uhlenbeck process

TL;DR

This paper analytically investigates two phase definitions for a two-dimensional Ornstein-Uhlenbeck process, revealing their geometric structure, dependence on system parameters, and differences in temporal progression.

Contribution

It provides explicit formulas and geometric insights for asymptotic phases in a tractable stochastic oscillator model, clarifying their properties and differences.

Findings

Isochrons form spokes of a wheel shape.

Backward phase is determined by the deterministic vector field.

Forward phase depends on both the deterministic field and noise matrix.

Abstract

Stochastic oscillators play a prominent role in different fields of science. Their simplified description in terms of a phase has been advocated by different authors using distinct phase definitions in the stochastic case. One notion of phase that we put forward previously, the \emph{asymptotic phase of a stochastic oscillator}, is based on the eigenfunction expansion of its probability density. More specifically, it is given by the complex argument of the eigenfunction of the backward operator corresponding to the least negative eigenvalue. Formally, besides the `backward' phase, one can also define the `forward' phase as the complex argument of the eigenfunction of the forward Kolomogorov operator corresponding to the least negative eigenvalue. Until now, the intuition about these phase descriptions has been limited. Here we study these definitions for a process that is analytically tractable, the two-dimensional Ornstein-Uhlenbeck process with complex eigenvalues. For this process, (i) we give explicit expressions for the two phases; (ii) we demonstrate that the isochrons are always the spokes of a wheel, but that (iii) the spacing of these isochrons (their angular density) is different for backward and forward phases; (iv) we show that the isochrons of the backward phase are completely determined by the deterministic part of the vector field, whereas the forward phase also depends on the noise matrix; and (v) we demonstrate that the mean progression of the backward phase in time is always uniform, whereas this is not true for the forward phase except in the rotationally symmetric case. We illustrate our analytical results for a number of qualitatively different cases.

Tables

Table 1. Table 1: Notation for different phase notions used in the paper. Limit cycle period is T 𝑇 T .

Symbol	Meaning
$\tau \in [0,T)$	Asymptotic phase of deterministic LC (timelike).
$\theta \in [0,2\pi )$	Asymptotic phase of deterministic LC (circular).
$\vartheta \in [0,2\pi )$	Geometric phase (standard polar coordinates).
$\Psi \in [0,2\pi )$	Backward phase for stochastic oscillator.
$\Phi \in [0,2\pi )$	Forward phase for stochastic oscillator.

Equations

$$\frac{d{x}_1}{dt}$$

$$\frac{d{x}_2}{dt}$$

$$P(\mathbf{x},t|{\mathbf{x}}_0,t_0)=\frac{1}{|d\mathbf{x}|}\mathrm{Pr}\{\mathbf{X}(t)\in [\mathbf{x},\mathbf{x}+d\mathbf{x})|\mathbf{X}(t_0)={\mathbf{x}}_0\}$$

$$\frac{\partial }{\partial t}P(\mathbf{x},t|{\mathbf{x}}_0,t_0)={\mathcal{L}}_{\mathbf{x}}[P(\mathbf{x},t|{\mathbf{x}}_0,t_0)]$$

$$-\frac{\partial }{\partial t_0}P(\mathbf{x},t|{\mathbf{x}}_0,t_0)={\mathcal{L}}_{{\mathbf{x}}_0}^{†}[P(\mathbf{x},t|{\mathbf{x}}_0,t_0)].$$

$$\mathcal{L}[P_{\lambda }]$$

$${\mathcal{L}}^{†}[Q_{\lambda }^{\ast }]$$

$$⟨Q_{\lambda }|P_{{\lambda }^{\prime }}⟩$$

$$P(\mathbf{x},t|{\mathbf{x}}_0,t_0)=P_0(\mathbf{x})+\sum _{\lambda }{e}^{\lambda (t-t_0)}{P}_{\lambda }(\mathbf{x})Q_{\lambda }^{\ast }({\mathbf{x}}_0).$$

$$P_1(\mathbf{x})=v(\mathbf{x})e^{-i\Phi (\mathbf{x})},Q_1^{\ast }({\mathbf{x}}_0)=u({\mathbf{x}}_0)e^{i\Psi ({\mathbf{x}}_0)},$$

$$\Psi (\mathbf{x})=\arg \left[Q_1^{\ast }(\mathbf{x})\right],\Phi (\mathbf{x})=-\arg \left[P_1(\mathbf{x})\right],$$

$$\frac{P(\mathbf{x},t|{\mathbf{x}}_0,t_0)-P_0(\mathbf{x})}{2u({\mathbf{x}}_0)v(\mathbf{x})}\simeq$$

$$e^{\mu (t-t_0)}\cos \left[\omega (t-t_0)+\Psi ({\mathbf{x}}_0)-\Phi (\mathbf{x})\right],\text{ as }t\to \infty .$$

$$\Psi (t)\equiv \Psi (\mathbf{X}(t)),\Phi (t)\equiv \Phi (\mathbf{X}(t)),$$

$$\frac{d\mathbf{x}}{dt}=\mathbf{F}(\mathbf{x}),\mathbf{x}(0)={\mathbf{x}}_0.$$

$$\frac{d\tau }{dt}=1,$$

$$\underset{t\to \infty }{\lim }\left|\mathbf{x}(t)-\gamma (t+\tau ({\mathbf{x}}_0))\right|=0,\text{ as }t\to \infty ,$$

$$\frac{d\mathcal{T}(\mathbf{x}(t))}{dt}=\nabla \mathcal{T}(\mathbf{x})\cdot \mathbf{F}(\mathbf{x})=1,$$

$$\frac{d\mathbf{x}}{dt}=\left(\begin{array}[]{cc}\mu&-\omega\\ \omega&\mu\end{array}\right)\mathbf{x},\quad\text{or,}\quad\dot{r}=\mu r,\quad\dot{\vartheta}=\omega$$

$$\mathcal{T}(r,\vartheta )=k\vartheta +\frac{1-k\omega }{|\mu |}\ln r,$$

$$\dot{\mathbf{x}}=A\mathbf{x}+B\xi ,⟨{\xi }_i(t){\xi }_j(t^{\prime })⟩={\delta }_{i,j}\delta (t-t^{\prime }),i,j=1,2.$$

$$A$$

$${\lambda }_{\pm }$$

$$A{\mathbf{v}}_{\pm }={\lambda }_{\pm }{\mathbf{v}}_{\pm },{\mathbf{w}}_{\pm }^{\ast }A={\lambda }_{\pm }{\mathbf{w}}_{\pm }^{\ast },$$

$$\mathbf{v}_{\pm}=Ne^{i\eta}\left(\begin{array}[]{c}\alpha_{\omega}-\omega_{0}\\ \pm i\omega-\alpha_{\mu}\end{array}\right),\;\;\mathbf{w}_{\pm}=\left(\begin{array}[]{c}\alpha_{\omega}+\omega_{0}\\ \mp i\omega-\alpha_{\mu}\end{array}\right),$$

$$D$$

$$=$$

$$\mathcal{L}[P(\mathbf{x})]$$

$${\mathcal{L}}^{†}[Q^{\ast }(\mathbf{x})]$$

$$A{\Sigma }_x+{\Sigma }_x{A}^T+D=0.$$

$${\Sigma }_x$$

$${\Sigma }_{x,1,1}$$

$${\Sigma }_{x,2,2}$$

$${\Sigma }_{x,1,2}={\Sigma }_{x,2,1}$$

$${\Sigma }_{x,1,1}^{-1}$$

Full text

I Introduction

Both oscillations and noise are ubiquitous in many systems of interest in physics Hempstead and Lax (1967); Stratonovich (1967), biology Martin et al. (2003); Kummer et al. (2005); Elowitz and Leibler (2000) and neuroscience Ermentrout (2014); Laing and Lord (2010). Synchronization and entrainment of deterministic oscillators may be analyzed by finding a coordinate transformation to a one-dimensional phase variable defined in terms of the asymptotic phase of a limit cycle Brown et al. (2004); Guckenheimer (1975); Kuramoto (1975); Schwemmer and Lewis . The classical approach breaks down in several cases of interest including: (a) limit-cycle systems endowed with noise Hempstead and Lax (1967); Ebeling et al. (1986); Lindner (2002); Ushakov et al. (2005); Gleeson and O’Doherty (2006); Jülicher et al. (2009); Stiefel et al. (2010); Ditlevsen and Greenwood (2013), (b) spiral-sink systems with noise-sustained oscillations, also known as quasicycles Uhlenbeck and Ornstein (1930); Schimansky-Geier and Zülicke (1990); Lugo and McKane (2008); Wallace et al. (2011); Brooks and Bressloff (2015), and (c) systems with an attracting heteroclinic cycle Shaw et al. (2012); Thomas and Lindner (2014); Giner-Baldo et al. (2017).

Consider as an example for (b), the two-dimensional Ornstein-Uhlenbeck process in the form of two linear stochastic differential equations

[TABLE]

In these equations the $\xi_{i}$ represent independent delta-correlated Gaussian white noise sources with $\langle\xi_{i}(s)^{T}\xi_{j}(t)\rangle=\delta_{ij}\delta(s-t)$. A trajectory for a particular choice of parameters, corresponding to a spiral sink, is shown in Fig. 1. This system resembles a well behaved oscillator in many respects: it displays a noisy rotation around the origin in phase space (A), shows roughly oscillatory behavior also in the single components (B), and exhibits a pronounced peak at a non-vanishing frequency in the power spectrum of one of the components (C). The quality factor of the oscillation (frequency of the spectral peak divided by its width at half maximum height), is rather high in the chosen example ($Q_{\text{F}}\approx 10$). Despite the apparent oscillations, however, considering the noiseless system (setting $B_{ij}=0$ in eq. (I)), we cannot define a phase of this spiral-sink system (see below subsec. II.2).

In Thomas and Lindner (2014) we introduced a new definition of the asymptotic phase for robustly oscillatory stochastic systems, based on a spectral decomposition of the backward (or adjoint) Kolmogorov density operator. We note that an equivalent decomposition was independently introduced in the context of dephasing of genetic oscillators in Potoyan and Wolynes (2014) and an alternative definition of the phase for stochastic oscillators, based on a mean first passage time construction, was put forward in Schwabedal and Pikovsky (2013); in this paper we exclusively focus on the asymptotic phases based on the eigenfunction expansion.

Our asymptotic phase is well defined for noisy systems whether the underlying mean-field dynamics exhibits a stable limit cycle, spiral sink or a stable heteroclinic orbit, provided the eigenvalues and eigenfunctions of the Kolmogorov operator and its adjoint operator satisfy a set of natural conditions Thomas and Lindner (2014), detailed below. Moreover, in the case of a stable limit cycle system, when the classical asymptotic phase is well defined, the isochrons for the stochastic system with small noise levels closely resemble the isochrons of the deterministic system.

Many questions regarding the asymptotic phase and the forward phase are still open. We lack intuition about the shape of the isochrons, their dependence on the noise strength, and the difference between the backward and forward phases. The purpose of this paper is to explore these issues for a case where we can obtain analytical insights, that is, the two-dimensional Ornstein-Uhlenbeck process with complex eigenvalues of the drift matrix (a special case of eq. (I)). For this system, we derive explicit expressions for the asymptotic (backward) phase and for the forward phase. We hope to convince the reader that this example is nontrivial. Inter alia, we find qualitative differences in the forward and backward phases. For instance, the backward phase isochrons do not depend at all on the noise properties of the system, i.e. on $B_{ij}$, whereas the forward phase’s isochrons do. Furthermore, if we plot the asymptotic phase $\Psi(t)=\Psi(x_{1}(t),x_{2}(t))$ as a function of time (see Fig. 1D), it progresses more steadily than the geometric phase $\vartheta(t)=\arctan(x_{2}(t)/x_{1}(t))$. In fact, we will show that the average rate of increase for the backward phase does not depend on time, nor on the geometric phase.

Our paper is organized as follows. In the next section we recall the definitions of the backward and forward phases for stochastic systems, as well as the definition of the deterministic phase. In sec. III we introduce the two-dimensional Ornstein-Uhlenbeck model and our specific parametrization of it. In sec. IV we derive the two eigenfunctions of interest, extract the dependence of backward and forward phases on the geometric phase, and discuss their general properties. In sec. V we look at a number of examples that give some insight into the differences between the two phase types. We conclude in sec. VI with a brief discussion and outlook.

Table 1 provides a list of notation.

II Review of the phase definitions in deterministic and stochastic systems

In this section we recall notions of the asymptotic phase of an oscillator for deterministic and stochastic dynamical systems. For an $n$-dimensional dynamical system, the reduction to a $1$-dimensional phase description assigns a scalar variable (the phase) to each point in the underlying space. For a deterministic limit cycle oscillator, the phase variable $\theta$ should satisfy $d\theta/dt=\text{const}$ along trajectories. In the setting of stochastic dynamical systems, we can obtain a reduced description of an oscillatory Markov process in terms of the eigenfunctions of the generator of the process.

II.1 Stochastic setting

As in Thomas and Lindner (2014), we assume the state of the system is given by a point $\mathbf{X}$ in a (discrete or continuous) finite dimensional space, $\mathbf{X}\in\mathcal{X}$, and evolves in continuous time $t\in\mathbb{R}$ as a Markov process with transition density

[TABLE]

where $t>t_{0}$. We consider processes that are homogeneous in time, meaning $P(\mathbf{x},t\>|\>\mathbf{x}_{0},t_{0})=P(\mathbf{x},t-t_{0}\>|\>\mathbf{x}_{0},0)$.

We assume the density evolves according to a differential operator $\mathcal{L}$ with formal adjoint $\mathcal{L}^{\dagger}$. That is, the evolution with respect to the latter time $t$ is given by a forward Kolmogorov equation (or forward Fokker-Planck equation, in the case of a system driven by additive Gaussian white noise)

[TABLE]

while the evolution with respect to the earlier time $t_{0}$ is given by a backward or adjoint Kolmogorov equation (or backward Fokker-Planck equation, in the Gaussian case),

[TABLE]

The operator $\mathcal{L}^{\dagger}$ is also known as the generator of the Markov process Ethier and Kurtz (2009). We will assume the forward and backward operators $\mathcal{L}$ and $\mathcal{L}^{\dagger}$ have a biorthogonal set of eigenfunctions $P_{\lambda}$, $Q^{*}_{\lambda}$ satisfying

[TABLE]

where $\delta_{\lambda,\lambda^{\prime}}$ is the Kronecker delta. We consider the case in which this biorthogonal system is complete, in the sense that for any $\mathbf{x}\in\mathcal{X}$, the conditional density may be written in terms of the spectral decomposition

[TABLE]

We assume the system has a unique invariant probability distribution corresponding to the forward eigenfunction $P_{0}$ of the trivial eigenvalue $\lambda_{0}=0$; the corresponding backward eigenfunction is $Q_{0}\equiv 1$. We moreover assume the remaining eigenvalues, which may be real or complex, have negative real part.

As in Thomas and Lindner (2014), we will define a system to be “robustly oscillatory” if, in addition, it fulfills the following conditions:

The eigenvalue with least negative real part forms a complex conjugate pair $(\lambda,\lambda^{*})$, with $\lambda=\mu+i\omega$ and $\omega>0$. 2. 2.

The pair $(\lambda,\lambda^{*})$ has real part sufficiently far from the rest of the eigenvalue spectrum, specifically $\Re\lambda^{\prime}\leq 2\mu$ for all $\lambda^{\prime}\not=\mu\pm i\omega$. 3. 3.

The relaxation rate of the probability density is significantly slower than the oscillation, in the sense that $|\mu|\ll\omega$.

Under these assumptions, we may write the eigenfunction pair corresponding to the slowest decaying eigenvalue, $\lambda$, in complex notation as

[TABLE]

which can be regarded as the definitions of the forward ($\Phi$) and backward ($\Psi$) asymptotic phases of the system. That is, we may define the backward and forward phases for a point $\mathbf{x}$ as functions of the coordinates $\mathbf{x}$ as

[TABLE]

where $\arg(a+ib)=\tan^{-1}(b/a)$ is the complex argument of $z=a+ib$. As shown in Thomas and Lindner (2014), the long-term behavior of the probability density, as it approaches the steady-state distribution $P_{0}(\mathbf{x})$, will be dominated by the difference between the backward phase at the initial point $\Psi(\mathbf{x}_{0})$ and the forward phase at the later point $\Phi(\mathbf{x})$:

[TABLE]

This asymptotic behavior allows one to extract the forward and backward phase from numerical simulations of the stochastic process, and in principle also from data.

When considering a stochastic trajectory $\mathbf{X}(t)$ we may alternatively represent the phases as functions of time

[TABLE]

i.e. as stochastic processes in their own right.

II.2 Deterministic setting

Here, instead of an evolving density $P(\mathbf{x},t\>|\>\mathbf{x}_{0},t_{0})$, we consider a trajectory $\mathbf{x}(t)$ satisfying an ordinary differential equation with initial condition $\mathbf{x}_{0}$:

[TABLE]

If the system has a periodic solution with period $T$, that is, $\mathbf{x}(t)=\mathbf{x}(t+T)$, then we may define a time-like phase $\tau$ as a map from the orbit $\Gamma=\{\mathbf{x}(t)\>|\>0\leq t<T\}$ to the circle $\tau\in[0,T)$ satisfying the scalar differential equation

[TABLE]

with initial condition $\tau(0)=\tau_{0}$ Ermentrout and Terman (2010). The phase $\tau_{0}$ associated with a reference point $\gamma_{0}\in\Gamma$ may be chosen arbitrarily, and $\tau$ is interpreted mod $T$. If all trajectories with initial conditions close to the periodic orbit converge to $\Gamma$ then we have a stable limit cycle $\mathbf{x}=\gamma(t)$. The set of initial conditions converging to $\Gamma$ is its basin of attraction, $\mathcal{B}$. For any initial condition $\mathbf{x}_{0}\in\mathcal{B}$, one defines its asymptotic phase as the scalar $\tau(\mathbf{x}_{0})$ such that

[TABLE]

for the trajectory $\mathbf{x}(t)$ with initial condition $\mathbf{x}(0)=\mathbf{x}_{0}$. This phase reduction, introduced in Guckenheimer (1975), has proven invaluable in the study of weakly coupled and weakly driven oscillators Kuramoto (1975); Ermentrout and Kopell (1984); Brown et al. (2004); Hoppensteadt and Izhikevich (1997); Schwemmer and Lewis . As originally discussed by Guckenheimer Guckenheimer (1975), the asymptotic phase may equivalently be described in terms of isochrons, which are the level sets of a differentiable function $\mathcal{T}(\mathbf{x})$ satisfying

[TABLE]

for all points $\mathbf{x}\in\mathcal{B}$. Under the natural boundary conditions imposed by continuity of $\mathcal{T}$ at the limit cycle, we may identify $\tau=\mathcal{T}$, up to an additive constant (compare eq. (14)). For a limit cycle with period $T$, the normalization of the phase variable to the interval $[0,T)$ is conventional. In order to emphasize the phase as a map to the circle, the phase normalization $\theta\equiv\left(2\pi\tau/T\right)\in[0,2\pi)$ may be used instead. In this case eq. (16) is changed to $d\theta/dt=2\pi/T$.

Transient oscillatory activity may arise in deterministic models that do not possess limit cycles, for example spiral sink systems and stable heteroclinic cycles. A deterministic dynamical systems has a stable heteroclinic cycle if there is a closed attracting set $\Gamma_{\text{het}}$ composed of trajectories connecting a repeating sequence of saddle equilibrium points Holmes and Stone (1992); Krupa and Melbourne (1995); Rabinovich et al. (2006); Shaw et al. (2012); Horchler et al. (2015). In this situation, there is no periodic trajectory with a finite period. Instead, trajectories near $\Gamma_{\text{het}}$ traverse the same neighborhood of phase space with progressively longer and longer intervals required to pass each saddle point in turn. Because there is no finite period, the phase and the asymptotic phase cannot be defined; see Shaw et al. (2012) for a discussion of the phase reduction problem for the deterministic case, and Thomas and Lindner (2014) for the stochastic case.

A spiral sink system possesses a stable equilibrium point for which the Jacobian matrix has a complex conjugate pair of eigenvalues. As a simple example, consider the rotationally symmetric system

[TABLE]

(in standard polar coordinates) with $\mu<0$ and $\omega>0$. It is well known that one cannot assign an asymptotic phase to points in the basin of attraction of a spiral sink fixed point such as eq. (17), because (unlike for a limit cycle, cf. eq. (15)) all initial points converge to the same trajectory at long times, namely $\mathbf{x}(t)\to\mathbf{0}$ as $t\to\infty$. If we seek solutions of eq. (16) for the vector field eq. (17), we find a one-parameter family of solutions

[TABLE]

for arbitrary constant $k\in\mathbb{R}$. The corresponding “isochrons” are given, in polar coordinates, as logarithmic spirals $\vartheta(r)=\vartheta_{0}+\frac{k\omega-1}{|\mu|}\ln r$. Setting $k=1/\omega$ gives evenly spaced “spokes of a wheel” isochrons, $\vartheta=\text{const}$, but this is only one choice in an infinite collection of solutions consistent with constant “phase” progression.

In contrast, we will show in the remainder of the paper that when a system with spiral sink dynamics is perturbed by noise, the isochrons of the forward and backward phase may be defined uniquely and unambiguously. This surprising result is a major contribution of our paper.

III The model: two-dimensional spiral sink with white noise

We consider a two-dimensional Ornstein-Uhlenbeck process in a setting such that the origin becomes a stable sink. This is of course a classical stochastic process that has been well studied Uhlenbeck and Ornstein (1930); much of the information we review in the following and adapt to our special notation can be found in standard text books, e.g. Risken (1984) and Gardiner (1985).

The general Langevin equation is

[TABLE]

We first discuss suitable choices and the properties of the two matrices $A$ and $B$.

We assume that the two eigenvalues of $A$ are complex (they are then a complex conjugate pair). As we will show below, for the stochastic spiral sink, the eigenvalues of $A$ coincide with the slowest decaying eigenvalues of the forward and backward operators, discussed in subsec. II.1. Therefore, to avoid introducing superfluous notation, we will write the eigenvalues of $A$ as $\lambda_{\pm}=\mu\pm i\omega$; recall that $\mu<0$ and $\omega>0$. We will use a specific notation for the matrix $A$ by which we can write $\lambda_{\pm}$ in a convenient way:

[TABLE]

where the last equation defines the frequency of rotation $\omega$. Throughout, we require $\mu<0$ (stability of the fixed point at the origin) and $\omega_{0}>0$ (counter-clockwise rotation). Setting $\alpha_{\mu}=\alpha_{\omega}=0$, we would obtain a rotationally symmetric stable sink, thus, $\{\alpha_{\mu},\alpha_{\omega}\}$ quantify the deviation from this special case. Obviously, to keep complex eigenvalues, we have to require that $\omega_{0}^{2}>\alpha_{\mu}^{2}+\alpha_{\omega}^{2}$ which we will assume to hold true throughout the following. We note that the phases will not change when we use a different unit of time and, hence, we could set either the real or the imaginary part of the eigenvalue to unity, without loss of generality. However, for the sake of broader applicability, we keep $\mu$ and $\omega_{0}$ in the following.

The complex-valued left and right eigenvectors of $A$ satisfy

[TABLE]

(the asterisk denotes complex-conjugate transpose of a vector). They are given in terms of our parameters by

[TABLE]

where the normalization factors $N>0$ and $e^{i\eta}$ are chosen so that $\mathbf{w}_{\pm}^{*}\mathbf{v}_{\pm}=1$, from which follows $\eta=\pi-\mbox{atan}\left(\frac{\alpha_{\mu}}{\omega}\right)$.

Turning to the noise matrix $B$, we note that what enters the theory is only the symmetric diffusion matrix

[TABLE]

Here we have expressed $D$ by deviations from the isotropic-noise case ($\beta_{D}=\beta_{c}=0$) that can occur when the noise is of different strength in the two variables (quantified by $-1\leq\beta_{D}\leq 1$) or when the noise in both variables is correlated ($\beta_{c}>0$) or anticorrelated ($\beta_{c}<0$), with the constraint $\beta_{c}^{2}+\beta_{D}^{2}\leq 1$. We note that because the system is linear there is no qualitative change to be expected when we turn up the overall noise intensity $\varepsilon$. Without the fixed values of the decay rate $-\mu$ and noise intensity $\varepsilon$ we have the five parameters left: $\omega_{0},\alpha_{\mu},\alpha_{\omega},\beta_{D},$ and $\beta_{c}$.

The forward and backward operators corresponding to eq. (19) read (using Einstein’s summation convention)

[TABLE]

The covariance matrix $\Sigma_{x}$, corresponding to the stationary solution $\mathcal{L}[P]=0$, satisfies the Lyapunov equation

[TABLE]

Moreover, it can be expressed as Gardiner (1985)

[TABLE]

In what follows, we will assume $\Sigma_{x}$ to be invertible.

Inserting all the different terms, we arrive at

[TABLE]

By means of the inverse of the covariance matrix

[TABLE]

we can express the stationary probability density as follows:

[TABLE]

We note that the case of a symmetric probability density is more general than the special case of the rotationally symmetric drift matrix plus an isotropic noise. Formally, we may ask under what conditions the stationary variance will be a multiple of the identity, $\Sigma_{x}\propto I$. Setting $\Sigma_{x,1,1}=\Sigma_{x,2,2}$ and $\Sigma_{x,1,2}=0$ in eq. (36) yields a unique solution

[TABLE]

Thus a necessary condition for an isotropic distribution is that any asymmetry in the individual dissipation rates is compensated by a matched asymmetry in the individual noise terms ($\alpha_{\mu}=\mu\beta_{D}$), and any asymmetry in the coupling between the variables is compensated by an anticorrelation in the driving noise ($\alpha_{\omega}=\mu\beta_{c}$).

In the absence of a driving noise, the stationary distribution collapses to a delta distribution at the origin, and the covariance matrix reverts to the zero matrix. As shown in subsec. II.2, in this case one can no longer uniquely define an asymptotic phase function. However, as we will show in sec. IV, as long as the noise has finite amplitude and $\Sigma_{x}$ has full rank, the asymptotic phase obtained from the backward equation, as well as the phase from the forward equation, is well defined.

The cross-spectral matrix with the power spectra of $x_{1}(t)$ and $x_{2}(t)$ on the diagonal can be obtained from the well-known expression Gardiner (1985):

[TABLE]

Specifically, for the first element of the matrix, the power spectrum of the first variable, we obtain explicitly in terms of our parameters

[TABLE]

We note that the quality factor (the ratio between the peak frequency $f_{\textrm{peak}}$ and the full width at half maximum, $\Delta f$) of a robustly oscillatory OU process is very close to half the ratio of the imaginary and real parts of the system’s eigenvalue, i.e.

[TABLE]

IV Analytical expressions for the asymptotic phases

We aim to find the eigenfunctions with the smallest negative real part for the forward and backward operators. In particular, we seek the complex arguments of these functions, which define our forward and backward phases. For the general system eq. (19), no expresssions of the eigenfunctions in terms of elementary functions is known. However, Leen et al. Leen et al. (2016) recently derived expressions for a simpler system that we will use. This system is given by

[TABLE]

and is simpler because its covariance matrix is a multiple of the identity matrix:

[TABLE]

We can transform our original model, eq. (19), to eq. (43) obeying eq. (44) as follows:

[TABLE]

where $\mathbf{u}_{k}$ and $\gamma_{k}$ with $k=1,2$ are the normalized eigenvectors and eigenvalues of the original (symmetric) covariance matrix $\Sigma_{x}$, respectively. Indeed $\Sigma_{y}=\left\langle\mathbf{y}\mathbf{y}^{T}\right\rangle=\left\langle C\mathbf{x}\mathbf{x}^{T}C^{T}\right\rangle=C\Sigma_{x}C^{T}$ and using the property $u^{T}_{i}u_{j}=\delta_{ij}$, we obtain eq. (44). Furthermore, with this transformation follow relations between drift and diffusion matrices of the original and the transformed system:

[TABLE]

If, for long times, the conditional probability density in $\mathbf{y}$ is given by

[TABLE]

where $P_{\pm}(\mathbf{y})$ and $Q_{\pm}(\mathbf{y}_{0})$ are the corresponding eigenfunctions of forward and backward operators with smallest real part of the corresponding eigenvalues $\bar{\lambda}_{\pm}$. Applying eqs. (7), (9), and (12) in Leen et al. Leen et al. (2016), we obtain these eigenfunctions in terms of $\mathbf{v}_{y,\pm}=C\mathbf{v}_{\pm}$ and $\mathbf{w}^{*}_{y,\pm}=\mathbf{w}^{*}_{\pm}C^{-1}$, the corresponding right and left eigenvectors of the transformed drift matrix $A_{y}=CAC^{-1}$. They read:

[TABLE]

Inserting these expressions into eq. (32)-eq. (33), together with eq. (34) and the fact that $\Sigma_{y}=\frac{1}{2}I$, it can be checked that these functions satisfy the forward and backward equations, respectively, with eigenvalues $\lambda_{0}=0$ (for $Q_{0}$ and $P_{0}$) and $\bar{\lambda}_{\pm}=\lambda_{\pm}$ (for $P_{\pm}$ and $Q^{*}_{\pm}$, respectively) identical to the eigenvalues of $A$ given in eq. (22).

From the above asymptotic expansion of the density $P(\mathbf{y},t|\mathbf{y}_{0},0)$ and its transformation $P(\mathbf{x},t|\mathbf{x}_{0},0)=P(\mathbf{y}(\mathbf{x}),t|\mathbf{y}_{0}(\mathbf{x}_{0}),0)|d\mathbf{y}/d\mathbf{x}|$, we obtain

[TABLE]

Noting that for a linear transformation the Jacobian $|d\mathbf{y}/d\mathbf{x}|=const$, we can conclude that the first eigenfunctions in the original variables read

[TABLE]

The forward and backward phases can now be extracted from these expressions taking into account that the only complex-valued objects in them are $\mathbf{v}_{+}$ and $\mathbf{w}^{*}_{+}$. Specifically, for the backward phase we obtain the concise expression in terms of the geometric phase $\vartheta$:

[TABLE]

while the expression for the forward phase, involving the inverse of the covariance matrix, is more lengthy

[TABLE]

where

[TABLE]

It will prove useful to have the expressions for the derivatives, that can be interpreted as densities of the respective isochrons. For the backward phase, we obtain

[TABLE]

whereas for the forward phase, this derivative reads

[TABLE]

A number of conclusions can be drawn analytically from eqs. (54-57):

The asymptotic phase of a stochastic oscillator is uniquely defined in terms of the backward phase eq. (54). Hence, in marked contrast to the deterministic case, where we obtain an entire family of possible phase definitions (cf. eq. (18) and the surrounding discussion), we do not have any ambiguity in the phase definition in the stochastic case, except for an arbitrary additive off-set. As we will see, our asymptotic phase defines the phase uniquely even in the case of vanishing noise and extracts the most simple, i.e. the non-logarithmic, definition of phase (see points 2 and 3 below). 2. 2.

The isochrons (the lines of equal phase) are rays starting from the origin: a rescaling of the vector $\mathbf{x}_{0}$ (which describes a straight line starting in the origin), for instance, cannot change the complex phase of the complex number $\mathbf{w}_{+}^{*}\mathbf{x}_{0}$ (this factor just rescales the complex number itself); the same argument holds true for the forward phase, i.e. the complex argument of $\mathbf{x}^{T}\Sigma_{x}^{-1}\mathbf{v}_{+}$ upon varying $\mathbf{x}$. Another manifestation of this simple geometry is that both forward and backward phases are only functions of the geometric phase but not of the radial coordinate. 3. 3.

Changing the geometric phase $\vartheta$ of the vector $\mathbf{x}_{0}$, however, we do not necessarily get a proportional change in the complex argument of $Q^{*}_{+}$. This means that although the isochrons are the spokes of a wheel, these spokes are sometimes closer together and sometimes further apart, i.e. the functions $\Psi(\vartheta)$ and $\Phi(\vartheta)$ are only in exceptional cases linear (the phase densities $d\Psi(\vartheta)/d\vartheta$ and $d\Phi(\vartheta)/d\vartheta$ are then constant). 4. 4.

The backward phase cannot depend on the noise properties of the system, because it is entirely determined by the drift matrix $A$, more specifically, by the left eigenvector of $A$. Furthermore, the backward phase is also independent of the decay rate $-\mu$ (which does not enter the eigenvectors) but depends only on $\omega_{0},\alpha_{\omega},$ and $\alpha_{\mu}$. 5. 5.

The geometric phase will not progress at a steady pace, even in the deterministic case. To see this, transform the system to polar coordinates to find in the deterministic case ($\varepsilon=0$)

[TABLE]

Because of the nonlinearities on the right hand side, progress in the geometric phase will not be constant and this will not be different in general for a finite noise intensity $\varepsilon>0$. In contrast, taking the time derivative of the complex-valued backward function, that we write as a product of an amplitude $u$ and the complex exponential of the phase $\Psi$, we obtain

[TABLE]

Focusing on the underlined parts, dividing by the complex exponential, and averaging over the noise, we arrive at two simple equations for the time-dependent mean values

[TABLE]

Hence, the mean amplitude decays exponentially and the average backward phase evolves at a constant velocity. 6. 6.

The forward phase does in general depend on the noise matrix $D$, as well as on the drift matrix $A$, because it involves the inverse of the covariance matrix that is shaped by the coefficients in the matrix $D$. Hence, in general forward and backward phases are distinct. 7. 7.

With a few algebraic manipulations and a comparison of coefficients of trigonometric functions in the derivatives eq. (57) and eq. (56), it can be shown that for an isotropic noise ($\beta_{D}=\beta_{c}=0$), forward and backward phases (apart from their off-set) are shifted versions of each other with respect to the geometric phase, i.e. the derivatives are

[TABLE] 8. 8.

In the limit as $\mu\to 0$ (hence, quality factor $Q_{\text{F}}\to\infty$) both forward and backward densities are equal

[TABLE]

If $\mu$ is small compared to the effective frequency $\omega$, i.e. for large quality factor, the two phases will be close to each other. 9. 9.

For the case of a rotationally symmetric probability density, i.e. when eq. (39) holds true, it can be shown that the densities of the phase are shifted versions of each other and the shift is $\pi/2$:

[TABLE]

V Illustration of specific cases

Here we give several numerical examples of isochrons, their density, and the stationary probability density of the system. It is useful to classify the distinct cases of anisotropy that can be caused by anisotropic terms in the drift matrix ($\alpha_{\mu},\alpha_{\omega}$) or in the diffusion matrix ($\beta_{D},\beta_{c}$) by the two parameter vectors

[TABLE]

In the following, we demonstrate that i) for the complete isotropic case ($\mathbf{\alpha}=0,\mathbf{\beta}=0$) both phase isochrons are uniformly distributed (as can be expected); ii) for a uniform density of isochrons of the backward phase, the forward phase’s isochrons can be non-uniformly distributed (e.g. for ($\mathbf{\alpha}=0,\mathbf{\beta}\neq 0$); iii) for isotropic noise ($\mathbf{\beta}=0$) but anisotropic drift matrix ($\mathbf{\alpha}\neq 0$) both isochron densities are non-uniform, but are shifted versions of each other; iv) for an isotropic stationary distribution $P_{0}$ (which happens for $\mathbf{\alpha}=\mu\mathbf{\beta}$) there can be non-uniform densities of backward and forward phases; v) backward and forward phase are particulary different for $\mathbf{\alpha}=-\mu\mathbf{\beta}$.

Starting with the first case of a completely isotropic system ($\mathbf{\alpha}=\mathbf{\beta}=0$), we show in the first row of Fig. 2 the isochrons of the backward (A1) and forward (A3) phases in the $(x_{1},x_{2})$ plane, their density (i.e. their derivative with respect to the geometric phase $\vartheta$) as a function of $\vartheta$ (A2), and the stationary probability density (A4) vs. $(x_{1},x_{2})$. In this case, not surprisingly, we find a rotationally symmetric probability density and uniformly spaced spokes of a wheel as the isochrons of both forward and backward phase; the densities for both phases are flat. This is a direct demonstration that our asymptotic phase singles out the simplest of the possible solutions from the family of functions, eq. (18), possible for the deterministic system.

For an isotropic drift matrix ($\mathbf{\alpha}=0$) but anisotropic diffusion matrix ($\mathbf{\beta}\neq 0$), the asymptotic phase (backward phase) will remain uniform because it does not depend on $\beta$ (cf. conclusion 4 in sec. IV) and this we can indeed see in Fig. 2B1 and B2. The forward phase, in contrast, depends on $\beta$ and is now nonuniform as becomes apparent in Fig. 2B2 and B3. Also with an anisotropic noise we obtain in general an anisotropic probability density (cf. Fig. 2B4) as can be expected due to the correlations and unequal intensities of the noise sources driving the two components of the Ornstein-Uhlenbeck process.

Next, we choose an isotropic diffusion matrix ($\mathbf{\beta}=0$) but anisotropic drift matrix ($\mathbf{\alpha}\neq 0$). In this case, in accordance with conclusion 7 in sec. IV, eq. (60), both phase densities are shifted versions of each other as becomes evident in Fig. 2C1-3. The probability density is again anisotropic (Fig. 2C4).

As we have observed in sec. III, for a specific relation between drift and diffusion anisotropies, eq. (39), or in terms of our vector notation, for

[TABLE]

the stationary probability density is isotropic despite the possible anisotropies of the single $A$ and $D$ matrices. Indeed, if we combine the two anisotropies from Fig. 2B and C, we obtain parameters that obey this condition, and consequently, the stationary probability density is rotationally symmetric (Fig. 3A4). In this case, the isochrons of the two phases are still nonuniformly distributed (Fig. 3A2) and the directions of maximum density are orthogonal to each other (cf. Fig. 3A1 and A3). Indeed, as we have demonstrated in conclusion 9 in sec. IV, the two isochron densities are identical after a shift by $\pm\pi/2$ in the argument $\vartheta$, eq. (62).

For comparison, Fig. 3B shows what happens if we rotate the $\alpha$ and $\beta$ vectors by $\pi/2$ (with $\alpha$ and $\beta$ still obeying the condition eq. (64)): evidently, the probability density is still isotropic (Fig. 3B4) while the isochrons have been rotated by $\pi/4$ (see Fig. 3B1 and B3) and the isochron densities have merely shifted by the same angle.

We finally look at a case that seems to be far away from matching eq. (64) by switching the sign of $\mathbf{\beta}$ compared to Fig. 3. Then we get a pronounced anisotropy in the stationary probability density (cf. Fig. 4A4) and the densities of the isochrons of backward and forward phases look now qualitatively different (Fig. 4A2). These densities are not merely shifted versions of each other, but the forward phase displays much stronger variations in its density, which also can be seen on comparison of the spacing of the isochrons shown in Fig. 4A1 and A3. As before, an equal rotation of the vectors $\mathbf{\alpha}$ and $\mathbf{\beta}$ (as done in (Fig. 4B) does not change the picture qualitatively, but everything (isochrons, stationary density) is rotated.

VI Conclusions

In this paper we have derived and studied the asymptotic phase for a simple stochastic oscillator, the two-dimensional Ornstein-Uhlenbeck process. Difficulties that are pertinent to the deterministic version of this system (the linearized version of a spiral sink) can be resolved by introducing our unambiguous asymptotic phase and taking the limit of a vanishing noise intensity. Our phase, defined by the complex argument of the eigenfunction of the backward Kolmogorov operator, singles out the simplest solution (spokes of a wheel) from the entire family of possible phases, eq. (18). To obtain the deterministic phase we do not have to let the noise intensity go to zero because (for the specific system considered here, the two-dimensional Ornstein-Uhlenbeck process) the asymptotic phase does not depend on the intensity and correlation properties of the driving white Gaussian noise at all.

Although the asymptotic phase is independent of the radial variable leading to the spokes of a wheel for the isochrons, there is a nontrivial dependence of the density of isochrons as a function of the geometric phase angle. The density change is such that the mean progression in the asymptotic phase is always uniform - in contrast to the progression of the geometric phase if the drift function deviates from the isotropic case (i.e. when $\alpha_{\mu}\neq 0$ or $\alpha_{\omega}\neq 0$). This can be regarded as the defining feature of the asymptotic phase.

We also studied the phase related to the forward Kolmogorov operator’s eigenfunction with smallest real part of the eigenvalue, or the forward phase, for short. The isochrons for this phase turned out to be the spokes of a wheel as well, but the density of the spokes as a function of geometric phase was shown to be generally different from that of the backward phase and to be dependent on the diffusion matrix (i.e. on the intensities and correlations among the noise forces driving the two components of the system). This marked difference between forward and backward phase can be understood based on the role of forward and backward eigenfunctions in the evolution of the system, eq. (11). The forward eigenfunction which characterizes the evolution at long times forward in time is shaped by how correlated the noise may be in the two components. The backward phase (the true asymptotic phase) characterizes the initial conditions that brought the system to the considered point in time. In our examples we also demonstrated that even for a system with perfectly isotropic stationary probability density, the density of backward and forward phase may depend non-uniformly on the geometric phase.

We identified several cases where the densities of backward and forward phases are merely shifted versions of each other: for a rotationally symmetric stationary probability density (with a shift of $\pi/2$), for an isotropic noise matrix (with a shift of $\arctan(\mu/\omega_{0})$), and in the limit of infinite quality factor by $\mu\to 0$ (with a vanishing phase shift).

A simple extension of the results considered here would be an $n$-dimensional Ornstein-Uhlenbeck process with $n>2$ when two of the eigenvalues of the Kolmogorov operator form a complex conjugate pair with real part that is much smaller in absolute value than the real parts of all other eigenvalues. In this case, the system would rapidly collapse to a two-dimensional manifold (the slow fiber, see Wiggins (1994)), the phase reduction can be accomplished by projection onto this plane, and the isochrons of the asymptotic phase will be spokes of a wheel in this plane.

Our results may also prove useful in nonlinear stochastic systems where a spiral sink may coexist with other dynamics, e.g. with a limit cycle. For example, the subthreshold dynamics of the stochastic planar Morris-Lecar model can be described by a two-dimensional Ornstein-Uhlenbeck process Ditlevsen and Greenwood (2013). Subthreshold oscillations have been experimentally manipulated by Stiefel et al. Stiefel et al. (2010) as a way of using phase resetting (by inhibitory stimuli) as a mechanism for extended temporal integration. Phase resetting presumes the existence of a phase even though until now there was no well-defined notion of phase related to subthreshold oscillations (essentially a noise-driven dynamics close to a spiral sink). Our results provide a framework for further investigations in this direction.

VII Acknowledgments

PJT was funded by NSF grant DMS-1413770.