Spontaneous collective synchronization in the Kuramoto model with additional non-local interactions

TL;DR

This paper explores the complex phase diagram of the Kuramoto model with non-local interactions and thermal noise, revealing both equilibrium and non-equilibrium phase transitions through exact analytical methods.

Contribution

It provides the first exact analytical results for the Kuramoto model with non-local interactions in both noise-free and identical-frequency limits.

Findings

Non-local interactions can both promote and hinder synchronization.

Exact solutions are obtained using Ott-Antonsen ansatz and transfer matrix methods.

The phase diagram includes rich equilibrium and non-equilibrium transitions.

Abstract

In the context of the celebrated Kuramoto model of globally-coupled phase oscillators of distributed natural frequencies, which serves as a paradigm to investigate spontaneous collective synchronization in many-body interacting systems, we report on a very rich phase diagram in presence of thermal noise and an additional non-local interaction on a one-dimensional periodic lattice. Remarkably, the phase diagram involves both equilibrium and non-equilibrium phase transitions. In two contrasting limits of the dynamics, we obtain exact analytical results for the phase transitions. These two limits correspond to (i) the absence of thermal noise, when the dynamics reduces to that of a non-linear dynamical system, and (ii) the oscillators having the same natural frequency, when the dynamics becomes that of a statistical system in contact with a heat bath and relaxing to a statistical equilibrium state. In the former case, our exact analysis is based on the use of the so-called Ott-Antonsen ansatz to derive a reduced set of nonlinear partial differential equations for the macroscopic evolution of the system. Our results for the case of statistical equilibrium are on the other hand obtained by extending the well-known transfer matrix approach for nearest-neighbor Ising model to consider non-local interactions. The work offers a case study of exact analysis in many-body interacting systems. The results obtained underline the crucial role of additional non-local interactions in either destroying or enhancing the possibility of observing synchrony in mean-field systems exhibiting spontaneous synchronization.

Equations

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}={\omega }_j+\frac{J}{N}\left(\sum _{k=1}^{N}X({\theta }_k)\right)Z({\theta }_j),$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}={\omega }_j+\frac{J}{N}\sum _{k=1}^N\sin ({\theta }_k-{\theta }_j).$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=\Delta {\omega }_j+\frac{J}{N}\sum _{k=1}^N\sin ({\theta }_k-{\theta }_j).$$

$$r(t)\equiv \frac{1}{N}\sum _{k=1}^N{e}^{i{\theta }_k(t)}.$$

$$r^{\mathrm{st}}\equiv |r(t\to \infty )|$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=\Delta {\omega }_j+\frac{J}{N}\sum _{k=1}^N\sin ({\theta }_k-{\theta }_j)+\frac{K}{2M}\sum _{k=-M}^M\sin ({\theta }_{j+k}-{\theta }_j),$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=\omega +\frac{1}{2M}\sum _{k=-M}^M\sin ({\theta }_{j+k}-{\theta }_j-\alpha ),$$

$$g(\omega )=\frac{1}{\pi }\frac{1}{{\omega }^2+1}.$$

$$\overline{t}\equiv Jt,\overline{\Delta }\equiv \frac{\Delta }{J},\overline{K}\equiv \frac{K}{J};$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}\overline{t}}=\overline{\Delta }{\omega }_j+\frac{1}{N}\sum _{k=1}^N\sin ({\theta }_k-{\theta }_j)+\frac{\overline{K}}{2M}\sum _{k=-M}^M\sin ({\theta }_{j+k}-{\theta }_j).$$

$$H=\sum _{j=1}^N\frac{p_j^2}{2m}+\frac{\tilde{J}}{2N}\sum _{j,k=1}^N\left[1-\cos ({\theta }_j-{\theta }_k)\right]-\frac{\tilde{K}}{4M}\sum _{j=1}^N\sum _{k=-M}^M\cos ({\theta }_{j+k}-{\theta }_j),$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=v_j,$$

$$m\frac{\mathrm{d}{v}_j}{\mathrm{d}t}=\gamma \tilde{\Delta }{\omega }_j-\gamma v_j+\frac{\tilde{J}}{N}\sum _{k=1}^N\sin ({\theta }_k-{\theta }_j)+\frac{\tilde{K}}{2M}\sum _{k=-M}^M\sin ({\theta }_{j+k}-{\theta }_j)+\sqrt{\gamma }{\eta }_j(t),$$

$$⟨{\eta }_j(t)⟩=0,⟨{\eta }_j(t){\eta }_k(t^{\prime })⟩=2T{\delta }_{jk}\delta (t-t^{\prime }).$$

$$\overline{t}\equiv t\frac{\tilde{J}}{\gamma },\overline{\Delta }\equiv \gamma \frac{\tilde{\Delta }}{\tilde{J}},{\overline{\eta }}_i(\overline{t})\equiv {\eta }_i(t)\sqrt{\frac{\gamma }{\tilde{J}}},\overline{K}\equiv \frac{\tilde{K}}{\tilde{J}},\overline{T}\equiv \frac{T}{\tilde{J}},$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}\overline{t}}=\overline{\Delta }{\omega }_j+\frac{1}{N}\sum _{k=1}^N\sin ({\theta }_k-{\theta }_j)+\frac{\overline{K}}{2M}\sum _{k=-M}^M\sin ({\theta }_{j+k}-{\theta }_j)+{\overline{\eta }}_j(\overline{t}),$$

$$⟨{\overline{\eta }}_j(\overline{t})⟩=0,⟨{\overline{\eta }}_j(\overline{t}){\overline{\eta }}_k({\overline{t}}^{\prime })⟩=2\overline{T}{\delta }_{jk}\delta (\overline{t}-{\overline{t}}^{\prime }).$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=\Delta {\omega }_j+\frac{1}{N}\sum _{k=1}^N\sin ({\theta }_k-{\theta }_j)+\frac{K}{2M}\sum _{k=-M}^M\sin ({\theta }_{j+k}-{\theta }_j)+{\eta }_j(t);$$

$$⟨{\eta }_j(t)⟩=0,⟨{\eta }_j(t){\eta }_k(t^{\prime })⟩=2T{\delta }_{jk}\delta (t-t^{\prime }),$$

$$P_{\mathrm{eq}}(\{{\theta }_j\})\propto \exp [-\mathcal{V}(\{{\theta }_j\})/T],$$

$$\mathcal{V}(\{{\theta }_j\})\equiv \frac{1}{2N}\sum _{j,k=1}^N\left[1-\cos ({\theta }_j-{\theta }_k)\right]-\frac{K}{4M}\sum _{j=1}^N\sum _{k=-M}^M\cos ({\theta }_{j+k}-{\theta }_j);$$

$$2={\int }_{-\infty }^{\infty }\mathrm{d}\omega \frac{Tg(\omega )}{T^2+{\omega }^2{\Delta }_c^2(K=0,T)}.$$

$${\Delta }_c(K=0,T)=\frac{1}{2}-T.$$

$$Z_j(t)\equiv \frac{1}{2M}\sum _{k=-M}^M{e}^{i{\theta }_{j+k}(t)}.$$

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=\Delta {\omega }_j+\frac{1}{2i}[r(t)e^{-i{\theta }_j}-r^{\star }(t)e^{i{\theta }_j}]+\frac{K}{2i}[Z_j(t)e^{-i{\theta }_j}-Z_j^{\star }(t)e^{i{\theta }_j}],$$

$${\int }_0^{2\pi }\mathrm{d}\theta f(\theta ,\omega ,x,t)=g(\omega )\forall x,t.$$

$$Z(x,t)={\int }_0^{2\pi }\mathrm{d}yG(x-y){\int }_{-\infty }^{\infty }\mathrm{d}\omega {\int }_0^{2\pi }\mathrm{d}\theta e^{i\theta }f(\theta ,\omega ,y,t),$$

