Hysteretic behavior of spatially coupled phase-oscillators

TL;DR

This paper studies how spatially coupled phase oscillators with inertia exhibit hysteretic synchronization behavior, influenced by local interactions and kernel radii, with implications for biological systems like firefly swarms.

Contribution

It introduces a generalized Kuramoto model with inertia and spatially dependent coupling, analyzing how local interactions affect hysteresis in synchronization.

Findings

Locally locked clusters form at low coupling strengths.

Hysteresis disappears with small kernel radii.

Kernel radius relative to maximal oscillator distance determines network behavior.

Abstract

Motivated by phenomena related to biological systems such as the synchronously flashing swarms of fireflies, we investigate a network of phase oscillators evolving under the generalized Kuramoto model with inertia. A distance-dependent, spatial coupling between the oscillators is considered. Zeroth and first order kernel functions with finite kernel radii were chosen to investigate the effect of local interactions. The hysteretic dynamics of the synchronization depending on the coupling parameter was analyzed for different kernel radii. Numerical investigations demonstrate that (1) locally locked clusters develop for small coupling strength values, (2) the hysteretic behavior vanishes for small kernel radii, (3) the ratio of the kernel radius and the maximal distance between the oscillators characterizes the behavior of the network.

Figures

Equations

$$R{e}^{\mathrm{i}\mathrm{\Phi }}=\sum _{j=1}^N{e}^{\mathrm{i}{\theta }_{\mathrm{j}}},$$

$$\frac{K}{N}\sum _{j=1}^N{A}_{ij}\sin ({\theta }_j-{\theta }_i),$$

$${\theta }_i=\theta ({\text{r}}_i).$$

$$K\sum _{j=1}^{n_i}\sin ({\theta }_j-{\theta }_i){\hat{W}}_{ij},$$

$${\hat{W}}_{ij}=\frac{W_{ij}}{{\sum }_{j=1}^{n_i}{W}_{ij}},$$

$$\frac{{\ddot{\theta }}_i}{m}+{\dot{\theta }}_i={\Omega }_i+K\sum _{j=1}^{n_i}\sin ({\theta }_j-{\theta }_i){\hat{W}}_{ij}i=1,..,N,$$

$${\Omega }_i=-{\Omega }_M+\frac{2(i-1){\Omega }_M}{N-1},$$

$$\begin{array}{rl}& {\theta }_i(t=0)=2\pi \frac{i}{N},\\ & {\dot{\theta }}_i(t=0)={\Omega }_i,\end{array}$$

$$\begin{array}{rl}& {\theta }_i(t=0)=0,\\ & {\dot{\theta }}_i(t=0)=0.\end{array}$$

