March 2012
EPL, 97 (2012) 50009
doi: 10.1209/0295-5075/97/50009

www.epljournal.org

Colored noise induces synchronization of limit cycle oscillators
Wataru Kurebayashi1(a) , Kantaro Fujiwara1 and Tohru Ikeguchi1,21Graduate School of Science and Engineering, Saitama University - 255 Shimo-ohkubo, Sakura-ku, Saitama-city,
Saitama, 338-8570 Japan
2
Saitama University Brain Science Institute - 255 Shimo-ohkubo, Sakura-ku, Saitama-city, Saitama, 338-8570 Japan
received 12 November 2011; accepted in final form 7 February 2012published online 5 March 201205.45.Xt – Synchronization; coupled oscillators
02.50.Ey – Stochastic processes
PACS 05.40.Ca – Noise
PACS
PACS

Abstract – Driven by various kinds of noise, ensembles of limit cycle oscillators can synchronize.
In this letter, we propose a general formulation of synchronization of the oscillator ensembles driven
by common colored noise with an arbitrary power spectrum. To explore statistical properties ofsuch colored noise-induced synchronization, we derive the stationary distribution of the phase
diﬀerence between two oscillators in the ensemble. This analytical result theoretically predicts
various synchronized and clustered states induced by colored noise and also clarifies that these
phenomena have a diﬀerent synchronization mechanism from the case of white noise.
c EPLA, 2012
Copyright ⃝Introduction. – Driven by common noise, many
nonlinear dynamical systems can synchronize. Thisphenomenon is called noise-induced synchronization,
which is observed in various kinds of the nonlinear
dynamical systems, for example, neural networks [1,2],
electric circuits [3], electronic devices [4], microbial
cells [5], lasers [6] and chaotic dynamical systems [7,8]. It
has been theoretically proven that limit cycle oscillators
can synchronize driven by common noise [9]. Many
studies have investigated the synchronization property
in case of various types of drive noises, for example,
Gaussian white noise [10–12] and Poisson impulses [13].
In ref. [10], using a formulation of limit cycle oscillators
driven by common and independent Gaussian white
noises, Nakao et al. analytically obtained the probability
density function (PDF) of phase diﬀerences between two
oscillators, which enables us to eﬀectively characterize thesynchronization property. However, although there aresome numerical studies [8,14,15], analytical conventional
studies are limited to the case that drive signals arewhite noise (temporally uncorrelated noise). If we canassume that the drive signal is white noise, we can use the
Fokker-Planck approximation [16] to explore statisticalproperties of oscillator ensembles. However, such an ideal
condition is rare in the real world. For example, in neuralcircuits, it is known that colored noise with negative
autocorrelation plays a key role to propagate synchronous(a) E-mail: kure@nls.ics.saitama-u.ac.jp

activities [17]. However, it still remains unclear how theoscillators behave if they are driven by common colorednoise.
Recently, it has been clarified how a limit cycle oscillator behaves if it is driven by colored non-Gaussian noise[18–20]. In this letter, utilizing eﬀective white-noise
Langevin description proposed in ref. [19], we extend
the formulation in ref. [10] to colored noise that has anarbitrary power spectrum. We then analytically derive
the PDF of the phase diﬀerence between the oscillators
if these oscillators are driven by common colored noise.
We also conducted numerical simulations to verify our
analytical results. The results show that the PDF of thephase diﬀerence explicitly depends on the power spectrumof the drive noise.
Model. – We used the following system that consists of
N identical limit cycle oscillators subject to common andindependent multiplicative colored noises. The dynamics
of the j -th oscillator is described by√Ẋ (j) = F (X (j) ) + DG(X (j) )ξ(t)√+ ϵH(X (j) )η (j) (t),(1)
for j = 1, . . . , N , where X (j) ∈ Rn is the n-dimensional
state variable of the j -th oscillator; F (X (j) ) ∈ Rn is anunperturbed vector field that has a stable T -periodic limit
cycle orbit S(t); ξ(t) ∈ Rm is the common noise, which
drives all of the oscillators; η (j) (t) ∈ Rm (j = 1, . . . , N ) is
the independent noise, which is received independently

50009-p1

