Solitary States and Partial Synchrony in Oscillatory Ensembles with Attractive and Repulsive Interactions

TL;DR

This paper investigates how networks of coupled oscillators transition between various synchronized states, including solitary and partially synchronous states, as the balance shifts from attraction to repulsion among elements.

Contribution

It provides both numerical and analytical insights into the transitions between synchronous, solitary, and quasiperiodic states in oscillatory ensembles with mixed interactions.

Findings

Identification of solitary states in oscillator networks

Analysis of transitions driven by increasing repulsion

Characterization of partially synchronous quasiperiodic dynamics

Abstract

We numerically and analytically analyze transitions between different synchronous states in a network of globally coupled phase oscillators with attractive and repulsive interactions. The elements within the attractive or repulsive group are identical, but natural frequencies of the groups differ. In addition to a synchronous two-cluster state, the system exhibits a solitary state, when a single oscillator leaves the cluster of repulsive elements, as well as partially synchronous quasiperiodic dynamics. We demonstrate how the transitions between these states occur when the repulsion starts to prevail over attraction.

Equations

$${\dot{\theta }}_j^{\sigma }={\omega }_{\sigma }+\sum _{{\sigma }^{\prime }=1}^M\frac{K_{\sigma {\sigma }^{\prime }}}{N}\sum _{k=1}^{N_{{\sigma }^{\prime }}}\sin ({\theta }_k^{{\sigma }^{\prime }}-{\theta }_j^{\sigma }+{\alpha }_{\sigma {\sigma }^{\prime }}),$$

$${\dot{\phi }}_j$$

$${\dot{\psi }}_j$$

$$H=h{e}^{i\Phi }=\frac{N_a}{N}{e}^{i{\alpha }_a}{Z}_a-\frac{N_r}{N}(1+\epsilon )e^{i{\alpha }_r}{Z}_r,$$

$${\dot{\phi }}_j$$

$${\dot{\psi }}_j$$

$$H=h{e}^{i\Phi }=\frac{1}{2}\left[Z_a-(1+\epsilon )Z_r\right].$$

$$h\cos {\Phi }_0$$

$$h\sin {\Phi }_0$$

$$\omega =h[\sin {\Phi }_0-\sin ({\Phi }_0-{\psi }_0)].$$

$$\sin {\psi }_0=-\frac{2\omega }{\epsilon }.$$

$$|\omega |\le -\epsilon /2.$$

$$\nu =\frac{1+\epsilon }{\epsilon }\omega .$$

$${\dot{\psi }}_{\pm }=\omega +h\sin ({\Phi }_0-{\psi }_0\mp \alpha ).$$

$$\dot{\alpha }=-\alpha h\cos ({\Phi }_0-{\psi }_0).$$

$$\cos {\psi }_0-(1+\epsilon )>0.$$

$$\omega =\pm \sqrt{-\frac{{\epsilon }^3}{2}-\frac{{\epsilon }^4}{4}}.$$

$$\cos {\psi }_0<\frac{1}{1+\epsilon }.$$

$$\dot{\phi }$$

$$\dot{\psi }$$

$$\dot{\delta }=\omega +\frac{\epsilon }{2}\sin \delta .$$

$$h{e}^{i{\Phi }_0}$$

$$\nu$$

$$\text{mbox}Im[e^{i({\Phi }_0-{\delta }_1)}]=\frac{1}{1+\epsilon }\text{mbox}Im(e^{i{\Phi }_0}).$$

$$h=\frac{1}{2}\sqrt{1-{\text{mbox}Im(e^{i{\Phi }_0})}^2}-\frac{N_r-2}{2N_r}\sqrt{{(1+\epsilon )}^2-{\text{mbox}Im(e^{i{\Phi }_0})}^2}.$$

$$0=x\sqrt{1-x^2}-\frac{N_r-2}{N_r}x\sqrt{{(1+\epsilon )}^2-x^2}-2\frac{1+\epsilon }{\epsilon }\omega .$$

$$f(x,N_r,\epsilon )=x\sqrt{1-x^2}-\frac{N_r-2}{N_r}x\sqrt{{(1+\epsilon )}^2-x^2}.$$

$$(1-2x^2)N_r\sqrt{{(1+\epsilon )}^2-x^2}=[{(1+\epsilon )}^2-2x^2](N_r-2)\sqrt{1-x^2}.$$

$${\dot{\delta }}_2=\frac{1}{2}[\sin {\delta }_2((1+\epsilon )-\cos {\delta }_1)-\cos {\delta }_2\sin {\delta }_1+\sin {\delta }_1].$$

$${\dot{\delta }}_2=A[\sin ({\delta }_2-{\delta }_2^{\ast })+\sin {\delta }_2^{\ast }],$$

$$\sin {\delta }_2^{\ast }=\frac{1}{1+\epsilon }\sin {\Phi }_0,$$

