Synchronization Dynamics in the Presence of Coupling Delays and Phase Shifts

In systems of coupled oscillators, the effects of complex signaling can be captured by time delays and phase shifts. Here, we show how time delays and phase shifts lead to different oscillator dynamics and how synchronization rates can be regulated by substituting time delays by phase shifts at constant collective frequency. For spatially extended systems with time delays, we show that fastest synchronization can occur for intermediate wavelengths, giving rise to novel synchronization scenarios.

In systems of coupled oscillators, the effects of complex signaling can be captured by time delays and phase shifts. Here, we show how time delays and phase shifts lead to different oscillator dynamics and how synchronization rates can be regulated by substituting time delays by phase shifts at constant collective frequency. For spatially extended systems with time delays, we show that fastest synchronization can occur for intermediate wavelengths, giving rise to novel synchronization scenarios.

.Ks
It has been known for over three-hundred years that interacting dynamic oscillators generally tend to synchronize, even if interactions are weak [1]. This synchronization occurs robustly and independent of the details of the interaction mechanism. A simple model for the generic features of synchronization is the Kuramoto model [2,3]. It describes the phase dynamics of instantaneously coupled oscillators using a periodic coupling function. The instantaneous frequency of a given oscillator is influenced by the phase received from other oscillators. In general, coupling tends to keep the phase difference between the oscillator and the received signal at a constant value α. For α = 0, two oscillators tend to synchronize, for α = π, they tend to lock in anti-phase.
In many systems, signaling processes are complex and signaling times cannot be ignored. For example, in biological systems, dynamic oscillators are often coupled via complex molecular signaling processes [4,5]. If the processes involved take a time comparable to the oscillation cycle, these time delays in the coupling can play a significant role for the dynamics of the system and the properties of synchronized states [6]. In principle, effects of signaling times could be captured either by modifying the phase shifts α or by introducing an explicit time delay τ [7]. The effects of time delays have been studied extensively. In particular, coupling delays can lead to multistability of synchronized states and affect their collective frequency [6,8,9]. While often considered as an undesired but inevitable feature, constructive roles of coupling delays on synchronization have been reported [10,11]. It has been shown for systems with both phase shifts and time delays, that in the synchronized state, time delays effectively induce an additional effective phase shift between coupled oscillators [6], suggesting that phase shifts alone may capture the essential effects of delays. This raises the question whether phase shifts and time delays play a similar role in networks of coupled oscillators. In this Letter, we show that grad-ually substituting time delays by phase shifts, keeping the collective frequency constant, there exists a specific combination of time delay and phase shift for which the rate of synchronization is fastest. This applies both to globally coupled oscillators as well as different coupling topologies. In spatially extended systems, substituting time delays by phase shifts can regulate the length scale at which synchronization is fastest. Our results demonstrate how the phase shift α and the delay τ account for different physical effects of complex oscillator coupling.
We obtain our results using a Kuramoto model for a network of identical coupled oscillators, which takes into account the time delay τ and phase shift α [6,12], Here, θ i (t) is the phase of oscillator i, N is the total number of oscillators, ω is the intrinsic frequency of the oscillators. Oscillator coupling of strength K is described by the 2π-periodic function Γ(ϑ). The adjacency matrix a ij with a ij ≥ 0 defines the coupling topology and ρ i = j a ij is the total weight of links of oscillator i. Because of the normalization of the coupling strength by ρ i in Eq. (1), in-phase synchronized states with θ i (t) = Ωt always exist. The collective frequency Ω obeys the equation [6,13] The collective frequency thus depends on α and τ . For τ > 0, several synchronized states with different collective frequency can coexist. By simultaneously changing α and τ , it is possible to keep the collective frequency constant. For any synchronized state with collective frequency Ω = Ω 0 obeying Eq. (2), the transformation τ → τ , α → α + Ω 0 (τ − τ ) preserves the existence of the synchronized state with Ω = Ω 0 . This transformation implies that for a given value of Ω 0 τ + α, there is a one-parameter family of systems in the (τ, α)-plane that can exhibit the same collective frequency. These systems can be parameterized as where ψ is a constant that sets the collective frequency Ω 0 as Varying τ and using the phase shift α(τ ), the collective frequency does not change. However, here we will show that the synchronization dynamics does change. To study the dependence of the synchronization dynamics on time delays and phase shifts at constant frequency, we introduce a small perturbation to the synchronized state and determine its exponential relaxation rate r 0 . We compute r 0 below both analytically and from numerical simulations. Numerically, r 0 can be determined from the exponential relaxation time of the perturbation to perfect synchrony, monitored by the Kuramoto order parameter [14] [15].
As a first example, we consider globally coupled oscillators with Γ(ϑ) = sin ϑ. Fig. 1 displays the relaxation rate of the system as a function of the time delay τ , obtained by numerical integration of Eq. (1) (circles). This result shows that in this system, there is a characteristic value of the coupling delay for which the synchronization rate is maximal. The analytical solution for the synchronization rate (solid line), derived below, displays a characteristic cusp where this maximum is attained.
As a second example, we consider a system with nearest-neighbor coupling in one dimension with periodic boundary conditions, see Fig. 2. In this case, as shown below, we find that spatial Fourier modes of the oscillator lattice relax independently, each with a relaxation rate r 0 (k) that depends on the wavevector k = 2πp/N of the Fourier mode. There exists a discrete set of wavevec- considered to be even, is the system size. Note that because of the delays, there is in fact a discrete set of relaxation rates r n (k) for a given wavevector. However, synchronization is governed by the slowest rate r 0 . The relaxation rate r 0 of long-wavelength modes, |k| < π/2, decreases with increasing wavelength and time delay, see Fig. 2 (dashed red lines). Fourier modes with short wavelengths, |k| > π/2, display a cusp-like maximum, a behavior that was already observed in a different system as shown in Fig. 1. The delay τ corresponding to this maximum and the corresponding relaxation rate r 0 depend on the wavevector k.
In order to better understand these examples and to obtain basic insights in the behavior of the large class of systems described by Eq. (1), we perform a general study of the relaxation rate as a function of coupling delay and phase shift. We consider the linear dynamics near the synchronized state, θ i (t) = Ω 0 t + εξ i (t) with ε 1. For simplicity, we focus on coupling topologies, for which the normalized adjacency matrix b ij = a ij /ρ i is symmetric and b ii = 0 for all i. The two examples introduced above fall within this class of systems. Following [16], the inphase synchronized state of Eq. (1) is stable if and only if γ > 0, where γ ≡ K d dϑ Γ(ϑ)| ϑ=−ψ . We only consider these stable cases.
To first order in ε, the time evolution of the perturbation is given by We introduce the collective relaxation modes φ i (t) ≡ where d ij is defined by jk d −1 ij b jk d kl = u i δ il and u i are the N eigenvalues of the matrix b ij . The eigenvalues u i are real and obey |u i | ≤ 1 [16]. To compute the relaxation of collective modes, we take the time derivative of φ i (t) and use Eq. (5) to replace ξ j (t). Inserting the identity δ ij = k d ik d −1 kj enables to express the result in terms of the collective modes φ i and the eigenvalues u i . The collective modes relax independently according to The ansatz φ(t) = e −λt yields the characteristic equation for the relaxation rates λ of the collective mode φ [17,18], where we have dropped the index i for notational simplicity. Solutions to Eq. (7) can be expressed in terms of the Lambert W function [19], defined by the relation W (z)e W (z) = z for z ∈ C. This function has discrete branches W n (z) separated by branch cuts, where n is the branch index [20]. Each branch n of W corresponds to one relaxation rate r n = Re λ n . Note that our sign convention for λ implies that for stable states, all r n are positive. Here, we focus on the slowest relaxation rate r 0 for a given collective mode φ, which corresponds to the long time behavior of φ. Solving Eq. (7) for λ, the solution λ 0 is where z τ ≡ uγτ e γτ , since the principal branch W 0 has the property Re W 0 ≥ Re W n [21]. The dependence of λ 0 on the coupling delay τ thus depends on the properties of the principal branch W 0 of the Lambert function.
To discuss the properties of the slowest relaxation rates r 0 , we consider separately collective modes with u > 0 and u < 0. In nearest-neighbor coupled systems, collective modes with u > 0 are the Fourier modes with long wavelengths, as shown below. For these modes, r 0 decreases monotonically and converges to zero for τ → ∞ (dashed lines in Fig. 2). This can be shown using Eq. (7) and writing γ − r = uγe τ r cos(τ ν) , −ν = uγe τ r sin(τ ν) , where ν = Im λ. The smallest value of r = r 0 corresponds to cos(τ ν 0 ) = 1. From Eq. (10), it then follows that ν 0 = 0. Using Eq. (7), we find Furthermore, Eq. (7) implies that τ = ln(u[1 − r 0 /γ])/r 0 , which reveals that r 0 → 0 corresponds to τ → ∞. Therefore, r 0 vanishes for large τ . Hence, the collective modes corresponding to u > 0 become stationary for large time delay. Eq. (6) furthermore implies that given two eigenvalues u 1 ≥ u 2 , the respective exponents satisfy for all τ by an argument similar to the one leading to Eq. (11). In Fig. 2, this is illustrated by the fact that the dashed lines never cross.
In the case of collective modes with u < 0 in Eq. (7), it can be shown that r 0 displays a cusp at τ = τ * , where At τ = τ * , r 0 is not analytic and its first derivative has a jump which stems from the definition of the principal branch W 0 of the Lambert function. We now show that dr0 dτ has opposite sign in the two regions τ ≤ τ * and τ > τ * . Therefore, the cusp is located at the maximum of r 0 , as suggested by Figs. 1 and 2. Differentiating Eq. (8) with respect to τ and using the defining relation of the Lambert W function to compute its derivative, we obtain, with W 0 (z τ ) ≡ U + iV , For τ ≤ τ * , we find V = 0. This follows from the properties of the principal branch W 0 , in particular Im For τ > τ * , we show that dr0 dτ is negative. In Eq. (13), the factor U − γτ is negative. This can be seen by taking the real part of Eq. (8) and using the fact that r 0 > 0. Furthermore, V = 0 in this region. The factor U + U 2 + V 2 is positive: For τ > τ * , we have U = −V cot V and the numerator of the second factor in (13) can be rewritten as Since V ∈ [0, π], the above expression is positive. Altogether, we conclude that for u < 0, d dτ r 0 ≤ 0 for τ > τ * . The corresponding collective modes therefore resynchronize slower as the time delay increases. The maximal resynchronization rate r * 0 at τ = τ * is given by The behavior of r 0 in the limit of large τ can be obtained from an expansion of r 0 in powers of τ −1 , r 0 = − ln |u|/τ + O(τ −2 ), which reveals that for collective modes with eigenvalues u and −u, the synchronization rate r 0 approaches the same asymptotic behavior for large τ . The inset of Fig. 2 reflects this property for the case of nearest-neighbor coupling.
The mode structure shown in Fig. 2 can be understood as follows.
We illustrate this result for the case of a twodimensional system, see Fig. 3. The synchronization rate r 0 is displayed as a function of the wavevector (k x , k y ) in Fig. 3A and 3D for two systems with no delay and finite delay τ , respectively, and α chosen according to Eq. (3) to impose the same collective frequency. For no delay, Eq. (8) leads to the classical scenario where short-wavelength collective modes decay quickly while long-wavelength modes decay slowly (dark corners in Fig. 3A). Interestingly, for long delays collective modes decay fastest at intermediate wavelengths ( Fig. 2 and dark diamond in Fig. 3D). This remarkable behavior is confirmed by full simulations of Eq. (1), Figs. 3B,C,E and F. The inset of Fig. 3E shows partially synchronized clusters on intermediate length scales with persisting phase differences on nearest-neighbor scale. This behavior reflects the fact that the curves for short-wavelength collective modes in Fig. 2 reverse their ordering as τ increases. A similar mode reversal has been observed in small systems of chaotic oscillators as a function of the coupling strength [22].
The behavior of a globally coupled system, Fig. 1, can be understood as follows. The normalized adjacency matrix is given by b ij = (1 − δ ij )(N − 1) −1 . The largest eigenvalue, u = 1, corresponds to the neutrally stable global phase shift. All other collective modes have eigenvalue u = (1 − N ) −1 . These modes therefore exhibit the same synchronization rate whose τ -dependence is nonmonotonic. According to Eq. (12), the maximal synchronization rate of a system with global coupling is located at τ * = W 0 (e −1 [N − 1])/γ and depends on the system size and properties of the coupling.
In this work we have shown how coupling delays and phase shifts play a different role in regulating synchronization in systems of coupled phase oscillators. Our results show that synchronization rates can exhibit maxima as a function of time delay when the collective frequency is kept constant by adjusting phase shifts. Interestingly, in spatially extended systems with time delays, the relaxation rate does not always decrease with increasing wavelength but intermediate wavelengths may relax faster than short ones, giving rise to novel relaxation scenarios. Phase shifts alone cannot give rise to this behavior.
Fast synchronization improves the resilience of the synchronized state in the presence of fluctuations or diversity [23][24][25]. Here we have considered identical oscillators, but in natural systems diversity of oscillators can introduce a distribution of frequencies. For a narrow frequency distribution, we have confirmed numerically (data not shown) that a maximum of the synchronization rate still occurs for a non zero coupling delay. If the function of a system demands the collective frequency to be in a specific small range, the possibility to regulate synchronization rates using phase shifts and time delays at constant frequency might be important. Examples for such systems are the core pacemaker of the circadian clock, regulating metabolism in higher organisms with a period of about 24 hours [26][27][28], the segmentation clock of vertebrates [29], whose collective frequency determines the length of body segments [8,9,30,31], and engineered systems of coupled lasers or electronic oscillators [32,33]. Our work shows that together with phase shifts, coupling delays can play an important role for the regulation of dynamic behaviors and the resilience of synchronized oscillator networks.
We thank Lucas Wetzel for many fruitful discussions and Douglas B. Staple for critical comments on the manuscript. We thank Andy Oates and members of