# Synchronous solutions and their stability in nonlocally coupled phase oscillators with propagation delays

## Abstract

We study the existence and stability of synchronous solutions in a continuum field of non-locally coupled identical phase oscillators with distance-dependent propagation delays. We present a comprehensive stability diagram in the parameter space of the system. From the numerical results a heuristic synchronization condition is suggested, and an analytic relation for the marginal stability curve is obtained. We also provide an expression in the form of a scaling relation that closely follows the marginal stability curve over the complete range of the non-locality parameter.

###### pacs:

05.45.Ra, 05.45.Xt, 89.75.-k## I Introduction

Phase oscillators are simple models of oscillatory activity where the
essential variable is an angular position along a limit cycle. Systems of
coupled phase oscillators provide a convenient paradigm for studying the
collective behavior of weakly interacting limit cycle oscillators
(1); (2); (3)
and are in turn useful for modeling a wide variety
of cooperative phenomena
observed in biological, physical and chemical systems (4); (5); (6); (7); (8).
One of the most well known models of coupled phase oscillators is the Kuramoto
model (2) which considers a system of globally
(all-to-all) coupled set of oscillators with a spread of intrinsic frequencies.
The oscillators synchronize to a common frequency as soon as the coupling strength
goes beyond a threshold value. For a system of identical oscillators one achieves phase
synchrony where, starting from an initial
random distribution of phases, all the oscillators acquire the same phase within a time
that is inversely proportional to the coupling strength. When the members of the set of
identical oscillators have a
spatial identity (e.g. when the coupling is of the nearest-neighbour kind or varies
in strength over space) the collective states can exhibit a broader class of phase-locking
behaviour including acquiring the same phase (synchrony) or a fixed amount of phase difference between
neighbouring oscillators (traveling waves). For a nonlocally coupled continuous system of this kind it
is also possible to sustain a more exotic state whereby phase-locked and incoherent activity
can simultaneously exist at different spatial locations (9),
giving rise to a
spatio-temporal pattern that has been termed as a *chimera state *
(10).

While the original Kuramoto model is based on the idealized assumption of instantaneous coupling between the oscillators, there has been a growing interest in extending the model to include finite time delay in the coupling and to study the existence and stability of phase-locked states in such a situation (11); (12); (13); (14); (15); (16); (17). Time delayed coupling is physically relevant in many real life systems and accounts for finite propagation times of signals, reaction times of chemicals, individual neuron firing periods in neural networks etc. One of the earliest investigations in this direction was made by Schuster and Wagner (11) who considered the collective states of a pair of delay-coupled phase oscillators and observed changes in the onset conditions for synchrony as well as the existence of additional (higher) branches of synchronized frequency states. Niebur et al. (12) considered a two-dimensional array of phase oscillators that interact via time delayed nearest neighbour coupling and found that time delay led to a reduction in the collective frequency of the system. Their numerical studies also found higher frequency states which were metastable and decayed to the lowest branch in the presence of noise. Kim et al (13) studied multistability of synchronized and desynchronized states in a time-delayed Kuramoto system in the presence of a pinning force. More recently Yeung and Strogatz (15) have carried out a systematic stability study of the synchronous states of a delayed Kuramoto system and have obtained analytic stability curves for the limiting case of coupled identical oscillators. Since the original Kuramoto model adopts a global coupling between the oscillators, its time-delayed version considered in all the above mentioned studies have assumed a single fixed amount of delay in the coupling. In a more realistic situation where the oscillators are spatially distributed and are coupled in a nonlocal fashion, such an assumption is not appropriate and one needs to consider distance-dependent delays.