$${\rho }_r={\rho }_a/(1+\epsilon ).$$

Full text

Networks of coupled oscillators are a popular model for many engineered or natural systems. The main effect – emergence of a collective mode via synchronization – is now well-understood and therefore focus of research shifted recently to analysis of different complex states. These states include chimeras, when a population of identical units splits into a synchronous and asynchronous part, quasiperiodic partially synchronous states, characterized by the difference of frequencies of individual units and of the collective mode, and clusters and heteroclinic cycles, to name just a few. Of particular interest are ensembles where some elements have only attractive connections while others have only repulsive ones. This model is motivated by studies of neuronal networks that are built from excitatory and inhibitory neurons. In this paper we analyze how the state of such a setup changes with the interplay of attraction and repulsion. We demonstrate that if the frequency mismatch between attractive and repulsive units is smaller than some critical value then desynchronization occurs via appearance of the solitary state. With the further increase of repulsion the system undergoes a transition to quasiperiodic partial synchrony. In the latter state the attractive units remain synchronized, while the repulsive group settles between synchrony and asynchrony so that the mean fields of both groups remain locked, but the frequency of the repulsive elements is larger than that of their mean field. For a large frequency mismatch of attractive and repulsive groups desynchronization immediately leads to partial synchrony.

I Introduction

Investigation of coordinated dynamics of many interactive oscillatory elements is relevant for the understanding of various phenomena from different branches of science. Probably, the most important and also mostly studied effect is the emergence of a collective mode, observed in populations of flashing fireflies Kaempfer (1906), groups of pedestrians on footbridges Strogatz et al. (2005) or metronomes placed on a common support Martens et al. (2013), electronic circuits Watanabe and Strogatz (1994), populations of cells Richard et al. (1996), synthetic genetic oscillators Prindle et al. (2011), etc. Besides of collective synchrony, oscillatory networks exhibit many other interesting dynamical states like clusters and heteroclinic switching Hansel, Mato, and Meunier (1993), chimeras Kuramoto and Battogtokh (2002), collective chaos Hakim and Rappel (1992), traveling waves Hooper and Grimshaw (1988), quasiperiodic partial synchrony Van Vreeswijk (1996); Rosenblum and Pikovsky (2007); Pikovsky and Rosenblum (2009); Clusella, Politi, and Rosenblum (2016), solitary states Maistrenko, Penkovsky, and Rosenblum (2014), and so on. Analysis of such states and transitions between them is in the focus of current research.

Some of mentioned effects can be studied within the framework of the famous Kuramoto model Kuramoto (1984) and of its immediate extension, the Kuramoto-Sakaguchi model Sakaguchi and Kuramoto (1986), that treat phase oscillators with the sine-coupling. Though this is a rather simplistic description of real-world oscillators, these models became extremely popular due to the possibility of analytical treatment Acebrón et al. (2005); Pikovsky and Rosenblum (2015). For example, they allow for theoretical description of synchronization transitions (that, in dependence on the distribution of oscillatory frequencies, can be alike second- or first-order Pazó (2005) phase transitions). Due to their specific mathematical properties, sine-coupled phase oscillators also often admit a low-dimensional description via the Watanabe-Strogatz (WS) Watanabe and Strogatz (1993, 1994) and Ott-Antonsen (OA) Ott and Antonsen (2008, 2009) theories. All this explains why the Kuramoto-Sakaguchi model became a paradigmatic one, with applications ranging from explanation of social effects Kaempfer (1906); Strogatz et al. (2005) to neuroscience Breakspear, Heitmann, and Daffertshofer (2010).

In most variants of the Kuramoto-Sakaguchi model researchers treat networks with attractive interactions and the existing literature extensively covers this case Montbrió, Kurths, and Blasius (2004); Abrams et al. (2008); Barreto et al. (2008). Networks of repulsive elements attract much less attention, although they show interesting effects Van Vreeswijk, Abbott, and Bard Ermentrout (1994); Tsimring et al. (2005); Pimenova et al. (2016). Not much attention is also paid to mixed networks Hong and Strogatz (2011a, b); Anderson et al. (2012); Iatsenko et al. (2013); Vlasov, Macau, and Pikovsky (2014); Qiu et al. (2016), consisting of both attractive and repulsive elements, though systems of this type are common in neuroscience, because real neurons interact via excitatory and inhibitory connections Wilson and Cowan (1972); Van Vreeswijk and Sompolinsky (1996); Peyrache et al. (2012); Dehghani et al. (2016).