$$G(x)=\left\{\begin{array}[]{ll}\frac{1}{4\pi\sigma}~{}{\rm if}~{}|x|<2\pi\sigma,\\ 0~{}{\rm otherwise},\\ \end{array}\right.\\$$

$$r(t)=\frac{1}{2\pi }{\int }_0^{2\pi }\mathrm{d}x{\int }_{-\infty }^{\infty }\mathrm{d}\omega {\int }_0^{2\pi }\mathrm{d}\theta e^{i\theta }f(\theta ,\omega ,x,t).$$

$$\displaystyle\frac{\partial f}{\partial t}+\frac{\partial}{\partial\theta}\Big{[}\Big{(}\Delta~{}\omega+\frac{1}{2i}[r(t)e^{-i\theta}-r^{\star}(t)e^{i\theta}]+\frac{K}{2i}[Z(x,t)e^{-i\theta}-Z^{\star}(x,t)e^{i\theta}]\Big{)}f\Big{]}=0.$$

$$f(\theta,\omega,x,t)=\frac{g(\omega)}{2\pi}\left[1+\sum_{n=1}^{\infty}\Big{(}\widetilde{f}_{n}(\omega,x,t)e^{in\theta}+[\widetilde{f}_{n}(\omega,x,t)]^{\star}e^{-in\theta}\Big{)}\right],$$

Full text

• Invited contribution to the Journal of Physics A: Math. Theor. Special Issue “Emerging Talents” being published on the occasion of the 50th anniversary celebrations in 2017 of the Journal of Physics series

Keywords: Nonlinear dynamics, Statistical mechanics, Synchronization, Phase transitions

1. 1 Introduction

2. 2 Model

3. 2.1 Relation to long-range interacting systems

4. 2.2 Summary of known results and our queries

5. 3 Analysis for the $(\Delta,K)$-plane with $T=0$

6. 3.1 Uniformly incoherent state: Stability

7. 3.2 Synchronized twisted state

8. 3.2.1 Existence

9. 3.2.2 Stability

10. 3.2.3 Stability of the zero-twist state

11. 3.3 Phase transition in the Kuramoto order parameter

12. 4 Analysis for the $(K,T)$-plane with $\Delta=0$

13. 4.1 $K=0$

14. 4.2 $K\neq 0$

15. 4.2.1 An equivalent Ising problem

16. 5 Simulation results for a general point in the $(\Delta,K,T)$-space

17. 6 Conclusions and perspectives

18. 7 Acknowledgments

19. A Proof that the dynamics (17) does not satisfy detailed balance unless $\Delta=0$

20. B Numerical scheme to integrate the equation of motion (17)

1 Introduction

Dynamical systems comprising a large number of interacting constituents with a non-linear evolution in time generically exhibit a variety of rich emergent behaviors that go well beyond the behaviour of the constituent elements [1, 2]. Perhaps the most fascinating one is that of collective synchronization, in which a large population of oscillating units that have diverse frequencies and are interacting weakly with one another adjust their individual rhythms to spontaneously evolve to a state in which the units operate in unison [3, 4]. Spontaneous synchronization lies at the heart of many physical phenomena in nature. Even sustenance of life requires the heart to beat as a result of harmonized contractions of the cardiac cells in a synchronous wave. A highly non-linear cooperative effect, synchronization is observed in yeast cell suspensions [5], among flashing fireflies [6, 7], in crickets chirping in synchrony [8], in an audience clapping in unison for a laudable performance in a concert hall [9, 10], among pedestrians on footbridges [11], and in a variety of experiments involving electrochemical [12] and electronic [13] oscillators, metronomes [14], Josephson junctions [15], laser arrays [16], etc. The reader may refer to Ref. [17] for an engrossing exposition of synchronization in biological systems.

An early theoretical approach to synchrony is due to Winfree, who couched the problem in the framework of $N\gg 1$ oscillators of nearly-identical frequencies that are weakly coupled to one another. The weak coupling leads to a fast relaxation of the oscillators to their limit cycles [1, 2], so that for subsequent times, they may be characterized solely by their phases. In the following, the word “phase” would be used to also refer to a thermodynamic phase of a macroscopic system. In order to avoid any possible confusion between the two different usages of the word “phase”, we will from now on use the term “angle” to mean oscillator phase, and the term “phase” to exclusively mean a thermodynamic phase. On timescales longer than the one over which the oscillators are characterized by their angles alone, the latter would evolve in time due to the coupling and the frequency differences between the oscillators. To avoid complications arising from a spatial distribution of the oscillators, Winfree endowed the system with a mean-field geometry: every oscillator responds to the collective effect of the whole population. The angles $\theta_{j}\in[0,2\pi);~{}j=1,2,\ldots,N$ evolve in time as [18, 19]

[TABLE]

where $\omega_{j}$ is the natural oscillation frequency of the $j$-th oscillator. The natural frequencies of all the oscillators form a set of quenched disordered random variables with a common probability distribution $G(\omega)$. In Eq. (1), $X(\theta_{k})$ denotes the influence of the $k$-th oscillator on the $j$-th one, which responds to the collective influence $\left(\sum_{k=1}^{N}X(\theta_{k})\right)$ through the sensitivity function or the so-called phase response curve $Z(\theta_{j})$. The model (1) assumes that the functions $Z$ and $X$ are identical for all the oscillators, and that the coupling $J>0$ is the same for every pair of oscillators. The scaling of $J$ by $N$ in Eq. (1) ensures that the model is well-behaved in the thermodynamic limit $N\to\infty$. Winfree showed in simulations that for a given $J$ and for sufficiently small diversity in the natural frequencies, the system (1) exhibits a synchronized state at long times. However, analytical results had to wait for Kuramoto, who came up with a simplification of the Winfree model that made it amenable to an exact treatment in the thermodynamic limit $N\to\infty$ (see however Ref. [20] and references therein). He assumed the weak-coupling condition $|\omega_{j}|\gg J~{}\forall~{}j$, so that Eq. (1) for every oscillator may be averaged over its oscillation frequency, thus obtaining for $Z(\theta_{j})X(\theta_{k})$ a function solely of the angle difference given by $d(\theta_{j}-\theta_{k})$ (see, e.g., Ref. [21] for details). He made a simple choice, $d(\theta)=\sin\theta$, thus obtaining the dynamics [22, 23, 24, 25, 26, 27, 28]

[TABLE]

We show later that Eq. (2) corresponds to the driven overdamped dynamics of globally-coupled $XY$ spins. However, we emphasize that there is no a priori justification to treat interacting limit-cycle oscillators as $XY$ spins. We note in passing that just as the Kuramoto model is obtained as the weak-coupling limit of the Winfree model, one may also obtain (see, e.g., Refs. [29, 30]) a Winfree-type ensemble of oscillators by considering a suitable weak-coupling limit of the dynamics of the so-called pulse-coupled leaky integrate-and-fire neuron system [31] that is extensively employed in the field of computational neuroscience to study neuronal dynamics.

Let us consider for $G(\omega)$ a unimodal distribution, i.e., one which is symmetric about the mean $\Omega_{0}$ and decreases monotonically to zero with increasing $|\omega-\Omega_{0}|$. The effect of $\Omega_{0}$ can be gotten rid of from Eq. (2) by viewing the dynamics in a frame rotating uniformly with frequency $\Omega_{0}$ with respect to an inertial frame; this tantamounts to implementing the Galilean shift $\theta_{j}\to\theta_{j}+\Omega_{0}t~{}\forall~{}j$ that leaves Eq. (2) invariant. Denoting the half-width-at-half-maximum (HWHM) of $G(\omega)$ by $\Delta>0$ 111For a Gaussian distribution with standard deviation equal to $\sigma_{\rm Gaussian}$, the HWHM is given by $\Delta=\sigma_{\rm Gaussian}\sqrt{2\ln 2}$., we may put in evidence the dependence of the dynamics (2) on $\Delta$ by replacing the term $\omega_{j}$ by $\Delta~{}\omega_{j}$, and concomitantly, consider from now on the $\omega_{j}$’s as dimensionless random numbers with a common distribution $g(\omega)$ that has zero mean and unit width and the normalization $\int_{-\infty}^{\infty}{\rm d}\omega~{}g(\omega)=1$. We thus obtain the dynamics

[TABLE]

From Eq. (3), we note that for a given $g(\omega)$, the frequency term alone induces independent oscillations of every oscillator at its own natural frequency, a tendency that is opposed by the global coupling that favors equal angles for all the oscillators, thereby promoting global synchrony. It is convenient to visualize the angles as points moving on a unit circle under the dynamics (3). A synchronized or a clustered state then corresponds to a macroscopic cluster of these points that is immobile in time 222This is because the dynamics (3) refers to a frame that is rotating uniformly with angular frequency $\Omega_{0}$ with respect to an inertial frame. When viewed in the latter, however, the cluster moves around the circle at angular frequency $\Omega_{0}$., while an unsynchronized state has points randomly distributed over the circle. For low $J$ and initial $\theta_{j}$’s that are all equal, the points on the circle while starting bunched together spread out on a timescale $\sim 1/\Delta$ due to the diversity in the $\omega_{j}$’s. By contrast, for sufficiently high $J$ and an initial state with small bunching, the interaction term in (3) grows in time by pulling in more and more oscillators towards the bunch, thus inducing a relaxation to a synchronized state. For a given initial condition (or an ensemble of initial conditions), whether synchrony is sustained at long times and the amount of it is determined by an interplay of the interaction with the diversity in the natural frequencies of the oscillators.

To characterize quantitatively the amount of synchrony in the system, Kuramoto introduced the (complex) synchronization order parameter $r(t)$ defined as [23]

[TABLE]

Correspondingly, one has a vector in the complex-$r$ plane with $x$ and $y$ components $(r_{x},r_{y})\equiv\left(\frac{1}{N}\sum_{k=1}^{N}\cos\theta_{k},\frac{1}{N}\sum_{k=1}^{N}\sin\theta_{k}\right)$; the length $|r|\equiv\sqrt{r_{x}^{2}+r_{y}^{2}}$ measures the amount of synchrony, while $\tan^{-1}(r_{y}/r_{x})$ gives the average angle. When the oscillators are unsynchronized so that over a stretch of time or in an ensemble of configurations at a given time, one has with equal probabilities $e^{i\theta}$ equal to any complex number with modulus unity, $|r|$ averages to zero. On the other hand, $|r|$ has a non-zero average when a finite fraction of oscillators have angle differences that are constant in time. We will denote by

[TABLE]

the stationary value of the synchronization order parameter, which may be obtained by averaging $r_{x}^{2}$ and $r_{y}^{2}$ over the stationary ensemble of configurations. Based on the discussions above, we expect $r^{\rm st}$ to exhibit qualitatively different behaviors as $\Delta$ is tuned from low (thus favoring $r^{\rm st}\neq 0$) to high (favoring $r^{\rm st}=0$) values at a fixed $J$. Indeed, it has been rigorously established that under such a tuning of $\Delta$, the system (3) in the thermodynamic limit undergoes in the stationary state a continuous phase transition, from a low-$\Delta$ synchronized phase ($r^{\rm st}\neq 0$) to a high-$\Delta$ incoherent phase ($r^{\rm st}=0$), at the critical threshold $\Delta_{c}=\pi Jg(0)/2$ [23, 32, 26]. In the thermodynamic limit, the system is well characterized by the probability density function $f(\theta,\omega,t)$, defined such that $f(\theta,\omega,t){\rm d}\theta{\rm d}\omega$ gives out of all oscillators with frequency in $[\omega,\omega+{\rm d}\omega]$ the fraction at time $t$ that have their angle in $[\theta,\theta+{\rm d}\theta]$. While the incoherent state with $f(\theta,\omega,t)=1/(2\pi)$ is linearly neutrally stable at all $\Delta$’s, a stable branch corresponding to a synchronized state bifurcates continuously for $\Delta\leq\Delta_{c}$ [24].

Over the years, the Kuramoto model has served as a paradigm to study spontaneous collective synchronization in many-body interacting systems, and has moreover initiated a wide variety of studies criss-crossing several disciplines and involving physicists, mathematicians, and applied scientists. For an overview of recent progresses and perspectives on the model and its many variants, see Ref. [33]. Results from extensive studies of the model have found numerous applications in areas ranging from bridge engineering and social sciences to neuroscience, and have even led to the introduction of novel theoretical concepts in nonlinear science such as the chimera states [34, 35], see Ref. [36] for a recent review. Chimeras are broken-symmetry states occurring in identical, symmetrically-coupled oscillator ensembles in which synchronized and desynchronized sub-populations coexist. These states have been observed in a variety of experimental situations involving, e.g., chemical and mechanical oscillators and photoelectrochemical devices [36], and also in many theoretical frameworks besides the Kuramoto setting, e.g., in a system of globally-coupled complex Ginzburg-Landau oscillators [37] and in a network of coupled-map lattices [38].

In this work, we investigate as to how the stationary behavior of the mean-field Kuramoto model (3), summarized above, gets modified by the inclusion of competing interactions that are non-local in space. Specifically, in the setting of a one-dimensional periodic lattice with sites occupied by limit-cycle oscillators, we consider in addition to a global coupling of the form in Eq. (3) a non-local coupling of strength $K$ between the angles of oscillators on one site with those of $M$ nearest-neighbor oscillators to the left and to the right. Moreover, we consider the dynamical evolution to take place in presence of a stochastic noise, modelled as a Gaussian, white noise with strength characterized by an effective temperature $T$. We take the coupling $K$ to be either positive or negative. In the former case, the non-local interaction acts in conjunction with the one due to the global coupling in inducing synchrony in the system, thereby leading to $r^{\rm st}\neq 0$. For $K<0$, on the other hand, the non-local interaction competes with the global coupling and may thus destroy the possibility of observing synchrony in the system at long times. The dynamics of our model is characterized by three parameters, namely, the HWHM $\Delta$ of the frequency distribution, the non-local coupling $K$, and the temperature $T$. Interestingly, for $\Delta=0$, when all the oscillators have the same natural frequency of oscillation, the dynamics may be reduced to that of a Hamiltonian system in contact with a heat bath for which the stationary state has the usual Gibbs-Boltzmann form [39] of phase-space distribution $\sim\exp(-H/T)$, with $H$ being the underlying Hamiltonian. For $\Delta\neq 0$, however, the dynamics relaxes at long times to a nonequilibrium stationary state [40] (in technical terms, unlike its equilibrium counterpart, the corresponding phase-space distribution does not satisfy detailed balance). Thus, for general non-zero values of $\Delta,K,T$, the dynamics of our model is dissipative, noisy, and is moreover out of equilibrium. A combination of all of these factors, together with the non-linear nature of the dynamical equations, offers a rather rich playground to observe interesting collective effects, while rendering at the same time the task of pursuing an exact analytical treatment of the system one of great difficulty.

A powerful exact method that has been recently developed to study non-noisy dynamics of coupled oscillator ensembles is the proposition and the implementation of the so-called Ott-Antonsen (OA) Ansatz [41, 42], which allows to rewrite in the thermodynamic limit the dynamics of coupled networks of phase oscillators in terms of a few collective variables. Specifically, in the context of the Kuramoto model (3) with a Lorentzian distribution of the oscillator frequencies, the ansatz studies the evolution in phase space by considering in the space ${\cal D}$ of all possible phase-space distributions $f(\theta,\omega,t)$ a particular class defined on and remaining confined to a manifold ${\cal M}$ in ${\cal D}$ under the time evolution of the angles. As a result of the choice of the particular class of $f(\theta,\omega,t)$, one obtains a single first-order ordinary differential equation for the evolution of the synchronization order parameter $r(t)$. The power and the usefulness of the ansatz lies in its remarkable ability to capture precisely and quantitatively through this single equation all, and not just some, of the order parameter attractors and bifurcations of the dynamics (3) (which may be obtained by performing numerical integration of the $N$ coupled non-linear equations (3) for $N\gg 1$ and evaluating $r(t)$ in numerics), for a Lorentzian $g(\omega)$. The success of the approach has led to hundreds of publications in applied mathematics and physics; A few recent ones are Refs. [43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55].

Within the dynamical setting of our model, we show that the system in the stationary state exhibits in the $(\Delta,K,T)$-space a very rich phase diagram exhibiting regions of synchronized and unsynchronized phases, and lines and surfaces of continuous transitions between them. The schematic phase diagram is shown in Fig. 1. In the backdrop of the highly non-trivial nature of the dynamics, as highlighted in previous paragraphs, exact results for the phase transitions are obtained in two contrasting limits, namely, (i) the limit $T\to 0$, and (ii) the limit $\Delta\to 0$. For case (i), the dynamics reduces to that of a non-linear dynamical system, and our exact analysis is based on the use of the Ott-Antonsen ansatz to derive a reduced set of nonlinear partial differential equations for the macroscopic evolution of the system. On the other hand, in the case of (ii), when the dynamics becomes that of a statistical system in contact with a heat bath and relaxing to a statistical equilibrium state, we derive our results by invoking the transfer matrix approach of the nearest-neighbor Ising model well known from theories of equilibrium statistical mechanics, and by extending it to consider non-local interactions. Besides offering a case study of exact analysis in many-body interacting systems, our work underlines the crucial role that non-local interactions may play in synchronizing systems in either destroying or enhancing the possibility of observing global synchrony in the system.

The paper is organized as follows. In Section 2, we give a precise definition of our model and summarize known results on stationary-state phase transitions observed in specific limits of the model. In Section 3, we give a detailed derivation of the phase diagram of the model in the $(\Delta,K)$-plane, while the same in the $(K,T)$-plane is discussed in Section 4. Simulation results on phase transitions for a general point in the $(\Delta,K,T)$-space are discussed in Section 5. The paper ends with conclusions and perspectives in Section 6. A rigorous proof that the dynamics of our model in the general case does not verify detailed balance is given in A, while B contains details of a numerical scheme to integrate the equations of motion of our model system.

2 Model

Our model given by Eqs. (17) and (18) is a variant of the Kuramoto model (3). To derive it, let us first consider a one-dimensional ($1d$) lattice of $N$ sites with periodic boundary conditions (site $j+N\equiv$ site $j$, with $j=1,2,\ldots,N$), where each site is occupied by a limit-cycle oscillator that is characterized completely by its angle $\theta_{j}\in[0,2\pi)$ and its natural frequency $\omega_{j}\in[-\infty,\infty]$. As in the Kuramoto model (3), the $\omega_{j}$’s are dimensionless numbers distributed according to a common unimodal distribution $g(\omega)$ that has zero mean and unit width, with the normalization $\int_{-\infty}^{\infty}{\rm d}\omega~{}g(\omega)=1$. The angles evolve in time according to the dynamics

[TABLE]

where the second term on the right hand side is the usual global coupling of the Kuramoto model, while the additional term that we introduce is the non-local ($M$-neighbor) interaction between the oscillators represented by the third term on the right hand side. Here, $K$ stands for the strength of coupling between an oscillator on a site with each of $M$ neighboring oscillators to the left and to the right, with $K>0$ (respectively, $K<0$) implying attractive (respectively, repulsive) interaction. Note that for $M=N/2$, the model (6) reduces to the Kuramoto model with a global coupling constant equal to $J+K$. Our aim in this work is to study the modification to the Kuramoto behavior due to non-local interactions, Hence, we consider the allowed range of values of $M$ to be satisfying $M<N/2$. Setting $K$ to zero allows to recover the Kuramoto model (3). In passing, we note that a dynamics similar to Eq. (6), but without the global coupling term and with a phase lag in the interaction (that is, having the $M$-neighbor interaction to be of the form $\sin(\theta_{j+k}-\theta_{j}-\alpha)$, with $\alpha\in(0,\pi/2)$ being the phase lag) has been analyzed in Ref. [43]; as is well known, the presence of a phase lag has important consequences on the behavior of the Kuramoto model [56]. Reference [52] considered a model similar to Eq. (6), but without the inclusion of the global coupling term. We note that as regards observing chimera states, the essential dynamical setup introduced, e.g., in Ref. [57], is given by the equation of motion

[TABLE]

with $\omega$ being the common natural frequency of the oscillators, and with $\alpha\in(0,\pi/2)$ being the phase-lag parameter. We thus observe an important difference in the form of the non-local interaction in (7) with respect to the dynamics (6), namely, the presence of the phase lag $\alpha$, whose inclusion has been argued to be crucial for observing the chimeras [36]. The issue of whether chimeras are observed in the dynamics (6) on including a phase lag in the non-local interaction and on making all the natural frequencies to be identical is an important question that is relegated to future studies.

A representative example of $g(\omega)$ that we specifically consider in this work to demonstrate our results and for which we obtain exact analytical results for relevant macroscopic properties is that of a Lorentzian distribution:

[TABLE]

Before proceeding, we rewrite Eq. (6) in a dimensionless form, by defining dimensionless quantities

[TABLE]

we get the dimensionless equation

[TABLE]

2.1 Relation to long-range interacting systems

We now establish a relation of the dynamics (10) to a specific limit of a certain Hamiltonian dynamics, which would prove quite useful later in the paper in studying the model. To this end, let us consider the Hamiltonian of a mean-field (classical) XY model in presence of additional non-local interactions on a $1d$ periodic lattice:

[TABLE]

where $p_{j}\equiv mv_{j}$ is the momentum conjugate to $\theta_{j}$ ($p_{j}$ and $\theta_{j}$ together constitute the set of canonically conjugate dynamical variables associated with the site $j$), $m$ is the mass, $v_{j}$ is the velocity, while $\widetilde{J}>0$ and $\widetilde{K}$ are respectively the global and the non-local ($M$-neighbor) coupling constant.

Note that in the Hamiltonian (11), the global-coupling term involves every spin interacting with every other with the same strength $\widetilde{J}$. Such an interaction is the extreme form (the mean-field limit) of the so-called long-range interactions exhibited by physical systems. Long-range interacting (LRI) systems are those in which the constituent particles interact with each other with a strength that decays slowly with their separation $r$ as $r^{-\alpha}$ for large $r$, with $0\leq\alpha\leq d$ in $d$ spatial dimensions [58, 59, 60, 61, 62]. LRI systems are encountered across disciplines, in astrophysics, hydrodynamics, plasmas, atomic and nuclear physics, and condensed matter physics. These systems are intrinsically non-additive so that they cannot be trivially divided into independent macroscopic sub-parts, a feature that leads to many fascinating phenomena not observed with short-range interactions, e.g., inequivalence of statistical ensembles [59, 60]. Other striking effects are breaking of ergodicity: the phase space is broken up into subspaces not connected by local dynamics. A very interesting dynamical feature of LRI systems is the occurrence of quasistationary states during relaxation to equilibrium. These states involve a slow relaxation of macroscopic observables over times that diverge algebraically with the system size, so that in the thermodynamic limit, the system remains trapped in them and never attains the Boltzmann-Gibbs equilibrium [63, 61, 60].

The Hamiltonian (11) has in addition to a long-range interaction a short-range one described by the coupling among $M$ nearest-neighbors. For small $M$, we expect the long-range behavior to dominate, and indeed, the equilibrium properties of the Hamiltonian (11) with $M=1$ and within microcanonical and canonical ensembles have demonstrated the feature of ensemble inequivalence emerging as a consequence of long-range interactions [64, 65].

In contact with a heat bath that induces noise into the system, and in presence of a friction constant $\gamma>0$, the dynamics derived from the Hamiltonian (11) and with additional external drives in the form of quenched disordered external toques $\widetilde{\Delta}~{}\omega_{j}$ acting on the individual spins is given by the set of equations

[TABLE]

where $\widetilde{\Delta}>0$ is a given parameter characterizing the strength of the external torques, while $\eta_{j}(t)$ is a Gaussian, white noise with

[TABLE]

Here, $T$ is the temperature of the heat bath in units of the Boltzmann constant, while angular brackets denote averaging over noise realizations. Let us now define the following dimensionless quantities:

[TABLE]

where note that for given values of $\widetilde{J}$ and $\widetilde{K}$, the ratio $\widetilde{K}/\widetilde{J}$ may not equal the quantity $\overline{K}$ as defined in the paragraph preceding Eq. (10); if this is the case, the equality may be achieved by multiplying both $\widetilde{J}$ and $\widetilde{K}$ by the same factor. Similarly, for given values of $\gamma,\widetilde{J}$ and $\widetilde{\Delta}$, the ratio $\gamma\widetilde{\Delta}/\widetilde{J}$ may not equal the quantity $\overline{\Delta}$ as defined in the paragraph preceding Eq. (10), and when this is the case, the equality may be achieved by multiplying both $\widetilde{J}$ and $\widetilde{\Delta}$ by the same factor. Using the definitions in Eq. (14), we obtain from Eq. (12) and in the limit $m/\gamma\ll 1$ the overdamped dynamics

[TABLE]

with

[TABLE]

From Eqs. (15) and (16), it is evident that as $\overline{T}\to 0$, the overdamped dynamics (15) reduces to Eq. (10).

On the basis of the foregoing discussions, we conclude that the general first-order dynamics that incorporates in specific limit the dynamics (10) is given by

[TABLE]

where all the quantities are dimensionless, and we have dropped the overbars in order not to overload our notation. Equations (17) and (18) define our model of interest in this work.

2.2 Summary of known results and our queries

We may ask: what is the nature of the stationary state that the dynamics (17) relaxes to at long times (i.e., in the limit $t\to\infty$)? For $\Delta=0$, the dynamics (17) relaxes to a Boltzmann-Gibbs (BG) equilibrium state, so that the probability distribution of the angles $\{\theta_{j}\}_{1\leq j\leq N}$ has the usual form

[TABLE]

with ${\cal V}$ being a potential energy function:

[TABLE]

for a proof, see A. For $\Delta\neq 0$, the dynamics (17) does not correspond to a Hamiltonian system because of the natural frequency term that cannot be derived from any potential $V_{\rm pot}(\{\theta_{j}\})$ that satisfies the periodicity of the system, namely, $V_{\rm pot}(\{\theta_{j}+2\pi\})=V_{\rm pot}(\{\theta_{j}\})$, and thus be intrinsic to the system. The external drives in the form of the natural frequencies continuously pump energy into the system. In this case, the dynamics at long times relaxes to a nonequilibrium stationary state [40], which does not have the BG form of angle distribution, and which violates detailed balance. The latter property is proven in A.

Note that the dynamics (17) is characterized by three dimensionless parameters $(\Delta,K,T)$, see Fig. 1. In this work, we consider the dynamics in the thermodynamic limit, and obtain its stationary-state phase diagram in the $(\Delta,K,T)$-space. It is pertinent to discuss the range of values of $M$ our results for the phase diagram apply to. To this end, let us define an interaction radius $\sigma$ as $\sigma\equiv M/N$. Note that our model (17) is to be considered for $M<N/2$, that is, for $\sigma<1/2$. Suppose one takes first the thermodynamic limit $N\to\infty$, and then gradually increase $M$ to larger and larger values while keeping $\sigma<1/2$. The limiting phase diagram that one gets as $M$ approaches infinity is what we obtain in this work, with the corresponding analysis for the $(\Delta,K)$-plane described in Section 3 and that for the $(K,T)$-plane described in Section 4. We will see that the actual value of $\sigma$ has only the role of a parameter that characterizes the dynamics of the system.

We note that, to the best of our knowledge, our quest for the complete phase diagram of the model (17) has not been addressed before. Only certain limits of the dynamics and the associated phase diagrams have been considered by two different communities of physicists, namely, the dynamical physicists and the statistical physicists. We now summarize these contributions.

Before proceeding, let us discuss qualitatively some general features of the dynamics. Considering Eq. (17) and a given frequency distribution $g(\omega)$, we note that the effect of the frequency term and the noise term is to induce every oscillator to oscillate at its natural frequency on an average, with thermal fluctuations superimposed on the average behavior. In contrast to the Kuramoto model (3), this tendency is now opposed by both the global and the non-local coupling among the oscillators. While the former favors equal angles for the oscillators, thereby promoting global synchrony among all the oscillators, the latter induces a local (that is, among $M$ neighboring oscillators to the left and to the right of a given oscillator) order. The latter can be either ferromagnetic or antiferromagnetic (borrowing terminologies from spin systems), depending respectively on whether $K$ is positive or negative. The interplay of these various tendencies ultimately determines whether synchrony among the oscillators is sustained in the stationary state and the amount of it. In this backdrop, we now summarize the known phase transitions exhibited by the dynamics (17).

• •

The case $K=T=0$ corresponds to the Kuramoto model, which is thus confined to the $\Delta$-axis, see Fig. 1. As already discussed, the system in the stationary state undergoes a continuous phase transition as a function of $\Delta$, from a low-$\Delta$ synchronized phase to a high-$\Delta$ incoherent phase at the critical threshold $\Delta_{c}(K=0,T=0)=\pi g(0)/2$ [23, 32, 26]. For the Lorentzian distribution (8), we find that $\Delta_{c}(K=0,T=0)=1/2$.

• •

The case $K=\Delta=0$ corresponds to the so-called Brownian mean-field (BMF) model [66], a set-up to study statics and dynamics of LRI systems in contact with an external heat bath. The underlying Hamiltonian is obtained from Eq. (11) by setting $\widetilde{K}$ to zero. The Hamiltonian describes a system of globally-coupled (classical) $XY$ spins, which allows to draw analogies with magnetic systems, and to refer to the corresponding stationary phases, the synchronized and the incoherent phase, as the magnetized and the unmagnetized phase, respectively. In equilibrium, the system exhibits a continuous transition between the two phases at the critical temperature $T_{c}(\Delta=0,K=0)=1/2$.

• •

The case $K=0,\Delta\neq 0,T\neq 0$ corresponds to the Kuramoto dynamics in presence of Gaussian, white noise, which was studied to account for stochastic fluctuations of the $\omega_{j}$’s in time [67]. In the stationary state, the transition point $\Delta_{c}(K=0,T=0)$, mentioned above, goes over to become a line of continuous transition between the synchronized and the incoherent phase, whose equation $\Delta_{c}=\Delta_{c}(K=0,T)$ is obtained by solving [67, 32]

[TABLE]

For the Lorentzian distribution (8), one may evaluate the integral on the right hand side of Eq. (21) by converting it to a complex integral, and then choosing a contour consisting of the real-$\omega$ axis closed on the lower-half complex-$\omega$ plane by an infinite semicircle on which the integral gives zero contribution. Evaluating the integral and using Eq. (21), one obtains

[TABLE]

From Eq. (22), it is easily checked that the line $\Delta_{c}=\Delta_{c}(K=0,T)$ has an intercept on the $T$-axis equal to $1/2$, in agreement with the phase transition point for the BMF model mentioned above.

The aforementioned points and lines of continuous transitions, with $r^{\rm st}$ as the order parameter, are indicated schematically in the $(\Delta,T)$-plane in Fig. 1. Note that these transitions all refer to the mean-field limit of the dynamics (17). In this paper, our primary objective is to investigate as to how this mean-field behavior is modified by the inclusion of the $M$-neighbor interaction. In other words, referring to Fig. 1, we ask: how does the phase diagram in the $(\Delta,T)$-plane extend to the whole of the $(\Delta,K,T)$-space ? In the rest of the paper, we use interchangeably the terms “magnetized” and “synchronized” to describe the clustered phase, and the terms “unmagnetized” and “homogeneous” for the unsynchronized/incoherent phase in model (17).

3 Analysis for the $(\Delta,K)$-plane with $T=0$

In this section, we consider the dynamics (17) with $T=0$, and turn to a discussion of its stationary state properties. To this end, let us introduce in the spirit of $r(t)$ a local synchronization order parameter $Z_{j}(t)$ defined at the $j$-th site by the equation [43]

[TABLE]

In terms of $Z_{j}(t)$ and $r(t)$, the equation of motion (17) with $T=0$ takes a form convenient for further analysis in the present section:

[TABLE]

where $\star$ denotes complex conjugation. From Eq. (24), it is evident that $r(t)$ (respectively, $Z_{j}(t)$) plays the role of a complex global (respectively, local) mean field, with both driving the dynamical evolution of the angles.

As mentioned in Section 2, in order to obtain our desired phase diagram, we consider the limits $N\to\infty$ and $M\to\infty$, keeping $\sigma=M/N<1/2$. In such a situation, it is reasonable and convenient to invoke a continuum limit of the dynamics in order to pursue its analytical treatment. The continuum limit corresponds to fixing the total length of the periodic lattice to be $2\pi$ and denoting the spatial location of the $j$-th site by $x_{j}\equiv 2\pi j/N$, so that as $N\to\infty$, the variable $x_{j}$ turns into a continuous variable $x\in[0,2\pi)$. In the continuum limit, the system is characterized by the probability density function $f(\theta,\omega,x,t)$, defined such that $f(\theta,\omega,x,t){\rm d}\theta{\rm d}\omega{\rm d}x$ gives the probability at time $t$ that an oscillator in position $[x,x+{\rm d}x]$ and with its natural frequency in $[\omega,\omega+{\rm d}\omega]$ has its angle in $[\theta,\theta+{\rm d}\theta]$. The density function satisfies $f(\theta+2\pi,\omega,x,t)=f(\theta,\omega,x,t)$, and the normalization

[TABLE]

In the continuum limit, the local mean field becomes

[TABLE]

with

[TABLE]

while the Kuramoto order parameter becomes

[TABLE]

In view of having a periodic spatial domain, all spatial integrals are to be evaluated using periodic boundary conditions. Note that one has $\int_{0}^{2\pi}{\rm d}x~{}Z(x,t)=2\pi~{}r(t)$.

The density $f(\theta,\omega,x,t)$ evolves in time according to the continuity equation that follows from the conservation of the total number of oscillators under the dynamics (24):

[TABLE]

Being $2\pi$-periodic in $\theta$, we expand the density $f(\theta,\omega,x,t)$ in a Fourier series in $\theta$:

[TABLE]

where $\widetilde{f}_{n}(\omega,x,t)$ is the $n$-th Fourier coefficient. Using $\int_{0}^{2\pi}{\rm d}\theta~{}e^{in\theta}=2\pi\delta_{n,0}$, we check that the above expansion satisfies Eq. (25).

We now implement the Ott-Antonsen (OA) ansatz that consists in restricting to the class of Fourier coefficients [41, 42]

[TABLE]

with $\alpha(\omega,x,t)$ an arbitrary function, and with $|\alpha(\omega,x,t)|<1$, so that the infinite series in Eq. (30) is converging. The OA ansatz also assumes that $\alpha(\omega,x,t)$ may be analytically continued to the whole of the complex-$\omega$ plane, that it has no singularities in the lower-half complex-$\omega$ plane, and that $|\alpha(\omega,x,t)|\to 0$ as ${\rm Im}(\omega)\to-\infty$ [41, 42].

Using Eqs. (30) and (31) in Eqs. (26) and (28), we obtain

[TABLE]

Equations (30) and (31) and the above expressions for $r(t)$ and $Z(x,t)$ on substituting in Eq. (29), and then on collecting and equating the coefficient of $e^{in\theta}$ to zero yield

[TABLE]

For the Lorentzian $g(\omega)$, Eq. (8), one may evaluate $r(t)$ and $Z(x,t)$ by using Eq. (8) in Eqs. (32) and (33) to get

[TABLE]

and

[TABLE]

where the contour $C$ is shown in Fig. 2, and we have used the fact that the contribution to the integral from the semicircular part of the contour vanishes in view of $|\alpha(\omega,x,t)|\to 0$ as ${\rm Im}(\omega)\to-\infty$. Evaluating the integrals in Eqs. (35) and (36) by the residue theorem, we get

[TABLE]

Calling $\alpha(-i,x,t)=u(x,t)$, Eq. (34) then gives

[TABLE]

The stationary solution $u_{\rm st}(x)$ of the above equation satisfies

[TABLE]

with

[TABLE]

3.1 Uniformly incoherent state: Stability

The uniformly incoherent state $u^{\rm inc}_{\rm st}(x)=0~{}\forall~{}x$, yielding $r_{\rm st}=0$ and $Z_{\rm st}(x)=0~{}\forall~{}x$, evidently satisfies Eq. (40), and is thus a stationary solution of Eq. (39). Let us study the linear stability of such a state by linearizing Eq. (39) about the state. To this end, we write

[TABLE]

which on using in Eq. (39) yields to leading order in $\delta u$ the equation

[TABLE]

with

[TABLE]

Let us use the expansion $\delta u(x,t)=a(q)e^{iqx}e^{\lambda t}$, with real $\lambda$, and with the wave number $q$ being an integer (in view of having a periodic spatial domain). Substituting in Eq. (43), and using $r(t)=1/(2\pi)\int_{0}^{2\pi}{\rm d}x~{}u^{\star}(x,t)$, see Eq. (37), we get $r(t)=1/(2\pi)\int_{0}^{2\pi}{\rm d}x~{}a^{\star}(q)e^{iqx}e^{\lambda t}=a^{\star}(q=0)e^{\lambda t}$, so that $r(t)$ being real implies that $a(q)$ has to be real. On using $\delta u(x,t)=a(q)e^{iqx}e^{\lambda t}$ in Eq. (44), we get the spectral equation determining the parameter $\lambda$:

[TABLE]

where

[TABLE]

is the Fourier transform of the coupling function $G$; note that we have $\widetilde{G}(0)=1$. In arriving at Eq. (48), we have used the fact that $G$ is an even function of its argument, $G(x)=G(-x)$, see Eq. (27).

Now, note that $\lambda(q)$ in Eq. (47) is real. Depending on whether $\lambda(q)$ is larger or smaller than zero makes the perturbation $\sim e^{iqx}$ grow or decay in time, respectively; the threshold between the two behaviors is obtained by setting $\lambda(q)$ to zero in Eq. (47), thereby obtaining for a fixed $\Delta$ the threshold

[TABLE]

In particular, one has $K_{c,{\rm inc}}^{(0)}(\Delta)=2\Delta-1$. Equation (47) may be rewritten as

[TABLE]

which implies that for values of $q$ such that $\widetilde{G}(q)>0$, the perturbation $\sim e^{iqx}$ grows and is thus sustained (respectively, decays, and is thus non-sustained) in time for $K>K_{c,{\rm inc}}^{(q)}(\Delta)$ (respectively, $K<K_{c,{\rm inc}}^{(q)}(\Delta)$). On the other hand, for values of $q$ such that $\widetilde{G}(q)<0$, the perturbation $\sim e^{iqx}$ grows and is thus sustained (respectively, decays, and is thus non-sustained) in time for $K<K_{c,{\rm inc}}^{(q)}(\Delta)$ (respectively, $K>K_{c,{\rm inc}}^{(q)}(\Delta)$). For values of $q$ such that $\widetilde{G}(q)=0$, Eq. (49) yields infinite $K_{c,{\rm inc}}^{(q)}$; these modes therefore do not exist for any finite $K$. The boundaries $K_{c,{\rm inc}}^{(q)}(\Delta)$ are shown in Fig. 3. From the figure, it is evident that the uniformly incoherent state is stable with respect to perturbations $\sim e^{iqx}~{}\forall~{}q$ so long as $|K|$ is small and lies in the central region around zero that has no overlap with any of the bounded regions. In fact for $K>0$, the incoherent state destabilizes as a whole when it destabilizes with respect to perturbations with wave number $q=0$. On increasing $K$, the uniformly incoherent state becomes unstable with respect to a long-wavelength perturbation ($q=0$) at $K=2\Delta-1$. On decreasing $K$, however, the instability is with respect to a perturbation with a shorter wavelength ($q\neq 0$).

3.2 Synchronized twisted state

3.2.1 Existence

Let us look for solutions of Eq. (39) given by plane waves, or, the so-called (partially) synchronized uniformly twisted states [43], which has the form $u(x,t)=ae^{i(qx+\nu t)}$, with real $a,\nu$, and integer $q$. Here, the wave number $q$ characterizes the “twist” of the state, giving the rate of change of angle with $x$ at a fixed $t$, while $\nu$ measures the temporal rate of rotation of the twisted state. The parameter $a$ with $0<a<1$ measures the level of coherence between the oscillator angles. We now obtain the conditions on the parameters $a,q,\nu$ for such a twisted state to be a solution of Eq. (39). Substituting $u(x,t)=ae^{i(qx+\nu t)}$ in Eq. (39) gives

[TABLE]

Equating for real and imaginary parts from both sides of Eq. (51), and using Eq. (48) and $\int_{0}^{2\pi}{\rm d}x^{\prime}~{}e^{iq(x^{\prime}-x)}=2\pi\delta_{q,0}$, we get

[TABLE]

The second equation implies that the synchronized uniformly twisted state is actually a stationary solution of Eq. (39). On the other hand, the first equation gives

[TABLE]

and

[TABLE]

Requiring $a(q=0)$ to be real implies that a zero-twist state is a stationary solution of Eq. (39) provided that for fixed $\Delta$, the parameter $K$ is larger than the critical value $K_{c}^{(0)}\equiv 2\Delta-1$. On the other hand, Eq. (54) implies that a $q\neq 0$-twist state is stationary for $K>K_{c}^{(q)}\equiv 2\Delta/\widetilde{G}(q)$ for $\widetilde{G}(q)>0$ and for $K<K_{c}^{(q)}$ for $\widetilde{G}(q)<0$. Comparing with Eq. (49), we see that $K_{c}^{(q)}=K_{c,{\rm inc}}^{(q)}$, and thus, the instability boundaries for the uniformly incoherent state shown in Fig. 3 are also the existence boundaries for the twisted states. In other words, the twisted state of wave number $q$ emerges as the incoherent state destabilizes with respect to perturbations of the same wave number. When exists, the synchronized uniformly twisted stationary state is represented as $u^{\rm syn}_{\rm st}(x)=a(q)e^{iqx}$, with $a(q)$ given by Eqs. (54) and (55).

Corresponding to the zero-twist stationary state, one has the stationary value $r^{\rm st}=a(q=0)$, while a $q\neq 0$-twist state yields the stationary value $r^{\rm st}=0$, where $a(q=0)$ is given by Eq. (55). On the other hand, Eq. (38) yields the stationary value

[TABLE]

with $a(q)$ given by Eq. (54).

3.2.2 Stability

In this section, we study the linear stability of the twisted stationary state, whose existence has been considered in Section 3.2.1. To this end, consider a twisted state with a given wave number $q_{0}$: $u^{\rm syn}_{\rm st}(x)=a(q_{0})e^{iq_{0}x}$. Using $r^{\rm st}=a(q_{0}=0)\delta_{q_{0},0}$, $Z^{\rm st}(x)=a(q_{0})e^{-iq_{0}x}\widetilde{G}(q_{0})$, with $a(q_{0})$ and $a(q_{0}=0)$ given respectively by Eqs. (54) and (55), and writing $u(x,t)$ as

[TABLE]

Eq. (39) yields to leading order in $\delta u$ the equation

[TABLE]

where we have

[TABLE]

Substituting $u^{\rm syn}_{\rm st}(x)=a(q_{0})e^{iq_{0}x}$, and introducing $\delta U(x,t)\equiv\delta u(x,t)e^{-iq_{0}x}$, Eq. (58) yields

[TABLE]

Let us introduce the column matrix

[TABLE]

in terms of which Eq. (61) may be rewritten as

[TABLE]

with

[TABLE]

Next, we seek solutions to Eq. (74) of the form

[TABLE]

Substituting in Eq. (74), and using the identities

[TABLE]

we get

[TABLE]

with $\mathbb{I}$ being the $2\times 2$ identity matrix, and

[TABLE]

On requiring the existence of the solution (79), we obtain from Eq. (95) the secular equation determining $\lambda$:

[TABLE]

where we have

[TABLE]

The solutions of Eq. (100) may be written as

[TABLE]

From Eq. (112), it follows that

[TABLE]

and hence, we get

[TABLE]

Choosing $K>K_{c,{\rm inc}}^{(q_{0})}=2\Delta/\widetilde{G}(q_{0})$ for $\widetilde{G}(q_{0})>0$ or $K<K_{c,{\rm inc}}^{(q_{0})}=2\Delta/\widetilde{G}(q_{0})$ for $\widetilde{G}(q_{0})<0$, so that the twisted state $\sim e^{iq_{0}x}$ under consideration exists in the first place, we may evaluate $\lambda_{\pm}(q)$ by using Eq. (113); having a positive (respectively, negative) $\lambda_{\pm}(q)$ implies that the perturbation $\sim e^{iqx}$ to the twisted state grows (respectively, decays) and is thus sustained (respectively, non-sustained) in time. Note that we have $\lambda_{+}(q)>\lambda_{-}(q)$.

3.2.3 Stability of the zero-twist state

For the particular case $q_{0}=0,q\neq 0$, we need to choose $K>K_{c,{\rm inc}}^{(0)}=2\Delta-1$. In this case, Eq. (120) reduces to

[TABLE]

so that

[TABLE]

and hence,

[TABLE]

Using Eq. (55) and the fact that $K>2\Delta-1$, and that $\widetilde{G}(q\neq 0)<1$, we may write $\lambda_{+}(q\neq 0)<\Lambda\equiv-\Delta+K\frac{\Delta}{1+K}$; simplifying, we get $\Lambda=-\frac{\Delta}{1+K}$. Using again $K>2\Delta-1$, and noting that $\Delta$ is a positive quantity, it follows that $1+K$ is also positive, and hence, $\Lambda$ is negative, implying that $\lambda_{+}(q\neq 0)<0$.

On the other hand, for $q=0$, we have

[TABLE]

and hence, we get

[TABLE]

where we have used Eq. (55) and the fact that $K>K_{c,{\rm inc}}^{(0)}=2\Delta-1$; we thus have

[TABLE]

Using the facts $\lambda_{+}(q\neq 0)<\lambda_{+}(0)$ and $\lambda_{+}(q)>\lambda_{-}(q)$, we conclude that the stability of the zero-twist state is determined by the behavior of $\lambda_{+}(0)$. Namely, at a fixed $\Delta$, and for $K>2\Delta-1$, we have $\lambda_{+}(0)=0$, so that referring to Fig. 3, we see that the zero-twist state stabilizes as soon as the incoherent one destabilizes. At a fixed $\Delta$, the transition between the two states takes place at $K=K_{c,{\rm inc}}^{(0)}=2\Delta-1$.

3.3 Phase transition in the Kuramoto order parameter

On the basis of our discussions in Sections 3.1 and 3.2, we conclude that the stationary state Kuramoto order parameter $r^{\rm st}$ undergoes a continuous transition at a fixed $\Delta$, from a low-$K$ zero value, corresponding to an incoherent state, to a high-$K$ non-zero value $r^{\rm st}=\sqrt{1-\frac{2\Delta}{1+K}}$, corresponding to a zero-twist synchronized state, at the critical threshold

[TABLE]

The line of transition $K_{c}(T=0,\Delta)$ is shown in the phase diagram, Fig. 1. The line intercepts the $\Delta$-axis at the point $\Delta=1/2$, which matches with the prediction for this point made on the basis of the analysis of the bare Kuramoto model, that is, in the absence of any non-local interactions, see the discussions in Section 2.2. Note that we may write $r^{\rm st}=\sqrt{\frac{K-K_{c}(T=0,\Delta)}{1+K}}$, from which we obtain as $K\to K_{c}^{+}(T=0,\Delta)$ the scaling $r^{\rm st}\sim\left(K-K_{c}(T=0,\Delta)\right)^{\delta}$, where the critical exponent $\delta$ has the value $\delta=1/2$.

4 Analysis for the $(K,T)$-plane with $\Delta=0$

In this section, we discuss the phase diagram of the model (17) in the $(K,T)$-plane, i.e., for $\Delta=0$. We first consider the case of finite $M$, present an analysis of the phase diagram as a function of $M$, and, in the end, consider the limit $M\to\infty$ of the results.

As discussed in Section 2, the dynamics (17) for $\Delta=0$ relaxes at long times to an equilibrium stationary state with the BG distribution for the angles, Eq. (19). The canonical partition function for a $1d$ periodic chain of $N$ sites is thus given by

[TABLE]

with $\beta\equiv 1/T$.

Next, using the Hubbard-Stratonovich transformation,

[TABLE]

in Eq. (137), and introducing auxiliary fields $\widetilde{z}_{1}\equiv\beta z_{1}$ and $\widetilde{z}_{2}\equiv\beta z_{2}$, we obtain

[TABLE]

Writing $z_{1}=z\cos\phi,z_{2}=z\sin\phi$, with real $z=(z_{1}^{2}+z_{2}^{2})^{1/2}>0$ and $\phi\in[0,2\pi)$ given by $\phi=\tan^{-1}(z_{2}/z_{1})$, we get

[TABLE]

In view of the invariance of the potential (20) under rotation by an equal amount of all the $\theta_{j}$’s, we get

[TABLE]

Below, we consider separately the cases $K=0$ and $K\neq 0$.

4.1 $K=0$

For $K=0$, Eq. (141) yields

[TABLE]

In the thermodynamic limit, $Z_{N}$ may be approximated by invoking the saddle-point method to perform the integration in $z$ on the right hand side; one gets

[TABLE]

where the saddle-point value $z_{s}$ solves the equation

[TABLE]

where $I_{n}(x)=(1/(2\pi))\int_{0}^{2\pi}{\rm d}\theta~{}\exp(x\cos\theta)\cos(n\theta)$ is the modified Bessel function of first kind and order $n$. It may be shown that $z_{s}$ is nothing but the synchronization order parameter $r^{\rm st}$, see Section 4.2 below. Equation (144) has a trivial solution $r^{\rm st}=0$ valid at all temperatures, while a non-zero solution exists for $\beta\geq\beta_{c}=2$ [60]. In fact, the system shows a continuous transition, from a synchronized/magnetized phase ($r^{\rm st}\neq 0$) at low temperatures to an incoherent/unmagnetized phase ($r^{\rm st}=0$) at high temperatures at the critical temperature $T_{c}=1/2$. The latter point coincides with the BMF phase transition point $T_{c}(\Delta=0,K=0)$ indicated in the phase diagram in Fig. 1; this is a consequence of the fact that the BMF model and our model with $\Delta=K=0$ have the same distribution of the angles in equilibrium given by Eqs. (19) and (20) with $K=0$ [66].

4.2 $K\neq 0$

For $K\neq 0$, Eq. (141) gives

[TABLE]

Here, we may identify the factor ${\cal Z}_{N}$ with the canonical partition function of a $1d$ periodic chain of $N$ phase-only oscillators (equivalently, classical $XY$-spins), where each oscillator interacts with strength $K/(4M)$ with $M$ neighboring oscillators to the left and to the right, and also with an external field of strength $z$ along the $x$ direction. Note from Eq. (146) that under $z\to-z$, one has

[TABLE]

so that on using the invariance of the potential (20) under rotation of all the $\theta_{j}$’s by an amount equal to $\pi$, we obtain that ${\cal Z}^{\prime}_{N}={\cal Z}_{N}$. We thus conclude that the factor ${\cal Z}_{N}$ is an even function of $z$.

Invoking the above mentioned analogy with the $1d$ periodic chain of oscillators, we now proceed to compute the factor ${\cal Z}_{N}$ for large $N$. Our approach is based on a combination of a matrix formulation that was developed to study a general spin model in $1d$ with an $n$-neighbor interaction, with $n$ arbitrary and finite [68, 69], and a transfer operator method that generalizes the well-known transfer matrix approach for Ising spins [39] to the case of continuous spins [70].

The starting point is to consider a $1d$ periodic chain of total number of sites equal to ${\cal N}M$, and then to divide it into ${\cal N}$ blocks of $M$ sites. Let us relabel the oscillators (equivalently, the sites accommodating them), such that $\theta_{j}^{(\alpha)}$ refers to the angle of the $j$-th oscillator within the $\alpha$-th block, with $\alpha=1,2,\ldots,{\cal N}$, and $j=1,2,\ldots,M$. By virtue of such a construction, an oscillator in the $\alpha$-th block interacts with the oscillators in the same block and with those in the $(\alpha-1)$-th and $(\alpha+1)$-th blocks. As a result, the total energy of any configuration of the oscillator angles may be expressed as a sum of (i) energies due to interaction of the oscillators with the external field of strength $z$, (ii) interaction energies of oscillators within the same block, and (iii) interaction energies of oscillators from adjacent blocks. Next, let us denote the configuration of the $\alpha$-th block by $C_{\alpha}\equiv\{\theta_{j}^{(\alpha)}\}_{1\leq j\leq M}$. With this notation, we may express the energy of the system in configuration $C\equiv(C_{1},C_{2},\ldots,C_{\cal N})$ as

[TABLE]

where $X_{C_{\alpha}}$ denotes the energy contribution due to interaction of types (i) and (ii), and $Y_{C_{\alpha},C_{\alpha+1}}$ denotes the energy due to interaction of type (iii) contributed by the oscillators in the $\alpha$-th and $(\alpha+1)$-th blocks:

[TABLE]

The transfer operator method [70] introduces an operator ${\cal T}(C,C^{\prime})$ as

[TABLE]

Let $\{\lambda_{q}\}$ be the set of eigenvalues 333Note that the operator ${\cal T}(C,C^{\prime})$ is not symmetric in $(C,C^{\prime})$, so that one has to distinguish between its left and right eigenvalues and eigenvectors. Here, $\lambda_{q}$’s refer to the set of the right eigenvalues of ${\cal T}(C,C^{\prime})$. of the transfer operator ${\cal T}(C,C^{\prime})$. In other words, denoting the eigenfunctions of ${\cal T}(C,C^{\prime})$ as $f_{q}(C)$, we have

[TABLE]

In terms of $\{\lambda_{q}\}$, we obtain the canonical partition for a $1d$ ring of ${\cal N}M$ sites as

[TABLE]

where we have $\int{\rm d}C_{\alpha}\equiv\int(\prod_{j=1}^{M}){\rm d}\theta^{(\alpha)}_{j}$. For large ${\cal N}$, the sum in Eq. (152) is dominated by the largest eigenvalue $\lambda_{\rm max}=\lambda_{\rm max}\left(\beta z,\frac{\beta K}{M}\right)$, yielding 444For any infinitesimal discretization of the $\theta_{j}$’s, the operator ${\cal T}(C,C^{\prime})$ becomes a finite-dimensional real square matrix with positive entries, so that the application of the Perron-Frobenius theorem [71] implies the existence of the largest eigenvalue that is real and non-degenerate.

[TABLE]

For our system of interest, Eq. (17), we have ${\cal N}M=N$, giving

[TABLE]

Equation (154) when combined with the fact shown earlier that ${\cal Z}_{N}$ is an even function of $z$ implies that $\lambda_{\rm max}$ is an even function of $z$. Substituting Eq. (154) in Eq. (145), we obtain in the thermodynamic limit the result

[TABLE]

In the same limit, one may further approximate $Z_{N}$ by invoking the saddle-point method to perform the integration in $z$; one gets

[TABLE]

where $z_{s}$ solves the saddle-point equation

[TABLE]

with $\widetilde{\phi}(\beta,z)$ being the free-energy function:

[TABLE]

The saddle-point equation may thus be written as

[TABLE]

From Eq. (156), one obtains the dimensionless free energy per oscillator, $\phi(\beta)\equiv-\lim_{N\to\infty}(\ln Z_{N})/N$, as

[TABLE]

where we have suppressed the dependence of $\phi(\beta)$ on $K$. We thus have

[TABLE]

Note that the free energy at a given temperature has a definite value given by Eq. (161), and is obtained by substituting the saddle-point solution $z_{s}$ into the expression for the free-energy function $\widetilde{\phi}(\beta,z)$.

As it turns out, the quantity $z_{s}$ in Eq. (159) is nothing but the stationary Kuramoto order parameter $r^{\rm st}$. To demonstrate that this is the case, consider the dynamics (17) with $\Delta=0$ and in presence of an additional potential $V_{\rm ext}(\{\theta_{j}\})\equiv-h\sum_{j=1}^{N}\cos\theta_{j}$ due to an external field of strength $h$ along the $x$-direction, so that the partition function (137) is modified to $Z_{N}^{(h)}\equiv\int\left(\prod_{j=1}^{N}{\rm d}\theta_{j}\right)\exp[-\beta\{{\cal V}(\{\theta_{j}\})+V_{\rm ext}(\{\theta_{j}\})\}]$. In this case, one obtains in the same way as one arrives at Eqs. (159) and (161) the following analogous equations

[TABLE]

where note that $z_{s}$ in Eq. (163) is a function of $\beta h$ and $\beta K/(4M)$ by virtue of Eq. (162). On the other hand, the stationary Kuramoto order parameter in presence of the field $h$ has values $r_{x}^{\rm st}(h)\neq 0$ and $r_{y}^{\rm st}(h)=0$, so that one obtains for $r^{\rm st}(h)=r_{x}^{\rm st}(h)$ in the thermodynamic limit

[TABLE]

where in obtaining the third equality, we have used the result that $\phi(\beta,h)=-\lim_{N\to\infty}(\ln Z_{N}^{(h)}))/N$, while in obtaining the last two equalities, we have used Eqs. (162) and (163). Comparing Eqs. (162) and (164), we conclude that $r^{\rm st}(h)=z_{s}$; It is evident from the derivation of this result that it holds for all values of $h$, including $h=0$. We have thus established the assertion made above that the quantities $z_{s}$ and $r^{\rm st}$ are identical, so that we may rewrite Eq. (159) as