Models with distance-dependent delays have been proposed in the past and analyzed under simplifying assumptions or in convenient limits. Crook et al. (14) considered a continuum system of coupled identical oscillators with a spatially decaying interaction kernel and modeled the space dependent time delay contribution through an effective phase shift term in the interaction. They also took the size of the system to be infinite for mathematical convenience but thereby effectively reduced the nature of the mutual coupling to a “local” one (since the interaction length is always much smaller than the system length). Zanette (16) studied another simplified version of the generalized system where he adopted a distance independent (global) coupling between the oscillators but introduced a distance dependent time delay in the interaction. He obtained numerical results on the stability of synchronous and propagating traveling waves and also some analytic results in the limit of small delay. Ko & Ermentrout (18) recently investigated the effects of distance dependent delays in sparsely connected oscillator systems and found that a small fraction of connections with time delay can destabilize the synchronous states. We note that distance-dependent delays have also been considered in other related contexts, e.g., in the stability of equilibria of continuum models of neural tissue with axonal propagation delays (19); (20). To the best of our knowledge, however, there has been no systematic study of the synchronization properties of the full nonlocal system of phase oscillators with distance-dependent delays. Our present work is motivated by a desire to address this problem. The presence of both non-locality and time delay in the coupling pose a formidable mathematical challenge for an analytic solution of the problem and one has to perforce adopt a numerical approach. In the recent past we have studied the generalized system to obtain chimera states in the presence of delay (21). In the present work we carry out a detailed study on the existence and stability properties of synchronous solutions of the system. Through detailed stability curves in the parameter space of the mean time delay () and the non-locality parameter () (for a fixed intrinsic oscillator frequency ) we delineate the relative influence of and on the stability of the synchronous state whose frequency is . From an extensive scan over several frequencies , we further obtain a comprehensive stability diagram in the parameter space of and that is independent of . Our numerical results further suggest that the instability mechanism corresponds to a saddle-node bifurcation and the most unstable perturbation has the lowest mode number. These results help us obtain an analytic expression for the synchronization condition. The synchronization condition can be closely approximated by an offset exponential power law relationship between and , which is valid over a wide range of coupling regimes and provides a useful criterion for determining the critical value of the delay-induced maximum phase shift below which the synchronous state is stable.

## Ii Model system and its synchronous states

We consider a continuum of identical phase oscillators, arranged on a circular ring and labeled by with periodic boundary conditions, whose dynamics is governed by

(1) |

where is the phase of the oscillator at location and time , whose intrinsic oscillation frequency is , is the coupling strength and is an even function describing the coupling kernel. The quantity denotes the signal propagation speed which gives rise to a time delay of for distance from the locations . As the oscillators are located on a ring with circumference , the distance between any two oscillators is given by the shorter of the Euclidean distance between them along the ring. In this configuration, the maximum distance between the coupled oscillators is and thus the maximum time delay would be .

We choose the coupling kernel to have an exponentially decaying nature and its normalized form is taken to be

(2) |

where is the inverse of the interaction scale length and is a measure of the nonlocality of the coupling. We make time and space dimensionless in Eqs.(1) and (2) by the transformations , , , , and and obtain

(3) |

(4) |

We look for synchronous solutions of Eq. (3) that have the form:

(5) |

The value of can be taken to be zero by a translation. Substituting Eq. (5) in Eq. (3) we get

(6) |

which is an implicit equation in . Being a transcendental equation its solution can in principle be multi-valued in for a given set of parameters and , and can lead to higher branches of collective frequencies as pointed out by Schuster and Wagner (11) for a system of two coupled oscillators.

We further define a mean delay parameter by

(7) |

which weights the individual delays with the corresponding connection weights. With the exponential connectivity given by Eq. (4), this translates into

(8) |

This gives values for the limiting cases: for local () and for global () coupling.

Fig. 1 plots the numerical solutions of Eq. (6) as a function of for and several values of . As the curve for shows, it is possible to have multiple solutions for a given value of . The stability of these higher states will be discussed in later sections of the paper. One also notes that the lowest branch shows frequency suppression as a function of the mean time delay .

We can also examine two simple limits of Eq. (6) which help us understand the influence of the non-locality parameter on the nature of the equilibrium solutions. For , corresponding to the global coupling limit, Eq. (6) reduces to,