In this paper we concentrate on emergence of solitary state and quasiperiodic partial synchrony in networks with attractive and repulsive connections. The solitary state, when a single repulsive unit leaves the synchronous cluster, was for the first time found and analyzed in Ref. Maistrenko, Penkovsky, and Rosenblum (2014) and later in Refs. Brezetskyi, Dudkowski, and Kapitaniak (2015); Jaros, Maistrenko, and Kapitaniak (2015); Chouzouris et al. (2018); Chen, Engelbrecht, and Mirollo (2019); Majhi, Kapitaniak, and Ghosh (2019). A generalized solitary state, where several oscillators exhibit dynamics different from that of the synchronous cluster received attention in Refs. Kapitaniak et al. (2014); Hizanidis et al. (2016); Rybalova et al. (2017); Semenova et al. (2017); Jaros et al. (2018); Semenova, Vadivasova, and Anishchenko (2018); Shepelev, Strelkova, and Anishchenko (2018); Rybalova, Strelkova, and Anishchenko (2018); Mikhaylenko et al. (2019); Sathiyadevi et al. (2019). This state appears at the border between synchrony and asynchrony, as soon as repulsion starts to prevail over attraction. Our setup is an extension of the finite-size two-group Kuramoto model treated in Ref. Maistrenko, Penkovsky, and Rosenblum (2014), where all oscillators were identical.

We demonstrate that for small frequency mismatches between the groups and a weak repulsion, there appears a small region, where the attractive units build a synchronous cluster, while the repulsive oscillators exhibit quasiperiodic partially synchronous dynamics. Slightly stronger repulsion leads to the solitary state, which is replaced by quasiperiodic dynamics again for bigger repulsion. For large mismatches in the frequency the solitary state is not observed, but only quasiperiodic dynamics.

II The Model

A popular version of the standard Kuramoto-Sakaguchi model is a system of $M$ interacting groups of identical units, described by the following equations:

[TABLE]

where $\theta^{\sigma}_{j}$ is the phase of the $i$th oscillator in the group $\sigma$ and $\sigma=1,\ldots,M$. Here $\omega_{\sigma}$ and $N_{\sigma}$ are the natural frequency and the number of oscillators in the group $\sigma$, $N=\sum_{\sigma}N_{\sigma}$, and $K_{\sigma\sigma^{\prime}}$ and $\alpha_{\sigma\sigma^{\prime}}$ are respectively the strength of the coupling and the phase shift characterizing interaction between groups $\sigma$ and $\sigma^{\prime}$.

In the following we analyze a two-group Kuramoto-Sakaguchi model wherein the coupling coefficients and the phase shift parameters depend on the acting group only, i.e. $K_{\sigma\sigma^{\prime}}=K_{\sigma^{\prime}}$ and $\alpha_{\sigma\sigma^{\prime}}=\alpha_{\sigma^{\prime}}$. We concentrate on a particular case, motivated by neuroscience applications, when the coupling within the first group is attractive while in the second group it is repulsive. We denote phases of the units in these groups by $\varphi$ and $\psi$, respectively. By re-scaling the time and performing a transformation to a reference frame co-rotating with the frequency of the attractive group, we write the model as

[TABLE]

where subscripts $a$ and $r$ stand for “attractive” and “repulsive”, respectively. Quantification of coupling has been reduced to a single parameter $K_{r}/K_{a}=-(1+\varepsilon)$, with $\varepsilon$ being the excess of repulsive coupling. An $\varepsilon<-1$ indicates that interaction within both groups is attractive and, trivially, the whole system synchronizes. For $\varepsilon=-1$ the second group is uncoupled and in the range $-1<\varepsilon<0$ the repulsive coupling is weaker than the attractive coupling. For $\varepsilon=0$ their magnitudes are identical and for $\varepsilon>0$ the repulsive coupling dominates.

Introducing the Kuramoto mean fields for both groups, $Z_{a}=\rho_{a}e^{i\Theta_{a}}=1/N_{a}\sum e^{i\varphi_{j}}$, $Z_{r}=\rho_{r}e^{i\Theta_{r}}=1/N_{r}\sum e^{i\psi_{j}}$, and the common forcing

[TABLE]

we re-write the model in a compact form as

[TABLE]

For the further analysis we restrict ourselves to the case of equally sized groups $N_{r}=N_{a}=N/2$ and $\alpha_{a}=\alpha_{r}=0$. Equation (3) then reduces to

[TABLE]

We notice that according to the Watanabe-Strogatz (WS) theory Watanabe and Strogatz (1993, 1994) the dynamical description of $n>3$ identical oscillators subject to a common force can be reduced to equations for three global variables and $n-3$ constants of motion. Thus, for $N_{a,r}>3$ and $\omega\neq 0$ the model (4,5) is in fact 6-dimensional and can be described by two coupled systems of WS equations, see Ref. Pikovsky and Rosenblum (2008). For $\omega=0$ all oscillators become identical and the whole ensemble can be described by three WS equations.

III Synchronous state

