## A PARTIAL DIFFERENTIAL EQUATION FOR THE MEAN-RETURN-TIME PHASE OF PLANAR STOCHASTIC OSCILLATORS ∗

ALEXANDER CAO , BENJAMIN LINDNER , † ‡ AND PETER J. THOMAS §

Abstract. Stochastic oscillations are ubiquitous in many systems. For deterministic systems, the oscillator's phase has been widely used as an effective one-dimensional description of a higher dimensional dynamics, particularly for driven or coupled systems. Hence, efforts have been made to generalize the phase concept to the stochastic framework. One notion of phase due to Schwabedal and Pikovsky is based on the mean-return time (MRT) of the oscillator but has so far been described only in terms of a numerical algorithm. Here we develop the boundary condition by which the partial differential equation for the MRT has to be complemented in order to obtain the isochrons (lines of equal phase) of a two-dimensional stochastic oscillator, and rigorously establish the existence and uniqueness of the MRT isochron function (up to an additive constant). We illustrate the method with a number of examples: the stochastic heteroclinic oscillator (which would not oscillate in the absence of noise); the isotropic Stuart-Landau oscillator, the Newby-Schwemmer oscillator, and the Stuart-Landau oscillator with polarized noise. For selected isochrons we confirm by extensive stochastic simulations that the return time from an isochron to the same isochron (after one complete rotation) is always the mean first-passage time (irrespective of the initial position on the isochron). Put differently, we verify that Schwabedal and Pikovsky's criterion for an isochron is satisfied. In addition, we discuss how to extend the construction to arbitrary finite dimensions. Our results will enable development of analytical tools to study and compare different notions of phase for stochastic oscillators. 1

Key words. stochastic oscillator; noise; first-passage-time problems; phase description

AMS subject classifications. 60G99, 60J70, 60J25, 92B25

- 1. Introduction. Stochastic oscillations are an important phenomenon in many areas of science and technology. The activity of certain brain areas [4, 23], repetitive motor activities [13, 32], the motion of small organelles called hair bundle in our inner ear organs [2, 20], and the lasing intensity of lasers under certain conditions [35] are but a few examples that can be modeled by a noisy oscillator.

The dynamic and stochastic mechanisms for generating noisy oscillations are diverse as well. The classical linear example system is the underdamped harmonic oscillator driven by white Gaussian noise [31, 34], showing a finite phase coherence because of amplitude and phase fluctuations (similar systems appear in ecological and neural models and are referred to as quasicycles [18, 33]).

More difficult to analyze are noisy oscillations occurring in nonlinear stochastic systems. A prominent example is a limit cycle perturbed by noise [9, 28]. Qualitatively different is the role of noise in excitable systems [17], in heteroclinic oscillators [26, 29], or in certain spatially extended systems [5, 16] - in these cases the driving fluctuations

∗

Funding: B.L. was funded through Deutsche Forschungsgemeinschaft Grant LI 1046/2-1 and BMBFGrant 01GQ1001A. A.C. was supported by a grant from the Arthur I. Mendolia '41 Fellowship from the Case Alumni Foundation. A.C. and P.J.T. were funded by National Science Foundation grant DMS-1413770.

† Department of Mathematics, Applied Mathematics, and Statistics. Case Western Reserve University, Cleveland, Ohio, 44106, USA (axc487@case.edu).

‡ Bernstein Center for Computational Neuroscience and Department of Physics. Humboldt University, 10115 Berlin, Germany (benjamin.lindner@physik.hu-berlin.de).

§ Department of Mathematics, Applied Mathematics, and Statistics. versity, Cleveland, Ohio, 44106, USA (pjthomas@case.edu).

Case Western Reserve Uni-

1 Some results presented here appeared previously in the Master's thesis of the first author.

are required to observe an oscillation (although a noisy one).

In the description of deterministic oscillations and weakly coupled deterministic oscillators a phase description has turned out to be extremely useful in many situations but in particular in neuroscience [6, 7, 15, 25]. Hence, more recently, different approaches have been developed to define a phase also for stochastic oscillations [24, 29].

Schwabedal and Pikovsky (henceforth referred to as S&amp;P) introduced a notion of phase that is based on an algorithmic numerical procedure involving the meanfirst-return-time; for brevity, we will refer to the mean-return-time (MRT) , described below. Alternatively, two of the present authors used the complex phase of the eigenfunction of the backward Kolmogorov operator (the generator of the Markov process) to introduce an asymptotic phase for stochastic oscillators [29, 31]. The two approaches can be regarded as generalizations of the two (deterministically equivalent) notions of phase, one based on return times among a system of Poincar´ sections, the e other on the asymptotic convergence of trajectories [12]. The question of whether one of these phase notions is superior to the other one was subject to some debate [22, 30] but will not be pursued here. What will be addressed in this paper is the question of how the algorithmic procedure introduced by S&amp;P in [24] can be reformulated as the solution of a partial differential equation. We mainly restrict attention to planar Markovian systems that display stochastic oscillations. This is certainly more limited in scope than the algorithm suggested in [24], which has been also applied to non-Markovian examples. Deriving a partial differential equation to determine the isochrons for the two-dimensional system (the lines that define equal phase), however, opens opportunities of more rigorous mathematical explorations of the MRT phase.

In [24] S&amp;P proposed a definition for the phase of a stochastic oscillator in terms of a system of Poincar´ sections e { glyph[lscript] MRT ( φ , ) 0 ≤ φ ≤ 2 π , } foliating a domain R⊂ R 2 and possessing the MRT property. Intuitively, a section glyph[lscript] MRT satisfies the MRT property if for all points x ∈ glyph[lscript] MRT , the mean return time from x back to glyph[lscript] MRT is constant. S&amp;P write:

'For a noisy system we define the isophase surface J as a Poincar´ e surface of section, for which the mean first return time J → J , after performing one full oscillation, is a constant T, which can be interpreted as the average oscillation period. In order for isophases to be well defined, oscillations have to be well defined as well: for example in polar coordinates, the radius variable must never become zero, so that one can reliably recognize each 'oscillation.' Random processes for which this is not the case should be treated with care.'

S&amp;P describe an iterative method for obtaining such a surface from an ensemble of two-dimensional trajectories (as obtained, for example, by simulating the Langevin eq. (2.1)). Clearly, when naively interpreted, the MRT corresponds to a mean first passage time from x ∈ glyph[lscript] MRT to glyph[lscript] MRT , and is therefore identically zero; to avoid a trivial definition, S&amp;P require that the mean be calculated for all trajectories 'after performing one full oscillation.'

While the intended meaning is intuitively clear, a mathematical formulation of the problem has not previously been given. It is furthermore not clear whether the method proposed by S&amp;P will converge, nor whether the isochrons (lines sharing a common MRT phase) are uniquely determined. Last but not least, it is unclear whether the isochrons could also be obtained by the solution of the partial differential equation (PDE) that governs the mean-first-passage time to an absorbing boundary [10].

S&amp;P apply their iterative method to a wide range of systems (including examples of the kind above, eq. (2.1), but also three-dimensional Markovian systems, e.g. twodimensional dynamics with colored noise). We show here for the special class of two-dimensional stochastic oscillators driven by white Gaussian noise, under natural assumptions given below, that isochrons with the MRT property are uniquely defined, and that they correspond to the level curves of the solution of the meanfirst-passage-time PDE with a periodic-plus-jump condition on an arbitrary simple curve connecting the inner and outer boundaries.

Our paper is organized as follows. In § 2 we introduce a class of two-dimensional Markovian models and state our assumptions. Furthermore, we develop the MRT isochrons as the level sets of a function T , unique up to an additive constant; we formulate the boundary conditions complementing the standard partial differential equation governing this function. We provide two distinct derivations of the main result. In § 2.2 we utilize the unwrapping construction and argue that the MRT function should assume a certain limiting form far upstream of any given transverse boundary. In § 3 we leverage properties of strongly elliptic second order differential operators to give a mathematically precise statement and rigorous proof of the same result. Specifically, through an application of Fredholm theory for strongly elliptic PDEs, we establish existence and uniqueness of the system of mean-return-time isochrons under mild conditions. In addition, in § 4 we illustrate our results by a number of examples, for which we solve the corresponding PDE with the periodic-plus-jump condition numerically. In § 5 we briefly address the extension of our results to stochastic oscillators in n &gt; 2 dimensions. We conclude with a short discussion of our results and possible extensions to more general systems. In the supplementary materials we discuss details of the numerical method and provide an additional example.

- 2. Model and derivation of the main result. We first present the model of a two-dimensional stochastic oscillator driven by white Gaussian noise, and derive a PDE and its boundary conditions, the solution of which provides the isochrons of the MRT phase.
- 2.1. The considered model and the MRT phase. We consider a planar Langevin system of the form

$$\dot { x } & = f _ { x } ( x, y ) + g _ { x 1 } ( x, y ) \xi _ { 1 } ( t ) + g _ { x 2 } ( x, y ) \xi _ { 2 } ( t ) \\ \dot { y } & = f _ { y } ( x, y ) + g _ { y 1 } ( x, y ) \xi _ { 1 } ( t ) + g _ { y 2 } ( x, y ) \xi _ { 2 } ( t )$$

where ξ 1 2 , ( ) is white Gaussian noise with t 〈 ξ i ( ) t ξ j ( t ′ ) 〉 = δ ij δ t ( -t ′ ); the multiplicative noise terms are interpreted in the sense of Ito. 2 The domain R has the topology of an annulus, formed by the complement of one smooth (piecewise C 1 ) star-shaped domain within another (with common center point), see Fig. 2.1a. Outer and inner boundaries will be referred to as R outer and R inner , respectively.

We assume that the total noise amplitudes do not vanish anywhere in the domain, i.e.

$$( 2. 2 ) \quad \ g _ { x 1 } ^ { 2 } ( x, y ) + g _ { x 2 } ^ { 2 } ( x, y ) > 0 \ \text{ and } \ g _ { y 1 } ^ { 2 } ( x, y ) + g _ { y 2 } ^ { 2 } ( x, y ) > 0 \ \forall ( x, y ) \in \mathcal { R },$$

2 We could equivalently write the system dx = f x ( x, y ) dt + g x 1 ( x, y ) dW 1 ( ) + t g x 2 ( x, y ) dW 2 ( ) t and dy = f y ( x, y ) dt + g y 1 ( x, y ) dW 1 ( ) + t g y 2 ( x, y ) dW 2 ( ), t with dW 1 and dW 2 are the increments of independent standard Wiener processes. See § 3 for a rigorous formulation. For completeness we repeat the equations in Stratonovich form in § 10.

4

Fig. 2.1: Original ring-like domain with a counterclockwise net probability flow and reflecting inner and outer boundaries (a). An isochron (sketched in red) can be defined by the property that the mean return time from any point of this line to itself (after one complete rotation) is given by the overall mean rotation time of the stochastic oscillator. Five trajectories of the noisy heteroclinic oscillator, eq. (4.5), on the original domain (b) and in unwrapped coordinates (c). The common initial condition is marked with a black dot. The inner boundary of the annulus (b) is the square max( | x , | | y | ) = glyph[epsilon1] = 0 05, which is mapped to the lower edge of the unwrapped . domain, β = -1. The outer boundary of the annulus is mapped to the upper edge of the unwrapped domain, β = +1. The transformation maps point ( x, y ) to ( α, β ) such that tan( α ) = y/x , with α continuous along trajectories. For the quarter wedge x ≥ | y | of the annular domain, we set β = -1+2( x -glyph[epsilon1] ) / ( π 2 -glyph[epsilon1] ). In the other quarter wedges of the annular domain, β is determined by a similar construction. Vertical dashed line marks one full rotation from the initial point. Compare Fig. 4.5.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-4ebd064a-5b49-46a7-8748-e5b80779d2a3.jpg)

that the noise be nonsingular, i.e. 3

glyph[negationslash]

$$( 2. 3 ) \quad \ g _ { x 1 } ( x, y ) g _ { y 2 } ( x, y ) - g _ { x 2 } ( x, y ) g _ { y 1 } ( x, y ) \neq 0 \ \forall ( x, y ) \in \mathcal { R },$$

and that all drift and diffusion functions are smooth ( C 2 ) and bounded; in particular we require that a unique stationary probability density exists and we have a nonvanishing net stationary probability current. Without loss of generality, our proof will assume the current circulates counterclockwise. In § 4 we consider specific examples, some of which rotate counterclockwise and others clockwise.

2.2. Disambiguation of the return-time problem by mapping to an infinite strip. The region R is diffeomorphic to the annulus 0 &lt; R -≤ r ≤ R + &lt; ∞ . The latter can be mapped to an angular variable α x, y ( ) ∈ [0 , 2 π ) and an amplitudelike variable β x, y ( ) ∈ -[ 1 1], corresponding to a rectangular domain. , One possible

3 This condition is required for strong ellipticity, see the proof of existence in § 3.

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

E

p

w

d

K

7

9

n

u

F

V

B

z

W

q

r

y

C

0

A

3

8

&gt;

H

c

N

S

J

X

O

k

j

G

Q

v

2

Y

T

f

M

D

5

+

P

g

/

m

I

Z

R

o

U

L

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

+

w

S

B

P

L

8

n

d

v

z

X

A

2

q

G

f

p

g

&gt;

7

c

V

N

E

J

3

U

r

/

9

0

m

o

M

F

y

Z

I

Q

C

j

W

Y

k

u

O

H

K

5

R

T

D

mapping from the original domain to the rectangle is

$$\alpha ( x, y ) = \begin{cases} \text{arctan} ( y / x ), & | x | > 0 \\ \pi / 2, & x = 0, y > 0 \\ 3 \pi / 2, & x = 0, y < 0 \end{cases}$$

$$\beta ( x, y ) = - 1 + 2 \frac { \sqrt { x ^ { 2 } + y ^ { 2 } } - R _ { \text{inner} } ( \alpha ( x, y ) ) } { | R _ { \text{outer} } ( \alpha ( x, y ) ) - R _ { \text{inner} } ( \alpha ( x, y ) ) | }$$

Fig. 2.2: The original ring-like domain (a) (Fig. 2.1) is mapped to a true annulus (b), then to the geometric phase and amplitude with periodic boundary conditions (c) and, finally, to a system with unwrapped phase and amplitude (d). In (a), the first return to a simple connection between the two boundaries, glyph[lscript] , after one rotation has been performed, is not well defined: we cannot simply turn the line glyph[lscript] into an absorbing boundary because then the condition of the performed rotation is not met and if we start on glyph[lscript] the mean return time will be zero. In the unwrapped phase and amplitude (d), we deal with infinitely many concatenated copies of the system and thus the return-time problem is now well posed: it is given by the first passage problem between two adjacent copies of the line (e.g. from glyph[lscript] to glyph[lscript] ′ or from glyph[lscript] ′ to glyph[lscript] ′′ ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-66ede540-e056-4ff8-99b8-0d4c6f80f954.jpg)