[TABLE]

In line with our set-out objective of obtaining the phase diagram, we now need to solve Eq. (165) for $r^{\rm st}$ as a function of $\beta$ and $K$, in the limit $M\to\infty$. One has to then first compute the largest eigenvalue $\lambda_{\rm max}$ for finite $M$, then solve Eq. (165) for $r^{\rm st}$, thereby obtaining the phase diagram in the $(K,T)$-plane for the given value of $M$, and, finally, take the limit $M\to\infty$ of the results. A roadblock in pursuing this program is the analytic computation of $\lambda_{\rm max}$ for general $M$. One may alternatively estimate $\lambda_{\rm max}$ numerically by discretizing the angles over the interval $[0,2\pi)$, for example, as $\theta_{j}^{(a_{j})}=a_{j}\Delta\theta$, with $a_{j}=1,2,\ldots,P$ and $\Delta\theta=2\pi/P$ for any large positive integer $P$; The operator ${\cal T}(C,C^{\prime})$ then takes the form of a matrix of size $P^{M}\times P^{M}$, which even for any reasonable value of $P$ and for not-so-large $M$ becomes numerically quite unmanageable in order that we may reliably estimate using Eq. (165) the phase diagram in the $(K,T)$-plane, leave alone the limit that is of interest to us, namely, the limit $M\to\infty$. In order to gain insights into the phase transitions for finite $M$ and their limiting behavior as $M\to\infty$, it proves insightful to consider an equivalent Ising system for which the operator ${\cal T}(C,C^{\prime})$ is much more manageable numerically, and, consequently, the estimation of $\lambda_{\rm max}$ is simpler and reliable, as we demonstrate in Section 4.2.1. Our subsequent analysis and line of argument will proceed along the following directions. We will first show for the equivalent Ising system that in the limit $M\to\infty$, the phase diagram in the $(K,T)$-plane obtained numerically using the approach of the transfer operator may be derived analytically by considering the model in the mean-field approximation. This observation hints at an apparent mean-field dominance in dictating the stationary properties of the equivalent Ising system in the limit $M\to\infty$. Assuming a similar mean-field dominance to also be at work for the oscillator problem at hand, Eq. (17) with $\Delta=0$, we then perform an explicit mean-field approximation of the model in equilibrium to determine its phase diagram in the $(K,T)$-plane, and show that the results are fully consistent with the phase diagram in the $(\Delta,K)$-plane derived in Section 3.