First we analyze conditions of existence and stability of a synchronous state, where $\varphi_{j}=\varphi$ and $\psi_{j}=\psi$ for all $j$ and observed frequencies are $\dot{\varphi}=\dot{\psi}=\nu$. Notice that generally $\varphi\neq\psi$, i.e. synchrony in this setup shall be understood as existence of a two-cluster state. Notice also that for $\varepsilon<-1$ both groups are attractive and synchronize regardless of $\omega$, therefore we are interested in the interval $\varepsilon>-1$. Let $\varphi=\nu t$, $\psi=\nu t+\psi_{0}$, and $\Phi=\nu t+\Phi_{0}$. Then real and imaginary parts of Eq. (6) provide

[TABLE]

Condition of existence.

Subtracting Eq. (4) from Eq. (5) and using $\dot{\psi}_{0}=0$ we find that

[TABLE]

Writing the second term as $\sin\Phi_{0}\cos\psi_{0}-\cos\Phi_{0}\sin\psi_{0}$ and excluding $\sin\Phi_{0}$ and $\cos\Phi_{0}$ using Eqs. (7) we obtain

[TABLE]

It follows, that synchrony does not exist for $\varepsilon=0$, when attraction and repulsion are balanced. For $\varepsilon>0$ the repulsion becomes stronger than attraction and therefore the synchronous two-cluster state cannot be expected either. This consideration yields the border of the synchronous domain for $\varepsilon<0$:

[TABLE]

In order to find the observed frequency $\nu$ we expand (5) and insert (7). Together with (9) this yields

[TABLE]

Notice that the ratio $(1+\varepsilon)/\varepsilon$ is negative in the region of existence, so that two synchronous clusters rotate in the direction, opposite to the one determined by $\omega$. (We remind that we consider the motion in a frame, co-rotating with the natural frequency of the attractive group.)

Condition of stability.

The next step is to determine stability of the two-cluster configuration. For this purpose we first consider the linear stability of the repulsive cluster with respect to a symmetric perturbation Yeldesbay, Pikovsky, and Rosenblum (2014). It means that phases of two perturbed oscillators become $\psi_{\pm}=\nu t+\psi_{0}\pm\alpha$, where $\alpha\ll 1$. This assures that the mean field $Z_{r}$ remains unchanged in the first-order approximation in $\alpha$. The perturbed oscillators then evolve according to

[TABLE]

In the first order in $\alpha$ we find

[TABLE]

Thus, the cluster is stable for $h\cos(\Phi_{0}-\psi_{0})>0$. With the help of Eqs. (7) this condition can be re-written as

[TABLE]

Hence, the border of stability is determined by the condition $\cos\psi_{0}=1+\varepsilon$. Now, using Eq. (9), we exclude $\psi_{0}$ and obtain the stability boundary as

[TABLE]

Using the same approach for the attractive group we find the condition for the stability to be

[TABLE]

In the domain where the synchronous state exists we have $\varepsilon<0$ and the latter condition is fulfilled.

Next, we have to consider the stability of the two-cluster configuration with respect to a shift of one of the clusters. For this purpose we re-write Eqs. (4,5) for the special case of $\varphi_{j}=\varphi$ and $\psi_{j}=\psi$. Using Eq. (3) we obtain

[TABLE]

which yields the Adler equation Adler (1946) for the distance between the clusters $\delta=\psi-\varphi$:

[TABLE]

This equation has a stable fixed point for $|\omega|<-\frac{\varepsilon}{2}$, i.e. in the whole domain of existence of the two-cluster solution.

The final conclusion is that the stability of the synchronous two-cluster state is given by Eq. (15). This result fits very well the numerical results shown in Fig. 1. As one can see, the stable domain is smaller than the region where full synchrony exists.

IV Nontrivial States Beyond the two-Cluster synchrony

IV.1 Solitary State

The next solution we observe is the three-cluster state. As has been shown in Ref. Maistrenko, Penkovsky, and Rosenblum (2014), the system (4,5) with $\omega=0$, exhibits, beyond the fully synchronous one-cluster solution, a peculiar solitary state, where a cluster of $N_{a}$ attractive and $N_{r}-1$ repulsive oscillators coexists with a phase-shifted solitary oscillator. This state is not of full measure, so that not every initial condition leads to it. The range of the coupling values, where this solution exists shrinks as $1/N$ for $N\to\infty$. This makes the solitary state reliably observable only for small system sizes. The picture we observe for $\omega\neq 0$ is slightly different. Though the loss of synchrony here also occurs via appearance of a solitary unit, now one finds a three-cluster state: a cluster of $N_{a}$ attractive oscillators, a cluster of $N_{r}-1$ repulsive oscillators, and a solitary repulsive unit. The phase shifts between clusters are constant, so that the whole configuration rotates with the same constant observed frequency $\nu$. An illustration of this can be found in Fig. 2.

Condition of existence.

