Stochastic switching in delay-coupled oscillators

A delay is known to induce multistability in periodic systems. Under influence of noise, coupled oscillators can switch between coexistent orbits with different frequencies and different oscillation patterns. For coupled phase oscillators we reduce the delay system to a non-delayed Langevin equation, which allows us to analytically compute the distribution of frequencies, and their corresponding residence times. The number of stable periodic orbits scales with the roundtrip delay time and coupling strength, but the noisy system visits only a fraction of the orbits, which scales with the square root of the delay time and is independent of the coupling strength. In contrast, the residence time in the different orbits is mainly determined by the coupling strength and the number of oscillators, and only weakly dependent on the coupling delay. Finally we investigate the effect of a detuning between the oscillators. We demonstrate the generality of our results with delay-coupled FitzHugh-Nagumo oscillators.


I. INTRODUCTION
In recent years dynamical systems with delays have evolved as a major topic in nonlinear sciences [1,2]. Time delays arise naturally, and might play a role in many areas of physics, biology and technology, such as nonlinear optics [3,4], gene regulatory circuits [5], population dynamics [6,7], traffic flows [8,9], neuroscience [10], and social or communication networks [11,12].
A well established effect of a delay in the dynamics is the possibility to induce multistability [13,14]. In oscillatory systems a delay gives rise to coexistent periodic orbits with different frequencies [15][16][17] and possibly different oscillation patterns [18][19][20][21][22]. Such coexistent patterns could be related to memory storage and temporal pattern recognition, especially in neural networks [23][24][25]. However, noise, which is unavoidable in real networks, can place important limitations to the capacity of a memory element, as it can induce mode hoppings between coexistent attractors.
We study the statistical properties of such mode hoppings in small networks of oscillators.
We consider a single oscillator with delayed feedback, two delay-coupled oscillators and a unidirectional ring, and we briefly discuss globally coupled elements. The number of possible frequencies scales with the roundtrip delay time, but the noisy system visits only a fraction of these frequencies, which scales with the square root of the delay time. While without noise the range of frequencies also scales with the coupling strength, we find that in the stochastic system it does not depend on the coupling strength. In contrast, the robustness of the orbits to noise, measured by the average residence time, is mainly determined by the coupling strength, while the delay has a minor effect. Complementary to local stability analysis, the study of coupled systems subject to noise also provides information about the robustness of certain oscillation patterns. We find that depending on network topology, an oscillation pattern might dominate: in unidirectional rings the oscillators spend equally much time in all the possible phase configurations, a globally coupled network shows a clear preference to in-phase synchrony. This paper is organized as follows. In section II we discuss stochastic switching for a single phase oscillator. We compare our numerical results to those obtained by a potential model [26], and discuss the model in the limit of strong coupling and large delay. We discuss stochastic switching of two coupled phase oscillators in section III, and extend the potential model. In section IV we extend our results to a unidirectional ring of delaycoupled oscillators. Finally, we demonstrate the generality of our results with delay-coupled FitzHugh-Nagumo oscillators in section V. We discuss our results in section VI.

