**1. Introduction**

In many complex systems of various nature a similar pattern of collective behaviour can be observed: The adjustment of rhythms of oscillating systems due to their interaction. This phenomenon is called synchronization and is ubiquitous in the real world with examples ranging from the synchronous firing of pacemaker cells in the heart and neurons in the brain to the synchronous rotation of electric generators in power grids [1–3]. For a long time it was a mystery how the synchronization can emerge despite the inevitable diversity in the natural frequencies of different units. It was Kuramoto who introduced a mathematical model of coupled oscillators that allowed this problem to be solved [4–6]. Motivated by the behavior of chemical and biological oscillators, this model later turned out to be quite general and applicable to such systems as coupled arrays of Josephson junctions [7] or populations of of biological neurons [8,9].

What is now known as the "Kuramoto model" consists of *N* phase oscillators with the harmonic interaction function:

$$\theta\_{\dot{\jmath}}(t) = \omega\_{\dot{\jmath}} - \frac{K}{N} \sum\_{k} \sin(\theta\_{\dot{\jmath}} - \theta\_k) \,\prime \tag{1}$$

where *<sup>j</sup>* = 1, ... , *<sup>N</sup>* is the unit number, *<sup>θ</sup><sup>j</sup>* ∈ *<sup>S</sup>*<sup>1</sup> are the phase variables, *<sup>ω</sup><sup>i</sup>* are the natural frequencies, and *K* is the global coupling strength. In his pioneering work Kuramoto showed this model to be mathematically tractable and to demonstrate the phase transition from an asynchronous state to synchronization at the critical coupling strength.

**Citation:** Klinshov, V.V.; Zlobin, A.A. . *Mathematics* **2023**, *11*, 2325. https:// doi.org/10.3390/math11102325

Academic Editor: Alexandra Kashchenko

Received: 18 April 2023 Revised: 7 May 2023 Accepted: 10 May 2023 Published: 16 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In subsequent decades the Kuramoto model became a classical paradigm for studying synchronization phenomena and it found numerous applications (see reviews [10,11] for the examples). From the other side, the simplicity and generality of the Kuramoto model makes it perfect for various kinds of generalization and modification. One of the most natural ways to make the model more realistic is to include coupling delays which are inevitably in real world due to the finite speed of signal propagation [12–14]. To take into account coupling delays the phases *θk*(*t*) in the sum from the r.h.s. of (1) are replaced by their delayed versions *θk*(*t* − *τ*). In the case of two coupled oscillators the delay causes multistabilty of synchronous solutions with different frequencies [15] (interestingly, multistability also emerges if the delay is introduced not into the coupling, but into the oscillators themselves [16]). The similar effect of multistability was found for globally coupled oscillators with time delays, which results in the possibility of discontinuous transitions between different regimes in addition to classical smooth transitions [17,18]. In the case of local coupling, the delays were shown to induce complex spatio-temporal patterns and clusters [19]. In the case of distant-dependant delays similar patterns emerge and become propagating structures [20,21]. The case of heterogeneous delays was also considered and the emergence of clusters with various phase relations was demonstrated [22]. In rings of oscillators with local coupling, the delay was shown to act differently depending on the network symmetry: It induces multistability for unidirectional coupling and stabilizes the most symmetric solution for bidirectional coupling [23]. In rings with nonlocal coupling, the delays give birth to travelling waves [24]. For networks of oscillators with complex connectivity, such as scale-free networks, coupling delays may induce explosive synchronization [25].

Another natural generalization of the Kuramoto model is the addition of external forces to the oscillator dynamics which results in the emergence of a phase-dependent and/or time-dependent term in the r.h.s. of (1). This term is often taken in the form *a* sin *θj*, in which case the phase oscillators turn into the so-called "active rotators". The dynamics of active rotators is much richer than that of phase oscillators since the former can demonstrate either oscillatory or excitatory behaviour depending on the relation between the individual frequency *ω<sup>j</sup>* and the parameter *a* which is often called a "non-isochronicity parameter". Ensembles of globally coupled active rotators were first considered in Refs. [26,27], and two different collective regimes were found: Entrainment with the external force and mutual entrainment. The introduction of heterogeneous external forces leads to similar results [28]. The transition between locked and unlocked states can be similar to either a super-critical Andronov-Hopf bifurcation or a saddle-node bifurcation [29]. In the case of heterogeneous assemblies made up of excitable and oscillatory units, it has been demonstrated that the transition to synchrony may be classical or re-entrant, depending on the particular form of the frequency distribution [30]. In addition to the parameter regions with synchronous and asynchronous regimes, the regions of bistability between different regimes were found as well [31].