After the change of variables, the SDE attains a new form

$$( 2. 6 )$$

$$\dot { \alpha } & = f _ { \alpha } ( \alpha, \beta ) + g _ { \alpha 1 } ( \alpha, \beta ) \bar { \xi } _ { 1 } ( t ) + g _ { \alpha 2 } ( \alpha, \beta ) \bar { \xi } _ { 2 } ( t ) \\ \dot { \beta } & = f _ { \beta } ( \alpha, \beta ) + g _ { \beta 1 } ( \alpha, \beta ) \bar { \xi } _ { 1 } ( t ) + g _ { \beta 2 } ( \alpha, \beta ) \bar { \xi } _ { 2 } ( t )$$

with drift and diffusion coefficients uniquely defined by the functions f x,y , g x 1 2 , , g y 1 2 , [10] and again we interpret the equations in the sense of Ito. Note that all coefficients and their derivatives are smooth and periodic in α , i.e. ∀ k ∈ Z and for j = 1 2, ,

$$f _ { \alpha, \beta } ( \alpha, \beta ) & = f _ { \alpha, \beta } ( \alpha + k \cdot 2 \pi, \beta ),$$

$$g _ { \alpha, j } ( \alpha, \beta ) & = g _ { \alpha, j } ( \alpha + k \cdot 2 \pi, \beta ),$$

$$g _ { \beta, j } ( \alpha, \beta ) & = g _ { \beta, j } ( \alpha + k \cdot 2 \pi, \beta ).$$

The exact form of the coefficients is not needed in the following, but we do require that also in the new variables the total noise intensity does not vanish anywhere and

the noise remain nonsingular glyph[negationslash]

$$g _ { \alpha 1 } ^ { 2 } ( \alpha, \beta ) & + g _ { \alpha 2 } ^ { 2 } ( \alpha, \beta ) > 0 \\ g _ { \beta 1 } ^ { 2 } ( \alpha, \beta ) & + g _ { \beta 2 } ^ { 2 } ( \alpha, \beta ) > 0 \\ g _ { \alpha 1 } ( \alpha, \beta ) g _ { \beta 2 } ( \alpha, \beta ) & - g _ { \alpha 2 } ( \alpha, \beta ) g _ { \beta 1 } ( \alpha, \beta ) \neq 0 \quad \forall ( \alpha, \beta ).$$

Now instead of using periodic boundary conditions with respect to α on the rectangular domain, we can extend the domain to the left and right by infinitely many copies of it, i.e. we can employ a well-known construction and use an unwrapped angle variable and in this way keep track of the rotations; see Fig. 2.2 for a sketch of the different mappings.

Specifically, the unwrapped phase now solves the problem of the ill-posed returntime problem in the original setup: it is unclear how to impose the condition of 'one performed rotation' in the original x -y space because with a stochastic driving we can never exclude that the trajectory might encircle the inner boundary in the opposite direction from the mean rotation. In the unwrapped-phase-amplitude domain, the return-time problem has a corresponding and well-posed first-passage-time problem from one copy of the line glyph[lscript] to the adjacent one on the right, e.g. from glyph[lscript] to glyph[lscript] ′ or from glyph[lscript] ′ to glyph[lscript] ′′ .

The problem of finding the phase implied by the mean-return or the mean-firstpassage time corresponds in this setting to the problem of finding the simple connecting curve glyph[lscript] such that the mean-first-passage time from any point on glyph[lscript] to the absorbing boundary at glyph[lscript] ′ (and reflecting boundaries at β = ± 1 and at α →-∞ ) is equal to a constant (the mean rotation time of the stochastic oscillator).

- 2.3. Forward and backward Kolmogorov operators, probability density and mean first-passage time function on the infinite strip. We can write down the two versions of the Kolmogorov equation, the first one of which is the FokkerPlanck equation

$$\partial _ { t } P ( \alpha, \beta, t ) = \mathcal { L } P ( \alpha, \beta, t ) = - \nabla \vec { J } ( \alpha, \beta, t ).$$

Here, P α,β,t ( ), if started with the initial condition P α,β, ( 0) = δ α ( -α 0 ) δ β ( -β 0 ), is the transition probability density for being at α and β at time t , if the system was at α 0 and β 0 at time t = 0. The vector glyph[vector] J = ( J α , J β ) denotes the probability current, the α component of which, J α ( α, β, t ) is related to the mean number of trajectories crossing a line α = const at a certain radial coordinate β . The probability current is defined in terms of the drift vector glyph[vector] f = ( f α , f β ) ᵀ and diffusion matrix G = gg ᵀ by

$$\vec { J } = \vec { f } P - \frac { 1 } { 2 } \left [ \begin{array} { c } \partial _ { \alpha } ( \mathcal { G } _ { \alpha \alpha } P ) + \partial _ { \beta } ( \mathcal { G } _ { \alpha \beta } P ) \\ \partial _ { \alpha } ( \mathcal { G } _ { \beta \alpha } P ) + \partial _ { \beta } ( \mathcal { G } _ { \beta \beta } P ) \end{array} \right ],$$

with g = ( g α 1 g α 2 g β 1 g β 2 ) . The reflecting boundary condition at β = ± 1 reads:

$$J _ { \beta } ( \alpha, \pm 1 ) \equiv 0 \ \forall \alpha,$$

where J β is the second (vertical) component of glyph[vector] J in the rectangular coordinates.

The time-dependent solution P α,β,t ( ) in the unwrapped phase does not approach a steady-state solution because the probability is constantly moving towards the right and diffusively spreading away from the center of mass. One can restore a steady-state

solution by lumping all probability and all currents into one period ( i.e. reverse the mapping from d to c in Fig. 2.2) by considering

$$( 2. 1 4 ) \ p ( \alpha, \beta, t ) = \sum _ { k = - \infty } ^ { \infty } P ( \alpha + 2 k \pi, \beta, t ) \xrightarrow { t \to \infty } \ p _ { 0 } ( \alpha, \beta )$$

$$( 2. 1 5 ) \ \vec { j } ( \alpha, \beta, t ) = \ \vec { \sum } \ \vec { J } ( \alpha + 2 k \pi, \beta, t ) \xrightarrow { t \to \infty } \ \vec { j _ { 0 } } ( \alpha, \beta ) = ( j _ { \alpha, 0 } ( \alpha, \beta ), j _ { \beta, 0 } ( \alpha, \beta ) ).$$

$$p ( \alpha, \beta, t ) & = \sum _ { k = - \infty } ^ { \infty } P ( \alpha + 2 k \pi, \beta, t ) \xrightarrow { t \to \infty } \quad p _ { 0 } ( \alpha, \beta ) \\ \vec { \jmath } ( \alpha, \beta, t ) & = \sum _ { k = - \infty } ^ { \infty } \vec { J } ( \alpha + 2 k \pi, \beta, t ) \xrightarrow { t \to \infty } \quad \vec { \jmath } _ { 0 } ( \alpha, \beta ) = ( j _ { \alpha, 0 } ( \alpha, \beta ), j _ { \beta, 0 } ( \alpha, \beta ) ).$$

The asymptotic (stationary) solution corresponds to the solution of eq. (2.11) with periodic boundary conditions over one period, e.g. p 0 (0 , β ) = p 0 (2 π, β ). This follows from the smoothness of the function P α,β,t ( ) as the solution of the Fokker-Planck equation under our assumptions on the coefficients. The corresponding stationary probability current glyph[vector] j 0 = ( j α, 0 , j β, 0 ) is related to the mean rotation rate as follows. If we integrate the α -component of the current over the entire range of the radial variable, this yields the mean rotation rate or, equivalently, the inverse of the mean rotation time:

$$\overline { J } = \underset { R _ { - } } { \overset { R _ { + } } { \int } } d \beta \, j _ { \alpha, 0 } ( \alpha, \beta ) = \frac { 1 } { \overline { \overline { T } } }.$$

Here α is an arbitrary but fixed value from the interval [0 , 2 π ). Hence, one way to determine the mean rotation time is to solve the stationary Fokker-Planck equation with periodic boundary conditions over the period and reflecting boundary conditions on the lower and upper boundaries in the radius and then to integrate the α -component of the current over the radius.

We may also write down the backward Kolmogorov equation for the probability density, which is a differential equation in terms of the initial position α , β 0 0 , but we use the operator of this equation, the backward operator, L † , in a different context. Let ˜ be a line glyph[lscript] parameterized as ˜ = glyph[lscript] { (˜( ˜ ) α β , β ˜ ) : -1 ≤ ˜ β ≤ - } 1 , with ˜ α a C 2 function of ˜ , β and α &lt; 2 π &lt; α β ˜( ), so that ( α, β ) is upstream of ˜ . glyph[lscript] In the system with the unwrapped phase we can formulate the equation for the mean-first-passage time ˜ ( T α,β ) from initial position ( α , β 0 0 ) to such a line, with adjoint reflecting boundary conditions at the interior and exterior boundaries of the cylindrical domain for β = ± 1 and for α →-∞ as follows:

$$\mathcal { L } ^ { \dagger } \tilde { T } ( \alpha, \beta ) = - 1,$$

where

(2.18)

$$\mathcal { L } ^ { \dagger } [ u ] = \left [ f _ { \alpha } \frac { \partial } { \partial \alpha } + f _ { \beta } \frac { \partial } { \partial \beta } + \frac { 1 } { 2 } \left ( \mathcal { G } _ { \alpha \alpha } \frac { \partial ^ { 2 } } { \partial \alpha ^ { 2 } } + \mathcal { G } _ { \alpha \beta } \frac { \partial ^ { 2 } } { \partial \alpha \beta } + \mathcal { G } _ { \beta \alpha } \frac { \partial ^ { 2 } } { \partial \beta \alpha } + \mathcal { G } _ { \beta \beta } \frac { \partial ^ { 2 } } { \partial \beta ^ { 2 } } \right ) \right ] [ u ].$$

See (3.3) for the adjoint reflecting (Neumann) boundary conditions. Equation (2.17) holds true in the infinite strip and the values ( α, β ) can be as far off to the left as we wish. In particular, because of the nonvanishing noise intensity, the solution ˜ ( T α,β ) is everywhere smooth and differentiable with respect to both variables [19].

Next we give an intuitive derivation of the jump-periodic boundary condition the MRT function should satisfy, which coincides with that established rigorously in

§ 3. We cannot simply lump everything back into one period as we did above with P α,β,t ( ) to obtain p α, β, t ( ), because ˜ T is a time and not a probability density; a direct naive summation would certainly not make sense. However, a kind of averaging summation of proper differences is meaningful as we will see in the following. We first note that since the coefficients in the Langevin equations are all periodic in 2 π , this differential equation is identical for ˜ ( T α,β ) and ˜ ( T α -2 kπ, β ) (i.e. for k = 1 2 3 , , . . . , L † [ ˜ ]( T α, β ) ≡ L † [ ˜ ]( T α -2 kπ, β ))). Furthermore, if the initial point is in our reference interval, i.e. 0 ≤ α ≤ 2 π , we expect that the time ˜ ( T α -2 kπ, β ) for k →∞ approaches the k -fold multiple of the mean rotation time T plus an order-one correction, which would imply that there exists a T max such that

$$\left | \tilde { T } ( \alpha - 2 k \pi, \beta ) - k \overline { T } \right | < T _ { \max }, \ \forall k \in \mathbb { Z }.$$

Then, the following sum of differences can be regarded as a kind of mean over the mean-first-passage times

$$T _ { N } ( \alpha, \beta ) = \frac { 1 } { N } \sum _ { k = 1 } ^ { N } ( \tilde { T } ( \alpha - 2 k \pi, \beta ) - k \overline { T } ).$$

The function T N ( α, β ) obeys

$$\mathcal { L } ^ { \dagger } T _ { N } ( \alpha, \beta ) = - 1.$$

because the constant terms ( kT ) drop out (every term of the backward operator starts with a derivative) and each of the N non-constant term contributes -1 /N to the right-hand side, i.e. in total the r.h.s. is -1 again.

We now define the function

$$T ( \alpha, \beta ) = \lim _ { N \to \infty } T _ { N } ( \alpha, \beta )$$

as the limit for infinite summation (assuming the limit exists) and argue that this function should obey a periodic-plus-jump condition. We first look at the boundary condition for T N for fixed values of α and β with the reference interval

$$& T _ { N } ( \alpha - 2 \pi, \beta ) = \frac { 1 } { N } \sum _ { k = 1 } ^ { N } ( \tilde { T } ( \alpha - 2 ( k + 1 ) \pi, \beta ) - k \bar { T } ) = \frac { 1 } { N } \sum _ { k = 2 } ^ { N + 1 } ( \tilde { T } ( \alpha - 2 k \pi, \beta ) - ( k - 1 ) \bar { T } ) \\ & = \sum _ { k = 1 } ^ { N } \frac { \tilde { T } ( \alpha - 2 k \pi, \beta ) - k \bar { T } } { N } + \bar { T } + \frac { \tilde { T } ( \alpha - 2 ( N + 1 ) \pi, \beta ) - ( N + 1 ) \bar { T } - ( \tilde { T } ( \alpha - 2 \pi, \beta ) - \bar { T } ) } { N }. \\ & - \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

The last term vanishes in the limit N →∞ and we obtain

$$T ( \alpha, \beta ) + \overline { T } = T ( \alpha - 2 \pi, \beta ).$$

There is no constraint on the value of α except that it be smaller than ˜( ˜ ) (paramα β eterizing the absorbing boundary ˜ ) glyph[lscript] and thus we can express the function T α,β ( ) by a steady decrease in α and a purely periodic function in α :

$$( 2. 2 4 ) \quad \ \ T ( \alpha, \beta ) = U ( \alpha, \beta ) - \frac { \overline { T } } { 2 \pi } \alpha \text{ with } U ( \alpha, \beta ) = U ( \alpha + 2 k \pi, \beta ), \ k \in \mathbb { Z }.$$

Two conclusions can be drawn. First, we can solve

$$\mathcal { L } ^ { \dagger } T ( \alpha, \beta ) = - 1$$

on the reference domain [0 , 2 π ] × -[ 1 1] with the periodic-plus-jump condition eq. (2.23) , taken at α = 0 and adjoint-reflecting boundary conditions on the lines β = ± 1:

$$T ( 2. 2 6 ) & & T ( 2 \pi, \beta ) + \overline { T } = & T ( 0, \beta ), \ \forall \beta \in [ - 1, 1 ],$$

$$\mathcal { L } _ { \text{ref} } T ( \alpha, \beta ) | _ { \beta = \pm 1 } = \ 0, \quad \forall \alpha \in [ 0, 2 \pi ].$$

Here L refl u = 0 represents the adjoint reflecting boundary conditions, namely L refl = G βα ∂ α + G ββ ∂ β . As we establish in § 3, under these conditions, the function is only determined up to a constant and this constant is indeed the only trace left of the original absorbing boundary ˜ . glyph[lscript] Secondly, because eq. (2.23) holds true for all values of α , the contour lines of T α,β ( ) = T 0 fulfill the condition for an isochron glyph[lscript] MFP by S&amp;P, because from every point of the contour line glyph[lscript] MFP the mean time to the 2 π shifted version of the contour line glyph[lscript] ′ MFP is equal to T by virtue of the periodicity condition eq. (2.23), i.e.

$$T ( \alpha, \beta ) \Big | _ { ( \alpha, \beta ) \in \ell _ { M F P } } = T ( \alpha, \beta ) \Big | _ { ( \alpha, \beta ) \in \ell ^ { \prime } _ { M F P } } + \overline { T }.$$

Finally, we note that the MRT phase Θ is given in terms of T and an arbitrary constant Θ 0 as

$$\Theta ( \alpha, \beta ) = \Theta _ { 0 } - T ( \alpha, \beta ) \frac { 2 \pi } { \overline { T } }.$$

With this definition, the MRT phase then obeys partial differential equation

$$\mathcal { L } ^ { \dagger } \Theta = \frac { 2 \pi } { \overline { \overline { T } } }$$

with periodic-plus-jump condition Θ → Θ+2 π along an arbitrary cut.

If the problem is given originally in polar coordinates (which will be the case for a few of our examples), i.e. if the starting point are the Langevin eq. (2.6), we have solved the problem of finding an equation for the MRT phase proposed by S&amp;P. If the problem is stated originally in ( x, y ), eq. (2.1), we reverse the mapping in Fig. 2.2, formulate the problem as a mean-return time problem in the original domain with a periodic-plus-jump condition and find that the contour lines of the function T x, y ( ) are the isochrons of the MRT phase in the sense of S&amp;P. Indeed, in accordance with equation eq. (2.29), the MRT phase in ( x, y ) coordinates is given by

$$\Theta ( x, y ) = \Theta _ { 0 } - T ( x, y ) \frac { 2 \pi } { \overline { T } }.$$

Briefly, we first have to solve for the stationary probability density

$$( 2. 3 2 ) \quad \mathcal { L } p _ { 0 } ( x, y ) = 0 \text{ with } \mathcal { R } _ { p } p _ { 0 } ( x, y ) | _ { ( x, y ) \in R _ { \pm } } = 0, \, \int \int d x \, d y \, p _ { 0 } ( x, y ) = 1$$

from which we can extract the stationary current and the mean rotation period

$$\vec { j } = \mathcal { J } p _ { 0 } ( x, y ) = ( j _ { x } ( x, y ), j _ { y } ( x, y ) ), \ \overline { T } = \frac { 2 \pi } { \int _ { R _ { - } ( y = 0 ) } ^ { R _ { + } ( y = 0 ) } d x \, j _ { y } ( x, 0 ) }.$$

Here R p p 0 = n ᵀ ± ∇ p 0 = 0 gives the reflecting boundary conditions for the (forward) Fokker-Planck operator, where n ± is the local unit normal vector for the boundary at ( x, y ) ∈ R ± , respectively. We then solve the equation for the mean return time in the original domain with a jump across the line R -≤ y ≤ R ,x + = 0,

$$( 2. 3 4 ) \ \mathcal { L } ^ { \dagger } T ( x, y ) = - 1 \text{ with } \mathcal { R } _ { T } T ( x, y ) | _ { ( x, y ) \in R _ { \pm } } = 0, \lim _ { \varepsilon \to 0 ^ { + } } \left ( T ( - \varepsilon, y ) - T ( \varepsilon, y ) \right ) \overline { T }.$$

The backward operator is written out in the supplement ( 10) for both the Ito and § Stratonovich interpretations of the original stochastic differential equations. Furthermore, in eq. (2.34) R T u = ∑ j =1 2 , n j ∑ k =1 2 , G jk ∂ u k = 0 gives the adjoint reflecting or Neumann boundary conditions. The isochrons of the MRT phase are then given by the contour lines of T x, y ( ) modulo T if the line crosses the boundary condition on y = 0. Fig. 8.1 (supplemental material) illustrates the construction for a square domain used to solve eq. (2.34).

We round out our discussion of the MRT phase derivation with several comments:

Remark 1. The jump condition can be imposed on any simple connection between the inner and outer boundaries. Furthermore, by the PDE and its boundary conditions, T x, y ( ) is only determined up to an additive constant.

Remark 2. It may appear strange that we have to solve two PDEs (one for p 0 ( x, y ) and one for T x, y ( )) to solve the problem, especially because the first equation serves only to determine T , the mean rotation time, which should be computable from the PDE for the MRT in the first place. In fact, it can be shown by Green's-function techniques, that the boundary conditions of the PDE for T x, y ( ) uniquely determine in a self-consistent manner T . Thus, it would in principle suffice to solve only the PDE for T x, y ( ). In all investigated examples we found, however, that the subsequent solution of the two equations is numerically more practical and efficient.

Remark 3. The mapping to the annulus assumes a smooth boundary in the ( x, y ) domain. However, we do not see any problem in principle to generalize the arguments used to a mapping that is only piecewise smooth. One of our examples, the heteroclinic oscillator on the square domain, is of such a type and does not seem to be problematic.

Remark 4. Although, as elaborated in § 3 we require the noise to be nonsingular in order to rigorously establish existence and uniqueness of solutions, it is possible that this condition may not be strictly necessary. For example, in the numerical example with polarized noise below ( 4.3) the nonsingular noise condition is violated, § nevertheless the numerical procedure based on the PDE appears to give the correct system of isochrons.

## 2.4. Relation to the asymptotic phase for deterministic systems. If a deterministic system of ordinary differential equations

$$\frac { d x } { d t } = f ( x ), \ \ x \in \mathbb { R } ^ { n }$$

has a hyperbolically stable limit cycle solution γ (a closed, isolated periodic orbit), with period τ , we may define the phase Θ ∈ [0 , 2 π ) for points in Γ = { γ t ( ) | 0 ≤ t &lt; τ } so that Θ = 0 for some reference point on Γ and d Θ /dt = 2 π/τ . We can extend this function to define the asymptotic phase Θ for any point x 0 in the basin of attraction (the stable manifold) of Γ, so that the trajectory starting at x 0 converges to the

periodic trajectory γ t ( + τ Θ( x 0 ) / 2 π ) as t →∞ [8, 12]. Applying the chain rule, the asymptotic phase function must satisfy

$$\frac { d \Theta ( { \mathbf x } ) } { d t } = { \mathbf f } ( { \mathbf x } ) ^ { \mathbf \tau } \nabla \Theta ( { \mathbf x } ) = \frac { 2 \pi } { \tau }$$

for all x in the basin of attraction of Γ. The boundary condition for this linear firstorder and nonhomogeneous partial differential equation is set by continuity of Θ at the limit cycle Γ.

Now consider a family of stochastic differential equations of the form eq. (2.1), but with the noise scaled by a small parameter √ glyph[epsilon1] , and the corresponding family of solutions to the backward eq. (2.25), which we may write (for the planar case x ∈ R 2 ) as

$$\frac { d \mathbf x } { d t } = \mathbf f ( \mathbf x ) + \sqrt { \epsilon } \left ( g _ { 1 } ( \mathbf x ) \xi _ { 1 } ( t ) + g _ { 2 } ( \mathbf x ) \xi _ { 2 } ( t ) \right ) \\ \cdot \quad \epsilon \longrightarrow \quad \cdot$$

$$( 2. 3 8 ) \quad \ \ - 1 = \mathcal { L } _ { \epsilon } ^ { \dagger } \left [ T _ { \epsilon } ( { \mathbf x } ) \right ] = f ( { \mathbf x } ) ^ { \mathbf r } \nabla \left [ T _ { \epsilon } ( { \mathbf x } ) \right ] + \frac { \epsilon } { 2 } \sum _ { i j } \mathcal { G } _ { i j } \partial _ { i j } ^ { 2 } T _ { \epsilon } ( { \mathbf x } ).$$

Suppose that, as glyph[epsilon1] → 0, T glyph[epsilon1] ( x ) converges uniformly on compact subsets of the domain to a C 2 function T 0 ( x ). Since for any glyph[epsilon1] , T glyph[epsilon1] is defined only up to an additive constant, we consider convergence in the sense that for arbitrary nonzero vectors v ∈ R 2 ,

$$\mathbf v ^ { \tau } \left ( \nabla T _ { 0 } - \nabla T _ { \epsilon } \right ) \rightarrow 0,$$

for all x in the domain. Fixing x and setting v = ( ), we see that for each f x x ,

$$\left ( 2. 4 0 \right ) \quad f ^ { \intercal } \nabla T _ { 0 } - f ^ { \intercal } \nabla T _ { \epsilon } = f ^ { \intercal } \nabla T _ { 0 } - \left ( - 1 - \frac { \epsilon } { 2 } \sum _ { i j } \mathcal { G } _ { i j } \mathcal { Q } _ { i j } ^ { 2 } T _ { \epsilon } \right ) \to 0, \text{ as } \epsilon \to 0,$$

where we have used eq. (2.38). Consequently if T glyph[epsilon1] converges to a well-behaved function T 0 in this way, it must satisfy

$$\mathcal { L } _ { 0 } ^ { \dagger } \left [ T _ { 0 } \right ] = \mathfrak { f } ( \mathbf x ) ^ { \intercal } \nabla \left [ T _ { 0 } ( \mathbf x ) \right ] = - 1.$$

Comparing eq. (2.36) and eq. (2.41), evidently if eq. (2.35) has a stable limit cycle, then the function T 0 must correspond with the deterministic asymptotic phase function Θ through the linear relation

$$\Theta ( { \mathbf x } ) & = \Theta _ { 0 } - T _ { 0 } ( { \mathbf x } ) \frac { 2 \pi } { \tau },$$

for arbitrary constant Θ 0 .

3. Proof of existence and uniqueness of the isochron function T . In this section we use Fredholm theory for strongly elliptic second order operators, and a maximum principle, to prove that the partial differential equation defining the meanreturn-time isochron function T has a solution, and that the solution is unique up to an additive constant. As in the previous section, we adopt coordinates α ∈ [0 , 2 π ) (the angular, periodic coordinate) and β ∈ [ -1 1] (the radial coordinate). , We will refer to the fundamental or local domain Ω = [0 , 2 π ) × [0 , 1] and the extended domain Ω ext = R × -[ 1 1]. , We assume the following:

- A1. Transformed into ( α, β ) coordinates, the trajectories ( α t , β ( ) ( )) of a strongly t Markovian time-homogeneous process obey an Ito equation

$$d \alpha & = f _ { 1 } ( \alpha, \beta ) \, d t + g _ { 1 1 } ( \alpha, \beta ) \, d W _ { 1 } ( t ) + g _ { 1 2 } ( \alpha, \beta ) \, d W _ { 2 } ( t ) \\ d \beta & = f _ { 2 } ( \alpha, \beta ) \, d t + g _ { 2 1 } ( \alpha, \beta ) \, d W _ { 1 } ( t ) + g _ { 2 2 } ( \alpha, \beta ) \, d W _ { 2 } ( t )$$

where f , g i ij are C 2 on Ω ext .

(Note: we will refer often below to the matrix G = gg ᵀ , defined in terms of the g ij .)

- A2. The functions f i and g ij are periodic in the first coordinate with period 2 π , i.e. ∀ α ∈ R and i = 1 2, , f i ( α +2 π, β ) = f i ( α, β ), and likewise for each g ij .
- A3. The second order differential operator P is strongly elliptic, where P is defined (following McLean)

$$\mathcal { P } u = - \sum _ { j = 1 } ^ { 2 } \sum _ { k = 1 } ^ { 2 } \partial _ { j } ( A _ { j k } \partial _ { k } u ) + \sum _ { j = 1 } ^ { 2 } A _ { j } \partial _ { j } u \text{ on } \Omega _ { \text{ext} }$$

and A jk = - G 1 2 jk and A j = -1 2 ∑ 2 k =1 ∂ k G jk + f j . Along the lines β = ± 1 we impose Neumann boundary conditions

$$0 = \sum _ { j = 1 } ^ { 2 } \nu _ { j } \sum _ { k = 1 } ^ { 2 } \mathcal { G } _ { j k } \partial _ { k } u$$

where ν = (0 , ± 1) is the outward unit normal at the respective boundary. We note that in the planar case ( n = 2) we consider here, strong ellipticity is guaranteed if the matrix

$$\mathcal { G } = g g ^ { \intercal } = \left ( \begin{array} { c c } g _ { 1 1 } ^ { 2 } + g _ { 1 2 } ^ { 2 } & g _ { 1 2 } g _ { 2 2 } + g _ { 2 1 } g _ { 1 1 } \\ g _ { 1 2 } g _ { 2 2 } + g _ { 2 1 } g _ { 1 1 } & g _ { 2 1 } ^ { 2 } + g _ { 2 2 } ^ { 2 } \end{array} \right )$$

satisfies the nondegeneracy conditions

$$\mathcal { G } _ { 1 1 } > 0, \ \mathcal { G } _ { 2 2 } > 0, \text{ and } \det \mathcal { G } \neq 0.$$

McLean's differential operator P corresponds to Kolmogorov's backward operator L † occurring in the first-passage/return time problem, also known as the generator of the Markov process. The adjoint P ∗ corresponds to Kolmogorov's forward operator (the Fokker-Planck operator) with reflecting boundary conditions at β = ± 1, describing the evolution of probability densities forward in time. We call (3.3) adjoint reflecting boundary conditions because P and P ∗ are adjoint operators on the appropriate function spaces. Specifically, P acts on L ∞ (Ω) ∩ C 2 (Ω) (bounded, twice differentiable functions) while P ∗ acts on L 1 (Ω) ∩ C 2 (Ω) (twice differentiable functions integrable on the local domain).