II. STOCHASTIC SWITCHING IN A SINGLE PHASE OSCILLATOR WITH FEEDBACK
The most basic delay network is a single oscillator with delayed feedback. We consider a Kuramoto oscillator, which describes the oscillating dynamics by a single phase variable.
It is a universal model, as many oscillators can be reduced to phase oscillators in the weak coupling regime [27][28][29]. Thanks to its simplicity, the Kuramoto model allows for analytical insights while still capturing many essential features of synchronization. A Kuramoto oscillator with delayed feedback and noise is modelled bẏ The oscillator has a natural frequency ω 0 , the other parameters are the coupling delay τ , the coupling strength κ > 0 and the coupling phase θ. The system is subject to additive Gaussian white noise ξ(t), with ξ(t)ξ(t 0 ) = 2Dδ(t − t 0 ). As the dynamics is invariant under a transformation φ(t) → φ(t) +ωt, ω 0 → ω 0 +ω, θ → θ −ωτ , we can omit the coupling phase θ without loss of generality.
We first briefly discuss the deterministic dynamics of this system [30,31]. Without noise, the oscillator resides in one of the frequenciesφ = ω k given by A graphical determination of the frequencies ω k is shown in Fig. 1. The orbits for which κτ cos(ω k τ ) > 1 holds, are stable. For large coupling or long feedback delay κτ ≫ 1, the stable frequencies close to ω 0 are approximated as ω k τ ≈ 2kπ, whereas the spacing is given by ω k+1 − ω k ≈ 2π/τ . As all solutions of Eq. (2) are limited by ω 0 − κ ≤ ω k ≤ ω 0 + κ, the number of coexistent stable orbits is estimated as κτ /π.
If we add noise to the system, the oscillator switches between these coexistent orbits. We  (2)). The intersections with the thick decreasing slopes of the sine function correspond to stable orbits, and are marked with a circle. The coloring of the circles relates to the probability distribution p(ω(t)) of the corresponding stochastic oscillator (shown in Fig. 3): the probability that the oscillator has a frequency ω(t) ≈ ω k is large for the most central frequencies ω k ≈ ω 0 , marked with a black circle, while the probability to find the system's frequency ω(t) close to the outer frequencies ω k ≈ ω 0 ± κ, marked with an empty circle, is negligible. Parameters are ω 0 = 6, κ = 2, τ = 10 and D = 0.5 ω 0 = 6, τ = 10), without noise, the oscillator has six stable periodic orbits, with respective frequencies ω 1 ≈ 4.48, ω 2 ≈ 5.07, ω 3 ≈ 5.67, ω 4 ≈ 6.27, ω 5 ≈ 6.87 and ω 6 ≈ 7.46, shown in Fig. 1. A typical timetrace of the phase evolution, with multiple mode hoppings between ω 3 and ω 4 , is shown in Fig. 2(a). As an indicator for mode hoppings we use the frequency measure ω(t) = (φ(t)−φ(t−τ ))/τ , which is the driving term of the dynamics, and corresponds to the average of the instantaneous frequencyφ(t) over the past delay interval.
Moreover, this definition of ω(t) respects the origin of the frequency locking, which lies in the auto-phase locking of the instantaneous phase φ(t) onto the delayed phase φ(t − τ ). The time evolution of ω(t) is shown in Fig. 2(b), exhibiting clear jumps between the deterministic frequencies ω k .
The distribution of frequencies p(ω(t)), with ω(t) defined as above, is shown in Fig. 3(a).
One can clearly distinghuish multiple maxima, corresponding to the deterministic frequencies ω 2 , ω 3 , ω 4 and ω 5 . The frequencies closest to the eigenfrequency ω 0 of the oscillator, i.e. To calculate the residence times of the orbits, we apply the following procedure: at the starting point t 0 the oscillator is considered to reside in the orbit with a frequency ω k for which the distance |ω(t 0 ) − ω k | is minimal, and it stays there as long as |ω(t) − ω k | < ǫ.
After a transition, we determine the new locking frequency again as the frequency at minimal distance. We chose ǫ = 2/3(ω k − ω k−1 ); for weak noise ω(t) does not show large fluctuations around the locking frequency ω k and the choice of ǫ does not largely affect the residence times of the orbits. In our simulations we obtained around 10 6 transitions. The residence time distributions of two of the orbits (ω 2 and ω 3 ) are shown in Fig. 4(a). The distribution is exponential. Upon the exponential decay there are signatures of the delay time; these are shown in the inset. The peaks can be understood as delay echoes which result from a known stochastic resonance effect in delay systems [13,32,33]: A mode hopping causes a perturbation, which increases the probability for a mode hopping at multiples of the feedback delay. Moreover, the average residences times, shown in Fig. 4(b), are largest for orbits ω 3 and ω 4 with a frequency close to the natural frequency ω 0 .
In order to interpret the mode hopping dynamics, we approximate the delay system by an undelayed system. It is then possible to define a Langevin equation and to compute the frequency distributions and average residence times of the different periodic orbits.
Such an approach is possible thanks to the simplicity of the Kuramoto oscillator, as the dynamics of the oscillator is only characterized by a frequency. A similar method has been suggested in the context of mode hopping between external cavity modes in a single laser with delayed feedback [26,34]. Using this approximation, we show analytically how the frequency distribution and average residence times scale with the feedback strength, delay, and frequency of the orbit. Thereby we focus on the regime κτ ≫ 1, in which a multitude of orbits coexists.
In order to simplify the system, we first rewrite the dynamics in terms of the delay phase difference The main step is the following: We approximate the instantaneous frequencyφ(t − τ ) by  the frequency averaged over the future delay interval plus its noise sourcė Such assumption is justified for weak noise, when the oscillator resides in one of the periodic orbits during a delay interval. But also in case of a random walk (κ = 0) it leads to the correct stationary distribution. In this way we obtain a closed equation without delay for the phase difference x(t), that can be written in terms of a potential [26], . As the noise sources ξ(t) and ξ(t − τ ) are uncorrelated, the simplified oscillator is effectively subject to a magnified noise strength of ξ 2 (t) = 4D. The approximation by white noise in Eq. (5) does not preserve correlations around multiples of τ , like those shown in Fig. 4(a). The potential V (x) is shown in Fig.   5. It is a parabolic potential modulated by a cosine function. The local minima x k = ω k τ correspond to the frequencies in the noise-free case D = 0. Our reduction procedure does not affect them and their calculation by the potential extrema reveals Eq. (2). The local maxima x m correspond to unstable solutions of the deterministic system. The phase difference x(t) relates in a simple way to the frequency measure x(t)/τ = ω(t).
Hence, the stationary distribution of frequencies p(ω) is given by a Boltzmann factor [35] We recognize a Gaussian envelope with mean ω 0 and variance σ 2 = 2D/τ . This envelope corresponds to the probability distribution of a random walk. Thus, while the total frequency range is given by 2κ, the range of visited frequencies scales with D/τ . As the spacing between the orbits scales inversely with the feedback delay, the number of attended orbits grows as √ Dτ .
The coupling function, which appears in the second factor of Eq. (6), determines the location and the shape of the different peaks. As the feedback strength increases, the peaks in the distribution become more pronounced. In Fig. 3 we compare our analytical result for the simplified system (Eq. (6)) with numerical simulations of the original delay system (Eq. (1)). We find that our theoretical results provide an excellent approximation for the distribution of frequencies and thus prove the validity of the applied reduction method.
Also the average residence times can be approximated by the potential model (5): In the limit of low noise, the escape rates from an orbit with frequency ω k are given by the Kramers rate [36,37] where the suffix denotes whether the oscillator hops to a mode with a higher or a lower frequency. The average residence time T 0 (ω k ) reads then .
For strong coupling and large feedback delay κτ ≫ 1, a multitude of orbits are stable, with ω k τ ≈ 2nπ and x m = (2n + 1)π. The average residence time is then further approximated We compared the average residence time of the different periodic orbits with our theoretical result (Eq. (7)) in Fig. 4, and the approximation gives good results. Consequently, the average residence time T 0 (ω k ) increases with the feedback strength κ, which determines the depth of the potential wells, and decreases with the noise strength D. For a fixed frequency ω k the feedback delay τ has a limited influence on the residence times, for long delays the delay dependency even vanishes.
Only the orbits with a frequency close to the natural frequency have a considerable average residence time, and are in this sense robust to noise. This range of these frequencies scales with D, and does not depend on the delay time or the coupling strength. Due to the frequency difference of 2π/τ , the number of orbits that is robust to noise scales approximately as Dτ . Moreover, there is a difference in the mode hopping behavior at long and short delay times: For long delays, the D/τ -range of attended orbits is much smaller than the range of robust frequencies, so that all visited orbits have a similar average residence time. If the delay is shorter, as it is the case for our choice of parameters, significant differences in the residence times of the orbits are observed.