$$W_{ij}^0=\{\begin{array}{ll}& 1\mathrm{if}{\mathrm{d}}_{\mathrm{ij}}\le \mathrm{\Delta },\\ & 0\mathrm{if}{\mathrm{d}}_{\mathrm{ij}}>\mathrm{\Delta },\end{array}$$

$$W_{ij}^1=\{\begin{array}{ll}& 1-\frac{d_{ij}}{\Delta },\mathrm{if}{\mathrm{d}}_{\mathrm{ij}}\le \mathrm{\Delta },\\ & 0,\mathrm{if}{\mathrm{d}}_{\mathrm{ij}}>\mathrm{\Delta },\end{array}$$

$$q:=\frac{\Delta }{d_{ij}^{max}}$$

Full text

1. Introduction

Synchronization is a collective behavior observed in many fields. It is a result of the interaction between oscillators capable of adjusting their rhythms/natural frequencies. To model synchronization of coupled phase oscillators the Kuramoto model proposed in 1975 [1] is often used, for example, to investigate the collective behavior of lasers [2], neurons [3, 4, 5], social groups [6] and even crickets [7]. This model was also used to describe the interesting phenomena related to Pteroptyx malaccae, a species of firefly capable of synchronous firing with almost no phase lag [8]. This can be attributed to this insect’s ability to alter its flashing frequency in response to external stimulus [9]. Motivated by behavior of fireflies, the original first order Kuramoto model was later extended with an inertial term by Tanaka et al. [10] which allows for the adaptation of the flashing frequency of one firefly. They showed that in a fully coupled system, the degree of synchrony depends on the coupling strength between the oscillators in a hysteretic manner. The critical coupling strength necessary for the system to transition from the incoherent state to the coherent state is larger than the critical coupling strength resulting in the breaking of synchrony. According to Tanaka et al., the coupling strength of a firefly depends on the ratio of the brightness of the firing and that of the environment. This can lead to the swarm dynamics to exhibit hysteretic behavior. As a result, if the brightness of the background is too large, the fireflies are unable to synchronize. However, the model suggests that by increasing the brightness of the background of a synchronously flashing swarm, they can also maintain synchrony in a much brighter environment.

The generalized model is used to describe the synchronization of Josephson junctions [11] and power-grids . The hysteretic behavior in the model is affected by many factors. In the original Kuramoto model, the coupling between the oscillators is usually considered to be undirected: they are either connected or not, which is a valid model for some applications such as power grids. Direction of current research include considering the topology of the network [18], the heterogeneity of the network connections [19], the effect of dilution [15, 20] and assortativity [21]. Heterogeneity can be also considered by assuming time-delay or frequency-weighted coupling [22, 23]. Spatial distribution of the oscillators can also result in the delay or weakening of the signal. As it was reported in [24], P. malaccae has a 3 feet range of vision which advocates the assumption of local interactions. The distance-dependency of the coupling strength was taken into account in the original Kuramoto model by using kernel functions [25, 26] .

In this work, we numerically analyzed the hysteretic behavior of spatially coupled phase oscillators described by the generalized Kuramoto-model with inertia. The paper is structured as follows. In the next Section, we briefly describe the model and our approach to include spatial coupling. We describe the simulations in Section 3, which is followed by the report of our main results in Section 4. Finally, in Section 5, we summarize our observations and discuss some possible applications.

2. The phase-oscillator model

2.1. Coupled oscillators with inertia

Also known as the damped driven pendulum model, the system of coupled oscillators with the inertial term extension has been introduced and numerically investigated by [10]. The system with $N$ mutually coupled oscillators reads as follows:

[TABLE]

where $\theta_{i}(t)$ and $\Omega_{i}$ are the phase and natural frequency of the $i$th oscillator respectively, $K$ is the coupling strength parameter expressing how quickly an oscillator can adapt to (the resultant of the) external stimuli, and $m$ is the inertial constant.

The global synchronization of the system of oscillators can be characterized by the complex order parameter

[TABLE]

where magnitude of the complex parameter $R$ – hereinafter referred to as order parameter – describes the level of synchronization, $\Phi=1/N\sum_{j}{\theta_{j}}$ is the arithmetic mean of the $\theta$’s. Using (2), $R=1$ corresponds to the complete synchronization ($\Phi=\theta_{j}$) and $R=0$ corresponds to an incoherent state of the oscillators.

Varying the coupling strength $K$ in the system leads to a hysteretic behavior in the model. For small $K$ values, the phases of the oscillators are incoherently distributed and $R\approx 0$. Increasing $K$ leads to the appearance of phase-locked oscillators and eventually the system reaches a coherent state. However, decreasing $K$ from a completely synchronized state results in a higher level of synchronization for smaller $K$ values.

2.2. Coupling with spatial collocation

The second term on the right-hand side of Eq. (1) represents the Kuramoto-synchronization term, which is often modified with the adjacency matrix $A_{ij}$ so that

[TABLE]

where the value of $A_{ij}$ is $1$ if and only if the $i$th and $j$th oscillators are coupled, and [math] otherwise [28]. In case $A_{ij}\neq 1\,\forall\,i,j$, the system is diluted. It is known, that some diluted systems exhibit a hysteretic behavior similar to the fully coupled case. The hysteretic dynamics of randomly diluted systems has recently been investigated [15]. However, even in case of the modified model Eq. (3), former investigations neglect the effect of spatial distribution on the hysteretic phase-synchronization dynamics.

As a generalization of the mean in Eqs. (1) and (3), we propose a spatial averaging technique for the computation of pairwise coupling strength as a function of internodal distances and local neighborhoods. Systems with spatially distributed phase-oscillators may require special treatment depending on the spatial distribution. Recently [25] and [26] investigated the phase-oscillator model without inertial term using wavelet-like and bell-shaped kernel functions, respectively. Both works consider the kernel as a function of the internodal distances for the weighting of the coupling strength.

We consider a set of $N$ spatially distributed nodes on the plane with positions $\textbf{r}_{i}$ with internodal Euclidean distances $d_{ij}=|\textbf{r}_{ij}|=|\textbf{r}_{i}-\textbf{r}_{j}|$. Using the nodal positions, we define the phase assigned to each node as

[TABLE]

In order for the pairwise coupling strength to be scaled as a function of the distances, we define a spherically symmetric kernel-function $W_{ij}=W(d_{ij},\Delta)$ with (finite or infinite) smoothing radius $\Delta$ and construct the weighted average for the Kuramoto phase-synchronizer term as

[TABLE]

where $n_{i}$ is the number of neighbors within the kernel radius $\Delta$ around the $i$th node. The $j$th oscillator is a neighbor of the $j$th oscillator if it is in the $\Delta$ neighborhood of the $i$th oscillator, i.e. if $d_{ij}\leq\Delta$. The normalization (6) is also known as Shepard’s correction of the weighted summation [29]. As a result, instead of Eq. (1) we consider the following equation

[TABLE]

where $K$ is the coupling parameter. Note, that in Eq. (7), the coupling parameter of the $i$th and $j$th oscillator is weighted with the kernel function, therefore the actual coupling strength between the oscillators is varying.

The model described by Eq. (7) is a generalization of Eq. (1). Choosing $W_{ij}=1$, in the $\Delta\rightarrow\infty$ limit we have $\sum_{j=1}^{n_{i}}{W_{ij}}=N$ and Eq. (7) reproduces the conventional Kuramoto model. Consequently, we keep the definition of the order parameter $R$ (Eq. (2)) to describe the level of synchronization of the spatially coupled system.

3. Simulations

We implemented the model in Nauticle, the general purpose particle-based simulation tool [30]. Facilitating the implementation and application of meshless numerical methods, Nauticle provides the sufficient flexibility in building arbitrary mathematical models with free-form definition of governing equations. Being a meshless open source simulation package, it is capable to solve large system of coupled ordinary differential equations, hence, the numerical model discussed in the present work can be configured easily in terms of both the equations and the geometrical layout.

Although the proposed model is suitable for arbitrary spatial distribution of the oscillators, as an initial study, unless stated otherwise, we investigate $N=n^{2}$ coupled phase oscillators placed on a two-dimensional equidistant square-grid with grid cell size $\Delta x$ and of edge length $L=(n-1)\Delta x$. Consequently, the position vector of an oscillator is $\mathbf{r}^{jk}=((j-1)\Delta x,(k-1)\Delta x)$, where $j,k=1,...,n$. Each oscillator is assigned with a scalar index parameter $i$, such that $i=(k-1)n+j$.

Following [10], we chose an evenly spaced natural frequency distribution on the interval of $[-\Omega_{M},\Omega_{M}]$ such that

[TABLE]

where $\Omega_{M}$ is a constant. We considered two types of initial conditions of Eq. (7), the uniformly diffused (IC 1) taking

[TABLE]

and the perfectly synchronized (IC 2), where

[TABLE]

In the present work we performed the simulations with two different, finite width spatial kernel-functions with a kernel radius $\Delta$. The zeroth order kernel-function is constant in the neighborhood of the $i$th oscillator

[TABLE]

and the first order kernel-function depends linearly on the distance between the $i$th and $j$th oscillators

[TABLE]

consequently it is maximal at the center and decreases toward the boundary of the neighborhood.

During the solution we applied the classic fourth order Runge-Kutta scheme for numerical integration of the system defined by Eq. (7). We fixed the time step size $\Delta t=0.1$ in all simulations. On the one hand, to let the value of $R$ reach a developed state, we run all cases for $500$ in simulation time and in order to eliminate the oscillations of $R$ in time, we computed the temporal average of the order parameter $\bar{R}$ of the last $2000$ steps. Also, an extended simulation of $10000$ steps had been performed to check if the results sufficiently converged but no significant changes were observed.

4. Results

Since the hysteretic behavior of the fully coupled system of oscillators with inertia was studied in [10], it is known that the temporal development of the order parameter $R$ may strongly depend on the initial conditions. In this section we present our numerical investigation in terms of the effects of the finite width spatial covering on the synchronization of the spatially distributed oscillators. We kept $\Omega_{M}=5$ and the inertia $m=0.85$ constant in all simulation cases and varied the kernel radii and the spatial distribution.

As a measure of the distance-based dilution we introduce a dimensionless parameter $q$, the relative kernel radius

[TABLE]

expressing the ratio between the kernel radius and $d_{ij}^{max}$ the maximal Euclidean distance between two oscillators of the system. In case of the zeroth order kernel-function $W_{ij}^{0}$ and $q\geq 1$, i.e. $\Delta>L$ the numerical setup reproduces the fully coupled second-order model (Eq. (1)).

4.1. The effect of spatial coupling

In order to present the local synchronization effects, we solved the equation for IC 1 and IC 2 for three different values of $q$ and both zeroth and first order kernel-functions (Eqs. (11) and (12), respectively).

If $q$ is large enough (Figs. 1A-D), the spatially coupled model exhibits similar hysteretic behavior as the fully-coupled model. However, spatial coupling leads to the appearance of a of the system (Fig. 1). Starting from an incoherent state (IC 1) and a small coupling parameter $K$, the average order parameter $\bar{R}$ stays near zero. As $K$ is increased, at a critical point $K_{C1}$ the synchronization begins, i.e. $\bar{R}$ start to increase. Further increasing $K$ leads to $\bar{R}$ reaching a maximum and it stays constant for a range of $K$. However it decreases before transitioning to the globally coherent state by starting to increase again at a critical value of the coupling parameter $K_{C2}$. Therefore in case of IC 1, the globally coherent state is reached through two critical points $K_{C1}$ and $K_{C2}$, the beginning of the local and global synchronization, respectively. In case of IC 2, i.e. starting from a coherent state and decreasing the value of $K$, the system either jumps into an incoherent state (Fig. 1A) or the transitional state corresponding to locally locked clusters (Fig. 1B-D) at a bifurcation point $K_{C3}$. For $K<K_{C3}$ there is almost no difference between IC 1 and IC 2. In particular, if the maximal coupling strength is small enough, the local behavior dominates.

Under a critical value of the relative kernel radius $q$, the local and global behavior are separated and the system exhibits no hysteresis (Figs. 1E-F), i.e. there is almost no difference between IC 1 and IC 2. Moreover, there is no constant part of the $\bar{R}-K$ diagram. Decreasing $q$ increases $K_{C2}$ and $K_{C3}$ and decreases $K_{C1}$. For $q\approx 0.14$ (Figs. 1E-F) the synchronization starts near $K\approx 0$.

There is no qualitative difference between the results calculated with the applied kernel-functions. However, for constant $\Delta$, $W_{ij}^{1}$ leads to a higher level of dilution compared to $W_{ij}^{0}$. As a result, having the same kernel radii, the range of the partial synchronization (i.e. $K_{C1}\leq K\geq K_{C2}$) is always larger when using a first order than a zeroth order kernel-function. Interestingly, the average order parameter always decreased before the global synchronization began at $K_{C2}$ in all of our simulations.

4.2. States of the system

We analyzed the states of the spatially coupled model depending on the coupling parameter $K$ in case $q$ is large enough for the system to exhibit the hysteretic behavior. Analysis of the solutions of Eq. 7 in time was carried out by taking uniformly diffused initial conditions (IC 1) and first order kernel-function $W_{ij}^{1}$ with fixed kernel radius $\Delta=0.6$ ($q\approx 0.46$) and $\Delta x=0.05$. By examining the time evolution of the order parameter $R$ (Fig. 2A) and the phases of the individual oscillators (Fig. 2B), three states of the system in the $\bar{R}-K$ diagram (Fig. 2C) was found. 2B is only a visualization of the data and the distances measured on the circle are independent of the strength of the interactions. .

If $0<K<K_{C1}$, the system is in the incoherent state. Since we keep the natural frequency distribution defined by Eq. (8) and it is a dominant part of the equation, there are two in the initial configuration ($\Omega_{i}<0$ and $\Omega_{i}\geq 0$). The order parameter $R$ varies in time in a random-like manner, but it remains close to zero (Fig. 2A, $K=2$) and the phases of the oscillators are heterogeneous (Fig. 2B $K=2$). The fluctuations of $R$ in time can be attributed to the finite number of oscillators in the system. Increasing $K$, the average order parameter $\bar{R}$ stays near zero until reaching a $K_{C1}$ value. Above $K_{C1}$, the degree of synchronization gradually increases.

Increasing $K$ from $K_{C1}$, Note, that the overlap of the clusters neither result in their interaction nor global synchronization. Further increasing $K$, the system rapidly jumps into a globally locked state in a bifurcation point, at a $K_{C2}$ value. At this point, the clusters merge as one of them take over the domination

If $K_{C2}<K$, there is only one cluster, the system is in a globally locked state (Fig. 2B, $K=20$) and $R$ is constant in time (Fig. 2A, $K=20$). As $K$ is further increased, the synchronization and $\bar{R}$ gradually increases (Fig. 2A-B, $K=200$) and the system seemingly reaches the globally coherent state at a large value of $K$. In contrast to the fully coupled model, where $R\approx 1$ is reached at moderate $K$, here the results imply, that $\bar{R}\rightarrow 1$ asymptotically. Accordingly, spatial coupling can hinder full synchronization for large system sizes in case the coupling strength of the individual oscillators is limited.

We also checked the solutions of the model for IC 2. For $0\leq K\geq K_{C3}$, the system has the same states as for IC 1. Above $K_{C3}$, there is one globally locked cluster with constant $R$ in time.

4.3. Effect of the relative kernel radius

In Fig. 1, $\Delta$ and $q$ were varied simultaneously in the different simulations. The diagrams show varying values of $K_{C1},K_{C2},K_{C3}$ and different $max(R_{C})$ in the locally locked state. To examine the effect of the parameters separately, we calculated $\bar{R}$ versus $K$ for both IC 1 and IC 2 keeping either the number of oscillators in a neighborhood of a general point constant (i.e. keeping $\Delta$ and $\Delta x$ constant) or the relative kernel radius $q$ fixed. Since, there is no qualitative difference between the first and zeroth order kernel-functions, we applied $W_{ij}^{1}$ in all cases.

As we can see in Fig. 7, although the number of oscillators in a general in-domain neighborhood is the same in all cases, both the local and global parts of the diagram are affected. In case of a smaller relative kernel radius, the local synchronization starts at smaller values of $K$ and a larger $K$ is necessary for the onset of the global synchronization. Furthermore, the maximal order $max(R_{C})$ of the system in the locally synchronized state is higher for $q_{2}$. Smaller relative kernel radius $q$ results in a larger range of $K$ corresponding to the locally synchronized states (i.e. $|K_{C1}-K_{C2}|$ increases), but it also requires a larger value of $K_{C3}$ for the global synchronization.

Fig. 8 shows, that there is only a small difference between the locations of $K_{C1}$ and $K_{C3}$, but the transitional state is not affected. The onset of the local synchronization requires larger $K$ as the system size is increased, but the global synchronization can be maintained for smaller values of $K$. Similar dependencies on the system size were reported in [10, 15] for the fully coupled and diluted systems, respectively. These results suggest, that the $\bar{R}-K$ diagram is characterized by $q$ and the local behavior is not affected by the system size. The maximal order $max({R_{C}})$ of the system in the locally synchronized state was also independent of the system size in our simulations.

Finally, we computed the maximal value of the order parameter in the locally locked state $\mathrm{max}(R_{C})$ for $K_{C1}<K<K_{C2}$ depending on $q$, keeping $L=1,\Delta x=0.05,N=400$ fixed and applying first order kernel-function. According to Fig. 9, the $max(R_{C})(q)$ function has a maximum. Consequently, there is an optimal value of the relative kernel radius in terms of the local synchronization. If $q$ is too large, the local behavior dominates and if the kernel radius is too small compared to the domain size, the local synchronization is hindered.

5. Conclusion

In summary, we carried out a numerical study to examine the effect of the spatial coupling of phase oscillators in the Kuramoto model with inertia. A finite size system was considered with compact zeroth and first order kernel functions to incorporate distance-dependent coupling strength between the oscillators. A dimensionless parameter, the relative kernel radius $q$ was introduced, expressing the ratio between the kernel radius and the maximal distance between the $i$th and $j$th oscillators. We examined the order of the global synchronization depending on the coupling parameter for different $q$ values. In the case of large kernel radii, the model gives back the original, fully coupled model.

We implemented the model in the Nauticle general purpose particle-based simulation tool [30] and examined the hysteretic behavior of the spatially coupled model. We fixed the domain size and investigated the behavior for different values of $q$. For a range of $q$, the system exhibits a hysteretic behavior depending on the initial conditions and the coupling parameter, but for a range of $K$ values locally locked clusters develop and dominate. In case of the presence of the local clusters, the order parameter $R$ is oscillates in time. We pointed out, that this oscillations can be attributed to locally locked clusters that develop and break up through some transitional states. If $q$ is small enough, the system exhibits no hysteresis, which agree with previous works showing that reducing the number of links between the oscillators lead to the vanish of the hysteretic behavior [15]. Analysis of the parameters in the problem showed, that the local behavior is characterized mainly by the relative kernel radius $q$, which has an optimal value leading to a maximal value of the synchronization in the locally locked state.

Our results suggest, that large systems of oscillators having a distance limited vision, such as the fireflies can develop locally synchronized clusters without transitioning into global synchronization. The proposed framework can be used to model more complex spatially coupled systems. Some bacteria such as Myxococcus xanthus [31] or systems of microgears [32] exhibit not only social interactions but also complex mechanical behavior. By the implementation of the model into Nauticle it is possible to consider a mathematical coupling of the Kuramoto model with other equations, such as mechanical models to examine complex systems of spatially moving oscillators.

Acknowledgement

The research reported in this paper was supported by the Higher Education Excellence Program of the Ministry of Human Capacities in the frame of Water Science & Disaster Prevention research area of Budapest University of Technology and Economics (BME FIKP-VÍZ). The research reported in this paper has been supported by the National Research, Development and Innovation Fund (TUDFO/51757/2019-ITM, Thematic Excellence Program). The research reported in this paper has been supported by the National Research, Development and Innovation Fund (TUDFO/51757/2019-ITM, Thematic Excellence Program).