Note that, although we consider two independent noise sources, dW 1 and dW 2 , we could also have k &gt; 2 independent noise sources (in this case g would be an n × k matrix) without a fundamental change in the results, provided G is nonsingular.

We further assume that:

- A4. The process viewed on Ω (taking α mod 2 π ) admits a density ρ α, β, t ( ) evolving according to

$$( 3. 6 ) \quad \frac { \partial \rho } { \partial t } = \mathcal { L } \rho = \mathcal { P } ^ { * } \rho = - \sum _ { j = 1 } ^ { 2 } \sum _ { k = 1 } ^ { 2 } \partial _ { j } ( A _ { k j } ^ { * } \partial _ { k } \rho ) - \sum _ { j = 1 } ^ { 2 } \partial _ { j } ( A _ { j } ^ { * } \rho )$$

where A ∗ kj = A jk and A ∗ j = A j , and we impose reflecting (Neumann) boundary conditions

$$0 = \sum _ { k = 1 } ^ { 2 } A _ { k j } ^ { * } \partial _ { k } \rho + A _ { j } ^ { * } \rho$$

at β = ± 1 (for all α ), and periodic boundaries in α , i.e.

$$\forall \beta \in [ - 1, 1 ], \ \rho ( 0, \beta ) = \rho ( 2 \pi, \beta ).$$

We assume that the system has a unique stationary distribution ρ ss ≥ 0, with 1 = ∫ Ω dαdβρ ss ( α, β ), i.e. satisfying the homogeneous equation

$$\mathcal { P } ^ { * } \rho _ { s s } = 0,$$

together with the boundary conditions (3.7)-(3.8).

We note the stationary flux vector J ss ( α, β ) corresponds componentwise to

$$J _ { \text{ss}, j } ( \alpha, \beta ) = f _ { j } \rho _ { \text{ss} } - \frac { 1 } { 2 } \sum _ { k = 1 } ^ { 2 } \partial _ { k } \left ( \mathcal { G } _ { j k } \rho _ { \text{ss} } \right ).$$

Lastly, we assume the mean drift is nonzero and oriented to the right:

A5. If γ : [ -1 1] , → [0 , 2 π ] is any C 1 function whose graph C γ = ( { α = γ β , β ( ) ) : -1 ≤ β ≤ 1 } connects the inner and outer domain, separating Ω ext into left and right connected components, with unit normal n ( β ) oriented into the right connected component, then the mean rightward flux through C γ is positive, i.e.

$$0 < \overline { J } \coloneqq \int _ { - 1 } ^ { 1 } d \beta \mathbf n ^ { \intercal } ( \beta ) J _ { \text{ss} } ( \gamma ( \beta ), \beta ).$$

The reciprocal of J is proportional to the mean period of the oscillator,

$$\overline { T } & = \left ( \overline { J } \right ) ^ { - 1 }.$$

Our goal is to establish the existence of a function T satisfying the inhomogeneous PDE with adjoint reflecting boundary conditions at β = ± 1 and jump-periodic boundary conditions, namely

$$\mathcal { P } T = - 1, \quad & \text{on} \ \Omega \\ \sum _ { k = 1 } ^ { 2 } \mathcal { G } _ { 2 k } \partial _ { k } T ( \alpha, \pm 1 ) = 0, \quad & \forall \alpha \in \mathbb { R } \\ T ( \alpha, \beta ) - T ( \alpha + 2 \pi, \beta ) = \overline { T }, \quad & \forall ( \alpha, \beta ) \in \Omega _ { \text{ext} }.$$

Theorem 3.1 (Existence and Uniqueness of Solutions). If Assumptions A1A5 hold, then the partial differential equation (3.13) with reflecting adjoint boundary conditions at β = ± 1 and jump-periodic boundary conditions in the α coordinate has a solution T α,β ( ) on Ω ext and, by restriction, on Ω . Moreover, the solution is unique up to an additive constant.

The proof relies on Theorem 4.10 of [19] (stated as Theorem 3.3 below), a version of the Fredholm alternative for strongly elliptic PDEs with general boundary conditions. This theorem employs the following notation. For a Lipschitz domain Ω, and s ∈ R , H s (Ω) denotes the Sobolev space of order s (based on the L 2 norm). H s D (Ω) denotes the subspace of functions in H s (Ω) that equal zero when restricted to the portion of the boundary of Ω where a Dirichlet condition is enforced. Since our boundaries do not include a Dirichlet component, for our problem H s D (Ω) ≡ H s (Ω). In addition, ˜ H s (Ω) denotes the closure of C ∞ comp (Ω) in H s ( R n ) (for the planar systems we consider, n = 2) where C ∞ comp (Ω) is the space of C ∞ functions with compact support in Ω. The space H -1 (Ω) is the set of distributions (or functions) h α, β ( ) that can be represented as h α, β ( ) = ∂ α h 1 ( α, β ) + ∂ β h 2 ( α, β ) for some functions h , h 1 2 ∈ L 2 (Ω). The inhomogeneities h α, β ( ) we consider will be C 2 functions of the coordinates, hence integrable over the compact domain Ω, and members of H -1 (Ω). For further details, see [19].

Theorem 4.10 will require that the operator P be coercive on H s D (Ω). Coercivity is a growth rate condition that holds for our system by virtue of the following

Theorem 3.2 (McLean 2000, Thm 4.7). Assume that P has scalar coefficients, and that P is strongly elliptic on Ω . If the leading coefficients satisfy

$$A _ { k j } = A _ { j k } \ \ o n \ \Omega, \ \ \ f o r \ a l l \ j \ a n d \ k,$$

then P is coercive on H 1 (Ω) .

The symmetry condition A kj = A jk holds for our system because G = gg ᵀ is a symmetric matrix. Strong ellipticity follows from our assumptions A1-A5, and will be established in the proof of our main theorem (3.1) below.

Finally, we will refer to the homogeneous problem

$$\mathcal { P } u & = 0, \quad \text{on} \ \Omega \\ \sum _ { k = 1 } ^ { 2 } \mathcal { G } _ { 2 k } \partial _ { k } u ( \alpha, \pm 1 ) & = 0, \quad \forall \alpha \in [ 0, 2 \pi ] \\ u ( 0, \beta ) - u ( 2 \pi, \beta ) & = 0, \quad \forall \beta \in [ - 1, 1 ]$$

associated with the general inhomogeneous problem (with periodic boundary conditions)

$$\mathcal { P } u = h, \ \text{ on } \Omega$$

$$\sum _ { k = 1 } ^ { 2 } \mathcal { G } _ { 2 k } \partial _ { k } u ( \alpha, \pm 1 ) & = g _ { N }, \quad \forall \alpha \in [ 0, 2 \pi ] \\ u ( 0, \beta ) - u ( 2 \pi, \beta ) & = 0, \quad \forall \beta \in [ - 1, 1 ]$$

McLean's theorem 4.10 is more general than our problem requires. We specialize to the case of scalar-valued functions ( m = 1, below), and our boundary contains no Dirichlet component (Γ D = ∅ ). Nevertheless for ease of comparison we state McLean's version of the theorem in full:

Theorem 3.3 (McLean 2000, Thm 4.10). Assume that Ω is a bounded Lipschitz domain, and that P is coercive on H 1 D (Ω) m . Let h ∈ ˜ H -1 (Ω) m , g D ∈ H 1 / 2 (Γ D ) m and g N ∈ H -1 / 2 (Γ N ) m , and let W denote the set of solutions in H 1 (Ω) m to the homogeneous problem (3.14) .There are two mutually exclusive possibilities:

- (i) The homogeneous problem has only the trivial solution, i.e., W = 0 { } . In this case, the homogeneous adjoint problem (3.9) , with boundary conditions (3.7) -(3.8) , also has only the trivial solution in H 1 (Ω) m , and for the inhomogeneous problem (3.15) we get a unique solution u ∈ H 1 (Ω) m . Moreover,

$$( 3. 1 6 ) \ | | u | | _ { H ^ { 1 } ( \Omega ) ^ { m } } \leq C | | h | | _ { \widetilde { H } ^ { - 1 } ( \Omega ) ^ { m } } + C | | g _ { D } | | _ { H ^ { 1 / 2 } ( \Gamma _ { D } ) ^ { m } } + C | | g _ { N } | | _ { H ^ { - 1 / 2 } ( \Gamma _ { N } ) ^ { m } }.$$

- (ii) The homogeneous problem has exactly p linearly independent solutions, i.e. dim W = p , for some finite p ≥ 1 . In this case, the homogeneous adjoint problem (3.9) , with boundary conditions (3.7) (3.8) -, also has exactly p linearly independent solutions, say v , . . . , v 1 p ∈ H 1 (Ω) m , and the inhomogeneous problem (3.15) is solvable in H 1 (Ω) if and only if

$$( 3. 1 7 ) \quad \ \ ( v _ { j }, h ) _ { \Omega } + ( \gamma v _ { j }, g _ { N } ) _ { \Gamma _ { N } } = \left ( \widetilde { B } _ { \nu } v _ { j }, g _ { D } \right ) _ { \Gamma _ { D } } \quad \text{for $1\leq j\leq p$.}$$

In (3.17) γ refers to the 'trace operator' that maps a function v to its value on the boundary.

Before proving Theorem 3.1 we state and prove a lemma concerning the solutions of the homogeneous problem, exploiting a maximum principle.

Lemma 3.4. If Assumptions A1-A5 hold, then every C 2 solution u of the homogeneous problem (3.14) is a constant function.

Proof of Lemma 3.4. Suppose a C 2 function u α, β ( ) satisfying (3.14) on the domain Ω = [0 2 , π ] × -[ 1 1] , attains a strict maximum u + at a point ( α , β 0 0 ) in the interior of Ω. At this point ∂ α u = ∂ β u = 0. Thus ∑ i A ∂ u i i ≡ 0 at ( α , β 0 0 ). Therefore (referring to (3.2) of A3)

$$( 3. 1 8 ) \quad & 0 = \sum _ { i, j } A _ { i j } \partial _ { i } \partial _ { j } u = \text{tr} ( \mathcal { A } ( \alpha _ { 0 }, \beta _ { 0 } ) H ( \alpha _ { 0 }, \beta _ { 0 } ) ), \quad \text{where} \\ & \mathcal { A } = ( A _ { i j } ) \text{ and} \\ & H = ( \partial _ { i } \partial _ { j } u ) \text{ is the Hessian matrix.}$$

By our strong ellipticity assumption (A3), det A &gt; 0. If α , β 0 0 is a strict local maximum then at this point det( H ) &gt; 0. Therefore det A H &gt; 0 as well. Both A and H are symmetric. By the trace identity, tr( A H ) = tr ( A 1 / 2 H A 1 / 2 ) , so tr( A H &gt; ) 0, contradicting tr( A H ) = 0. Therefore u cannot achieve a strict local maximum (or, by a parallel argument, a strict local minimum) anywhere in the interior of Ω, so u must be constant.

Rather than proving the existence of the function T α,β ( ) directly, we prove the existence of a function U α,β ( ) satisfying a related inhomogeneous PDE with periodic boundary conditions, and obtain T from the equation

$$T ( \alpha, \beta ) = U ( \alpha, \beta ) - \frac { \overline { T } } { 2 \pi } \alpha.$$

Proof of Theorem 3.1. Assuming A1-A5 hold, consider the inhomogeneous prob-

lem

$$\mathcal { P } U = - 1 - f _ { 1 } ( \alpha, \beta ) \frac { \overline { T } } { 2 \pi }, \quad \text{on } \Omega \\ ( 3. 2 0 ) \quad \sum _ { k = 1 } ^ { 2 } \mathcal { G } _ { 2 k } \partial _ { k } U ( \alpha, \pm 1 ) = \pm \mathcal { G } _ { 2 1 } ( \alpha, \pm 1 ) \frac { \overline { T } } { 2 \pi }, \quad \forall \alpha \in [ 0, 2 \pi ] \\ U ( 0, \beta ) - U ( 2 \pi, \beta ) = 0, \quad \forall \beta \in [ - 1, 1 ] \\ \dots \colon \colon \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \end{array}$$

which is equivalent to the system (3.15) (the object of Theorem 3.3) with inhomogeneities h = - -1 f T/ π 1 2 and g N = ±G 21 T/ π 2 in the interior and on the reflecting boundary of Ω, respectively. The ± in the inhomogeneous boundary condition reflects the direction of the outward unit normal at the upper boundary, ν = (0 1) at , β = 1, versus the lower boundary, ν = (0 , -1) at β = -1.

To see that Theorem 3.3 applies, note that

- 1. Ω = [0 2 , π ] × -[ 1 1] is a bounded Lipschitz domain, by construction; ,
- 2. the condition h ∈ ˜ H -1 (Ω) m holds since f 1 is assumed to be a C 2 function on a compact domain Ω;
- 3. the condition g D ∈ H 1 / 2 (Γ D ) m is satisfied trivially since Γ D (the subset of the boundary on which a Dirichlet condition is imposed) is the empty set;
- 4. the condition g N ∈ H -1 / 2 (Γ N ) m holds because G 21 , is assumed to be a C 2 function on the compact domain Ω, a fortiori on its boundary 0 ≤ α ≤ 2 π , β = ± 1; and
- 5. P is coercive, by Assumption A3 and Theorem 3.2.

Thus we may apply Theorem 3.3 to the system (3.20). By Lemma 3.4, all solutions of the homogeneous problem (3.14) are constant. On the other hand, any function u = const satisfies (3.14). Hence dim W = p = 1, and case (ii) of Theorem 3.3 applies. Therefore, the inhomogeneous problem (3.20) has solutions if and only if

$$0 & = ( v _ { 1 }, h ) _ { \Omega } + ( \gamma v _ { j }, g _ { N } ) _ { \Gamma _ { N } },$$

where v 1 is any solution of the adjoint homogeneous problem (3.9) with boundary conditions (3.7)-(3.8). By Assumption A4, all solutions of the adjoint homogeneous problem are multiples of the steady state density ρ ss , so the solvability condition (3.21) is equivalent to

$$0 & = \int _ { \Omega } d \alpha \, d \beta \, \rho _ { \text{ss} } ( \alpha, \beta ) \left ( - 1 + \frac { \overline { T } } { 2 \pi } f _ { 1 } ( \alpha, \beta ) \right ) \\ & \quad - \frac { 1 } { 2 } \int _ { \alpha = 0 } ^ { 2 \pi } d \alpha \, \left [ \rho _ { \text{ss} } ( \alpha, 1 ) \mathcal { G } _ { 2 1 } ( \alpha, 1 ) \frac { \overline { T } } { 2 \pi } - \rho _ { \text{ss} } ( \alpha, - 1 ) \mathcal { G } _ { 2 1 } ( \alpha, - 1 ) \frac { \overline { T } } { 2 \pi } \right ]. \\ \infty \, \cdots \, \cdots \, \cdots \, \cdots \, \alpha \, \quad \, \lambda \, \cdots \, \cdots \, \cdots \, \alpha \, \quad \, \lambda \, \cdots \, \cdots.$$