Wataru Kurebayashi et al.
by each oscillator; G(X (j) ) ∈ Rn×m and H(X (j) ) ∈ Rn×m
represent how the oscillators are coupled to the common
and independent noises; D and ϵ are parameters to control
the intensities of the common and independent noises. We
introduced the following three assumptions: i) ξ(t) ∈ Rm
and η (j) (t) ∈ Rm are independent, identically distributed
zero-mean colored noises, namely, ⟨ξ(t)⟩ = 0, ⟨η (j) (t)⟩ =0, ⟨ξ(t)η (j) (s)⊤ ⟩ = O, and ⟨η (j) (t)η (k) (s)⊤ ⟩ = O (j ̸= k),where ⊤ denotes the transpose and ⟨·⟩ represents the
temporal average; ii) ξ(t) and η (j) (t) can be approximatedas the convolution of an arbitrary filter function andwhite noise;and iii) ξ(t) and η (j) (t) have correlation
times shorter than the time scale of the phase diﬀusion
1
1(∼O(D− 2 , ϵ− 2 )).
To characterize the statistical properties of the drive
noises ξ(t) and η (j) (t), we define correlation matricesCξ (τ ) ∈ Rm×m and Cη (τ ) ∈ Rm×m as Cξ (τ ) = ⟨ξ(t)ξ(t −τ )⊤ ⟩ and Cη (τ ) = ⟨η (j) (t)η (j) (t − τ )⊤ ⟩ (j = 1, . . . , N ). For
the sake of simplicity, we assumed that all independent
noises η (j) (t) have the same statistical property characterized by Cη (τ ). The (i, j)-th element of Cξ (τ ) is the
cross-correlation function of the i-th and j -th elements ofthe common noise ξ(t). The diagonal elements of Cξ (τ ) areautocorrelation functions. In the same way, we can characterize the statistical property of η (j) (t) by using Cη (τ ).

f (θ, t) is described by the following eﬀective Fokker-Planck
equation:Eﬀective Langevin description. – To quantify thesynchronization property without loss of generality, we
consider the relationship of only two oscillators, thatis, the two-body problem of φ(1) (t) and φ(2) (t), anddefine the phase diﬀerence θ (:= φ(1) − φ(2) ). As we focus
on the stochastic dynamics of θ, we define f (θ, t) as the
PDF of the phase diﬀerence θ. Utilizing the eﬀective
white-noise Langevin description [19], the evolution of

×ZH (φ)⊤ Cη (τ )ZH (φ − ωτ ) + O(D 2 , ϵ 2 ). (7)

∂
1 ∂2
∂f
+ v (1) (θ)f − · 2 v (2) (θ)f = 0,
∂t ∂θ
2 ∂θ

(3)

where v (1) (θ) and v (2) (θ) are eﬀective drift and diﬀusion coeﬃcients. We have the drift coeﬃcient v (1) (θ) = 0
because ⟨θ̇⟩ = ⟨φ̇(1) − φ̇(2) ⟩ = 0. Meanwhile, the diﬀusion
coeﬃcient v (2) (θ) is obtained as
! +∞
"#
$#
$%
v (2) (θ) =dτ θ̇(t) − ⟨θ̇⟩ θ̇(t − τ ) − ⟨θ̇⟩
−∞

=

! +∞

dτ

−∞
# (1)

× φ̇

"# (1)$φ̇ (t) − φ̇(2) (t)$%(t − τ ) − φ̇(2) (t − τ ) ,(4)

where ⟨·⟩ represents the temporal average. For simplicity& +∞ "of notation, we define djk as djk = −∞ dτ [φ̇(j) (t) −%ω][φ̇(k) (t − τ ) − ω] . Then, we obtainv (2) (θ) = d11 + d22 − d12 − d21 = 2d11 − 2d12 .

(5)

The phase variable φ(j) (t) can be expanded as φ(j) (t) =