For a description of this state we write $\varphi=\nu t$, $\psi_{1,\ldots,N_{r}-1}=\nu t+\delta_{1}$, $\psi_{N_{r}}=\nu t+\delta_{1}+\delta_{2}$, and $\Phi=\nu t+\Phi_{0}$. This yields the equations

[TABLE]

From the last two equations it follows that $\mbox{Im}[e^{i(\Phi_{0}-\delta_{1})}]=\mbox{Im}[e^{i(\Phi_{0}-\delta_{1}-\delta_{2})}]$ and $\mbox{Re}[e^{i(\Phi_{0}-\delta_{1})}]=-\mbox{Re}[e^{i(\Phi_{0}-\delta_{1}-\delta_{2})}]$. This yields $2\delta_{1}+\delta_{2}=2\Phi_{0}-\pi$. Multiplying (20) with $e^{-i\Phi_{0}}$ and taking the imaginary part we find, by replacing $h$ with Eq. (21),

[TABLE]

By applying this relation to Eqs. (21,22) we find that observed frequency $\nu$ is described by the same Eq. (11) as in the synchronous state. However, while in the case of full synchrony $(1+\varepsilon)/\varepsilon$ was negative, here it is positive. Next, multiplying Eq. (20) by $e^{-i\Phi_{0}}$ and taking this time the real part we obtain, after replacing $\mbox{Re}(e^{i\Phi_{0}})=\sqrt{1-{\mbox{Im}(e^{i\Phi_{0}})}^{2}}$:

[TABLE]

Finally, replacing $h$ with the help of Eqs. (21,11) and introducing $x=\mbox{Im}(e^{i\Phi_{0}})$, we obtain

[TABLE]

To find the parameter domain of existence of the solitary state we need to find the range of $\omega$ so that Eq. (26) can be fulfilled for a given $\varepsilon$. First of all notice that Eq. (26) is invariant with respect to the transformation $x\to-x$ and $\omega\to-\omega$. The branch for $\omega>0$ is given by the solution for $x\in(0,1]$ and the other one can be inferred by using the transformation $\omega\to-\omega$. Consider the function $f$ consisting of the first two terms on the right hand side of Eq. (26):

[TABLE]

The border of the solitary state for $\omega>0$ can then be calculated as $\omega=\frac{\varepsilon}{2(1+\varepsilon)}f_{\max}(x,N_{r},\varepsilon)$. To find the maximum of $f$ we write $\partial{f}/\partial{x}=0$, which yields

[TABLE]

Squaring Eq. (28) and ordering it by powers of $x$ we get a cubic equation for $x^{2}$. The expression for the roots is too long to be shown here, but the calculated maximal $\omega$ for the solitary state fits the numerical results nicely, as shown in Fig. 3.

Phase shifts in the solitary state.

To determine the phase shifts $\delta_{1}$ and $\delta_{2}$, we first rewrite Eq. (2) in terms of $\delta_{1}$ and $\delta_{2}$:

[TABLE]

Next, similarly to the case of $\omega=0$ studied in Ref. Maistrenko, Penkovsky, and Rosenblum (2014), we write it as

[TABLE]

where $\tan\delta_{2}^{*}=\sin\delta_{1}/((1+\varepsilon)-\cos\delta_{1})$ and $A=\sin\delta_{1}/(2\sin\delta_{2}^{*})$. A stable state has the solution $\delta_{2}=0$ or $\delta_{2}=2\delta_{2}^{*}+\pi$. The first solution corresponds to the 2-cluster state and the second solution to the solitary state. As shown earlier the phase shifts in the solitary state are related via $2\Phi_{0}-\pi=2\delta_{1}+\delta_{2}$. This can also be expressed as $\Phi_{0}=\delta_{1}+\delta_{2}^{*}$. Equation (24) then allows one to write the relation between $\delta_{2}^{*}$ and $\Phi_{0}$ as

[TABLE]

and consequently allows for the calculation of $\delta_{2}$ and $\delta_{1}$ from $\Phi_{0}$. $\Phi_{0}$ can be calculated numerically from Eq. (26) and the resulting phase shifts coincide with the numerical results in Fig. 2.

Stability.

An analytical linear stability analysis shows that the value of $\delta_{2}$ is stable in the region of existence. Finding the stability for $\delta_{1}$ is not as simple and can only be done numerically. Still we find it to be stable in the whole region of existence for $N_{a}=N_{r}=5$. The stability analysis can be found in Appendix A.

Case $\omega=0$ vs. case $\omega\neq 0$.

Our numerical results indicate that for $\omega\neq 0$ in the parameter range where the solitary state exists, it is the only attractor. This is an essential difference with the previously studied case $\omega=0$, see Ref. Maistrenko, Penkovsky, and Rosenblum (2014), where the solitary state has not full measure. Indeed, for $\omega=0$ the system (4,5,6) admits splay state solutions $h=0$ with $\Theta_{a}=\Theta_{r}=\Phi$ and