4.2.1 An equivalent Ising problem

To define an equivalent Ising problem, consider a setting similar to the one for the dynamics (17), namely, a $1d$ lattice of $N$ sites with periodic boundary conditions, where we take each site to be occupied by an Ising spin $S_{j}=\pm 1$. There is an all-to-all ferromagnetic coupling between the spins. Additionally, each spin interacts with strength $K/(4M)$ with $M$ neighboring spins to the left and to the right, with $M<N/2$. The coupling $K$ can be of either sign, with $K>0$ (respectively, $K<0$) implying a ferromagnetic (respectively, an antiferromagnetic) $M$-neighbor interaction. The Hamiltonian of the Ising system comprises just the potential energy ${\cal V}_{\rm Ising}(\{S_{j}\})$ given by

[TABLE]

The two terms on the right hand side of Eq. (166) are just the Ising analog of the corresponding terms on the right hand side of Eq. (20), obtained by replacing the continuous variables $\theta_{j}$ in the latter with the discrete Ising variables $S_{j}$. For $M=1$, the model (166) reduces to the Ising model with nearest-neighbor and long-range interactions studied in Refs. [72, 73, 74, 75]. Moreover, in Eq. (166), setting $K$ to zero allows to recover the mean-field Ising model (the zero-field Curie-Weiss model of ferromagnet) [39, 76]. The magnetic order in the system is characterized by the magnetization