(9) |

This was the limiting case examined by Zanette (16), where the transcendental nature of Eq. (6) was retained, thereby permitting several multiple branches of solutions to exist, but the coupling was taken to be global. In the limit of , (9) gives,

(10) |

demonstrating frequency suppression of the collective mode due to finite time delay in this parameter regime. The other interesting scenario is that of , which corresponds to the situation of local coupling among the oscillators. In this case Eq. (6) becomes

(11) |

Eq. (11) is a third order polynomial equation in and can at most have three real solutions in contrast to the higher number of multiple roots of the transcendental Eq. (6). The local limit can also be approached by taking limit in the unscaled form of Eq. (6) in which is replaced by and by . The mean delay () in this system equals and Eq. (6) can be rewritten in terms of as

(12) |

Crook et al. (14) considered such a limit but with an additional simplifying assumption regarding the contribution of the delay term: Instead of explicitly treating the delay term in the argument of the phase function they modeled the delay contribution by a space dependent phase shift (implicit delay). Thus the basic model Eq. (1) was simplified in their case to be of the form

(13) |

with . The corresponding equation for the equilibrium solutions is given in this case by

(14) |

which shows only a single collective frequency that experiences time delay induced reduction for and a slow rise to for .

The degree to which the non-locality parameter can influence the frequency of the synchronous equilibrium solutions of the coupled oscillator system can be seen in Fig. 2. This figure provides a graphical comparison of this influence by displaying plots of vs obtained by solving Eq. (6) with different values of . The curve labeled NL represents a typical case of nonlocal coupling, here calculated for . The curve L is for local coupling with , whereas the global coupling curve G is calculated for . The and curves obtained from Eq. (9) and Eq. (11) respectively are not significantly different from what are shown in this figure. The curve plots Eq. (14) for comparison. We have used for the plots shown in this figure.

## Iii Stability of the synchronous solutions

We now examine the stability of the synchronous solutions obtained in the previous section.

### iii.1 Eigenvalue equation

The linear stability of solutions of Eq. (3) is determined by the variational equation

(15) |

where . With the ansatz , , , we obtain the eigenvalue equation

(16) |

Writing and separating Eq. (16) into its real and imaginary parts we get

(17) | ||||

(18) |

Since the perturbations corresponding to again yield synchronous solutions, the linear stability of the synchronous state requires that all solutions of Eq. (16) have for all non-zero integer values of . The marginal stability curve in the parameter space of is defined by and in principle can be obtained by setting in Eq. (17), solving it for and substituting it in Eq. (18). In practice it is not possible to carry out such a procedure analytically for the integral Eqs. and one needs to adopt a numerical approach, which is discussed in the next section.

### iii.2 Numerical determination of eigenvalues

To systematically determine the eigenvalues of Eq. (16) in a given region of the complex plane we use multiple methods in a complementary fashion. First, Eq. (6) is solved for for a given set of values of the parameters , and . Following that, we need to determine the complex zeros of the function defined as

which is equivalent to finding solutions of (16). To do this we have primarily relied on the numerical technique developed by Delves and Lyness (22) based on the Cauchy’s argument principle. By this principle the number of unstable roots of is given by

where the closed contour encloses a domain in the right half of the complex plane with the imaginary axis forming its left boundary. Once we get a finite number for we further trace the location of the roots by plotting the zero value contour lines of the real and imaginary parts of the function in a finite region of the complex plane . The intersections of the two sets of contours locate all the eigenvalues of Eq. (16) in the given region of the complex plane. The computations are done on a fine enough grid (typically ) to get a good resolution. A systematic scan for unstable roots is made by repeating the above procedure for many values of the perturbation number and by gradually extending the region of the complex plane. We have made extensive use of Mathematica in obtaining the numerical results on the stability of the synchronous states.

### iii.3 Results