[TABLE]

For $\omega\neq 0$ the state $h=0$ is not a solution and numerical studies indicate that the completely asynchronous case $\rho_{a}=\rho_{r}=h=0$ is unstable. Thus, the solitary state remains the only attractor.

Absence of other clustered states.

According to the WS theory Watanabe and Strogatz (1993, 1994); Pikovsky and Rosenblum (2011), the repulsive group can be fully described by two global angle variables $\Psi$ and $\Gamma$, global variable $0\leq\kappa\leq 1$, and $N_{r}$ constants $\chi_{k}$, $k=1,\ldots,N_{r}$. The latter depend on initial conditions and obey three additional constraints. The original phase variables can be obtained from the global ones with the help of the Möbius transformation Marvel, Mirollo, and Strogatz (2009); Pikovsky and Rosenblum (2015) as $e^{i\psi_{k}}=e^{i\Gamma}(\kappa+e^{i(\chi_{k}-\Psi)})/(\kappa e^{i(\chi_{k}-\Psi)}+1)$. For $\kappa<1$, general initial conditions, i.e. different $\chi_{k}$, yield different $\psi_{k}$ (for an example of such dynamics see the partially synchronous state described in the next Section). For $\kappa=1$ typically all $\psi_{k}=\psi$, i.e. one observes a one-cluster state. However, it is possible that $e^{i(\chi_{k}-\Psi)}=-1$ for some $k=n$ and then one phase $\psi_{n}$ differs from other clustered phases, i.e. the solitary state is observed Pikovsky ; Maistrenko, Penkovsky, and Rosenblum (2014). Other cluster states except for full synchrony and the $(N_{r}-1,1)$ configuration are therefore not allowed, see Ref. Engelbrecht and Mirollo (2014) for a rigorous proof. Certainly, similar consideration can be applied to the attractive group, but there the solitary state is unstable and only the trivial one-cluster state is observed.

IV.2 Self-Consistent Partial Synchronization

IV.2.1 Numerical analysis

Outside of the domains of full synchrony and solitary states we find a partially synchronized repulsive group, characterized by the order parameter $0<\rho_{r}<1$. As for the attractive group, we find that it remains synchronous even for such large values of $\varepsilon$ as 10. Though the condition of its full synchrony (16) can be easily extended for the general case of $\rho_{r}\leq 1$ to $\rho_{r}\cos(\Theta_{r}-\Theta_{a})<(1+\varepsilon)^{-1}$, we were not able to prove the synchrony analytically and only checked it numerically 111The attractive group remained fully synchronized even when the units were made non-identical by sampling the frequencies from a normal distribution with zero mean and standard deviation of $10^{-3}$. Hence, stability of the attractive group is not a numerical artifact.. A diagram of the states, including the domains of existence of full synchrony and of the solitary state, combined with the presentation of the time-averaged order parameter $\bar{\rho}_{r}$ 222In the following the time-averaged quantities are denoted by overlined letters. can be found in Fig. 4.

The observed partial synchrony can be seen as a self-organized quasiperiodic state, SOQ (or self-consistent partial synchrony, SCPS) Rosenblum and Pikovsky (2007); Pikovsky and Rosenblum (2009); Clusella, Politi, and Rosenblum (2016). The latter is characterized by the difference between the average frequency of the oscillators and their mean field. Indeed, in our setup the average frequency (observed frequency) $\bar{\nu}_{r}$ of repulsive units is larger than the average frequency $\bar{\Omega}_{r}$ of their mean field. (In fact, the instantaneous frequencies also differ nearly all the time.) Furthermore, the mismatch $\bar{\nu}-\bar{\Omega}_{r}$ increases with $\omega$. Nevertheless, both sub-populations remain synchronous on the macroscopic level, i.e. the average mean field frequencies coincide, $\bar{\Omega}_{r}=\bar{\Omega}_{a}$, see Fig. 5. We notice that close to the border of the solitary state these frequencies are not always well-defined, as indicated by small values of the minimal instantaneous order parameter. In this border domain we observe very long transients; precise identification of the dynamical states here requires a separate investigation.

Results for similar computations for a large range of $\omega$ are presented in Fig. 6. However, here the simulations were started from many different initial conditions. As one can see, partially synchronous states are characterized by a large degree of multistability: In fact, the whole range of SCPS is multistable, as can be seen in Fig. 6 as well as in Fig. 9 below. Different initial conditions result in different values of $\bar{\Omega}_{r}$ and $\bar{\nu}_{r}$ 333To obtain these quantities we have averaged the frequencies over the time interval of 500 units, after transient of 1000 units.. Interestingly, the variation of these quantities reduces with increasing $\omega$. For all these parameters the mean fields of both populations remain synchronized; we have also checked that their phases remain well-defined 444Even for such large values as $\omega=1$ and $\varepsilon=1$ the smallest observed order parameter over 100 different initial conditions was 0.08, with the average being 0.2..