[TABLE]

Similar to Eq. (19), the Ising system has the BG equilibrium distribution given by $P_{\rm eq}(\{S_{j}\})\propto\exp[-{\cal V}_{\rm Ising}(\{S_{j}\})/T]$, so that the canonical partition function reads

[TABLE]

Invoking the Hubbard-Stratonovich transformation (138), and proceeding similarly to the analysis presented in Sections 4 and 4.2, one obtains

[TABLE]

The above equations are equivalents of Eqs. (145) and (146). Similar to Eq. (146), one may interpret the factor ${\cal Z}_{N}^{\rm Ising}$ in Eq. (170) as the canonical partition function of a $1d$ periodic chain of $N$ Ising spins, with each spin interacting with strength $K/(4M)$ with $M$ neighboring spins to the left and to the right, and also with an external field of strength $z$. Employing such an analogy, one may proceed to evaluate the factor ${\cal Z}_{N}^{\rm Ising}$, by following the same line of analysis involving the transfer operator that was pursued in Section 4.2 to evaluate the factor ${\cal Z}_{N}$. A difference that arises in the present case of Ising spins with respect to the oscillator case is that the transfer operator now takes the form of a $2^{M}\times 2^{M}$ matrix ${\cal T}^{\rm Ising}$, with elements given by