III. TWO MUTUALLY COUPLED PHASE OSCILLATORS
More common than a single oscillator driven by its own delayed feedback are coupled oscillators. We consider here the simple case of two mutually delay-coupled oscillators with independent noise sources. This system is modelled bẏ with ω 01,02 = ω 0 ∓ ∆/2, and ∆ being the detuning between the oscillators. We repeat first the case of identical oscillators (ω 01 = ω 02 ≡ ω 0 ) without noise (D = 0) [18,31]. The The frequencies of the in-phase orbits are identical to the single feedback system; they are given by For the anti-phase orbits the frequencies can be found by solving ω k = ω 0 + κ sin(ω k τ ).
The coupled system thus has twice as many coexisting periodic orbits as the single system.
In-phase orbits are stable for cos(ω k τ ) > 0 and anti-phase orbits for cos(ω k τ ) < 0. A graphical determination of the frequencies is shown in Fig. 6(a): stable in-phase and antiphase frequencies alternate each other. For large κτ ≫ 1 the frequencies ω k close to the natural frequency ω 0 are approximated as ω k τ ≈ nπ, so that the separation between the frequencies approaches π/τ . Without noise, nonidentical oscillators still synchronize to a common frequency if the coupling is strong enough |∆| < 2κ [15]. Detuned oscillators, however, are no longer exactly in-phase or anti-phase with each other, but they exhibit a phase difference δ depending on the locking frequency and the detuning. We find for the frequencies ω k and the phase difference δ if the conditions cos(ω k τ + δ) > 0 and cos(ω k τ − δ) > 0 hold, an orbit is stable. We solve Eq.
(9) graphically in Fig. 6(b). For nonidentical oscillators the frequency range is reduced to 2κ − ∆; for large κτ however neither the locking frequencies ω k nor their stability is largely affected by the detuning. to the relative probability that the in-phase (anti-phase) orbit is visited by the corresponding stochastic system, darker labeling corresponds to a higher probability to find a frequency ω(t) ≈ ω k .
The corresponding probability distributions are shown in Fig. 3 (b,c). In panel (b) a stable orbit is labeled as in-phase if −π/2 < δ < π/2. Just like for a single feedback oscillator, the frequencies ω k ≈ ω 0 are most often visited, the width of the frequency distribution is however smaller. Parameters are ω 0 = 6, κ = 2, τ = 10 and (b) ∆ = 0.8.
We show the phase evolution of two identical delay-coupled oscillators in Fig. 7. Mode hopping happens in two stages: if one oscillator, the leader, changes its frequency, the other oscillator, the laggard, follows a delay time later. Looking at the evolution of the driving terms φ 1,2 (t) − φ 2,1 (t − τ ), shown in Fig. 7(b), it is clear that during a transition the driving term of the leader changes with 2π, while the laggard changes its frequency without a phase jump in its drive. As a frequency measure for the coupled system we use the mean frequency of the two oscillators averaged over the past delay interval ω(t) =  ism. We rewrite the system as a function of the driving terms x 1 (t) and x 2 (t), defined as x 1,2 (t) = (φ 1,2 (t) − φ 2,1 (t − τ )). We then assume that the oscillators are locked to the same fixed frequency over the delay interval, and as such, thatφ 1 (t − τ ) andφ 2 (t − τ ) only differ in the contribution of the noise. This leads to the main reductioṅ In this way we can rewrite the system as a function of a twodimensional potential: with x 0 = ω 0 τ andξ 1,2 (t) = ξ 1,2 (t) − ξ 2,1 (t − τ ). This potential is shown in Fig. 8. The wells are located at (x 1 , x 2 ) = (ω k τ + 2nπ − δ, ω k τ − 2nπ + δ). The frequency of the system is then given by the average frequency ω = (x 1 + x 2 )/(2τ ). As the phase difference between detuning. The arrows indicate the two pathways for a transition between two frequencies, thicker arrows correspond to more probable pathways. Parameters are ω 0 = 6, κ = 2, τ = 10 and (b) the oscillators is only determined upon a multiple of 2π, the potential is 4π-periodic with respect to x 1 − x 2 = x A . For identical oscillators there are thus two equally probable pathways for a transition: x 1 changes with almost 2π, while x 2 remains almost constant, and φ 1 (t) leads the dynamics, and vice versa. These pathways are indicated by arrows in Fig. 8. Transitions typically take place between orbits with a minimal frequency difference, and therefore with a different oscillation pattern. If the oscillators are identical, we obtain the frequency distribution p(ω) by integrating over the phase difference x A . We find with I 0 (y) being the modified Bessel function of the first kind, I 0 (y) = y 2n 2 2n (n!) 2 . Like for the single oscillator, the frequency distribution can be written as a Gaussian envelope multiplied with a factor determining the separate peaks. However, the variance of the envelope decreases with a factor 1/2 compared to the single feedback system. The Bessel function I 0 (y) is symmetric: we find alternating peaks corresponding to in-phase and anti-phase orbits, their height only depends on their respective frequencies ω k and not on the oscillation pattern. We compare our numerical and theoretical results for the frequency distribution in Fig. 3(b). The agreement is excellent. The oscillators are thus always synchronized (except during the delay interval following a transition), but they spend a proportion of time in in-phase as well as in anti-phase orbits. As a result, for long enough delays, the overall correlation between the oscillators vanishes at zero lag, but shows maxima at odd multiples of the coupling delay.
For nonzero detuning ∆ > 0, the potential (Eq. (11)) is tilted, as is shown in Fig.   8(b). Consequently the phase difference x A (t) between the oscillators preferentially increases during a mode hopping. The most probable, and the least probable transition pathway between two frequencies are also sketched on Fig. 8(b). The ratio between the transition rates is approximated by exp(−pi * Delta/2D), so that for large detunings it is reasonable to assume that all the transitions to a higher frequency are induced by the faster oscillator x 2 , and the transitions to a lower frequency by the slower one x 1 . For κτ sufficiently large, we can approximate the envelope by assuming detailed balance This corresponds to a Gaussian envelope of the frequency distribution with mean ω 0 = (ω 1 +ω 2 )/2 and variance σ 2 = D/(τ (1−ǫ)), with ǫ = 2 arcsin(∆/2κ)/π > 0. The distribution of frequencies thus becomes broader due to the detuning, in agreement with the numerical results for the full delay system. In Fig. 3(c) we show the approximated Gaussian envelope together with the simulated distribution of frequencies.
For identical oscillators, we approximate the residence times of the orbits by assuming all the transitions take place via the two optimal pathways. We obtain then for the mean residence time This corresponds to half of the lifetime of the orbits of a single oscillator with the a roundtrip delay 2τ . We compare numerical and theoretical results in Fig. 4(b).
It is possible to extend these results to a unidirectional ring of N oscillators. Such system is then modelled byφ with N + 1 ≡ 1. Without detuning, the coupling topology allows for in-phase oscillations φ n (t) = ω k t and several out-of-phase oscillation patterns φ n (t) = ω k t + n∆φ, with ∆φ = 2mπ/N. The corresponding frequencies are given by ω k = ω 0 −κ sin(ω k τ −∆φ), and they are stable if cos(ω k τ − ∆φ) > 0 [18]. This results in alternating orbits with a different oscillation pattern. For strong coupling and long delay, the frequencies are separated by 2π/(Nτ ).
, and assuming that the instantaneous frequencies of the oscillators can be approximated by the mean frequency averaged over the delay interval and their noise source,φ we find an N-dimensional potential x l . The frequency of the system is then measured by ω(t) = x S (t)/τ . It is no longer possible to compute the frequency distribution p(ω) in terms of simple analytical expressions as above. However, for identical oscillators (Delta n = 0) it is straighforward to see that the parabolic term in Eq. (16) leads to a Gaussian envelope. The variance of this envelope is given by σ 2 = 2D/(Nτ ), and thus scales inversely with the total roundtrip delay Nτ . As the frequency difference between the orbits is approximated by 2π/(Nτ ), the number of attended orbits scales as √ Nτ . Moreover, the potential is symmetric with respect to the different oscillation patterns, so that each pattern is equally often visited in the long delay limit.
For low noise, zero detuning and large κτ , we find that the average residence times scale inversely with the number of oscillators in the ring, and depend weakly on the total roundtrip delay. They are approximated by We compared the frequency distributions and residences times of the simplified non-delay system with simulations of three, four and five delay-coupled oscillators, and the agreement is excellent (not shown).