Other types of generalized Kuramoto models were studied as well, including timedependent [8,32,33], adaptive [34–38] and memristive [39] coupling, including noise [26,27,40], non-harmonic coupling functions [41,42] and pulse coupling [43,44], systems on smooth manifolds [45–47], etc. Typically, the analytic study of the Kuramto model and its generalizations is performed in the thermodynamic limit when the number of units tends to infinity. In this case the microscopic equations (1) are replaced by the continuity equation, and the individual frequencies *ω<sup>j</sup>* are replaced by the probability density *g*(*ω*). For the classical Kuramoto model the particular form of this distribution does not play any important role. First, the system is invariant under the change of variables *ω<sup>j</sup>* → *ω<sup>j</sup>* + Ω for arbitrary Ω, which means that the mean frequency *ω* can be always set to zero. Second, for unimodal symmetric distributions, the critical coupling at which the transition takes place equals *Kc* = 2/(*πg*(0)), i.e., it depends only on the height of the distribution peak but not on its particular shape. However, when coupling delays or non-isochronicity are included, the shape of the distribution comes into play because it becomes important how the typical frequencies *ω* relate to the delay *τ* and non-isochronicity parameter *a*.

In the present paper we study how the synchronization transition in the Kuramoto model with delay depends on the shape of the frequency distribution *g*(*ω*). For this sake we consider a series of rational distributions which are all unimodal and symmetric but are different in the flatness of the peak and decay rate of the tails. Using the Ott-Antonsen approach [48,49] we derive low-dimensional dynamical systems governing the collective dynamics of the population and perform its bifurcation analysis. We show that the shape of the frequency distribution can play a significant role for the system properties. For widely used Lorentzian distribution, the delay always prevents synchronization by raising the critical coupling strength. However, for distributions with lighter tails, the delay can also promote synchronization by lowering the critical coupling.

### **2. Model**

We consider a heterogeneous assembly of *N* oscillators with delayed coupling

$$\dot{\theta}\_{\dot{j}}(t) = \omega\_{\dot{j}} - a\sin\theta\_{\dot{j}}(t) - \frac{K}{N} \sum\_{k} \sin\left(\theta\_{\dot{j}}(t) - \theta\_{k}(t-\tau)\right),\tag{2}$$

where *<sup>j</sup>* = 1, ... , *<sup>N</sup>* is the unit number, *<sup>θ</sup><sup>j</sup>* ∈ *<sup>S</sup>*<sup>1</sup> are the phase variables, *<sup>ω</sup><sup>j</sup>* are the natural frequencies, *a* is the non-isochronicity parameter, *K* is the global coupling strength and *τ* is the coupling delay. Strictly speaking, *ω<sup>i</sup>* can be called "natural frequencies" only for *<sup>a</sup>* <sup>=</sup> 0. For non-zero *<sup>a</sup>*, the rotation becomes non-uniform with the frequency - *ω*2 *<sup>j</sup>* − *a*<sup>2</sup> for 0 < *a* < |*ωj*|. At |*ωj*| = *a* an isolated unit undergoes a SNIPER bifurcation toward the excitatory regime. This very dependence of the local oscillatory dynamics on the parameter *a* allows us to call it a "non-isochronicity parameter".

In order to characterize the degree of synchrony in the population we introduce the Kuramoto complex order parameter

$$R(t) = \frac{1}{N} \sum\_{j} e^{i\theta\_j(t)}.\tag{3}$$

The absolute value of this parameter serves the main indicator of the system synchronization. When it is close to zero, the phases of the oscillators are not correlated, and the system is asynchronous. When the order parameter is sufficiently different from zero it indicates the emergence of a bunch of oscillators rotating with close phases, which means synchronization. Using the Kuramoto order parameter allows us to rewrite (2) as

$$\dot{\theta}\_{\dot{\jmath}} = \omega\_{\dot{\jmath}} - \frac{a}{2i} (e^{i\theta\_{\dot{\jmath}}} - e^{-i\theta\_{\dot{\jmath}}}) + \frac{K}{2i} (R\_{\tau}e^{-i\theta\_{\dot{\jmath}}} - R\_{\tau}^{\*}e^{i\theta\_{\dot{\jmath}}}),\tag{4}$$

where *R<sup>τ</sup>* ≡ *R*(*t* − *τ*) and the asterisk denotes the complex conjugate. In the thermodynamic limit *N* → ∞, the macroscopic state of the system is described by the probability density function *f*(*θ*, *ω*, *t*), which obeys the continuity equation

$$\frac{\partial f}{\partial t} + \frac{\partial}{\partial \theta}(fv) = 0,\tag{5}$$

with the velocity *v* being the r.h.s. of Equation (4). The Kuramoto parameter in this case is evaluated as

$$R(t) = \int\_{-\infty}^{\infty} d\omega \int\_{0}^{2\pi} f(\theta, \omega, t) e^{i\theta} d\theta \,\tag{6}$$