[TABLE]

Noting that ${\cal T}^{\rm Ising}_{C,C^{\prime}}$ is a finite-dimensional real square matrix with positive entries, the application of the Perron-Frobenius theorem [71] implies the existence of its largest eigenvalue $\lambda^{\rm Ising}_{\rm max}$ that is real and non-degenerate. Hence, similar to Eqs. (154) and (156), we obtain in the limit $N\to\infty$ the result ${\cal Z}_{N}^{\rm Ising}=\left(\lambda^{\rm Ising}_{\rm max}\right)^{N/M}$, and consequently,

[TABLE]

where $z_{s}$ solves a saddle-point equation that has the same form as in Eq. (159). The dimenionsionless free energy has the form of Eq. (158). Moreover, following the arguments given in Section 4.2 to show that the quantity $z_{s}$ in the saddle-point equation (159) is nothing but the stationary Kuramoto order parameter $r^{\rm st}$, it may be shown that the quantity $z_{s}$ in Eq. (175) coincides with the stationary magnetization $m^{\rm st}$.

The basic program to identify the phase transition point $T_{c}(K,M)$ for given values of $K$ and $M$ is as follows:

• •

For a given value of the temperature $T$, we first form the matrix ${\cal T}^{\rm Ising}$ in Eq. (171), and then compute numerically its largest eigenvalue $\lambda_{\rm max}$ by invoking the so-called power method [77]; to this end, we employ a numerically efficient code that implements the method 555A FORTRAN90 library that implements the power method and is distributed under the GNU LGPL license is available at http://people.sc.fsu.edu/~jburkardt/f_src/power_method/power_method.html.

