**3. Reduction of the Collective Dynamics**

Following the theory of Ott and Antonsen [48,49] we will look for the long-term dynamics of the continuity Equation (5) in the form

$$f(\theta,\omega,t) = \frac{g(\omega)}{2\pi} \left( 1 + \sum\_{n=1}^{\infty} \left[ (z^\*(\omega,t))^n e^{in\theta} + (z(\omega,t))^n e^{-in\theta} \right] \right), \tag{7}$$

where

$$\lg(\omega) = \int\_0^{2\pi} f(\theta, \omega, t) d\theta \tag{8}$$

is the probability density function of the natural frequencies and

$$z(\omega, t) = \int\_0^{2\pi} f(\theta, \omega, t) e^{i\theta} d\theta,\tag{9}$$

is the local complex order parameter of the subpopulation with the natural frequency *ω*. Obviously, the global and the local order parameter are connected by the selfconsistency condition

$$R = \int\_{-\infty}^{\infty} g(\omega) z(\omega) d\omega. \tag{10}$$

Substituting (7) into (5), one obtains the following equations for *z*(*ω*, *t*):

$$z(\omega, t) = i\omega z + \frac{a}{2}(1 - z^2) + \frac{K}{2}(R\_\tau - R\_\tau^\* z^2),\tag{11}$$

which together with (10) defines a delay integro-differential equation describing the collective dynamics of the population in the thermodynamic limit.

As the next step we consider a family of rational distributions *g*(*ω*), namely

$$g\_n(\omega) = \frac{c\_n}{(\omega - \Omega)^{2n} + \Delta^{2n}} \tag{12}$$

where *n* is natural, Ω is the mean frequency, Δ is the distribution half-width, and

$$\mathcal{L}\_n = \frac{1}{\pi} n \sin \frac{\pi}{2n} \Delta^{2n-1} \tag{13}$$

is the normalization constant. For *n* = 1 this distribution turns into a classical Cauchy distribution, and for *n* → ∞ it converges to a uniform distribution on the interval *ω* ∈ [Ω − Δ; Ω + Δ]. Assuming the rational function *g*(*ω*) allows us to evaluate the integral (10) using the residue theorem. A similar approach was recently used for populations of quadratic integrate-and-fire neurons [50,51]. Consider the analytic extension of function *z*(*ω*, *t*) to complex *ω*, then the integration contour can be closed by an infinitely large arc in the upper complex half-plane. In this half-plane the function (12) has *n* simple poles

$$
\omega\_k = \Omega + \Delta e^{i\alpha\_k},\tag{14}
$$