Notice that transition from the solitary state to partial synchrony is accompanied by change of the direction of rotation with respect to the considered coordinate frame 555We remind that we use the frame, co-rotating with the natural frequency of oscillators in the attractive group.. Indeed, before the transition all frequencies are positive, while immediately after it they are negative, see Fig. 5. With a further increase of the parameter $\omega$, the frequency of the repulsive units $\bar{\nu}_{r}$ becomes positive and then tends to $\omega$. In fact, for large $\omega$ or for strongly repulsive systems, the repulsive units tend to have a uniform distribution of phases. However, they remain perturbed by the field of the synchronous attractive cluster, so that the uniform distribution can be reached only asymptotically.

We illustrate partially synchronous dynamics of the repulsive group by several snapshots in Fig. 7, for an intermediate value $\omega=0.1$. We see that repulsive oscillators form a group (a loose cluster), then the first oscillator in the group accelerates, stays for some instant in anti-phase with respect to others, so that we can speak about transient solitary state, and then joins the group again, now becoming the last one in the group. Then the group dissolves again, and now the oscillator that was initially the third in the group stays for some time in anti-phase to the rest of the group, then the group recombines, and so on. Notice that only every second oscillator undergoes the transient solitary state. This dynamics seem to be independent of the initial condition and was observed both for even and odd $N_{a,r}$. This bears some resemblance to a phenomenon observed in an ensemble of attractive and repulsive active rotators, see Ref. Zaks and Tomov (2016).

To conclude the discussion of the multistability of the partially synchronous state, we analyze a large system. In Fig. 8 we show two distributions of phases $\psi$ for $N_{a,r}=1024$. These distributions have been obtained by simulation started from different initial conditions: in one case, illustrated in a), we use a perturbed cluster state, while the case in b) corresponds to random initial conditions. The distributions differ in their form, as well as in their dynamics. In the first case the distribution is bounded and bimodal; it moves with time and “breathes”, changing its width. Generally the phase differences between the mean fields and common force vary in time. In the case of random initial conditions phases spread around the unit circle and their distribution is unimodal and nearly stationary (small time fluctuations are probably due to finite size effect). The differences between the distributions also lead to slight differences in the average frequencies. For the perturbed cluster we find $\bar{\nu}_{r}=0.262$ and $\bar{\Omega}_{r}=-0.196$ and for the random initial conditions we obtain $\bar{\nu}_{r}=0.264$ and $\bar{\Omega}_{r}=-0.181$.

IV.2.2 Theoretical analysis

Here we provide some analytical estimates for the state of partial synchrony. As already mentioned, for the case $\omega=0$ and partial synchrony of the repulsive units, the relation between order parameters of two groups is given by Eq. (32). Since the attractive group is always synchronized, $\rho_{a}=1$, we obtain $\rho_{r}=1/(1+\varepsilon)$. We expect that this expression can be used as an estimation also for small $\omega$. We also expect that this expression yields the upper limit for $\rho_{r}$, since an increase in $\omega$ can only lead to a decrease in the level of synchrony.

Next, we recall that according to the WS theory the description of the system (4,5,6) can be reduced to six equations for collective variables. (Below we use the WS equations in the form, suggested in Ref. Pikovsky and Rosenblum (2008).) Furthermore, we restrict the consideration to the Ott-Antonsen (OA) manifold Ott and Antonsen (2008, 2009) that corresponds to uniform distribution of the constants of motion in the WS theory Pikovsky and Rosenblum (2008). In this case the system is further simplified, with four equations for $\rho_{a,r}$ and $\Theta_{a,r}$. Moreover, since $\rho_{a}=1$, we obtain a three-dimensional system. The final equations follow from the WS equations Pikovsky and Rosenblum (2008) and read

[TABLE]

Introduction of the phase shift between the mean fields $\delta=\Theta_{r}-\Theta_{a}$ leads to the two-dimensional system

[TABLE]

Notice that since we are very far from the thermodynamic limit, the OA Ansatz can be considered only as a rather crude approximation and, hence, Eqs. (36) provide only some estimates.

We are interested in states, where the mean fields are locked, and therefore $\delta$ is bounded. We consider a weaker condition $\dot{\delta}=0$ and also neglect time variability of the order parameter, taking $\dot{\rho_{r}}=0$. Applying this approximation to Eqs. (36) we obtain an estimation for the average order parameter $\bar{\rho}_{r}$:

[TABLE]