• •

We then compute the free-energy function $\widetilde{\phi}(\beta,z)$ as a function of $z$ by using Eq. (158).

• •

We repeat the last two steps for several values of $T$, locating numerically for each $T$ the value of $z$ at which $\widetilde{\phi}(\beta,z)$ is minimum. Because of the symmetry of our problem, non-zero minimizers of $\widetilde{\phi}(\beta,z)$, if and when they exist, always occur in pairs symmetrically disposed on either side of zero: $z_{s}=\pm A$, with $0<A<1$.

• •

A continuous phase transition point is given by the value of $T$ at which the two non-zero minimizers occurring at lower temperatures merge with each other for the first time, so that the only minimizer at higher temperatures is at $z_{s}=0$, see Fig. 4.

• •

A first-order phase transition point is given by the value of $T$ at which there are three minimizers at $z_{s}=\pm A$ and at $z_{s}=0$, with $0<A<1$, such that the values of $\widetilde{\phi}(\beta,z)$ at these three minima coincide, see Fig. 4.

Following the above program, the phase diagram in the $(K,T)$-plane is reported in Fig. 5. At a fixed $K$, as $T$ is increased, the system undergoes a phase transition from a low-$T$ synchronized phase ($r^{\rm st}\neq 0$) to a high-$T$ unsynchronized phase ($r^{\rm st}=0$). The phase transition is continuous for smaller (in magnitude) values of $K$, and is of first order for larger (in magnitude) values of $K$, with the two separated by a tricritical point indicated in the figure. For $M=1$, one may solve exactly for the line of continuous transition and the tricritical point, which are respective given by $1/T=\exp\left(-K/T\right)$ and $K_{\rm CTP}=-(\ln 3)/(2\sqrt{3})$ [72, 73]; Note that for $K=0$, the equation $1/T=\exp\left(-K/(2T)\right)$ gives the phase transition point of the mean-field Ising model as $T_{c}=1$ [76]. From Fig. 5, it is evident that with increase of $M$, (i) the tricritical point approaches the $K$-axis, and (ii) the phase boundary approaches the line $T=1+K$, the latter being therefore the $M\to\infty$ limit of the phase boundary. We may thus conclude that in the limit $M\to\infty$, the model (166) exhibits a continuous transition at all finite temperatures and a first-order transition only at zero temperature, with the phase boundary given by

[TABLE]

The above equation yields the correct value for the phase transition point for $K=0$, and also gives a phase transition point at zero temperature for $K=-1$.

The result (176) may be understood physically as follows. On utilizing the translational invariance of the system, and on assuming pair factorization $\langle S_{j}S_{k}\rangle_{\rm eq}=\langle S_{j}\rangle_{\rm eq}\langle S_{k}\rangle_{\rm eq}$ in the joint limit $N\to\infty,M\to\infty$ 666The pair factorization is expected to hold exactly for a purely mean-field model [39]., Eq. (166) gives the average energy density in equilibrium as

[TABLE]

The above expression for the equilibrium energy density coincides (up to an inconsequential constant term) with that for the mean-field Ising model with coupling constant $J_{\rm eff}^{\rm mean-field~{}Ising}\equiv 1+K$. At any non-zero temperature, the mean-field Ising model exhibits a continuous transition as a function of temperature, from a low-$T$ magnetized phase to a high-$T$ unmagnetized phase at the critical temperature $T_{c}^{\rm mean-field~{}Ising}\equiv J_{\rm eff}^{\rm mean-field~{}Ising}=1+K$ [39, 76], from which one readily obtains Eq. (176). At zero temperature, one has a first-order transition between a magnetized phase in which all the spins are aligned parallel to each other, and a non-magnetized phase in which neighboring spins point in opposite directions.

On the basis of the above analysis for the equivalent Ising problem, (166), we may anticipate for the model (17) with $\Delta=0$ that in the limit $M\to\infty$, the model exhibits a continuous transition at all temperatures $T>0$ and a first-order transition at $T=0$, with the phase boundary given by

[TABLE]

Similar to the Ising case considered in the preceding paragraph, the above equation may actually be derived by considering the equilibrium average of Eq. (20), and by utilizing translational invariance and assuming pair factorization $\langle\sin\theta_{j}\sin\theta_{k}\rangle_{\rm eq}=\langle\sin\theta_{j}\rangle_{\rm eq}\langle\sin\theta_{k}\rangle_{\rm eq}$ and $\langle\cos\theta_{j}\cos\theta_{k}\rangle_{\rm eq}=\langle\cos\theta_{j}\rangle_{\rm eq}\langle\cos\theta_{k}\rangle_{\rm eq}$, to obtain the average energy density in equilibrium as

[TABLE]

Here, we have used $(r^{\rm st})^{2}\equiv(r^{\rm eq})^{2}=\langle r_{x}\rangle_{\rm eq}^{2}+\langle r_{y}\rangle_{\rm eq}^{2}$, as follows from Eq. (4). Up to an irrelevant constant term, Eq. (179) is the same as the equilibrium energy density of the mean-field XY model with effective coupling constant $J_{\rm eff}^{\rm mean-field~{}XY}\equiv 1+K$, and which exhibits a continuous phase transition at the critical temperature $T_{c}^{\rm mean-field~{}XY}\equiv J_{\rm eff}^{\rm mean-field~{}XY}=(1+K)/2$. Equation (178) yields consistently and correctly the phase transition point for $K=0$, namely, $T_{c}(\Delta=0,K=0)=1/2$, the phase transition point of the BMF model, see Fig. 1. The difference between the Ising and the Kuramoto model in that the former deals with discrete variables while the latter with continuous variables is reflected in the appearance of an extra factor of two in Eq. (178) with respect to Eq. (176). Let us note in passing that considering an Ising ferromagnet with pair-wise interactions and the classical XY model with the same couplings, the critical inverse temperatures in the two cases have been proved to satisfy $\beta_{c}^{\rm XY}\geq\beta_{c}^{\rm Ising}$ [78]; for the mean-field case of interaction that we study here, we indeed find that the equality holds. We checked the result (178) in direct simulations of the dynamics (17) by performing numerical integration of Eq. (17) with $\Delta=0$, by using the scheme detailed in B. The results, presented in Fig. 6, show the absence of hysteresis loops and abrupt jumps characteristic of a first-order transition, but rather a smooth variation of $r^{\rm st}$ with $T$ consistent with a continuous transition and in agreement with the analysis in the foregoing paragraphs.

5 Simulation results for a general point in the

$(\Delta,K,T)$-space

In the absence of analytical results, in this section, we report on simulation results on stationary state phase transitions for a general point in the $(\Delta,K,T)$-space in Fig. 1. To obtain the results, we performed numerical integration of Eq. (17) for the Lorentzian $g(\omega)$, Eq. (8). For details on the numerical scheme, see B. For given values of $K$ and $T$, and an initial configuration with oscillators at $\theta=0$, we let the system equilibrate at $\Delta=0$. Subsequently, we tune $\Delta$ adiabatically to high values and back in a cycle. Note that the tuning of $\Delta$ is performed for a fixed realization of the frequencies $\omega_{j}$’s; Referring to Eq. (17), we see that tuning of $\Delta$ is equivalent to changing the factor multiplying the frequency term in the equation of motion (17). Adiabatic tuning of $\Delta$ ensures that the system has sufficient time to attain stationarity before the value of $\Delta$ changes significantly. Figure 7 shows the behavior of the synchronization order parameter $r^{\rm st}$ for several values of $K$ at representative temperatures and for a fixed realization of the natural frequencies $\omega_{j}$’s. We have checked that up to numerical precision, the results do not change on changing the realization of the $\omega_{j}$’s. From the figure, we observe the absence of sharp jumps and hysteresis behavior characteristic of a first-order transition, but rather a continuous variation of $r^{\rm st}$ with $\Delta$ expected of a continuous phase transition. These features are consistent with the phase diagram shown in Fig. 1, and lend credence to the analysis of the model (17) presented in this work.

Referring to the phase diagram 1, we note that although the transition surface in the $(\Delta,K,T)$-space intersects with the $(\Delta,K)$-plane, the $(K,T)$-plane, and the $(\Delta,T)$-plane in straight lines, there is no a priori reason for the surface itself to be a plane. The simulation results in Fig. 7 are however consistent with the equation of a plane:

[TABLE]

the above equation may be solved for $\Delta_{c}(K,T)$, and its value is in good agreement with the transition point between the $r^{\rm st}\neq 0$ and $r^{\rm st}=0$ phase observed in Fig. 7. In this regard, let us recall that the phase diagram 1 is obtained for a Lorentzian distribution of the natural frequencies, Eq. (8). For a different unimodal distribution with a non-compact support (e.g., a Gaussian), the intersection of the transition surface with the $(K,T)$-plane, being applicable to the case when the natural frequency term is absent in the dynamics, remains a straight line, while its intersection with the $(\Delta,T)$-plane is obtained from Eq. (21). Let us choose a Gaussian $g(\omega)$ given by $g(\omega)=(1/\sqrt{2\pi})\exp(-\omega^{2}/2)$, for which the latter intersection is shown in Fig. 8, and is evidently not a straight line. The simulation results for a Gaussian $g(\omega)$ are shown in Fig. 9. Similar to the Lorentzian case, one finds the absence of sharp jumps and hysteresis behavior characteristic of a first-order transition, but rather a continuous variation of $r^{\rm st}$ with $\Delta$ that implies a continuous phase transition. From the figure, one may estimate the transition point between the $r^{\rm st}\neq 0$ and $r^{\rm st}=0$ phase, and find from the estimated values that the general transition surface in the $(\Delta,K,T)$-space for a Gaussian frequency distribution is not a plane (and thus its intersection with the $(\Delta,K)$-plane is not a straight line). Thus, we are led to conclude that having straight transition lines in the $(\Delta,T)$-plane and the $(\Delta,K)$ plane in Fig. 1 is typical to a Lorentzian frequency distribution and does not hold in general for other unimodal distributions with a non-compact support.

We observe a peculiar feature of the phase diagram 1: the phase transition line in the $(K,T)$-plane and the $(\Delta,K)$-plane suggests that the temperature $T$ in the dynamics corresponding to the former plane plays a role similar to $\Delta$ in the dynamics for the latter plane. Indeed, the transition line in the $(K,T)$-plane is $K_{c}(T,\Delta=0)=2T-1$, while the one in the $(\Delta,K)$-plane is $K_{c}(T=0,\Delta)=2\Delta-1$. This observation is somewhat counterintuitive, given that $\Delta$ corresponds to a quenched disordered (that is, a time-independent) noise in the dynamics, while $T$ signifies an annealed (that is, a time-dependent) noise, and that these two types of noise typically have very different consequences on the properties of many-body interacting systems, e.g., on surface growth dynamics [79].

Another point worth noting about the phase diagram 1 is that both the lines $K_{c}(T,\Delta=0)=2T-1$ and $K_{c}(T=0,\Delta)=2\Delta-1$ may be derived by considering the noisy Kuramoto model with an effective global coupling equal to $1+K$ (and thus only mean-field and no non-local interaction), as we demonstrate below. To this end, consider the equation of motion

[TABLE]

Then, in the $(\Delta,K)$-plane, the model reduces to the Kuramoto model with modified global coupling given by $J_{\rm eff}^{\rm Kuramoto}\equiv 1+K$. In terms of rescaled variable $\widetilde{t}\equiv t(1+K)$, the equation of motion has the form of Eq. (3) for the Kuramoto model:

[TABLE]

with

[TABLE]