where *k* = 1, *n* and *α<sup>k</sup>* = *π*(*k* − 0.5)/*n*. Thus, the integral (10) can be evaluated as

$$R(t) = -\frac{i}{\Delta} \sin \frac{\pi}{2n} \sum\_{k=1}^{n} (\omega\_k - \Omega) z(\omega\_{k'} t). \tag{15}$$

Writing Equation (11) for *ω*1, *ω*2,. . . , *ω<sup>n</sup>* allows us to obtain a closed set of *n* delay differential equations for complex variables

$$z\_k = i(\Omega + \Delta e^{ia\_k})z\_k + \frac{a}{2}(1 - z\_k^2) + \frac{K}{2}(R\_\pi - R\_\pi^\* z\_k^2),\tag{16a}$$

$$R = -i\sin\frac{\pi}{2n}\sum\_{k=1}^{n}e^{i\mathbf{a}\_k}z\_{k\prime} \tag{16b}$$

where *k* = 1, *n* and *zk*(*t*) ≡ *z*(*ωk*, *t*). Introducing *zk* = *xk* + *iyk* and *R* = *X* + *iY* allows us to rewrite these equations in the real form :

$$\begin{aligned} \mathbf{x}\_k^\prime &= -\Omega y\_k - \Delta (y\_k \cos \alpha\_k + \mathbf{x}\_k \sin \alpha\_k) + \dots \\ \frac{a}{2} (\mathbf{1} + y\_k^2 - \mathbf{x}\_k^2) &+ \frac{K}{2} (X\_\tau (\mathbf{1} + y\_k^2 - \mathbf{x}\_k^2) - 2Y\_\tau \mathbf{x}\_k y\_k)\_\prime \end{aligned} \tag{17a}$$

$$y\_k = \Omega x\_k + \Delta(x\_k \cos \alpha\_k - y\_k \sin \alpha\_k) - a x\_k y\_k + \dots$$

$$\frac{K}{2}(Y\_{\pi}(1+x\_k^2-y\_k^2)-2X\_{\pi}x\_ky\_k),\tag{17b}$$

$$X = \sin\frac{\pi}{2n} \sum\_{k=1}^{n} (x\_k \sin a\_k + y\_k \cos a\_k)\_{\prime} \tag{17c}$$

$$Y = \sin\frac{\pi}{2n} \sum\_{k=1}^{n} (y\_k \sin u\_k - x\_k \cos u\_k). \tag{17d}$$

This set of DDEs governs the collective dynamics of the population in the thermodynamic limit *N* → ∞. The following analysis is based on system (17).

## **4. Studying the Role of the Coupling Delay**

For the case of isochronous oscillations with *a* = 0, the system always has a trivial steady state *zk* = *R* = 0 corresponding to asynchronous dynamics of the oscillators. This state is stable for weak coupling *K* and destabilizes via an Andronov-Hopf bifurcation when the coupling becomes sufficiently strong, which constitutes a classical Kuramoto scenario [5,10]. The stable limit cycle born in this bifurcation corresponds to a partial synchronization of the oscillators.

The dynamics of the system is illustrated in Figure 1 where the dynamics of the individual phases and the Kuramoto order parameter are shown for the synchronous and asynchronous regimes. In both cases, the system starts from random initial conditions. In the case of asynchronous dynamics, all the phases rotate incoherently, and the order parameter remains close to zero. When the synchronization is achieved, a bunch of oscillators quickly emerge whose phases rotate with the same frequency, and the order parameter reaches a sufficiently non-zero value.

In order to illustrate the transition from the asynchronous to the synchronous state we plot the dependence of the Kuramoto order parameter on the coupling strength in Figure 2. The order parameter is small for weak coupling and rapidly grows as the coupling strength exceeds the critical value. Note that the results obtained by the simulation of the reduced model (17) coincide with those obtained for the microscopic system up to the fluctuations induced by the finite size effects. Note also that adding of the coupling delay might sufficiently influence the system dynamics and shift the critical value of the coupling strength. Further, we will analyze the role of the delay in detail with the help of the reduced system.

**Figure 1.** The dynamics of the system in the asynchronous (**a**,**b**) and synchronous (**c**,**d**) regimes. The top panels show the time traces of 10 randomly chosen phases *θj*, while the bottom panels show the time trace of the Kuramoto order parameter |*R*|. The coupling strength *K* = 1 for (**a**,**b**) and *K* = 4 for (**c**,**d**) The other parameters are *N* = 1000, *n* = 1, Ω = 1, Δ = 1, *τ* = 0.

**Figure 2.** The dependence of the Kuramoto order parameter on the coupling strength. Black circles: *τ* = 0, gray squares: *τ* = 1.5. The other parameters: *N* = 1000, *n* = 1, Ω = 3, Δ = 1. Solid lines indicate the results obtained by the simulation of the reduced system (17).

In the case of zero delay *τ* = 0 the critical coupling *Kc* depends on the distribution width Δ in a linear way. Indeed, according to the results of Kuramoto [52],

$$K\_{\mathbb{C}} = \frac{2}{\pi \mathfrak{g}(\Omega)} = \frac{2\Delta}{n \sin \frac{\mathfrak{T}}{2n}}.\tag{18}$$

For nonzero coupling delays the critical coupling becomes delay-dependent. In order to determine the bifurcation point it is necessary to write the characteristic equation for the trivial steady state which has the form |*D*(*λ*)| = 0, where *D* is (2*n*) × (2*n*) matrix

$$D = \left( \begin{array}{cc} D^{xx} & D^{xy} \\ D^{yx} & D^{yy} \end{array} \right) , \tag{19}$$

and *<sup>D</sup>xx*, *<sup>D</sup>xy*, *<sup>D</sup>yx*, *<sup>D</sup>yy* represent *<sup>n</sup>* × *<sup>n</sup>* matrices with the following elements:

$$D\_{km}^{xx} = \left(-\lambda - \Lambda \sin a\_k\right)\delta\_{km} + \frac{K}{2} \sin \frac{\pi}{2n} e^{-\lambda \tau} \sin a\_{k\tau} \tag{20}$$

$$D\_{\rm km}^{\rm xy} = -(-\Omega - \Delta \cos a\_k)\delta\_{\rm km} + \frac{K}{2}\sin \frac{\pi}{2n}e^{-\lambda \tau}\cos a\_{k\prime} \tag{21}$$

$$D\_{km}^{yx} = -\left(\Omega + \Lambda \cos a\_k\right)\delta\_{km} - \frac{K}{2}\sin\frac{\pi}{2\pi}e^{-\lambda\pi}\cos a\_k. \tag{22}$$

$$D\_{\rm km}^{yy} = -(-\lambda - \Delta \sin a\_k)\delta\_{\rm km} + \frac{K}{2}\sin \frac{\pi}{2n} e^{-\lambda \tau} \sin a\_{k\prime} \tag{23}$$

where *δkm* equals one for *k* = *m* and zero in other case.

At the Andronov-Hopf point a pair of roots *λ* = ±*iω* emerge, which allows us to determine the value of the coupling strength *Kb* at the bifurcation point by solving |*D*(*iω*)| = 0. For small delays, this equation can be solved numerically by taking (18) as the initial point, then the solution can be traced along the delay value as a parameter. The obtained dependence *Kb*(*τ*) is plotted in Figure 3a for the Lorentzian distribution of the oscillator frequencies (*n* = 1). The bifurcation coupling shows a minimum at *τ* = 0 and grows rapidly and monotonically for non-zero delays. Note that we have calculated the bifurcation curve for both positive and negative delays. Although negative time delays are not physical, we use them in the bifurcation analysis for a reason that will become clear later. Namely, they will help us to find other bifurcation curves existing for positive delays.

**Figure 3.** (**a**) The Andronov-Hopf bifurcation curve for system (17) with *n* = 1, *a* = 0, Ω = 3 and Δ = 0.1. (**b**) The same curve (black solid line) and its reappearing instances (black dashed lines). The gray thick line shows the synchronization border.

An important point is that the existence of one Andronov-Hopf bifurcation curve in a system with time delay implies the existence of other bifurcation curves at other delay values due to the so-called reappearance of periodic solutions [53]. Indeed, if the bifurcation takes place at the coupling strength *Kb* for the delay *τ*<sup>0</sup> this means the existence of a limit cycle with the period *T* = 2*π*/*ω* and vanishing amplitude. This implies the existence of the same limit cycle at the same coupling strength for the delays *τ<sup>k</sup>* = *τ*<sup>0</sup> + *kT*, where *k* is an integer, which means that each of these points are also Andronov-Hopf points. Note, however, that the stability of the emergent limit cycle can change, which means that the bifurcation can be either supercritical or subcritical.