In Fig. 1 the solid portion of the curve shows the stable synchronous states of Eq. (3) for and for various values of and . The terminal point on a given solid curve of vs marks the marginal stability point. The marginal stability point is seen to shift towards larger values of as one moves down to curves with lower values of . A more compact representation is obtained if one plots versus , since in this case the solutions corresponding to different values of for a given consolidate onto a single curve, as shown in Fig. 3 for and respectively.

It is seen from both figures that the stability domains of the synchronous solutions are restricted to the lowest branch where the curves are decreasing. This suggests a heuristic necessary (but not sufficient) condition for stability of the synchronous solutions: From Fig. 1 we have that for stable synchronous solutions, and from Eq. (6) we calculate

(19) |

where

(20) |

and

Since and is bounded, the denominator in Eq. (19) is positive for small values of . Hence, for positive , the requirement implies the condition , that is,

(21) |

An alternative approach to arrive at the necessary condition (Eq.21) would be to make use of the results presented in Fig.3. We recast the dispersion relation given by Eq.(6) in the form:

(22) |

where

It is seen from Fig.3 that the stability domain of the synchronous solutions is restricted to the lowest branch where the curves have a negative slope. This again suggests a heuristic necessary condition for stability of the synchronous solutions to be leading to the necessary condition given by Eq.(21). The prime indicates a derivative of .

As we will see later, the marginal stability curve obtained from or does lie above the true marginal stability curve (see Fig. 5), confirming that condition (21) is necessary but not sufficient for the stability of synchronous states.

The points where solid and dotted lines meet in the curves of Figs. 1 and 3 mark the marginal stability point for the respective values. These points are obtained for a range of values and are plotted in () space in Fig. 4 by filled points. They all lie on a single curve, which is analytically derived below. Our numerical results further reveal that for the marginal stability points, the imaginary part of the eigenvalue of the mode is zero—in other words the mode loses stability through a saddle-node bifurcation. It can easily be checked by inspection that is one of the solutions of Eq. (18) for any value of ; however, it is not evident analytically that this is the only possible solution for , and our numerical results have helped us confirm that this is indeed the case. Hence, putting in Eq. (17) we get the following integral relation between the parameters and .

(23) |

Further, we have also observed that the most unstable perturbation is the one with the lowest mode number, namely . Therefore Eq. (23) with defines the marginal stability curve, so the condition for synchronization takes the form

(24) |

The solid line in Fig. 4 is the analytical curve of marginal stability defined by setting the left side of Eq. (24) to zero, and it can be seen that the numerically calculated marginal values (represented by points) fit this curve perfectly. The figure also shows the stability curves obtained for the and perturbations (dashed and dotted lines, respectively) and these are seen to lie above the marginal stability curve. We have carried out a numerical check for a whole range of higher numbers and the results are consistent with the above findings.

For a system with constant delay (i.e. if is replaced by in Eq. (3)), the term can be taken outside of the integral in (24) and the remaining integrand is nonnegative. Hence, the synchronization condition in this case becomes simply

(25) |

This agrees with the results obtained previously for constant-delay systems (15); (17).Thus our result, as given by Eq. (23), generalizes the condition (25) to systems with space-dependent delays, and shows a nontrivial relation between the spatial connectivity and delays for the latter case.

In order to gain some intuition into the complex interaction between connectivity and delays, we have obtained an approximate expression for the marginal stability curve by a numerical fitting procedure, yielding the relation

(26) |

for the stability of synchronous oscillations. Here, the left side involves the temporal scales of the dynamics (namely, it is the average time delay normalized by the oscillation period of the synchronized solution) while the right hand side involves the spatial scales of connectivity. In this view, the synchronization condition is a balance between the temporal and spatial scales. For high connectivity (), the system can tolerate higher average delays in maintaining synchrony, and the largest allowable delays decrease roughly exponentially as the spatial connectivity is decreased. In the same figure we have also plotted with dotted curve the approximate condition (21), which is found to lie above the marginal stability curve in the entire range of . The disparity between the two curves becomes particularly noticeable at large values of .