V. GENERAL PERIODIC SYSTEMS WITH DELAYED COUPLING
In order to investigate whether our results are valid in a broader context, we compare the switching behavior of other nonlinear delay-coupled oscillators to our results for phase oscillators. The Kuramoto model is a weak-coupling limit, which only describes the phase dynamics, and does not take any influence on the amplitude into account; therefore we expect that our results mainly apply for weak coupling.
First, we sketch the deterministic periodic solutions in a general delay system. For a single oscillator, it is known that a feedback delay induces coexisting periodic orbits, with a frequency separation of 2π/τ [16]. We show here briefly that in a unidirectional ring of identical oscillators, a delay gives rise to alternating in-phase and out-of-phase orbits, in a similar way as for phase oscillators. For general limit cycle systems, unlike for phase oscillators, it is not so straightforward to determine the respective orbits and their stability properties.
Extending the approach of Yanchuk and Perlikowski for a single feedback system [16], we consider a set of N identical nonlinear systems coupled in a unidirectional ring with delaẏ where x N +1 ≡ x 1 . For the following we assume, that this network allows for an in-phase synchronized periodic solution x n (t) = x n−1 (t) = x n (t + T ) at a coupling delay τ = τ 0 .
Shifting the delay to τ 1 = τ 0 + T /N, the same periodic orbit is a solution of the system, the oscillators however exhibit a phase difference x n (t) = x n+1 (t − T /N) = x n (t + T ). Similarly, we find the same waveform appearing with all the other out-of-phase patterns that are allowed in the ring: a pattern corresponding to x n (t) = x n−1 (t − kT /N) = x n (t + T ) can be found at a delay τ k = τ 0 + kT /N. An orbit with a period T thus reappears when shifting the delay by an amount T /N.
The periodic solutions are organized in branches: as the delay increases, the period T of an orbit varies continuously between a minimal period T min and a maximal period T max . For a fixed delay τ , the number of coexistent orbits resulting from a single branch can then be estimated in the following way: we have τ ≈ nT max /N ≈ mT min /N. The number of periodic states is then estimated as m − n = Nτ (T −1 min − T −1 max ), with in-phase and out-of-phase orbits alternating each other. The frequency difference between two orbits is approximated by 2π/(Nτ ), just like for phase oscillators. It is possible to show that the stability of these orbits depends on their period, but not on the oscillation pattern. In the long delay limit the stability no longer depends on the number of oscillators in the ring, or the coupling delay.
As an examplary system, we investigate numerically stochastic switching between such coexistent orbits in FitzHugh-Nagumo oscillators. We simulated a single oscillator (N = 1) with delayed feedback, and two identical mutually delay-coupled oscillators (N = 2). Our oscillator is modelled by with (v N +1 , u N +1 ) ≡ (v 1 , u 1 ), and ξ n (t) being Gaussian white noise with a variance given by ξ 2 (t) = 2D. We choose our parameters so that without delayed coupling and without noise the oscillator(s) show periodic spiking dynamics. A typical timetrace of an oscillator with noise and feedback, which performs a mode hopping, is shown in Fig. 9(a).
We analyze the mode hopping in a similar way as for phase oscillators. We define the phase of the oscillators by the Hilbert-transform of the fast variable, φ n (t) = arg(H(v n (t)), but similar results were obtained by using the alternative definition φ n (t) = arctan(v n (t)/w n (t)).
In both cases the waveform is very different from sinusoidal, so the frequency shows large fluctuations within a period even without noise. The frequency measure is given by ω(t) = (φ n (t) − φ n (t − τ ))/(Nτ ), similar as for Kuramoto oscillators. We show the evolution of the frequency ω(t) in Fig. 9(b). Although the mode hopping event is hard to detect in the timetrace ( Fig. 9(a) asymmetric excursions from the deterministic frequency. Fig. 10 (a,b) compares the frequency distributions p(ω) of a single oscillator with feedback, and two mutually coupled oscillators. For the single oscillator, shown in Fig. 10(a), we find five different peaks, separated by a frequency difference of 2π/τ . The shape of the different peaks is asymmetric; this feature results from the asymmetric waveform of the spikes. The frequency distribution p(ω) for two mutually coupled oscillators is shown in Fig.   10(b): we find the peaks at the same frequency as for the single element; they correspond to in-phase orbits. Between the in-phase peaks we find maxima that can be associated to anti-phase orbits. Just like for phase oscillators, the frequency distribution for the coupled system has the same mean, and half of the variance as the single system. The corresponding average residence times of each orbit are shown in Fig. 10(c). The average residence times of the single oscillator (black dots) are larger than those of the coupled system (pink dots).
Moreover, the average residence times show the same trend as for coupled phase oscillators: the orbits with a central frequency are most robust against noise.
The agreement with the Kuramoto oscillators is even quantitative: In Fig. 10(a,b)   We also compared in Fig. 10(c) the Kuramoto residence times for a single oscillator (Eq. (7), pink dashed curve) and for two coupled oscillators (Eq. (14), blue dashed curve) to the residence times for FitzHugh-Nagumo elements. We thereby used the parameters D and ω 0 determined from the frequency distributions, the coupling strength κ can then be estimated from the average residence times. We find that, for the same coupling strength κ for the single and the coupled system, the residence times are well approximated by the phase model. Moreover, we find a single parameter set (ω 0 , κ, D, θ) which models the frequency distribution and the average residence times for both a single feedback and two coupled oscillators. Hence, the scaling properties of the stochastic periodic dynamics with the delay time and the oscillator number are reproduced.
We have studied the influence of additive noise on one, two and a ring of phase oscillators coupled with delay. In such systems multiple periodic orbits coexist, and under influence of noise the oscillators hop from one orbit to another. We approximated the system as a noisy particle in potential well; both for the distribution of frequencies as for residence times the approximation is excellent. Although our approximation only applies for weak noise, we obtain a good agreement for the distribution of frequencies for all noise strenghts.
However, it should be remarked that in the simplified model some dynamical phenomena are not reproduced. The most prominent example are the delay stochastic resonance peaks in the residence time distribution, shown in Fig. 4(a). Also frequency oscillations with a periodicity of a roundtrip delay time, which are typically present in the system, are no longer visible. Transient behavior is different as well: in the delay system a transient decays at a rate proportional to the inverse delay time, while in the reduced system the decay happens much faster.
We found that the oscillators only visit a fraction of the deterministic stable orbits: whereas the number of deterministic orbits scales with oscillator number N, delay time τ and coupling strength κ, the number of visited orbits scales as √ DNτ , and does not depend on the coupling strength κ. The orbit with a frequency closest to the natural frequency is the most probable, irrespective of its oscillation pattern.
Our results on the average residence times indicate the robustness of the orbits against weak noise. The most robust orbits are those with a frequency close to the natural frequency, also irrespective of the oscillation pattern. The sensitivity of an orbit against noise depends strongly on the coupling strength, the coupling delay plays only a minor role. The number of robust orbits scales as DNτ .
For two delay-coupled oscillators, and for unidirectional rings the systems does not show any preference for a particular oscillation pattern. The different oscillation patterns are equally often attended, in the long delay limit. However, this symmetry between in-phase and out-phase patterns depends on the coupling topology. We also simulated three, four and five delay-coupled Kuramoto oscillators in an all-to-all configuration. In this case however a description as a noisy particle in a potential is not accurate, as not only periodic dynamics is observed. The distribution of frequencies looks different: the peaks associated to in-phase orbits are considerably higher than those corresponding to out-of-phase dynamics. For long delays even only in-phase periodic orbits are visible in the frequency distributions. We find that the frequency distributions narrows with the number of oscillators: the variance of the frequency distribution scales as 1/(N − 1).
The Kuramoto model is a weak coupling approximation for limit cycle oscillators. Therefore we expect our results to apply for delay-coupled nonlinear systems showing stable periodic dynamics. In particular, the Kuramoto approximation applies when the coupling mainly influences the oscillation phase, while the waveform or oscillation amplitude is hardly affected. We found indeed a good correspondence between Kuramoto and FitzHugh-Nagumo oscillators in this case. However, we expect the approximation to break down as the coupling strength increases and amplitude instabilities play a role in the dynamics.
Not only in stable oscillatory systems, but also in a chaotic attractor a delay has the effect of inducing multiple periodic orbits. Hence, the chaotic attractor of two delay-coupled chaotic systems contains in-phase as well as anti-phase orbits, and they have similar stability properties (for long enough delay). Therefore, it is not surprising that we find the same correlation pattern, with a high correlation at the delay time, but no correlation at zero lag, for coupled noisy oscillators and chaotic systems with delay [38]. However, chaotic and stochastic systems show different scaling behavior with the delay time and the number of coupled elements.
It is worth noting that the envelopes of the frequency distributions are the same as those for a random walk. The delayed feedback only imposes restrictions on the distribution of the two point distribution of x(t) = φ(t) − φ(t − τ ), but it does not affect the envelope. On timescales much shorter than the delay t 0 ≪ τ , the influence of the feedback is even not visible: the two point distribution of φ(t) − φ(t − t 0 ) is identical to the one of a random walk.
A possible explanation of this surprising phenomenon lies in the fact, that the equations of motion do not impose any restrictions on this phase difference, as long as t 0 is different enough from τ . Hence, the random-walk can explore the possible range. However, on timescales equal or larger than the delay, the dynamics (i.e. the timetrace) of an oscillator with delayed feedback differs significantly from a random walk. Also the two point distributions show a clear fingerprint of the delay time, and for t 0 > τ , a larger variance than a random walk. We believe that the issue of two-/N-point distributions in delay systems is worth being studied in more detail.