But the steady-state probability flux (3.10) satisfies

$$\text{But the steady-state probability flux } ( 3. 1 0 ) \text{ satisfies} \\ ( 3. 2 3 ) \quad \frac { 2 \pi } { \overline { T } } = & \int _ { \Omega } d \alpha \, d \beta \, J _ { \text{ss}, 1 } \\ = & \int _ { \Omega } d \alpha \, d \beta \left [ f _ { 1 } \rho _ { \text{ss} } - \frac { 1 } { 2 } \left ( \frac { \partial } { \partial \alpha } \left ( \mathcal { G } _ { 1 1 } \rho _ { \text{ss} } \right ) + \frac { \partial } { \partial \beta } \left ( \mathcal { G } _ { 1 2 } \rho _ { \text{ss} } \right ) \right ) \right ] \\ ( 3. 2 4 ) = & \mathbb { E } [ f _ { 1 } ] - \frac { 1 } { 2 } \int _ { \beta = - 1 } ^ { 1 } d \beta \, \left ( \mathcal { G } _ { 1 1 } ( 2 \pi, \beta ) \rho _ { \text{ss} } ( 2 \pi, \beta ) - \mathcal { G } _ { 1 1 } ( 0, \beta ) \rho _ { \text{ss} } ( 0, \beta ) \right ) \\ ( 3. 2 5 ) \quad & - \frac { 1 } { 2 } \int _ { \alpha = 0 } ^ { 2 \pi } d \alpha \, \left ( \mathcal { G } _ { 1 2 } ( \alpha, 1 ) \rho _ { \text{ss} } ( \alpha, 1 ) - \mathcal { G } _ { 1 2 } ( \alpha, - 1 ) \rho _ { \text{ss} } ( \alpha, - 1 ) \right ).$$

$$( 3. 2 5 ) \quad \ \ & - \frac { 1 } { 2 } \int _ { \alpha = 0 } ^ { 2 \pi } d \alpha \left ( \mathcal { G } _ { 1 2 } ( \alpha, 1 ) \rho _ { s s } ( \alpha, 1 ) - \mathcal { G } _ { 1 2 } ( \alpha, - 1 ) \rho _ { s s } ( \alpha, - 1 ) \right ).$$

The first and last terms in (3.24)-(3.25) match corresponding terms in (3.22), and the middle term in (3.24) vanishes because both G and ρ ss are periodic in α . Therefore we may rewrite (3.22) in terms of the steady-state flux:

$$\overline { T } = \frac { 2 \pi } { \int _ { \Omega } d \alpha \, d \beta \, J _ { s s, 1 } }.$$

Thus the mean period T emerges as a parameter required to have a specific value in order that the solvability condition should be satisfied. But this condition is already given in Assumption A5. (Compare Remark 2 on page 10.)

We conclude that the solutions of the inhomogeneous problem for U (3.20) have the form U α,β ( ) + c for arbitrary c ∈ R . Having established the existence of the function U , it remains to show that the function T defined by (3.19) solves the inhomogeneous problem (3.13) with jump-periodic boundary conditions. From (3.2) it is clear that P ( c ) = 0 and P ( α ) = f 1 ( α, β ). By inspection of (3.20), it is evident that P ( T ) = -1, as required. The adjoint reflecting boundary conditions are satisfied, since clearly

$$\left [ \mathcal { G } _ { 2 1 } \frac { \partial } { \partial \alpha } + \mathcal { G } _ { 2 2 } \frac { \partial } { \partial \beta } \right ] ( \alpha ) = \mathcal { G } _ { 2 1 }.$$

The jump-periodic boundary condition on the local domain, T (0 , β ) = T (2 π, β ) + T , holds by construction. On the extended domain T α,β ( ) = T α ( +2 π, β ) + T for all ( α, β ) ∈ Ω ext .

Finally, T inherits uniqueness from U , up to an additive constant.

Remark: The properties of strong ellipticity and coercivity are preserved under C 1 diffeomorphisms ([19], p. 156). If we begin with a star-shaped domain and make the transformation described in § 2.2 we thus establish the existence of a unique T -function (up to an additive constant) the level curves of which satisfy the MRT property postulated by S&amp;P.

For an extension of theorem 3.1 to oscillators in n &gt; 2 dimensions, see § 5.

- 4. Detailed examples. We illustrate the solution of the backward equation and the resulting MRT isochrons for several example systems (one additional example may be found in supplemental section § 9). A detailed discussion of the numerical method may be found in supplemental section § 8.
- 4.1. An analytically solvable test case -the isotropic noisy StuartLandau oscillator. We consider a solvable case with a trivial solution for the phase, here simply given by the angle of the common polar coordinates (the simple geometric phase). The system of stochastic differential equations is in this case given in polar coordinates by

$$( 4. 1 ) \quad \ \dot { \theta } = \omega + \sqrt { 2 D } \xi _ { 1 } ( t ), \ \dot { r } = - \gamma r ( r ^ { 2 } - 1 ) + \sqrt { 2 D } \xi _ { 2 } ( t ).$$

Here ξ 1 2 , ( ) are independent Gaussian white noise sources with t 〈 ξ i ( ) t ξ j ( t ′ ) 〉 = δ i,j δ t ( -t ′ ). With γ &gt; 0 the noiseless system has a stable fixed point in the radius variable at r = 1; with noise the radius fluctuates around this value and there are also fluctuations in the phase velocity (see Fig. 4.1a for a trajectory plotted in cartesian coordinates) but there is no coupling between phase and amplitude. We solve the MRT equation numerically for a bounded annular region cutting out an inner circle and restricting the motion to a disk with a maximum radius.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-44c298d5-a41e-449b-989b-f456324aceb1.jpg)

x

Fig. 4.1: Isochronal Stuart-Landau oscillator clock model, eq. (4.1), with ω = 1, γ = 1, and D = 0 01. . Sample trajectory (with timestep ∆ t = 0 001) with initial . condition at (0 , 0) (a), MRT phase Θ( x, y ) obtained by a finite difference scheme (b) and its contour lines (c), which are the MRT isochrons. As expected, the isochrons are spokes-of-the-wheel and each 'spoke' is constant up to the order of 10 -12 . For the parameters, used T = 2 π .

If the drift term of the angular variable (in polar coordinates) does not depend on the radius, it is easy to show that the MRT function T θ, r ( ) is linear in θ , i.e.

$$( 4. 2 ) \quad \ \ T ( \theta, r ) = T _ { 0 } - \frac { \overline { T } } { 2 \pi } \theta = T _ { 0 } - \frac { \theta } { \omega }, \quad \text{therefore } \Theta ( \theta, r ) = \theta + \Theta _ { 0 }$$

(where T 0 and Θ 0 are arbitrary constants) is a solution of the PDE that satisfies the periodic-plus-jump condition: the radial part of the backward operator yields immediately zero, the linear increase in θ results in the inhomogeneity on the right hand side. The contour lines of this function, and hence also the MRT phase function Θ( x, y ), are simply the spokes of a wheel and this is also returned by our numerical routine, which solves the PDE numerically (see Fig. 4.1, middle and right). These would be also the isochrons for the deterministic oscillator. Our numerical solution, obtained for finite inner and outer boundaries, reproduces this trivial analytical solution, which is a first indication of the robustness of the solution method.

4.2. Stuart-Landau oscillator with antirotating phase. The next example is due to Newby and Schwemmer [21] and displays an interesting amplitude-phase coupling:

$$( 4. 3 ) \quad \dot { \theta } = \omega + \omega \gamma c ( 1 - r ) ^ { 2 } + \sqrt { 2 D } \xi _ { 1 } ( t ), \quad \dot { r } = - \gamma r ( r ^ { 2 } - 1 ) + \sqrt { 2 D } \xi _ { 2 } ( t )$$

With parameters γ &gt; 0 and c &lt; 0, the radius-dependent drift in the θ dynamics will cause an anti-rotation compared to the movement on the limit cycle at r = 1 whenever the radius sufficiently deviates from the limit cycle (both for smaller and larger values). The exact condition for this change in the deterministic force field (cf. also Fig. 4.2a) is (1 -r ) 2 &gt; ( -γc ) -1 . The shape of the isochrons reflects the reversed deterministic velocity away from the limit cycle. Both positions inside and outside the limit cycle have mean angular velocity larger than, and in the same direction as, the mean angular rotation, while points near the limit cycle have mean angular velocity that is smaller and in the reverse direction. Consequently, the isochrons, resulting as contour lines from the PDE solution, attain a hook-like shape similar to what is observed for the deterministic system. As the noise becomes stronger, we expect that

Fig. 4.2: Stuart-Landau oscillator with antirotating phase, eq. (4.3), with ω = 1, γ = 15, c = -15, and D = 0 198. . Sample trajectory (∆ t = 10 -4 ) for an initial condition at (1 , 0), limit cycle (grey bold line) and vector field shown by arrows (a), MRTphase Θ( x, y ) from the finite-difference-scheme solution (b) and its contour lines, i.e. the isochrons of the MRT phase (c). For the parameters used, T = 11 89. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-d78fb487-1d91-4064-9517-07073a20f707.jpg)

the isochrons of the MRT phase become more straight because the trajectory diffuses faster in the radial direction, which makes all the initial points for the race around the circle more equal. Put differently, with more noise in the system, less head start is required for extreme radial positions to complete one round in the same mean time.

4.3. Stuart-Landau oscillator with y -polarized noise. The previous examples shared a basic rotational symmetry: the dynamics in polar coordinates was independent of the geometric angle θ . As a consequence, all the isochrons can be mapped onto each other by a simple rotation.

Fig. 4.3: Noisy Stuart-Landau oscillator with y -polarized noise given by eq. (4.4) with ω = 1 99, . κ = 1, and σ = 0 2. . Example trajectory (a) (∆ t = 0 01, trajectories gen-. erated in polar ( r, θ ) coordinates and then converted to Cartesian ( x, y ) for analysis), MRT phase Θ( x, y ) (b) and its contour lines, i.e. isochrons of the MRT phase (c). For the parameters used, T = 20 93. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-080f59bb-e5f9-402c-ad86-044fb7c5e295.jpg)

We consider now a case that lacks this symmetry, the Stuart-Landau oscillator

with y -polarized noise 4 :

$$\text{with $y$-polarized noise$} \colon & \quad \begin{matrix} \dot { \theta } = f _ { 1 } ( \theta, r ) + \sqrt { 2 D } \sin ( \theta ) \xi ( t ) \\ \dot { r } = f _ { 2 } ( \theta, r ) + \sqrt { 2 D } r \cos ( \theta ) \xi ( t ) \\ f _ { 1 } ( \theta, r ) = \omega + r \cos \theta - \kappa r ^ { 2 } + D \cos \theta \sin \theta \\ f _ { 2 } ( \theta, r ) = r ( 1 - r ^ { 2 } ) + r D \left ( \cos ^ { 2 } \theta - \sin ^ { 2 } \theta \right ) \\ \quad \end{matrix} \quad \begin{matrix} \dot { \theta } = f _ { 1 } ( \theta, r ) + \sqrt { 2 D } \sin ( \theta ) \xi ( t ) \\ \quad \end{matrix} \quad \begin{matrix} \dot { \theta } = f _ { 2 } ( \theta, r ) + \sqrt { 2 D } r \cos ( \theta ) \xi ( t ) \\ \quad \end{matrix} \quad \begin{matrix} \dot { \theta } = f _ { 2 } ( \theta, r ) + \sqrt { 2 D } \sin ( \theta ) \xi ( t ) \\ \quad \end{matrix} \quad \begin{matrix} \dot { \theta } = f _ { 2 } ( \theta, r ) + \sqrt { 2 D } \sin ( \theta ) \xi ( t ) \\ \quad \end{matrix} \quad \begin{matrix} \dot { \theta } = f _ { 2 } ( \theta, r ) + \sqrt { 2 D } \sin ( \theta ) \xi ( t ) \\ \quad \end{matrix} \quad \begin{matrix} \dot { \theta } = f _ { 2 } ( \theta, r ) + \sqrt { 2 D } \sin ( \theta ) \xi ( t ) \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddot { \phi } = f _ { 1 } ( \theta, r ) + \sqrt { 2 D } \sin ( \theta ) \xi ( t ) \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddkappa \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots } \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \end{matrix} \quad \begin{matrix} \ddots \\ \ \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \ \end{matrix} \quad \begin{matrix} \ddots \\ \quad \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \begin{matrix} \ddots \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \quad \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ { \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \. \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ } \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ { \ \. \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \. \ \. \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Note that the θ drift is now also modified such that the angular velocity is sped up for θ ≈ π but slowed down around θ ≈ 0. In addition (4.4) is an example of an excitable system: in the absence of noise the dynamics have a stable fixed point rather than a limit cycle, and a circulating current appears only in the presence of noise.

Note also that the same noise sample ξ t ( ) drives both the θ and r coordinate. Thus in contrast to (4.1) and (4.3) the noise is singular, i.e. the nondegeneracy condition det G glyph[negationslash] = 0 (eq. (3.5)) is violated. The algorithm nevertheless produces the MRT phase (Fig. 4.3).

The MRT phase resulting from the PDE solution displays a rather heterogeneous structure. First of all, there is a clustering of isochrons on the left, where the angular velocity is systematically increased as explained above. Secondly, there is little noise in the θ variable when the trajectory is close to the x axis; at the same time, there is little noise in the radial dynamics close to the y-axis, because of the sine and cosine prefactors of the noise terms.

Consequently, the deviation of the isochrons from the spokes-of-a-wheel shape changes significantly: on the right hand side, isochrons cover a larger range of the angular variable than on the left hand side, resulting from an interplay between the angle dependence of phase velocity and noise level.

We take this most involved example also as a test case to check whether the algorithmic definition of the MRT phase is matched by our PDE result (Fig. 4.4). We pick out three contour lines of the PDE solution (Fig. 4.4a-c), distribute a number of initial points along this line (cyan dots), and start an ensemble of trajectories on each of the points (Fig. 4.4b). We simulate one rotation until the respective trajectory hits the back of the isochron; the corresponding mean return times are plotted in Fig. 4.4d-f and reveal a very good agreement for all initial points - the picked contour lines satisfy S&amp;P's algorithmic definition of the MRT property. For comparison we also show the same measurement in each case for another line, the spoke of a wheel passing through the same part of the noisy limit cycle (red lines in Fig. 4.4a-c). Points on these lines have different mean return times to the same line, hence, the spoke is not an isochron of the MRT phase in any of the cases considered.

## 4.4. Noisy heteroclinic oscillator. Our last example is given by the stochastic differential equations

$$\begin{array} { c c c c c } \mu & \mu & \mu & \mu & \mu & \mu & \mu \\ & \dot { x } = \cos ( x ) \sin ( y ) + \alpha \sin ( 2 x ) + \sqrt { 2 D } \xi _ { 1 } ( t ), \\ & \dot { y } = - \sin ( x ) \cos ( y ) + \alpha \sin ( 2 y ) + \sqrt { 2 D } \xi _ { 2 } ( t ). \end{array}$$

This model has been used in Ref. [29] to study another notion of phase, the asymptotic phase of a stochastic oscillator, and is an example of a noisy heteroclinic network ([1, 3, 26, 27], see also [14]). The model's deterministic dynamics ( D = 0) is characterized

4 S&amp;P study this system using the equivalent Stratonovich formulation, compare equation (5) of [24], ˙ θ = ω + r cos θ -κr 2 + √ 2 D sin( θ ) ◦ ξ t ( ), ˙ r = (1 r -r 2 ) + √ 2 Dr cos( θ ) ◦ ξ t ( ).

Fig. 4.4: Testing the MRT property for the noisy Stuart-Landau oscillator with y -polarized noise (4.4). Sample trajectory (a-c, black trace) and parameters are as in Fig. 4.3. Cyan curves (a-c) show isochronal sections obtained from the PDE solution; red curves are 'spoke' sections passing through the isochronal section. For each of 20 points per curve, equally spaced in the radial coordinate, we generate approximately 600,000 independent trajectories, using the Euler-Maruyama method (trajectories generated in polar ( r, θ ) coordinates with ∆ t = 0 001) and computing . the MRT from a starting point to the same section after completing one rotation (d-f). For initial conditions on the PDE-derived isochrons, mean return times (d-f, cyan dots) are within 0.5% of the mean period (d-f, black horizontal line). For initial conditions on a spoke, MRTs deviate significantly from the mean period (d-f, red dots). Thus the isochronal sections obtained as level curves of T x, y ( ) solved with the periodic-plus-jump condition, is a MRT isochron in the sense of S&amp;P, whereas the spoke shows a strong increase of the mean return time with growing radius and is thus not such an isochron.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-de967429-6636-4a92-89ac-0c3c0efa736f.jpg)