Using the known results about the phase transition in the Kuramoto model, we conclude that the dynamics (182) exhibits a phase transition between a low-$\Delta_{\rm eff}$ synchronized phase and a high-$\Delta_{\rm eff}$ unsynchronized phase at a critical value given for the Lorentzian distribution (8) by $\left[\Delta_{\rm eff}\right]_{c}=\pi g(0)/2=1/2$. Combining this result with Eq. (183), we obtain the result we had set out to demonstrate, namely, $K_{c}(T=0,\Delta)=2\Delta-1$. To obtain the other result, namely, $K_{c}(T,\Delta=0)=2T-1$, we start with Eq. (181) with $\Delta=0$. Noting that in this case, the dynamics relaxes to a Boltzmann-Gibbs (BG) equilibrium state with the probability distribution of the angles $P_{\rm eq}(\{\theta_{j}\})\propto\exp[-\widetilde{\cal V}(\{\theta_{j}\})/T];~{}\widetilde{\cal V}(\{\theta_{j}\})\equiv(1+K)/(2N)\sum_{j,k=1}^{N}\left[1-\cos(\theta_{j}-\theta_{k})\right]$, one may perform an analysis similar to that in Section 4.1 to arrive at the result $T_{c}(\Delta=0,K)=(1+K)/2$, which then yields $K_{c}(T,\Delta=0)=2T-1$. The above derivation of the results $K_{c}(T=0,\Delta)=2\Delta-1$ and $K_{c}(T,\Delta=0)=2T-1$ based on only mean-field interaction points towards an apparent mean-field dominance in the stationary state of the dynamics (17) in the $(K,T)$-plane, where one has an equilibrium dynamics, and in the $(\Delta,K)$-plane for a Lorentzian distribution, where one has a non-equilibrium dynamics. Mean-field dominance in the stationary state due to an equilibrium [80] and a nonequilibrium [81] dynamics has been recently observed in a variant of the Kuramoto model that comprises oscillators interacting with one another with a strength that decays as a power-law of their separation on a $1d$ lattice [82]. The origin of the features of the phase diagram mentioned in the present and the two preceding paragraphs and the necessary and sufficient conditions for their validity are open issues left for future studies.

6 Conclusions and perspectives

In this work, we addressed the issue of spontaneous collective synchronization in many-body interacting systems within the ambit of the paradigmatic Kuramoto model of globally-coupled phase oscillators of distributed natural frequencies. The model is known to exhibit as a function of the diversity of the natural frequencies a transition between a synchronized and an unsynchronized phase. Our objective in this work was to investigate the robustness of such a behavior with respect to additional interactions. Specifically, we considered the effect of including a non-local $M$-neighbor interaction between the oscillators residing on the sites of a one-dimensional periodic lattice of $N$ sites. Here, we dealt with the case of a unimodal frequency distribution, and in particular, a Lorentzian distribution. In presence of thermal noise of strength proportional to a temperature, the resulting dynamics is effectively characterized by three parameters: the width $\Delta$ of the frequency distribution, the strength $K$ of the $M$-neighbor interaction, and the temperature $T$. In obtaining our results, we considered the simultaneous limits $M\to\infty,N\to\infty$, while keeping the interaction radius $\sigma\equiv M/N$ to satisfy $\sigma<1/2$. The latter condition allows to have distinct forms of mean-field and non-mean-field interactions.

The analysis presented in this work revealed that in the stationary state, the dynamics of our model in the $(\Delta,K,T)$-space exhibits a very rich phase diagram that involves both equilibrium and non-equilibrium phase transitions. In two contrasting limits of the dynamics, namely, (i) the limit $T\to 0$, and (ii) the limit $\Delta\to 0$, we obtained exact analytical results for the phase transitions, by borrowing tools from diverse disciplines, namely, the ones of non-linear dynamics and statistical mechanics. For case (i), when the dynamics reduces to that of a non-linear dynamical system, our exact analysis is based on the use of the so-called Ott-Antonsen (OA) ansatz to derive a reduced set of nonlinear partial differential equations for the macroscopic evolution of the system. On the other hand, in the case of (ii), the dynamics becomes that of a statistical system in contact with a heat bath and relaxing to a statistical equilibrium state, and our analytical results are derived by extending the transfer matrix approach of the nearest-neighbor Ising model to consider non-local interactions. Referring to Fig. 1, the line $K_{c}(T=0,\Delta)$ is obtained by the OA ansatz, while the line $K_{c}(T,\Delta=0)$ is obtained by the transfer matrix method. Being an ansatz, it is remarkable that the OA approach is able to correctly predict the phase transition point $K_{c}(T=0,\Delta=0)=-1$, as is checked by obtaining the same point from the analysis for the $(K,T)$-plane by employing the well-established and exact transfer matrix approach [39]. It remains an outstanding problem to obtain analytical results for the phase transition at a general point in the $(\Delta,K,T)$-space.

While the analytical results presented using the Ott-Antonsen ansatz applies to a Lorentzian distribution of the natural frequencies, qualitatively similar results are expected to hold for other unimodal frequency distributions that have a non-compact support similar to the Lorentzian distribution. Referring to the phase diagram 1, the phase boundary in the $(K,T)$-plane is obviously independent of the choice of the frequency distribution, while the boundaries in the $(\Delta,K)$- and the $(\Delta,T)$-plane and, more generally, the transition surface in the $(\Delta,K,T)$-space would depend on the specific form of the frequency distribution. Nevertheless, we expect the general features of the phase diagram to hold for other non-compact unimodal distributions, but these would certainly change if one considers distributions that are unimodal with compact support or those that are not unimodal, e.g., a bimodal distribution [83, 84, 85]. Resolution of this issue is under investigation. Another open issue is to consider finite values of $M$. In this case, the phase diagram in the $(K,T)$-plane shows a tricritical point [65], similar to the one observed for the Ising case in Fig. 5, and we may expect the tricritical point to extend to a tricritical line in the $(\Delta,K,T)$-space. Our preliminary results indeed point in that direction, and a detailed analysis will be reported elsewhere [86].

We may mention other immediate and physically relevant offshoots of our work, for example, considering in the dynamics (10) the global coupling term to also include a second harmonic $\sim\sin\left(2(\theta_{k}-\theta_{j})\right)$ [87, 88], the presence of a local potential [89], a phase-lag parameter [56, 57], a time delay in the interaction between the oscillators [41, 52], and/or considering in place of the first-order dynamics investigated in this work the case of a second-order dynamics that includes the effect of a finite inertia of the oscillators, and which is known to alter significantly the behavior of the bare Kuramoto model [32, 90, 91, 92, 93].

7 Acknowledgments

This paper is dedicated to my beloved father, the most important person in my life. SG is grateful to Alessandro Campa for fruitful discussions and suggestions, and for a careful reading of the manuscript. SG is also grateful to the Laboratoire de Physique, École Normale Supérieure de Lyon for support and warm hospitality during his stay as Professeur Invité in June 2017 when this manuscript was being finalized.

Appendix A

Proof that the dynamics (17) does not satisfy detailed balance unless $\Delta=0$

In this Appendix, we give a formal proof that the dynamics (17) does not satisfy detailed balance in the stationary state unless one has $\Delta=0$. To this end, we first consider for simplicity of discussion the case of a bimodal $g(\omega)$, and then generalize our discussion to a general $g(\omega)$. Consider a given realization of $g(\omega)$, in which there are $N_{1}$ oscillators with natural frequencies equal to $\omega_{1}$ and $N_{2}$ oscillators with frequencies equal to $\omega_{2}$, with $N_{1}+N_{2}=N$. Note that we need to consider at least two values of the natural frequencies in order to have a non-zero $\Delta$. Let us define the $N$-oscillator distribution function $f_{N}(\theta_{1},\dots,\theta_{N_{1}},\theta_{N_{1}+1},\dots,\theta_{N},t)$ as the probability density at time $t$ to observe the system around the values $\{\theta_{j}\}_{1\leq j\leq N}$, with the normalization $\int\left(\prod_{j=1}^{N}{\rm d}\theta_{j}\right)f_{N}(\{\theta_{j}\},t)=1$. The time evolution of $f_{N}$ follows the $N$-dimensional Fokker-Planck equation that may be written down from the equation of motion (17) by following standard prescription [94]:

[TABLE]

where we have defined the $N\times 1$ column vector $\Omega$ with first $N_{1}$ entries equal to $\omega_{1}$ and the following $N_{2}$ entries equal to $\omega_{2}$, and where the superscript $T$ denotes matrix transpose operation: $\Omega^{T}\equiv\left[\omega_{1}~{}\omega_{1}\dots~{}\omega_{1}~{}\omega_{2}\dots~{}\omega_{2}\right]$.

The Fokker-Planck equation (184) may be rewritten as

[TABLE]

where we have defined

[TABLE]

In Eq. (185), the drift vector $A_{j}({\bf x})$ is given by

[TABLE]

while the diffusion matrix is

[TABLE]

The dynamics described by the Fokker-Planck equation of the form (185) satisfies detailed balance if and only if the following conditions are satisfied [94]:

[TABLE]

where $f_{N}^{s}({\bf x})$ is the stationary solution of Eq. (185). Here, $\epsilon_{j}\equiv\pm 1$ denotes the parity of $x_{j}$’s with respect to time reversal $t\to-t$: Under time reversal, the $x_{j}$’s transform as $x_{j}\rightarrow\epsilon_{j}x_{j}$, with $\epsilon_{j}=-1$ or $+1$ depending on whether $x_{j}$ is odd or even under time reversal. In our case, $\theta_{j}$’s are even variables, so that we consider $\epsilon_{j}=+1~{}\forall~{}j$ in the following discussion.

Using Eq. (188), we see that the condition (189) is trivially satisfied for our model. To check the other condition, namely, Eq. (190), let us formally solve the equation to obtain $f_{N}^{s}({\bf x})$, and check whether the solution solves Eq. (185) in the stationary state. Using $\epsilon_{j}=+1~{}\forall~{}j$, we see that the condition reduces to

[TABLE]

solving which yields

[TABLE]

Substituting Eq. (192) into Eq. (185), and requiring that $f_{N}^{s}({\bf x})$ given by the former is a stationary solution of the latter, we obtain the condition that $\Delta$ has to be equal to zero. We thus conclude from the foregoing discussions that the dynamics (17) does not satisfy detailed balance unless $\Delta=0$. The foregoing discussions for a bimodal $g(\omega)$, establishing the lack of detailed balance for $\Delta\neq 0$, obviously extend to any choice of $g(\omega)$. For $\Delta=0$, we get the stationary solution as

[TABLE]

where the function ${\cal V}(\{\theta_{j}\})$ is given by Eq. (20). We thus obtain for $\Delta=0$ the equilibrium distribution for the angles as

[TABLE]

Appendix B

Numerical scheme to integrate the equation of motion (17)

Here we give details of the numerical scheme to integrate the equation of motion (17). To this end, rewriting the mean-field term in Eq. (17) in terms of the quantities $(r_{x},r_{y})$ defined in the paragraph following Eq. (4), we obtain the net torque acting on the $j$-th oscillator as

[TABLE]

To simulate the dynamics (17) over a time interval $[0:{\cal T}]$ and for given values of $K,M$ and $T$, we first choose a time step size $\Delta t\ll 1$, and set $t_{n}=n\Delta t$ as the $n$-th time step of the dynamics, where $n=0,1,2,\ldots,{\cal N}_{t}$, and ${\cal N}_{t}={\cal T}/\Delta t$. One step of the numerical scheme then involves the following update of the $\theta_{j}$’s for $j=1,2,\ldots,N$:

[TABLE]

with $X_{n}$ a Gaussian-distributed random number with zero mean and variance given by

[TABLE]

and

[TABLE]

The above numerical scheme neglects terms of order higher than $\Delta t$, and may be derived by following the general techniques discussed in Ref. [95]. For the numerical results reported in this work, we have chosen a fixed time step $\Delta t=0.01$. We have checked that up to numerical precision, the obtained numerical results do not depend on $\Delta t$, so long as the latter is small.