√ (j)
√ (j)
(j)
(j)
(j)
φ0 (t) + DφD,1 (t)+ ϵφϵ,1 (t)+ DφD,2 (t)+ ϵφϵ,2 (t) + · · ·√√by using D and ϵ as expansion parameters, where
(j)
(j)
(j)
Phase reduction. – Under the assumption that the φ0 (t), φD,k (t) and φϵ,k (t) (k = 1, 2, . . .) are approxi(j)noise intensity is suﬃciently weak (D ≪ 1 and ϵ ≪ 1), we mate perturbed solutions of φ(j) (t). We have φ0 (t) =
(j)
(j)
(j)
can apply the phase reduction method [20,21] to eq. (1). φ(j)0 (0) + ωt, φ̇D,1 (t) = ZG (φ0 (t)) · ξ(t) and φ̇ϵ,1 (t) =By introducing a phase variable φ(j) , eq. (1) is reduced to(j)
ZH (φ0 (t)) · η (j) (t). Using these perturbed solutions,√ (j) √ (j)the following phase equation:eq. (2) can be written as φ̇(j) = ω + Dφ̇D,1 + ϵφ̇ϵ,1 +√O(D, ϵ). Using this approximation and the fact thatφ̇(j) = ω + DZG (φ(j) ) · ξ(t)
(j)
(k)
(j)
(k)
√
(j)
(j)
+ ϵZH (φ ) · η (t) + O(D, ϵ),(2) ⟨φD,1 (t)φϵ,1 (t − τ )⟩ = 0 and ⟨φϵ,1 (t)φD,1 (t − τ )⟩ = 0, we
obtain! +∞where φ(j) (t) ∈ [−π, +π] is a phase variable that corre" (j)
%
(k)
(j)
=
D
dτ φ̇D,1 (t)φ̇D,1 (t − τ )djksponds to the state of the j -th oscillator X , ω−∞(= 2πT −1 ) is the natural frequency, and ZG (φ(j) )
! +∞
" (j)
%
(j)
(k)
and ZH (φ ) are the phase sensitivity functions that+ϵdτ φ̇ϵ,1 (t)φ̇ϵ,1 (t − τ )represent the linear response of the phase variable−∞
3
3
φ(j) to the drive noises [20,21]. The phase sensitiv(6)+O(D 2 , ϵ 2 ).
(j)(j)
ity functions ZG (φ ) and ZH (φ ) are defined as
Thus, using eq. (6), we can calculate d11 as follows:
ZG (φ(j) ) = ∇X φ(j) |X=S(φ(j) ) · G(S(φ(j) )) and ZH (φ(j) )=
! +π
!
(j)
(j)
∇X φ |X=S(φ(j) ) · H(S(φ )). As discussed in ref. [20],
D +∞
d11 =
dτ
dφ
the O(D, ϵ) term is necessary to describe the exact phase
2π −∞−πdynamics, while the phase diﬀusion is not aﬀected by
×ZG (φ)⊤ Cξ (τ )ZG (φ − ωτ )the O(D, ϵ) term. As we will focus on the phase diﬀusion! +∞ ! +π
in the following sections, we do not take this term into
ϵ
+
dτ
dφ
account.
2π −∞
−π
3

3

In the same way, d12 is given by
!
! +π
D +∞
d12 =
dτ
dφ
2π −∞
−π

50009-p2

×ZG (φ)⊤ Cξ (τ )ZG (φ − θ − ωτ )3

3

+O(D 2 , ϵ 2 ).