by four saddle points in the corners ( ± π/ , 2 ± π/ 2) and an unstable focus at (0 0). , Without noise the system approaches the heteroclinic cycle connecting the saddles and the system does not sustain a finite period oscillation (the heteroclinic cycle has an infinite period). With a weak noise, the system displays pronounced oscillations of a finite phase coherence [29] and each of the components has pronounced peaks in its power spectrum [11]. In the plane these oscillations become manifest by a clockwise rotation around the origin (Fig. 4.5a).

Here we solve eq. (2.34) on the domain ( -π/ , π/ 2 2) × -( π/ , π/ 2 2) and cut out a small square in the middle. On these outer and inner boundaries reflecting boundary conditions are imposed while along the blue cut line in Fig. 8.1a, we apply the jump condition. The resulting mean return time (as a function of the starting position ( x, y )) is shown in Fig. 4.5, middle and the isoclines, which are the isochrons of the MRT phase, are displayed in Fig. 4.5, right.

The MRT isochrons wind inward around the origin in a counterclockwise direction.

Fig. 4.5: Noisy heteroclinic oscillator, eq. (4.5), with α = 0 1, . and D = 0 01125. . Sample trajectory (∆ t = 0 001) with initial condition at (0 . , 0) that moves clockwise (a). In the absence of noise, the system has a stable heteroclinic orbit with infinite period. With noise ( D &gt; 0), the small perturbations eventually knock the trajectory out of the corners to form a stochastic oscillator with finite mean period. MRT phase Θ( x, y ) finite-difference-scheme solution (b) and its contour lines, i.e. the isochrons of the MRT phase (c). The computational grid is 251 × 251 and the punctured square has side length of 0.05. For the parameters used, T = 16 23. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-f453354f-9d8c-45cd-a761-e698abaccf33.jpg)

In terms of the mean return time this geometry has a simple interpretation: we have to start at an advanced position if we start at the slow track on the outside (close to the square's sides), where the speed is reduced, compared to a starting position on the inside (close to the central square), where the speed is larger.

Fig. 4.6: (a) Heteroclinic trajectory in black for D = 0 01125 and . α = 0 1 with initial . condition at origin. The cyan curve is an isophase pulled from the contour plot in Fig. 4.5c. The cyan dots indicate 20 equally spaced initial conditions along this curve. The red curve is a 'spoke' running through the tenth isophase initial condition. The red dots indicate 20 equally spaced initial conditions. (b) The black line is the mean period calculated from the stationary flux (used in the finite difference scheme). Cyan dots are the MRTs back to the isophase after completing one full oscillation of the isophase initial conditions, from an ensemble of approximately 520,000 trajectories generated using the Euler-Maruyama method with ∆ t = 0 001. . The x -axis is the radial distance from the origin. Red dots are the MRTs back to the spoke after completing one full rotation from the spoke initial conditions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-c9198028-6bc0-48c1-a6f1-c0e34dfa8693.jpg)

As in the case of the Stuart-Landau oscillator with y -polarized noise, we can test

that the level curves of the solution of eq. (2.34) have the MRT property. Fig. 4.6 contrasts the mean return times for ensembles of trajectories with initial conditions located either along a level curve of Θ, the MRT phase, vs. along a simple radial section. Using an ensemble of independent trajectories for each initial condition, the mean return times for the spoke vary by a factor of two as the starting distance from the origin ranges from 0.5 to 1.5. In contrast, the MRTs for the isophase fall within 0 7% of the mean period regardless of starting radius, illustrating S&amp;P's MRT . criterion for an isochron.

- 5. Extension to n -dimensional oscillators. For continuity of presentation and ease of illustration, we gave our main theorem 3.1 for an oscillator in two dimensions. But the proof based on McLean's Theorem 4.10 carries over to oscillators in arbitrary finite dimension with little modification, so it is worth stating the general case here. As before, we assume that our noisy oscillator has been mapped through a coordinate transformation to a system with one periodic coordinate α ∈ S 1 ≡ [0 , 2 π ] and the remaining n -1 coordinates lying within a simply connected compact subset of R n -1 , with piecewise smooth boundary; call this set B n -1 1 ⊂ R n -1 . Thus the region is assumed to have the form of an ( n -1)-dimensional cylinder embedded in R n . We define the local domain Ω = S 1 × B n -1 1 and the extended domain Ω ext = R × B n -1 1 . The driving noise vector dW will have K ≥ n components dW ,... , dW 1 K . Recapitulating assumptions A1-A5 in this context, we assume the following:
- A1'. Transformed into ( α, glyph[vector] β ) coordinates, the trajectories ( α t , glyph[vector] ( ) β t ( )) of a strongly Markovian time-homogeneous process obey an Ito equation

$$d \alpha & = f _ { 1 } ( \alpha, \vec { \beta } ) \, d t + \sum _ { k = 1 } ^ { K } g _ { 1 k } ( \alpha, \vec { \beta } ) \, d W _ { k } ( t ) \\ d \beta _ { i } & = f _ { i } ( \alpha, \vec { \beta } ) \, d t + \sum _ { k = 1 } ^ { K } g _ { i k } ( \alpha, \vec { \beta } ) \, d W _ { k } ( t ), \text{ for } 2 \leq i \leq n, \\ \dots \dots \, \tau \quad - \dots \, \varkappa ^ { 2 } \, \dots \, \varkappa \, \colon \, \dots \, \tau \, \colon \, \tau \, \colon$$

where f , g i ik are C 2 on Ω ext for 1 ≤ i ≤ n and 1 ≤ k ≤ K . As before, we will refer to the n × n matrix G = gg ᵀ .

A2'. The functions f i and g ik are periodic in the first coordinate with period 2 π , i.e. ∀ α ∈ R and i = 1 2, , f i ( α +2 π, glyph[vector] β ) = f i ( α, glyph[vector] β ), and likewise for each g ik .

A3'. The second order differential operator P is strongly elliptic, where P is defined

$$\mathcal { P } u = - \sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } \partial _ { i } ( A _ { i j } \partial _ { j } u ) + \sum _ { i = 1 } ^ { n } A _ { i } \partial _ { i } u \text{ on } \Omega _ { \text{ext} }$$

and A ij = - G 1 2 ij and A i = -1 2 ∑ 2 j =1 ∂ j G ij + f i . At the boundary | glyph[vector] β | = 1 we impose Neumann boundary conditions

$$0 = \sum _ { i = 1 } ^ { 2 } \nu _ { i } \sum _ { j = 1 } ^ { 2 } \mathcal { G } _ { i j } \partial _ { j } u$$

where ν is the outward unit normal at the boundary.

A4'. The process viewed on Ω (taking α mod 2 π ) admits a density ρ α, glyph[vector] ( β, t ) evolving according to

$$( 5. 4 ) \quad \ \frac { \partial \rho } { \partial t } = \mathcal { L } \rho = \mathcal { P } ^ { * } \rho = - \sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } \partial _ { i } ( A _ { j i } ^ { * } \partial _ { j } \rho ) - \sum _ { i = 1 } ^ { 2 } \partial _ { i } ( A _ { i } ^ { * } \rho )$$

where A ∗ ji = A ij and A ∗ i = A i , and we impose reflecting (Neumann) boundary conditions

$$0 = \sum _ { j = 1 } ^ { 2 } A _ { j i } ^ { * } \partial _ { j } \rho + A _ { i } ^ { * } \rho$$

at | glyph[vector] β | = 1 (for all α ), and periodic boundaries in α , i.e.

$$( 5. 6 ) & & \forall \vec { \beta } \in B _ { 1 } ^ { n - 1 }, \quad \rho ( 0, \vec { \beta } ) = \rho ( 2 \pi, \vec { \beta } ).$$

We assume that the system has a unique stationary distribution ρ ss ≥ 0, with 1 = ∫ Ω dαd glyph[vector] βρ ss ( α, glyph[vector] β ), i.e. satisfying the homogeneous equation

$$\mathcal { P } ^ { * } \rho _ { s s } = 0,$$

together with the boundary conditions (3.7)-(3.8).

As before, the stationary flux vector J ss ( α, glyph[vector] β ) corresponds componentwise to

$$J _ { \text{ss}, i } ( \alpha, \vec { \beta } ) = - \frac { 1 } { 2 } \sum _ { j = 1 } ^ { n } \partial _ { j } \left ( \mathcal { G } _ { i j } \rho _ { \text{ss} } \right ) + f _ { i } \rho _ { \text{ss} }.$$

- A5'. We assume the mean drift is nonzero and oriented in the direction of increasing α . That is, if γ : B n -1 1 → [0 , 2 π ] is any C 1 function whose graph C γ = { ( α = γ β , glyph[vector] ( glyph[vector] ) β ) : glyph[vector] β ∈ B n -1 1 } cuts transversely through the extended domain, separating Ω ext into left and right connected components, with unit normal n ( glyph[vector] β ) oriented into the downstream connected component, then the mean rightward flux through C γ is positive, i.e.

$$0 < \overline { J } \coloneqq \int _ { - 1 } ^ { 1 } d \vec { \beta } \mathbf n ^ { \intercal } ( \vec { \beta } ) J _ { \text{ss} } ( \gamma ( \vec { \beta } ), \vec { \beta } ).$$

These assumptions suffice to establish the existence and uniqueness (up to an additive constant) of a MRT function T α, glyph[vector] ( β ) satisfying L † [ T ] = -1 subject to adjointreflecting boundary conditions ∑ n i =1 ν i ( α, glyph[vector] β ) G ij ∂ T j = 0 at the boundary ∂ Ω, and jump-periodic boundary conditions ∀ glyph[vector] β ∈ B n -1 1 , T α, glyph[vector] ( β ) = T α ( +2 π, glyph[vector] β ) + T , where T is the mean period of the oscillator, T = ( J ) -1 . As in the planar case, the proof would involve the auxiliary function U α, glyph[vector] ( β ) satisfying

$$\mathcal { P } U = - 1 - f _ { 1 } ( \alpha, \vec { \beta } ) \frac { \overline { T } } { 2 \pi }, \quad \text{on } \Omega \\ ( 5. 1 0 ) \quad \sum _ { i = 2 } ^ { n } \nu _ { i } \sum _ { j = 1 } ^ { n } \mathcal { G } _ { i j } \partial _ { j } U ( \alpha, \pm 1 ) = \mathcal { G } _ { i 1 } ( \alpha, \vec { \beta } ) \frac { \overline { T } } { 2 \pi }, \quad \forall \alpha \in [ 0, 2 \pi ] \text{ and } \vec { \beta } \in \partial B _ { 1 } ^ { n - 1 } \\ U ( 0, \vec { \beta } ) - U ( 2 \pi, \vec { \beta } ) = 0, \quad \forall \vec { \beta } \in B _ { 1 } ^ { n - 1 }, \\.. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \end{cases}$$

with ν the outward unit normal vector at the boundary of Ω.

- 6. Summary and conclusions. In this paper we have found an analytic way to define the mean-return-time (MRT) phase, originally proposed by Schwabedal and Pikovsky in terms of an algorithm, for the important class of smooth two-dimensional

stochastic oscillators that are driven by white Gaussian noise. We showed that the defining isochrons are given as the contour lines of the solution of the conventional PDE for the mean-first-passage-time function but with an uncommon periodic-plusjump condition. We illustrated this construction in a number of stochastic oscillator models and verified the algorithmic MRT-property of our PDE solution for the most involved examples (Stuart-Landau oscillator with y-polarized noise and heteroclinic oscillator). Some open questions remain.

It would be of interest to study the effect of increasing levels of noise on the shape of the isochrons. Preliminary results indicate that in many cases the shape of the isochrons change towards spokes of a wheel with increasing noise level, i.e. the isochrons become less curved. A similar phenomenon was observed in a version of the asymptotic phase for stochastic oscillators that we introduced previously, based on a spectral decomposition of the generator [29]. However, it is unclear what are the exact conditions for the stochastic oscillator ( i.e. for the functions appearing in the Langevin eq. (2.1)) leading to isochrons of greater or lesser curvature. Also, in cases in which no deterministic phase exists as for the heteroclinic oscillator, one should also consider the opposite limit and let the noise level shrink.