Eliminating $\delta$ by squaring the equations and reordering terms, we obtain a cubic equation for $\bar{\rho}_{r}^{2}$. The expression for the roots are too lengthy and therefore not shown; the results for the average order parameter $\bar{\rho}_{r}$ can be seen in Fig. 9. We see that for large $\omega$ the estimation of $\bar{\rho}_{r}$ is quite good.

Given $\bar{\rho}_{r}$ we find $\delta$ from Eqs. (37). In its turn, this yields the estimation of the average frequency of the repulsive mean field $\bar{\Omega}_{r}$ from Eq. (34) as $\bar{\Omega}_{r}=\omega+(1+\bar{\rho}_{r}^{2})\sin(\delta)/4\bar{\rho}_{r}$. For a known $\bar{\Omega}_{r}$ the average frequency of an oscillator can be calculated with the help of the WS theory Baibolatov et al. (2009). Using this we find the average frequency $\bar{\nu}_{r}$ of the repulsive oscillators (for the derivation see Appendix B) to be

[TABLE]

The estimated $\bar{\nu}_{r}$ fits the numerical results in Fig. 6 for large $\omega$ quite well; the estimate $\bar{\Omega}_{r}$ is not as good, but also corresponds to the numerics for large $\omega$.

V Conclusion

We have analyzed the interplay of attraction and repulsion in a two-group Kuramoto model. In the considered network each group consists of identical elements but the groups differ in their frequencies. We have found that if attraction is stronger than repulsion then there exist an interval of frequency mismatch $\omega$ where the system synchronizes, in the sense that each group forms a cluster. The stronger the repulsion, the smaller is this interval of two cluster synchrony. The shift between synchronous clusters is determined by $\omega$. A further increase of repulsion or of $|\omega|$ destroys the two-cluster synchrony. However, the attractive group remains synchronized while the repulsive one undergoes a transition to quasiperiodic partial synchrony. In this state the order parameter of the repulsive group is between zero and one, the mean field frequency remains locked to the frequency of the attractive group, but individual units have a different, generally incommensurate, frequency. For small $|\omega|$ the transition from two-cluster synchrony to partial synchrony occurs via formation of a solitary state. In this regime there exist two clusters (one with attractive units and one with all repulsive units but one) and one solitary repulsive oscillator. The borders of synchronous and solitary regimes have been obtained analytically. We notice that the domain of the solitary state solutions rapidly shrinks with the increase of ensemble size, whereas the partial synchrony persists for large ensembles as well. For large $|\omega|$ the frequencies of the individual units and of the mean field have been estimated with the help of the WS theory. We believe that our results can be useful for analysis of neuronal ensembles with excitatory and inhibitory connections.

Acknowledgements.

This paper was developed within the scope of the IRTG 1740 / TRP 2015/50122-0, funded by the DFG/ FAPESP. M. R. was supported by the Russian Science Foundation (Grant No. 17-12-01534). The authors thank A. Pikovsky and Y. Maistrenko for helpful discussions.

Appendix A Stability of the Solitary State

The linear stability analysis of Eq. (30) with a perturbation of strength $\alpha$ yields

[TABLE]

In the first order in $\alpha$ we find $\dot{\alpha}=-A\cos\delta_{2}^{*}\alpha$ and thus the condition for stability is $A\cos\delta_{2}^{*}>0$. Since $A=\sin\delta_{1}/(2\sin\delta_{2}^{*})$ this can also be written as

[TABLE]

With the help of the definition of $\delta_{2}^{*}$ as $\tan\delta_{2}^{*}=\sin\delta_{1}/((1+\varepsilon)-\cos\delta_{1})$ the condition for stability becomes

[TABLE]

Since the solitary states exists only for $\varepsilon>0$, this condition is always fulfilled.

To demonstrate the stability of $\delta_{1}$ is not that simple. The equation for $\delta_{1}$ has the form

[TABLE]

and cannot be reduced to a form similar to Eq. (30). So, we directly substitute he $\delta_{1\pm}=\delta_{1}\pm\alpha$ and find

[TABLE]

We analyze this equation numerically, by computing $\delta_{1}$ and $\delta_{2}$ with the help of Eq. (31); this analysis shows that $\delta_{1}$ is also stable.

Appendix B Oscillator frequency in the partially synchronous state

The WS theory operates with three collective variables; two of them are angles. The first one corresponds to the maximum of the distribution of individual phases. On the OA manifold this variable coincides with the phase of the mean field. The second angle variable (we denote it as $\chi$) determines phase shift of individual oscillators with respect to the mean field (or, generally, outside of the OA manifold, with respect to the first angle variable). Correspondingly, the average frequency of units can be obtained as (see Baibolatov et al. (2009) for details):

[TABLE]

where $\Theta,\chi$ obey the WS equations

[TABLE]

Expressing $\bar{\dot{\chi}}$ via $\bar{\dot{\Theta}}$ we obtain the individual frequency

[TABLE]