## Iv Conclusions and Discussion

We have investigated the existence and stability of the synchronous solutions of a continuum of nonlocally coupled phase oscillators with distance-dependent time delays. Our model system is a generalization of the original Kuramoto model by the inclusion of naturally occurring propagation delays. The equilibrium synchronous solutions of this system are shown to differ significantly from those of similar earlier models that had introduced simplifications either in the coupling or in the nature of the time delay. The equilibrium solutions of the lowest branch are seen to exhibit frequency suppression as a function of the mean time delay. We have carried out a linear stability analysis of the synchronous solutions and obtained a comprehensive marginal stability curve in the parameter domain of the system. Our numerical results show that the synchronous states become unstable via a saddle-node bifurcation process and the most unstable perturbation corresponds to an (or kink type) spatial perturbation on the ring of oscillators. These findings allow us to define an analytic relation, given by Eq. (24), as a condition for synchronization. We have also obtained approximate forms for the synchronization condition that provides a convenient means of assessing the stability of synchronous states. Our results indicate an intricate relation between synchronization and connectivity in spatially extended systems with time delays.

## Acknowledgments

GCS would like to acknowledge the hospitality of the Max-Planck Institute for the Physics of Complex Systems, Dresden, Germany, where part of the work was carried out during his sabbatical at BioND Group headed by Dr. Thilo Gross. GCS also acknowledges his discussions with Dr. Luis G. Morelli at MPI-PKS.

### References

- A. T. Winfree, J. Theor. Biol. 16, 15 (1967).
- Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics Vol. 39, edited by H. Araki (Springer-Verlag, Berlin, 1975); Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, Berlin, 1984).
- G. B. Ermentrout and N. Kopell, SIAM J. Math. Anal., 15, 215 (1984).
- T. J. Walker, Science 166, 891 (1969).
- J. Buck, Q. Rev. Biol. 63, 265 (1988).
- C. S. Peskin, Mathematical Aspects of Heart Physiology (Courant Institute of Mathematical Sciences, New York, 1975).
- S. H. Strogatz, R. E. Mirollo, and P. C. Matthews, Phys. Rev. Lett. 68, 2730 (1992).
- K. Wiesenfeld, P. Colet, and S. H. Strogatz, Phys. Rev. Lett. 76, 404 (1996); Phys. Rev. E 57, 1563 (1998).
- Y. Kuramoto and D. Battogtokh, Nonlin. Phenom. Compl. Sys., 5,380 (2002).
- D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett., 93, 174102 (2004); Int. J. Bif. Chaos, 16,21 (2006).
- H. G. Schuster and P. Wagner, Prog. Theor. Phys., 81, 929 (1989).
- E. Niebur, H. G. Schuster, and D. M. Kammen, Phys. Rev. Lett., 67, 2753 (1991).
- S. Kim, S. H. Park and C. S. Ryu, Phys. Rev. Lett. 79, 2911 (1997).
- S. M. Crook, G. B. Ermentrout, M. C. Vanier, and J. M. Bower, J. Comput. Neurosci., 4, 161 (1997).
- M. K. S. Yeung and S. H. Strogatz, Phys. Rev. Lett. 82, 648 (1999).
- D. H. Zanette, Phys. Rev. E, 62, 3167 (2000).
- M. G. Earl and S.H. Strogatz, Phys. Rev. E, 67,036204 (2003).
- T. K. Ko and G. B. Ermentrout , Phys. Rev. E, 76, 056206 (2007).
- A. Hutt and F. M. Atay, Physica D 203, 30 (2005); Phys. Rev. E 73, 021906 (2006).
- F. M. Atay and A. Hutt, SIAM J. Appl. Dynamical Systems 5, 670 (2006);
- G. C. Sethia, A. Sen and F. M. Atay, Phys. Rev. Lett., 100, 144103 (2008).
- L. M. Delves and J. N. Lyness, Math. Comput. 21, 543 (1967).