Finally, one should compare systematically to the asymptotic phase proposed in [29]. For the trivial case of the isotropic Stuart-Landau oscillator that lacks an angleamplitude coupling, both definitions yield the same stochastic phase which is the same as in the deterministic case [30]. Generally, in cases with an existing deterministic phase, both definitions of a stochastic phase yield this deterministic phase in the limit of vanishing noise and hence can be regarded as possible generalizations of the deterministic phase to the stochastic case. For a finite noise level, however, there seem to be small differences between the asymptotic phase and the MRT phase in most systems. The exact nature of these differences and their role in reduced descriptions of stochastic oscillators remain exciting topics of future research.

7. Acknowledgments. The authors thank J. Dobrosotskaya, W. McLean, and D. Gurarie for helpful discussions. PT acknowledges support from NSF grant DMS1413770, from Oberlin College, and from the Mathematical Biosciences Institute and the NSF under grant DMS-1440386. Large-scale Monte Carlo simulations made use of the High Performance Computing Resource in the Core Facility for Advanced Research Computing at Case Western Reserve University.

## REFERENCES

- [1] D. Armbruster, E. Stone, and V. Kirk , Noisy heteroclinic networks , Chaos, 13 (2003), pp. 71-79.
- [2] J. Ashmore, P. Avan, W. Brownell, P. Dallos, K. Dierkes, R. Fettiplace, K. Grosh, C. Hackney, A. Hudspeth, F. J¨licher, u et al. , The remarkable cochlear amplifier , Hearing research, 266 (2010), pp. 1-17.
- [3] Y. Bakhtin , Noisy heteroclinic networks , Probability Theory and Related Fields, (2010), pp. 142, http://dx.doi.org/10.1007/s00440-010-0264-0. 10.1007/s00440-010-0264-0.
- [4] D. Burke and A. M. de Paor , A stochastic limit cycle oscillator model of the eeg , Biol. Cybern., 91 (2004), p. 221.
- [5] K. Dierkes, F. J¨licher, and B. Lindner u , A mean-field approach to elastically coupled hair bundles , Eur. Phys. J. E, 35 (2012), p. 37.
- [6] B. Ermentrout , Type I membranes, phase resetting curves, and synchrony , Neural Comput., 8 (1996), p. 979.
- [7] G. Ermentrout and N. Kopell , Parabolic bursting in an excitable system coupled with a slow oscillation. , SIAM J. Appl. Math, 46 (1986), pp. 233-253.
- [8] G. B. Ermentrout and D. H. Terman , Foundations Of Mathematical Neuroscience , Springer,

- 2010. [9] R. Feistel and W. Ebeling , Deterministic and stochastic theory of sustained oscillations in autocatalytic reaction systems , Physica A: Statistical Mechanics and its Applications, 93 (1978), pp. 114-137.
- [10] C. W. Gardiner , Handbook of Stochastic Methods, 2nd Edition , Springer-Verlag, Berlin, 1997.
- [11] J. Giner-Baldo, P. Thomas, and B. Lindner , Power spectrum of a noisy system close to a heteroclinic orbit , J. Stat. Phys., 168 (2017), p. 447.
- [12] J. Guckenheimer , Isochrons and Phaseless Sets , J. Math. Biol., 1 (1975), p. 259.
- [13] J. Guckenheimer and A. Javeed , Locomotion: exploiting noise for state estimation , Biol. Cybern., (2018), p. 1.
- [14] M. W. Hirsch, S. Smale, and R. L. Devaney , Differential Equations, Dynamical Systems &amp; An Introduction to Chaos , Elsevier Academic Press, 2nd ed., 2004.
- [15] F. C. Hoppensteadt and E. M. Izhikevich , Weakly Connected Neural Networks , SpringerVerlag, New York, 1997.
- [16] R. Kawai, X. Sailer, L. Schimansky-Geier, and C. Van den Broeck , Macroscopic limit cycle via pure noise-induced phase transitions , Phys. Rev. E., 69 (2004), p. 051104.
- [17] B. Lindner, J. Garc´ ıa-Ojalvo, A. Neiman, and L. Schimansky-Geier , Effects of noise in excitable systems , Phys. Rep., 392 (2004), p. 321.
- [18] C.-A. Lugo and A. J. McKane , Quasicycles in a spatial predator-prey model , Phys. Rev. E, 78 (2008), p. 051911.
- [19] W. C. H. McLean , Strongly elliptic systems and boundary integral equations , Cambridge university press, 2000.
- [20] A. B. Neiman, K. Dierkes, B. Lindner, L. Han, and A. L. Shilnikov , Spontaneous voltage oscillations and response dynamics of a hodgkin-huxley type model of sensory hair cells , J. Math. Neurosci., 1 (2011), p. 11.
- [21] J. Newby and M. Schwemmer , Effects of Moderate Noise on a Limit Cycle Oscillator: Counterrotation and Bistability , Phys. Rev. Lett., 112 (2014), p. 114101.
- [22] A. Pikovsky , Comment on 'Asymptotic phase for stochastic oscillators' , Phys. Rev. Lett., 115 (2015), p. 069401.
- [23] O. V. Popovych, C. Hauptmann, and P. A. Tass , Control of neuronal synchrony by nonlinear delayed feedback , Biol. Cybern., 95 (2006), p. 69.
- [24] J. Schwabedal and A. Pikovsky , Phase Description of Stochastic Oscillations , Phys. Rev. Lett., 110 (2013), p. 204102.
- [25] K. M. Stiefel and G. B. Ermentrout , Neurons as oscillators , J. Neurophysiol., 116 (2016), p. 2950.
- [26] E. Stone and D. Armbruster , Noise and O(1) amplitude effects on heteroclinic cycles , Chaos, 9 (1999), p. 499.
- [27] E. Stone and P. Holmes , Random perturbations of heteroclinic attractors , SIAM J. Appl. Math, 50 (1990), pp. 726-743.
- [28] P. A. Tass , Phase resetting in medicine and biology: stochastic modelling and data analysis , Springer Science &amp; Business Media, 2007.
- [29] P. J. Thomas and B. Lindner , Asymptotic phase of stochastic oscillators , Phys. Rev. Lett., 113 (2014), p. 254101.
- [30] P. J. Thomas and B. Lindner , Thomas and Lindner reply: , Phys. Rev. Lett., 115 (2015), p. 069402.
- [31] P. J. Thomas and B. Lindner , Phase descriptions of a multidimensional Ornstein-Uhlenbeck process , Phys. Rev. E, 99 (2019), p. 062221.
- [32] A. M. van Mourik, A. Daffertshofer, and P. J. Beek , Deterministic and stochastic features of rhythmic human movement , Biological Cybern., 94 (2006), pp. 233-244.
- [33] E. Wallace, M. Benayoun, W. van Drongelen, and J. D. Cowan , Emergent Oscillations in Networks of Stochastic Spiking Neurons , Plos One, 6 (2011), p. e14804.
- [34] M. C. Wang and G. E. Uhlenbeck , On the Theory of the Brownian Motion-II , Rev. Mod. Phys., 17 (1945), p. 323.
- [35] S. Wieczorek , Stochastic bifurcation in noise-driven lasers and Hopf oscillators , Phys. Rev. E, 79 (2009), p. 036209.
- [36] A. Winfree , The Geometry of Biological Time , Springer-Verlag, 1980.

## SUPPLEMENTARY MATERIALS

8. Numerical Method for Obtaining the MRT Phase. In this section we describe in detail how we incorporate the jump condition into a finite difference numerical scheme to obtain the mean-return-time function T x, y ( ) for the heteroclinic oscillator example of § 4.4, and we briefly investigate the robustness of the isochrons to changes in the inner domain boundary size and specification of the mean period T .

## 8.1. Finite Difference Method. We solve equation (2.34):

$$\mathcal { L } ^ { \dagger } T ( x, y ) = - 1 \text{ with } \mathcal { R } _ { T } T ( x, y ) | _ { ( x, y ) \in R _ { \pm } } = 0, \lim _ { \varepsilon \to 0 ^ { + } } \left ( T ( - \varepsilon, y ) - T ( \varepsilon, y ) \right ) = \overline { T },$$

with adjoint reflecting boundary conditions, i.e.

$$\mathcal { R } _ { T } u = \sum _ { j = 1, 2 } n _ { j } \sum _ { k = 1, 2 } \mathcal { G } _ { j k } \partial _ { k } u = 0$$

along the inner and outer boundary of a square domain

$$\mathcal { D } = \left ( [ - \pi / 2, \pi / 2 ] \times [ - \pi / 2, \pi / 2 ] \right ) \, \langle \left ( [ - \epsilon / 2, \epsilon / 2 ] \times [ - \epsilon / 2, \epsilon / 2 ] \right ),$$

with a jump of T implemented along a vertical cut at x = 0 running from y = glyph[epsilon1]/ 2 to y = π/ 2, see Fig. 8.1A.

We employed a standard finite difference scheme to solve the MFPT-PDE with the mean period jump condition. The domain was uniformly discretized and the points within the punctured center removed. All derivatives were approximated using the 2nd order, centered finite differences at interior points of the domain (black grid points in Fig. 8.1A), excluding the cut line (blue points in Fig. 8.1A):

$$( 8. 1 ) \quad \ \ ( \partial _ { x } T ) _ { i, j } \to \frac { T _ { i + 1, j } - T _ { i - 1, j } } { 2 h } _ { \alpha \alpha \cdot \cdot \cdot \cdot }$$

$$\left ( \partial _ { x x } ^ { 2 } T \right ) _ { i, j } \rightarrow \frac { T _ { i + 1, j } \stackrel { \dots } { - 2 } T _ { i, j } + T _ { i - 1, j } } { h _ { \substack { n \\ n } } ^ { 2 } } \right ) _ { n }$$

$$( 8. 3 ) \quad \left ( \partial _ { x y } ^ { 2 } T \right ) _ { i, j } \rightarrow \frac { T _ { i + 1, j + 1 } \stackrel { \cdot \cdot } { - T _ { i + 1, j - 1 } } - T _ { i - 1, j + 1 } + T _ { i - 1, j - 1 } } { 4 h ^ { 2 } }$$

and similarly for the y derivatives.

At inner and outer boundary points (red points in Fig. 8.1A) the adjoint reflecting boundary conditions were implemented by introducing 'ghost points' one grid space beyond the edge of the domain. The normal vector ( n , n 1 2 ) ᵀ was taken to be ( ± 1 0) , ᵀ or (0 , ± 1) ᵀ respectively, as dictated by the geometry.

We implemented the mean period jump condition along a line segment connecting the outer and inner boundaries (blue dots in Fig. 8.1A). For the heteroclinic oscillator implementation, the cut was made along the vertical line segment from ( x = 0 , y = glyph[epsilon1]/ 2) to ( x = 0 , y = π/ 2). Denote the set of finite difference points on this line segment by writing ( i 0 , j ) ∈ glyph[lscript] . Then, for ( i 0 , j ) we implement the jump in T between the point ( i 0 , j ) and the immediately adjacent point to the left, which we denote ( i 0 +1 , j ). That is, we replace T i 0 +1 ,j with ( T + T i 0 +1 ,j ) . For these points we thus replace (8.1)-(8.2) with

$$\frac { T _ { i _ { 0 } + 1, j } - T _ { i _ { 0 } - 1, j } } { 2 h } \to \frac { ( \overline { T } + T _ { i _ { 0 } + 1, j } ) - T _ { i _ { 0 } - 1, j } } { 2 h } } { \sim \, \frac { 2 h } { \sim \, \frac { \sim \, \frac { \sim \, \sim } { \sim \, \frac { \sim \, \sim } { \sim \, \sim } } } }$$

$$( 8. 5 ) \quad \ \frac { T _ { i _ { 0 } + 1, j } - 2 T _ { i _ { 0 }, j } + T _ { i _ { 0 } - 1, j } } { h ^ { 2 } } = \frac { \binom { \overset { \dots } { T } + T _ { i _ { 0 } + 1, j } } { h ^ { 2 } } - 2 T _ { i _ { 0 }, j } + T _ { i _ { 0 } - 1, j } } { h ^ { 2 } }.$$

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

N

L

G

v

f

K

n

Q

d

V

J

p

9

8

y

o

&gt;

A

B

H

c

D

S

g

E

r

U

Y

O

7

T

0

X

M

F

2

k

m

z

w

u

W

P

I

j

q

/

5

+

R

C

Z

3

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

u

E

c

A

B

0

y

J

r

z

V

M

v

Y

7

8

X

k

&gt;

H

D

L

S

g

N

f

U

9

o

O

T

G

F

2

m

w

W

P

I

j

n

d

Q

q

/

5

+

R

p

C

K

Z

3

Fig. 8.1: Numerical scheme (heteroclinic oscillator example). (a) Finite difference grid for the heteroclinic oscillator. Interior points (black grid) use standard secondorder finite difference operators to approximate first and second derivatives. Adjoint reflecting boundary conditions are imposed at the inner and outer boundaries (red points). The size of the inner boundary is exaggerated for clarity. The jump condition is imposed along a vertical cut (blue points). The computational domain for Fig. 4.5 used a 251 × 251 grid with spacing ∆ x = ∆ y = π/ 250 and inner square side length glyph[epsilon1] = 2 π/ 125 ≈ 0 05. . (b) A 'time crystal', a la Winfree [36], providing a didactic illustration of the unwrapping of the domain to form a multi-valued surface, analogous to the unwrapped domain introduced in § 2.3 of the main text. Coordinates x and y correspond to ( x, y ) coordinates in panel A; coordinate z corresponds to the (multivalued) geometric phase z x, y ( ) = tan -1 ( y/x ). The jump condition eliminates the necessity of considering multiple copies of the domain.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-0364d065-ec47-4fb9-807b-940edc690a52.jpg)

A parallel procedure of subtracting T from T i 0 ,j is used for derivatives centered at location ( i 0 +1 , j ), whenever location ( i 0 , j ) ∈ glyph[lscript] .

The resulting system of equations is linear-affine, with all constants involving T written on the right-hand side.

As a technical point, the system is rank-deficient, leading to nonuniqueness of solutions. This point is readily addressed by specifying the value of T at a particular location ( i ∗ , j ∗ ). We prescribed the value of the northwest corner when solving for T for the heteroclinic oscillator:

$$T \left ( - \frac { \pi } { 2 }, \frac { \pi } { 2 } \right ) = T _ { 0, 0 } = 1 0 0$$

Fixing T at this point removes one unknown from the linear system, and the resulting reduced-dimensional system is uniquely solvable with standard linear algebraic methods.

The examples in § 4.1, § 4.2, § 4.3 and § 9 were implemented using standard polar coordinates rather than Cartesian coordinates. All examples in the main text used a 251 × 251 discretization of the computational domain. Code to reproduce the examples is available at https://github.com/pjthomas9/isophase .

In each of these examples, the annular computational domain contains the overwhelming bulk of the stationary probability distribution. As observed by S&amp;P, the definition of the MRT phase presumes that the oscillator trajectories avoid the region

near the 'center' of the domain. The inner boundary should be made small enough that the excised region at the center contains very little net probability. Similarly, if the underlying model has an unbounded domain, then the outer boundary should be large enough that the region excluded from the computational domain contains very little net probability.

## 8.2. Robustness of the Method.

8.2.1. Tolerance to the value of T . Although we did not systematically investigate parametric robustness of the method in detail, our heuristic observations suggest the method is robust against some variation in the T parameter. That is, we have observed that if the parameter representing T varies by up to 5% from the true value of T , the method nevertheless produces a reasonably smooth solution that is quantitatively close to the solution produced with the correct value of T .

- 8.2.2. Sensitivity to inner domain size. In order to test the sensitivity of the method for the heteroclinic oscillator system to changes in the size of the inner domain boundary, we ran the numerical algorithm using glyph[epsilon1] ∈ { 0 05 . , 0 1 . , 0 2 . } on a 201 × 201 finite difference map, with the other parameters the same as given in § 4.4. Fig. 8.2 compares the MRT function T x, y ( ) obtained when the interior boundary region [ -glyph[epsilon1]/ 2 , glyph[epsilon1]/ 2] × -[ glyph[epsilon1]/ 2 , glyph[epsilon1]/ 2] is imposed with different values of glyph[epsilon1] . Both the absolute difference in the value of T (panels A-B) and the shape of the isochrons (panels C-D) show modest sensitivity to glyph[epsilon1] restricted to a small (roughly O glyph[epsilon1] ( ) in size) neighborhood of the inner domain boundary. Outside a small neighborhood the values of T and the shape of the isochrons are insensitive to the precise value of glyph[epsilon1] dictating the inner boundary size.
- 9. A Further Example: Stuart-Landau oscillator with counterrotating phase. We discuss one additional nontrivial example which was suggested by Newby and Schwemmer, a Stuart-Landau oscillator with a true amplitude-phase dependence [21]. In polar coordinates, the model is given by

$$( 9. 1 ) \quad \dot { \theta } = \omega + \gamma c ( 1 - r ^ { 2 } ) + \sqrt { 2 D } \xi _ { 1 } ( t ), \quad \dot { r } = - \gamma r ( r ^ { 2 } - 1 ) + \sqrt { 2 D } \xi _ { 2 } ( t ).$$

On and inside the deterministic limit cycle ( r ≤ 1), the deterministic field still induces on average a counterclockwise (mathematically positive) rotation, while for points sufficiently outside the limit cycle ( r 2 &gt; 1 + ω/ γc ( )) the rotation is clockwise (cf. Fig. 9.1, left). The stochastic trajectory goes indeed both ways, depending on the value of the radius.

All the isochrons, here computed as the contour lines of the PDE solution, spiral inward, which can be understood as follows. Points with large radius have to get a head start because they go for a while in the wrong direction and need some time to turn (by a mean-driven or noise-induced reduction of the radius variable). Points on the inside of the circle, have a considerably higher rotation speed and have to be started with a certain delay, at an earlier geometric phase. Of course, these arguments also apply to the deterministic system and we see that noise sources with moderate intensity do not change this picture qualitatively.

- 10. Ito versus Stratonovich interpretation. For completeness, we provide here the form of the PDE governing the mean-return-time function T x, y ( ) when the underlying stochastic differential equations

$$( 1 0. 1 )$$

$$\dot { x } & = f _ { x } ( x, y ) + g _ { x 1 } ( x, y ) \xi _ { 1 } ( t ) + g _ { x 2 } ( x, y ) \xi _ { 2 } ( t ) \\ \dot { y } & = f _ { y } ( x, y ) + g _ { y 1 } ( x, y ) \xi _ { 1 } ( t ) + g _ { y 2 } ( x, y ) \xi _ { 2 } ( t )$$

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

N

L

G

v

f

K

n

Q

d

V

J

p

9

8

y

o

&gt;

A

B

H

c

D

S

g

E

r

U

Y

O

7

T

0

X

M

F

2

k

m

z

w

u

W

P

I

j

q

/

5

+

R

C

Z

3

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

M

3

7

v

j

J

k

I

K

w

V

2

5

&gt;

A

B

8

H

c

D

L

S

g

N

E

y

r

f

U

Y

9

o

O

T

0

G

X

F

m

z

u

W

P

n

d

Q

q

/

+

R

p

C

Z

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

u

E

c

A

B

0

y

J

r

z

V

M

v

Y

7

8

X

k

&gt;

H

D

L

S

g

N

f

U

9

o

O

T

G

F

2

m

w

W

P

I

j

n

d

Q

q

/

5

+

R

p

C

K

Z

3

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

Z

o

u

p

I

y

7

A

0

r

O

T

+

g

3

8

&gt;

B

H

c

V

N

S

w

M

E

W

q

/

9

J

v

X

Y

D

k

U

2

m

z

R

P

f

n

d

5

Q

K

C

F

G

j

L

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

C

8

z

W

H

X

D

R

v

U

f

Q

E

d

y

Z

w

&gt;

A

B

3

c

V

N

S

r

q

/

9

L

0

o

M

j

T

7

u

k

P

m

K

5

p

n

J

O

g

I

Y

G

+

F

2

Fig. 8.2: Sensitivity of the MRT function T to the size of the interior domain boundary glyph[epsilon1] . (a) Absolute difference | T x, y, glyph[epsilon1] ( 1 ) -T x, y, glyph[epsilon1] ( 2 ) | at each point in grid for glyph[epsilon1] 1 = 0 05 . and glyph[epsilon1] 2 = 0 1. . The inner square has side length 0.1. (b) Absolute difference | T x, y, glyph[epsilon1] ( 1 ) -T x, y, glyph[epsilon1] ( 2 ) | for glyph[epsilon1] 1 = 0 05 and . glyph[epsilon1] 2 = 0 2. . The inner square has side length 0.2. Mean period T = 16 23. . The relative shift in the time function due to changing the inner boundary size is less than one part in a thousand throughout most of the region. Immediately adjacent to the inner boundary, relative deviations occur of up to 8% of T . (c) Comparison of a single MRT isochron (chosen to have a common value of T at the outer boundary) for glyph[epsilon1] = 0 05 (solid black curve), . glyph[epsilon1] = 0 1 (red dashed . curve), and glyph[epsilon1] = 0 2 (cyan dashed curve). The inner square has side length 0.2. (d) As . in panel C, magnified to show detail. The isochrons obtained with different interior domain sizes glyph[epsilon1] differ appreciably only in a neighborhood of size O glyph[epsilon1] ( ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-f35d3844-7611-4726-bf7e-4e2f8b2d4f2c.jpg)

are understood in the Ito and the Stratonovich interpretations [10]. In both cases the mean-return-time function T x, y ( ) satisfies the inhomogeneous PDE L † [ T ] = -1 with adjoint reflecting boundary conditions on the inner and outer domain boundaries, and a jump of T across an arbitrary simple curve connecting the inner and outer boundaries. The form of the adjoint Kolmogorov operator (generator of the Markov process) is given below.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00487v1-figure-7eb64323-2f0e-423c-b435-682e048616d9.jpg)

x

Fig. 9.1: Stuart-Landau oscillator with counterrotating phase, eq. (9.1), with ω = 1, γ = 15, c = 4, and D = 0 18. . Sample trajectory (∆ t = 10 -4 ) with initial condition at (0 , -1) with the deterministic limit cycle shown by the bold grey line and the deterministic velocity field (a) [note that the field rotates clockwise for sufficiently large radius, which is the mentioned counterrotation], MRT phase Θ( x, y ) from the finite-difference-scheme solution (b) and its contour lines, i.e. the isochrons of the MRT phase (c). For the parameters used, T = 6 284. .

10.1. Ito interpretation. Under the Ito interpretation of the SDE eq. (10.1), the form of L † is the same as that given with coordinates ( α, β ) in the main text (cf. (2.18)).

$$( 1 0. 2 ) \ \mathcal { L } ^ { \dagger } [ u ] = \left [ f _ { x } \frac { \partial } { \partial x } + f _ { y } \frac { \partial } { \partial y } + \frac { 1 } { 2 } \left ( \mathcal { G } _ { x x } \frac { \partial ^ { 2 } } { \partial x ^ { 2 } } + \mathcal { G } _ { x y } \frac { \partial ^ { 2 } } { \partial x y } + \mathcal { G } _ { y x } \frac { \partial ^ { 2 } } { \partial y x } + \mathcal { G } _ { y y } \frac { \partial ^ { 2 } } { \partial y ^ { 2 } } \right ) \right ] [ u ],$$

where G = gg ᵀ and g = ( g x 1 g x 2 g y 1 g y 2 ) . The adjoint reflecting boundary conditions at an exterior boundary with normal vector ( n , n 1 2 ) ᵀ are

$$\sum _ { j } n _ { j } \sum _ { k } \mathcal { G } _ { j k } \partial _ { k } T = 0.$$

Writing out L † explicitly in terms of our drift and diffusion coefficients gives

(10.4)

$$\stackrel { \ll - \cdots \cdot \cdot } { \mathcal { L } ^ { \dagger } = f _ { x } } \frac { \partial } { \partial x } + f _ { y } \frac { \partial } { \partial y } + \frac { 1 } { 2 } \left ( ( g _ { x 1 } ^ { 2 } + g _ { x 2 } ^ { 2 } ) \frac { \partial ^ { 2 } } { \partial x ^ { 2 } } + 2 ( g _ { x 1 } g _ { y 1 } + g _ { x 2 } g _ { y 2 } ) \frac { \partial ^ { 2 } } { \partial x y } + ( g _ { y 1 } ^ { 2 } + g _ { y 2 } ^ { 2 } ) \frac { \partial ^ { 2 } } { \partial y ^ { 2 } } \right ) \cdot \psi$$

10.2. Stratonovich interpretation. If we interpret the SDE eq. (10.1) in the sense of Stratonovich, this is equivalent to an Ito interpretation of an SDE with modified drift coefficients ˆ f x , f ˆ : y

$$\dot { x } & = \hat { f } _ { x } ( x, y ) + g _ { x 1 } ( x, y ) \xi _ { 1 } ( t ) + g _ { x 2 } ( x, y ) \xi _ { 2 } ( t ) \\ \dot { y } & = \hat { f } _ { y } ( x, y ) + g _ { y 1 } ( x, y ) \xi _ { 1 } ( t ) + g _ { y 2 } ( x, y ) \xi _ { 2 } ( t ),$$

where

$$\hat { f } _ { x } \, = \, & f _ { x } + \frac { 1 } { 2 } \left ( g _ { x 1 } \frac { \partial } { \partial x } g _ { x 1 } + g _ { x 2 } \frac { \partial } { \partial x } g _ { x 2 } + g _ { y 1 } \frac { \partial } { \partial y } g _ { x 1 } + g _ { y 2 } \frac { \partial } { \partial y } g _ { x 2 } \right ) \\ \hat { f } _ { y } \, = \, & f _ { y } + \frac { 1 } { 2 } \left ( g _ { x 1 } \frac { \partial } { \partial x } g _ { y 1 } + g _ { x 2 } \frac { \partial } { \partial x } g _ { y 2 } + g _ { y 1 } \frac { \partial } { \partial y } g _ { y 1 } + g _ { y 2 } \frac { \partial } { \partial y } g _ { y 2 } \right ) ;$$

$$- \, _ { J \, y \, ^ { \top } } \, _ { 2 } \, \left ( ^ { y x _ { 1 } } \, _ { \partial x } ^ { y y _ { 1 } \, ^ { \top } \, ^ { y x _ { 2 } } } \, _ { \partial x } ^ { y y _ { 2 } \, ^ { \top } \, ^ { y y _ { 1 } } } \, _ { \partial y } ^ { y y _ { 1 } \, ^ { \top } \, ^ { y y _ { 2 } } } \, _ { \partial y } ^ { y y _ { 2 } } \right ) \,,$$

compare eqs. (4.3.40, 4.3.43) of [10].

Substituting the modified drift coefficients ˆ in to (10.2) or (10.4) gives the form of f the PDE satisfied by T . The change of interpretation from Ito to Stratonovich does not change the diffusion coefficients g ij , and thus the adjoint reflecting boundary conditions are identical in both cases.
Alexander Cao Benjamin Lindner Peter J. Thomas
A Partial Differential Equation for the Mean--Return-Time Phase of Planar Stochastic Oscillators
Stochastic oscillations are ubiquitous in many systems. For deterministic
systems, the oscillator's phase has been widely used as an effective
one-dimensional description of a higher dimensional dynamics, particularly for
driven or coupled systems. Hence, efforts have been made to generalize the
phase concept to the stochastic framework. One notion of phase due to
Schwabedal and Pikovsky is based on the mean-return time (MRT) of the
oscillator but has so far been described only in terms of a numerical
algorithm. Here we develop the boundary condition by which the partial
differential equation for the MRT has to be complemented in order to obtain the
isochrons (lines of equal phase) of a two-dimensional stochastic oscillator,
and rigorously establish the existence and uniqueness of the MRT isochron
function (up to an additive constant). We illustrate the method with a number
of examples: the stochastic heteroclinic oscillator (which would not oscillate
in the absence of noise); the isotropic Stuart-Landau oscillator, the
Newby-Schwemmer oscillator, and the Stuart-Landau oscillator with polarized
noise. For selected isochrons we confirm by extensive stochastic simulations
that the return time from an isochron to the same isochron (after one complete
rotation) is always the mean first-passage time (irrespective of the initial
position on the isochron). Put differently, we verify that Schwabedal and
Pikovsky's criterion for an isochron is satisfied. In addition, we discuss how
to extend the construction to arbitrary finite dimensions. Our results will
enable development of analytical tools to study and compare different notions
of phase for stochastic oscillators.