(8)Colored noise induces synchronization of limit cycle oscillators
The detailed derivations of eqs. (7) and (8) are shown
Using the Fourier representations defined above, wein appendix A.
can obtain the Fourier representations of the correlation
Finally, from eqs. (5), (7) and (8), we have the eﬃcient functions g(θ) and h(θ):
diﬀusion coeﬃcient v (2) (θ):
+∞
+∞
(
(
#
$
gl eilθ , h(θ) =hl eilθ ,(13)g(θ) =(2)v (θ) = 2D g(0) − g(θ) + 2ϵh(0),
(9)
l=−∞

where g(θ) and h(θ) are correlation functions defined as
1g(θ) =
2π

! +∞

dτ

! +∞

dτ

−∞

! +π

dφ

! +π

dφ

−π

×ZG (φ)⊤ Cξ (τ )ZG (φ − θ − ωτ ),h(θ) =

1
2π

−∞

−π

×ZH (φ)⊤ Cη (τ )ZH (φ − θ − ωτ ).

(10)

(11)

If we assume that the drive noise is white, namely, Cξ (τ ) =Cη (τ ) = δ(τ )Em , eqs. (10) and (11) are exactly equivalent
to eq. (6) in ref. [10], where Em is an m × m identity
matrix. The results show that eqs. (10) and (11) are anatural generalization of eq. (6) in ref. [10].
We obtain the explicit form of the Fokker-Planck equation of eq. (3) from eqs. (9)–(11). The stationary distribution of the phase diﬀerence f0 (θ) is given as the stationary
solution of eq. (3). Then, if we put ∂f /∂t = 0 in eq. (3),
we obtain
f0 (θ) =

ν
ν′
#
$
,
=
v (2) (θ) D g(0) − g(θ) + ϵh(0)(12)

where ν and ν ′ (= ν/2) are normalization constants.
Fourier representation. – To understand theresults obtained in the previous section, we rewrite
the correlation functions defined in eqs. (10) and (11)by using the Fourier representation. We introduced
the Fourier series expansion of the phase'+∞ sensitivity
functions ZG (φ) and ZH (φ) as ZG (φ) = l=−∞ YG,l eilφ'+∞and ZH (φ) = l=−∞ YH,l eilφ , where i denotes the& +π1
imaginary unit and YG,l ∈ Cm (= 2πdφZG (φ)e−ilφ )
−π
&
+π
1
and YH,l ∈ Cm (= 2πdφZH (φ)e−ilφ ) are Fourier−π
coeﬃcients (l = −∞, . . . , ∞).
andSubsequently,
wedefinePξ (Ω) ∈ Cm×m
m×m
as the Fourier transforms of Cξ (τ )Pη (Ω) ∈ C& +∞and Cη (τ ), that is, Pξ (Ω) = −∞ dtCξ (t)e−iΩt and& +∞Pη (Ω) = −∞ dtCη (t)e−iΩt . Let us note that Pξ (Ω) andPη (Ω) are Hermitian matrices, namely, Pξ (Ω) = Pξ (Ω)†
and Pη (Ω) = Pη (Ω)† because Cξ (τ ) = Cξ (−τ )⊤ andCη (τ ) = Cη (−τ )⊤ from their definitions, where † denotes
the adjoint. The (i, j)-th elements of Pξ (Ω) and Pη (Ω)
represent the cross-spectra of the i -th and j -th elementsof ξ(t) and η (j) (t). In particular, the diagonal elements ofPξ (Ω) and Pη (Ω) represent the power spectra.

l=−∞

†
†
where gl (= YG,lPξ (lω)YG,l ) and hl (= YH,lPη (lω)YH,l )are Fourier coeﬃcients (l = −∞, . . . , ∞). The derivationsof gl and hl will be shown in appendix B.
These expressions clearly suggest that the correlation
functions g(θ) and h(θ) only depend on Pξ (±lω) andPη (±lω) (l = 0, 1, 2, . . .), that is, the other frequency
components can be neglected. In the next section, we will
demonstrate that colored noise induces various synchronized and clustered states, which are clearly explained byeq. (13).

Numerical simulations. – To demonstrate the validity of our results, we perform numerical experiments for
two types of limit cycle oscillators. The first example
is the Stuart-Landau oscillator, which takes the normal
form of the supercritical Hopf bifurcation [21]: ẋ = x −c0 y − (x2 + y 2 )(x − c2 y), ẏ = y + c0 x − (x2 + y 2 )(y + c2 x),
where X = [x, y]⊤ is a state variable and c0 and c2
are parameters. In the simulation, we fixed c0 = 1, c2 =−2, G = H = diag(1, 1), D = 0.0095 and ϵ = 0.0005, wherediag(λ1 , . . . , λm ) denotes an m × m diagonal matrix thathas the diagonal elements λ1 , . . . , λm . This model isreduced to the phase equation that has the natural
frequency ω √= c0 − c2 = 3 and the phase sensitivity function Z(φ) = 2[sin(φ + 3π/4), sin(φ + π/4)]⊤ .
In the simulation, we use a two-dimensional drive
noise that has the correlation matrix Cex (τ ) ∈ R2×2
defined as Cex (τ ) = diag(Cex (τ ), Cex (τ )) and Cex (τ ) =γ −γ|τ |cos ω0 τ , where ω0 and γ are parameters that2e
represent the peak frequency and the characteristic decay
time. We define Pex (Ω), the Fourier transform of Cex (τ ),2
as Pex (Ω) = γ2 {[γ 2 + (Ω + ω0 )2 ]−1 + [γ 2 + (Ω − ω0 )2 ]−1 }.
A drive noise characterized by Cex (τ ) can be generated by
the damped noisy harmonic oscillator (see eqs. (43)–(49)in ref. [19] for details).
We use the common noises with (ω0 , γ) = (1, 1), (3, 1)
and (5, 1) and the independent noise with (ω0 , γ) = (0, 3).
The power spectra of these common noises are shown infig. 1(a). From eq. (13), the correlation functions g(θ) andh(θ) are given by g(θ) = {[1 + (ω0 + 3)2 ]−1 + [1 + (ω0 −3)2 ]−1 }cos θ and h(θ) = cos θ, for ω0 = 1, 3 and 5, whichcorrespond to the three types of the common noise. The
derivations of g(θ) and h(θ) will be shown in appendix C.
The correlation function g(θ) calculated above indicatesthat the eﬀective intensity of the common noise depends
on the peak frequency ω0 and is maximal at ω0 = ω. It
means that the synchronous degree is maximized at ω0 =
ω. In figs. 1(b)–(d), we compared the results of the direct
numerical simulation using the Stuart-Landau oscillator
and its corresponding phase oscillator with the analytical

50009-p3

Wataru Kurebayashi et al.

1.0
0.5
0.0

(a)

0

1

2

3

frequency Ω

0.5

FitzHugh-Nagumo

0.4
0.3
0.2
0.1
0.0

(c)

-π

-π/2

0

π/2

phase difference θπprobability density f (θ)

1.5

0.5

Stuart-Landau

0.4
0.3
0.2
0.1
0.0

(b)
probability density f (θ)

power spectrum P (Ω)

results. All PDFs are well fitted by the theoretical curves.
Our theory clearly predicts that the highest synchronous
degree is realized at ω0 = ω.
The second example is the FitzHugh-Nagumo oscillator [22,23]: v̇ = v − v 3 /3 − u + I0 , u̇ = µ(v + a − bu), whereX = [v, u]⊤ is a state variable and a, b, µ and I0 are
parameters. In the simulation, we fixed a = 0.7, b = 0.8,
µ = 0.08, I0 = 0.875, G = H = diag(1, 0), D = 0.045 andϵ = 0.005. For these parameters, this oscillator has thenatural frequency ω ≃ 0.1725. This oscillator models
bursting behavior of a neuron, and only the first variablev, which corresponds to the membrane potential of aneuron, is subject to noise.
In the simulation, we use the one-dimensional noise thathas the correlation function Cex (τ ). Diﬀerent from the firstexample, we use the same parameters (ω0 , γ) for both thecommon and independent noises. We used two parameter
sets (ω0 , γ) = (ω, 0.1) and (3ω, 0.1). The power spectra ofthese drive noises are shown in fig. 2(a). We obtain thecorrelation function g(θ) (= h(θ)) numerically as shown in
figs. 2(b).
In figs. 2(c) and (d), we compared the results of thedirect numerical simulation with the analytical results.
The numerical results are in good agreement with thetheoretical results. As theoretically predicted, a 3-cluster
state is realized as shown in fig. 2(d). If oscillators are
driven by white noise, clustered states are induced only bymultiplicative noise [10]. However, in case of colored noise,clustered states are induced not only by multiplicative
noise but also by additive noise.

Fig. 2: (Color online) Simulation results of the FitzHughNagumo oscillator. (a) Power spectra of the common noises areshown for ω0 = ω and 3ω. For these drive noises, g(θ) (= h(θ))
is shown (b). The PDFs of θ for ω0 = ω (synchronized state)(c) and for ω0 = 3ω (3-cluster state) (d) are shown.

probability density f (θ)

Fig. 1: (Color online) Simulation results of the Stuart-Landau
oscillator (crosses) and the corresponding phase oscillator(open circles). (a) Power spectra of the common noises areshown for ω0 = 1, 3 and 5. The PDFs of θ show the frequency
dependency of the synchronization property for ω0 = 1 (b),
ω0 = 3 (= ω) (c), and ω0 = 5 (d).

-π

-π/2

0

π/2

π

phase difference θ

0.5

Hodgkin-Huxley

0.4
0.3
0.2
0.1
0.0

(d)

-π

-π/2

0

π/2

π

phase difference θ

Fig. 3: (Color online) Simulation results of the limit cycle
oscillators subject to green noise. (a) Power spectra of thedrive noises. The PDFs of θ obtained by the theory (lines) andnumerical simulations (circles) for the Stuart-Landau oscillator(b), the FitzHugh-Nagumo oscillator (c) and the HodgkinHuxley oscillator (d).

In the third example, we used the Hodgkin-Huxley
oscillator [24], which enables us to demonstrate whether
the theory is applicable to higher-dimensional limit cycle
systems. We use green noise used in ref. [8], which is
generated by applying a high-pass filter to white noise.
The power spectrum is shown in fig. 3(a). Diﬀerent fromthe periodic noise characterized by Cex (τ ), the green noise
has a vanishing spectrum as Ω → 0. In the simulation, forthe sake of simplicity, we used the same type of drive noisefor the common and independent noises, and we set the50009-p4

Colored noise induces synchronization of limit cycle oscillators
noise intensities (D, ϵ) = (0.0002, 0.0001). Figures 3(b)–(d)We rewrite ZG (φ), ZH (φ), ξ(t) and η (1) (t) by using
compare the theoretical and numerical results, which show their elements and obtain! +∞ (
that our theory is also valid for these cases.
m (
m
"
(1)
ZH,k (φ0 (t))
=
D
dτ
d
11
Summary and discussions. – In this letter, we

extended a formulation to analyze various synchronized
and clustered states of uncoupled limit cycle oscillators
driven by common and independent colored noises. Using
this formulation, we derived the probability density function of the phase diﬀerence and rewrote it by the Fourier
representation. The obtained expressions clearly show thatthe synchronization property depends on the power spectrum of the drive noises. Such dependency has alreadybeen reported experimentally. For example, in ref. [25],the reliability, or synchronization across trials, is exploredin neuronal responses to periodic drive inputs with various frequencies. The reliability is maximized at a certain
frequency, which is similar to our results shown in fig. 1.
Our results in this letter theoretically support the results
in ref. [25], because a neuron in an oscillatory state can be
regarded as a noisy limit cycle oscillator.
Generally, noise in the real world often has a nonflat and characteristic power spectrum. In this sense, ourformulation is a useful tool to estimate the synchronization property for both theoretical and practical aspects.
Namely, the results obtained in this letter can be applied
to a wide range of purposes from mathematical modelings
to technological problems.

−∞

%
(1)
×ξk (t)ξl (t − τ )ZH,l (φ0 (t − τ ))
! +∞
m (
m
(
"
(1)
ZH,k (φ0 (t))
dτ
+ϵ

−∞
k=1 l=1
%
(1)
(1)
(1)
×ηk (t)ηl (t − τ )ZH,l (φ0 (t − τ ))3

d11 = D

+ϵ

3

3

3

+ O(D 2 , ϵ 2 ).

3

3

3

+O(D 2 , ϵ 2 ),

(A.1)

"
(1)
dτ ZG (φ0 (t))⊤ ξ(t)%(1)(t − τ )⊤ ZH (φ0 (t − τ ))(1)

ZH,k (φ0 (t))

k=1 l=1(A.4)

where Cξ,kl and Cη,kl are the (k, l)-th elements of Cξ (φ)
and Cη (φ). Finally, we rewrite eq. (A.4) by using ZG (φ),
ZH (φ), Cξ (τ ) and Cη (τ ) and obtain
! +π
!
D +∞
dτ
dφ
d11 =
2π −∞
−π
×ZG (φ)⊤ Cξ (τ )ZG (φ − ωτ )! +∞ ! +π
ϵ
+
dτ
dφ
2π −∞
−π

−∞

×η

m (
m
(
"

! +π
ϵ
+
dτ
dφ
2π −∞
−π
m (
m
(
ZH,k (φ)ZH,l (φ − ωτ )Cη,kl (τ )×%
(1)×ξ(t − τ )⊤ ZG (φ0 (t − τ ))
! +∞
"
(1)
dτ ZH (φ0 (t))⊤ η (1) (t)
+ϵ
(1)

dτ

k=1 l=1
! +∞

−∞

−∞

! +∞

+O(D 2 , ϵ 2 )
! +π
!
D +∞
=
dτ
dφ
2π −∞
−π
m (
m
(
ZG,k (φ)ZG,l (φ − ωτ )Cξ,kl (τ )
×

−∞

! +∞

(1)

ZG,k (φ0 (t))

k=1 l=1% " (1)%(1)(1)
×ZH,l (φ0 (t − τ )) φ ηk (t)ηl (t − τ ) t

(1)

=D

m (
m
(
"

k=1 l=1% "%(1)
×ZG,l (φ0 (t − τ )) φ ξk (t)ξl (t − τ ) t

ZH (φ0 (t)) · η (1) (t) into eq. (6), we obtain
! +∞
"
(1)
d11 = Ddτ [ZG (φ0 (t))⊤ ξ(t)]+ O(D , ϵ )dτ−∞Appendix A: derivations of eqs. (7) and (8).
(j)
(j)
(j)
– Substituting φ̇D,1 = ZG (φ0 (t)) · ξ(t) and φ̇ϵ,1 =

3
2

! +∞

−∞

%
(1)
×[ZG (φ0 (t − τ ))⊤ ξ(t − τ )]
! +∞
"
(1)
dτ [ZH (φ0 (t))⊤ η (1) (t)]
+ϵ

(A.3)

where ZG,l (φ) and ZH,l (φ) are the lth elements of ZG (φ)(1)
and ZH (φ), and ξl (t) and ηl (t) are the lth elements of(1)
ξ(t) and η (t).
We assume that the phase variable φ(1) and the drive
noises ξ(t) and η (1) (t) are approximately independent.
Under this assumption, the temporal average ⟨·⟩ can be& +πdivided into two parts; ⟨·⟩φ (:= (2π)−1 −π dφ · ) and ⟨·⟩t
& +s(:= lims→∞ (2s)−1 −s dt · ). Thus, we obtain

The authors would like to thank S. Ogawa and AGS
Corporation for their encouragement on this research
project.

3
2

3+ O(D 2 , ϵ 2 ),

∗∗∗

%
(1)
×[ZH (φ0 (t − τ ))⊤ η (1) (t − τ )]k=1 l=1(A.2)
50009-p5

×ZH (φ)⊤ Cη (τ )ZH (φ − ωτ )3

3

+O(D 2 , ϵ 2 ).

(A.5)

Wataru Kurebayashi et al.
In the same way, one can calculate d12 as follows. We is a parameter. In the same way, the Fourier coeﬃcient
†
1
(1)
(2)
use the fact that ⟨φϵ,1 (t)φϵ,1 (t − τ )⟩ = 0 and eliminate the hl is given by h±1 = YH,±l YH,±l Pex (ω)|ω0 =0,γ=3 = 2 and(2)phase variable of the second oscillator φ0 by substituting hl = 0 (l ̸= ±1). Substituting gl and hl into eq. (13), we
can obtain the explicit forms of g(θ) and h(θ).
(2)
(1)
(2)
φ0 = φ0 − θ into φ0 , and then, we obtain
d12 = D

! +∞
−∞

"
(1)
dτ [ZG (φ0 (t))⊤ ξ(t)]REFERENCES%(2)×[ZG (φ0 (t − τ ))⊤ ξ(t − τ )]3

3

+O(D 2 , ϵ 2 )
! +∞
"
(1)
=D
dτ ZG (φ0 (t))⊤ ξ(t)
−∞

%
(2)
×ξ(t − τ )⊤ ZG (φ0 (t − τ ))3

3

+O(D 2 , ϵ 2 )
!
! +π
D +∞
=
dτ
dφ
2π −∞
−π

×ZG (φ)⊤ Cξ (τ )ZG (φ − θ − ωτ )3

3

+O(D 2 , ϵ 2 ).

(A.6)

Appendix B: derivation of eq. (13). – Fromeq. (10), one can calculate the Fourier coeﬃcient gl as
follows. We introduce a new variable χ (:= φ − θ − ωτ ) and
use the fact that Pξ (Ω) is a Hermitian matrix. Then, we
obtain! +π1dθg(θ)e−ilθ
gl =
2π −π
! +π
! +∞ ! +π
1
1
=
dθ
dτ
dφZG (φ)⊤
2π −π
2π −∞
−π
×Cξ (τ )ZG (φ − θ − ωτ )e−ilθ
)
*
! +π
1
=
dφZG (φ)⊤ e−ilφ
2π −π
)! +∞
*)
*
! +π
1
dτ Cξ (τ )eilωτdχZG (χ)eilχ
×
2π −π
−∞
†
⊤
= YG,lPξ (lω) YG,l = YG,lPξ (lω)† YG,l†= YG,lPξ (lω)YG,l ,(B.1)

where · denotes the complex conjugate. From eq. (11),hl can be derived likewise.
Appendix C: derivations of the correlation functions g(θ) and h(θ). – For the StuartLandau oscillator we used in the simulations, we cancalculate the Fourier coeﬃcients YG,l and YH,l as
andYG,l = YH,l = 0
YG,±1 = YH,±1 = 12 [1 ± i, 1 ∓ i]⊤(l ̸= ±1). Thus, from eq. (13), the Fourier coeﬃcient gl†
is given by g±1 = YG,±l
YG,±l Pex (ω)|γ=1 = 12 {[1 + (ω0 +
2 −1
2 −1
3) ] + [1 + (ω0 − 3) ] } and gl = 0 (l ̸= ±1), where ω0[1] Mainen Z. F. and Sejnowski T. J., Science, 268 (1995)
1503.
[2] Galán R. F., Fourcaud-Trocmé N., Ermentrout
G. B. and Urban N. N., J. Neurosci., 26 (2006) 3646.
[3] Yoshida K., Sato K. and Sugamata A., J. Sound Vib.,
290 (2006) 34.
[4] Utagawa A., Asai T., Hirose T. and Amemiya Y.,
IEICE Trans. Fundam., 91 (2008) 2475.
[5] Zhou T., Chen L. and Aihara K., Phys. Rev. Lett., 95(2005) 178103.
[6] Uchida A., McAllister R. and Roy R., Phys. Rev.
Lett., 93 (2004) 244102.
[7] Zhou C. and Kurths J., Phys. Rev. Lett., 88 (2002)
230602.
[8] Wang Y., Lai Y.-C. and Zheng Z., Phys. Rev. E, 79(2009) 056210.
[9] Teramae J.-N. and Tanaka D., Phys. Rev. Lett., 93(2004) 204103.
[10] Nakao H., Arai K. and Kawamura Y., Phys. Rev.
Lett., 98 (2007) 184101.
[11] Yoshimura K., Davis P. and Uchida A., Prog. Theor.
Phys., 120 (2008) 621.
[12] Nagai K. H. and Kori H., Phys. Rev. E, 81 (2010)
065202.
[13] Nakao H., Arai K., Nagai K., Tsubo Y. andKuramoto Y., Phys. Rev. E, 72 (2005) 26220.
[14] Yoshimura K., Valiusaityte I. and Davis P., Phys.
Rev. E, 75 (2007) 026208.
[15] Hata S., Shimokawa T., Arai K. and Nakao H., Phys.
Rev. E, 82 (2010) 036206.
[16] Risken H., The Fokker-Planck Equation: Methods of
Solution and Applications (Springer-Verlag) 1996.
[17] Câteau H. and Reyes A. D., Phys. Rev. Lett., 96 (2006)
058101.
[18] Teramae J. and Tanaka D., Prog. Theor. Phys., 161(2006) 360.
[19] Nakao H., Teramae J.-N., Goldobin D. S. andKuramoto Y., Chaos, 20 (2010) 3126.
[20] Goldobin D. S., Teramae J.-N., Nakao H. and
Ermentrout G. B., Phys. Rev. Lett., 105 (2010)
154101.
[21] Kuramoto Y., Chemical Oscillations Waves, and Turbulence (Dover Publications) 2003.
[22] FitzHugh R., Biophys. J., 1 (1961) 445.
[23] Nagumo J., Arimoto S. and Yoshizawa S., Proc. IRE,
50 (1962) 2061.
[24] Hodgkin A. and Huxley A., J. Physiol., 117 (1952)500.
[25] Fellous J. M., Houweling A. R., Modi R. H., Rao
R. P. N., Tiesinga P. H. E. and Sejnowski T. J., J.
Neurophysiol., 85 (2001) 1782.

50009-p6