The bifurcation curve found by starting from the delay-less case reappears in the manner described above and leads to the emergence of other bifurcation curves as shown in Figure 3b. These curves demonstrate minimums at *τ* = *kT*0, where *T*<sup>0</sup> = 2*π*/*ω*<sup>0</sup> and *ω*<sup>0</sup> is the frequency at which the collective oscillations emerge for *τ* = 0. Obviously, it is the frequency of the distribution peak, i.e., *ω*<sup>0</sup> = Ω. For non-zero delays the frequency becomes different, therefore the different bifurcation curves do not exactly match each other in shape. However, they all still a single minimum at the multiples of *T*0. The trivial steady state, i.e., the asynchronous regime is stable when the coupling strength is below all the bifurcation curves. Thus, the synchronization border is defined by the lowest point of the curve and has a saw-like shape. This border coincides completely with that obtained in Ref. [17] for the same setting which corroborates the validity of our analysis.

The obtained results suggest that introduction of the coupling delay prevents the system synchronization: For non-zero delays, the critical coupling at which the oscillators start to synchronize increases with respect to the delay-free case (at best, the critical coupling does not change if the delay is a multiple of *T*0). This result seems to be obvious from the physical viewpoint: It is harder to adjust if one receives outdated information. However, it turns out that the coupling delay can in some cases promote synchronization. This surprising effect is observed when the distribution *g*(*η*) is different from Lorentzian, i.e., *n* > 1.

In order to illustrate this effect we calculated the synchronization borders on the *τ* − *K* plane for different values of *n*. We adjust the distribution half-width as

$$
\Delta = \frac{n}{2} \sin \frac{\pi}{2n'} \tag{24}
$$

so that the critical coupling for the delay-free case equals unity for all *n*. The results are plotted in Figure 4 and show a significant difference between the Lorentzian distribution and distributions with *n* > 1. For the Lorentzian distribution the delay always prevents synchronization and the critical coupling is always not less than unity. For the distributions of higher order *n* > 1, some delay values can promote synchronization so that the critical coupling becomes less than unity. Another feature of high-order distributions is the complex form of the synchronization border with many peaks and valleys of different shapes.

**Figure 4.** (**a**) The synchronization borders of system (17) for *n* = 1 (dotted line), *n* = 2 (dash-dotted line) and *n* = 5 (solid line). The mean frequency Ω = 3, the half-width Δ = *<sup>n</sup>* <sup>2</sup> sin *<sup>π</sup>* <sup>2</sup>*<sup>n</sup>* . (**b**) Enlarged part of the panel (**a**).
