Coupling disorder in a population of swarmalators

In recent studies [1, 2], oscillators that sync and swarm, called swarmalators, have been considered, and interesting long-term states have been found. The swarmalators can move in the space under the correlated dynamics of the phase and space variables of each oscillator. The mobile feature of the oscillators is found to induce nonstationary states such as the active-phase wave state and the splintered phase wave state. The collective property of finite number of swarmalators is studied [3], the implementation in robots is tried as an application [4], and various steady states by finite interaction-distance is reported [5]. Mathematical property [6], model extension [7, 8, 9], and minimal modeling [10] are of basic interest. Application to swarming and flocking of biological colony [11] is a primary interest. Control of engineering objects like drones or robots is also an application front [12].

In the study of swarmalators, the interesting incoherent states are usually the consequence of negative-definite coupling strength in phase dynamics (the positive-definite case is also in Ref. [2] where noise is added). On the other hand, the coupling characteristic among the constituents of many systems is rather complex and not so simple. For example, the interaction among the neurons in the neural network systems is given by the mixture with positive and negative ones [13, 14], in general. As another example, Japanese frogs’ calling behavior can be understood by considering the mixture of positive and negative interaction with each other [15]. Probably, the spin glass  [16, 17], to which the various interesting properties of condensed matter is attributed, can be the representative example for the non-definite sign of the couplings. Considering those features in nature, the mixed coupling with positive and negative strength deserves to be considered for understanding the collective behavior of the real systems, which motivates the present study.

In this paper, we consider a population of the oscillators that can sync and swarm, governed by the random interaction in the phase dynamics. We explore how the coupling-disorder affects the long-term states in the system. In particular, we pay attention to the possibility of the phase transition [18] in the system, and focus on whether the patterns observed in the absence of coupling disorder still appear. The effective annealing of the randomly quenched coupling strengths to their average is suggested in the mobility of swarmalators to understand the numerical results. Phase transition to the fully synchronized phase, reproducibility of the phases known in the original model of the no coupling disorder, and mixture of deformed long-term states by coupling disorders are explained from the viewpoint of the annealed couplings.

This paper is organized as follows: Section II introduces the model of coupled swarmalators with random couplings, and Sec. III shows the collective behavior and the phase transition in the system. In Sec. IV we derive the threshold of the transition by the linear stability analysis of the fully synchronized state. Various long-term states including the nonstationary ones are shown in Sec. V, and the system with couplings of more than one quenched randomness are understood as the combination of annealed-coupling systems in Sec. VI. A brief summary is given in Sec. VII.

The generalized model of $N$-coupled oscillators that we consider here is given by

for $i=1,\cdots,N$, where $\theta_{i}$ and ${\bf{r}}_{i}$ represent the phase $(0\leq\theta_{i}\leq 2\pi)$ and the position vector of the $i$th oscillator, respectively. $F(r_{ij})$ is a function for the spatial dynamics of the oscillators with $r_{ij}=|{\bf{r}}_{ij}|$ for ${\bf{r}}_{ij}\equiv{\bf{r}}_{j}-{\bf{r}}_{i}$, and $G(r_{ij},\theta_{ij})$ is a function for spatial dynamics with $\theta_{ij}\equiv\theta_{j}-\theta_{i}$. $K_{ij}$ denotes the random coupling strength between the oscillators $i$ and $j$, having the symmetric property $K_{ij}=K_{ji}$. It represents a sort of bond coupling between the oscillators, instead of the site one ($K_{i}$) that affects the oscillator itself. For simplicity, we randomly chose $K_{ij}$ from the two-peaks distribution

where $p$ is the probability of the positive coupling, and $K_{p}>0$ and $K_{n}<0$. For convenience, the ratio $Q\equiv-K_{n}/K_{p}(>0)$ has been chosen for the control parameter. With the functions $F$ and $G$ in the model, the phase variable $\theta$ and the spatial one $r$ become correlated, which means the oscillators can move around in the space. Such correlation in the dynamics of the space and the phase of the oscillators has been also studied in Ref. [1], where the authors found five long-term states including the two nonstationary ones. Note that the phase coupling in Ref. [1] does not have any disorder. We here consider coupling disorder that is given by the random values from the distribution in Eq. (3).

We notice that the special case with $F(r_{ij})=1$ and $G(r_{ij},\theta_{ij})=0$ is equivalent to the mean-field $XY$ model [17] with random coupling strength, governed by the Hamiltonian ${\cal{H}}=-\frac{1}{2N}\sum_{i\neq j}K_{ij}\cos(\theta_{j}-\theta_{i})$. Equation (1) with $F=1$ and $G=0$ leads to the overdamped version of the Hamiltonian dynamics at zero temperature, and interestingly it is found that the first-order phase transition occurs at the same threshold $p_{c}$ as that we study here [19].

In this paper, we consider $F=1/r_{ij}$. The function $G$ consists of the attraction and repulsion forces acting on each oscillator, where the force functions are taken as the algebraic ones with a power like $r_{ij}^{-a}$ following Ref. [1]. The model is then given by

where $A$ and $B$ are respectively the parameters for the attractive force and the repulsive force. Here, we choose $A=B=1$ for convenience. $J$ is the parameter which measures how the phase similarity enhances the spatial proximity. For example, a positive value of $J$ means like attracts like, i.e., the swarmalators tend to be near the other swarmalators with similar phases. On the other hand, a negative value of $J$ means the opposite case: The swarmalators prefer to be near the others with opposite phases. To keep the attractive force always positive, the values of $J$ are constrained like $-1\leq J\leq 1$, following the Ref. [1]. The difference from the model of Ref. [1] is the generalization of the constant $K$ to pair-dependent $K_{ij}$. For $K_{ij}>0$, the swarmalators prefer to have the similar phase, but for $K_{ij}<0$, the opposite tendency occurs. We are interested in the effect of random $K_{ij}$ in the formation of long-term states.

To see the collective states in the system and the effects of random coupling on the states we perform the numerical simulations on Eqs. (4) and (5). We first investigate the phase synchronization behavior of the oscillators, by measuring the complex order parameter defined by [20]

where $R$ measures the phase coherence, and $\Theta$ is the mean phase angle of the oscillators. For example, $R=0$ means the incoherent state where the all oscillators have the random phase $\theta_{i}\in[0,2\pi)$, and $R=1$ is the fully synchronized state where they all have the same one $\theta_{j}=\theta_{s}$ for all $j$’s.

Also, to see the correlation between the phase angle $\theta_{j}$ and the azimutal one $\phi_{j}=\tan^{-1}(y_{j}/x_{j})$ for the $j$th oscillator, which is induced by the mobile feature of the oscillators in the system, we measure another order parameter defined by [1]

where $S_{\pm}$ is the magnitude, and $\Psi_{\pm}$ is its mean phase. Here, the order parameter $S$ measures the “correlation” between the spatial information (by the $\phi_{j}$) and the phase information (by the $\theta_{j}$). The system exhibits “assortative and/or commutative” correlation between the two depending on the initial conditions, where the $S_{+}$ comes from the assortative correlation, and $S_{-}$ comes from the commutative correlation. We chose $S$ by taking the maximum value of $S_{+}$ and $S_{-}$, i.e., $S={\rm{max}(S_{\pm})}$ [1].

In order to examine $R$ and $S$, we basically performed the numerical simulations using the Python or C programs. The total $5\times 10^{4}$ time steps with the discrete time unit $dt=0.1$ have been considered, where the first $2.5\times 10^{4}$ steps were discarded for the equilibrium, and then the time average of the quantities has been done for the later steps. And, the ten samples have been used for the average. Initial phases and positions were, respectively, randomly sampled from the uniform distributions over $[0,2\pi)$ and $[-1,1)\times[-1,1)$.

Figure 1 shows the behavior of $R$ and $S$ as a function of $p$ for various values of $J$. The value of $Q=1.25$ was chosen, and $K_{p}=1$ and $K_{n}=-1.25$ were chosen for convenience. We find that the system shows the phase transition behavior from the incoherent state ($R=0$) to the fully synchronized one ($R=1$) at a finite value of $p$, i.e., $p_{c}~{}\sim 0.6$. Interestingly, we find that the behavior of $R$ and its threshold $p_{c}$ do not depend on $J$, which can be easily understood by the absence of $J$ in the phase dynamics given by Eq. (4). On the other hand, the behavior of $S$ is different from that of $R$ as shown in Fig. 1: In the regime of the fully synchronized state ($R=1$) for $p>p_{c}$, the order parameter $S$ is zero, which implies no correlation between the phase angle and the spatial angle. However, in the incoherent state with $R=0$ for $p<p_{c}$, $S$ shows finite value, which means there is some correlation between the phase angle and the spatial angle. Moreover, the value of $S$ for $p<p_{c}$ shows $S=1$ for a certain region of $p$, and the range of $p$ of this region depends on the value of $J$ as shown in Fig. 1 (b): The range gets wider as $J$ increases.

To pin down the threshold $p_{c}$, we now estimate it numerically. We measure the value of $p$ at which the first derivative $dR/dp$ shows the maximum. Figure 2 shows the behavior of $dR/dp$ as a function of $p$ for various values of $Q$, whereas the value of $p$ at which the $dR/dp$ reaches the maximum is denoted as $p^{*}$. The inset of Fig. 2 shows $p^{*}$ as a function of $Q$, which shows a good consistency with the theoretical prediction $p_{c}=Q/(1+Q)$ (see Sec. IV).

In this section, we investigate the linear stability of the fully synchronized state with $R=1$. Let a fully synchronized state have a common phase $\theta_{i}=\theta_{\rm s}$ for all $i$. In order to examine its linear stability, we consider a slightly perturbed situation such as $\theta_{i}=\theta_{\rm s}+\phi_{i}$ with small perturbation $\phi_{i}$ for all $i$. For this setting, Eq. (4) is linearized to

Separating the summation according to the values of $K_{ij}$, we rewrite Eq. (8) as

where $C_{p}(i)$ is the collection of $j$ having $K_{ij}=K_{p}$ and $C_{n}(i)$ is that for $K_{ij}=K_{n}$.

We here note that $K_{ij}$ is not included in the spatial dynamics [see Eq. (5)]. Thus, for a fixed $i$, the fraction of such $j$’s giving $K_{ij}=K_{p}$ in the total $N$ would be approximately $p$, where the deviation is expected to decrease as $N$ increases. Similarly, that for $K_{ij}=K_{n}$ is given by $(1-p)$. This consideration motivates us to replace $\sum_{j\in C_{p}(i)}$ and $\sum_{j\in C_{n}(i)}$ with $p\sum_{j(\neq i)}$ and $(1-p)\sum_{j(\neq i)}$, respectively, in Eq. (9) when $N$ is large enough. Therefore, for large $N$, Eq. (9) can be cast into

where $\langle K_{ij}\rangle\equiv pK_{p}+(1-p)K_{n}$ is the average of $K_{ij}$ for the two-peak distribution $h(K_{ij})$.

Introducing $m_{ij}=-1/(Nr_{ij})$ for $i\neq j$ and $m_{ii}=-\sum_{j\neq i}m_{ij}$, one writes Eq. (10) as

Interestingly, the matrix $m=\{m_{ij}\}$ is positive-definite for large $N$ as follows. It holds $|m_{ij}|\sim{\cal O}(1/\sqrt{N})$ for off-diagonal terms. This is based on the numerical observation that all swarmalators reside in the two-dimensional spatial region of the ${\cal O}(1)$-area and there is no concentration of them. Thus, $r_{ij}\sim 1/\sqrt{N}$ holds, from which $|m_{ij}|\sim{\cal O}(1/\sqrt{N})$ follows. It also follows that $m_{ii}=-\sum_{j\neq i}m_{ij}\sim{\cal O}(\sqrt{N})$. Then, as $N$ increases, $m=\{m_{ij}\}$ approaches $\bar{m}=\{\bar{m}_{ij}\}$, which is given by

where $\delta_{ij}$ is the Kronecker delta that assigns $1$ when it assigns $i=j$ while $0$ otherwise. The matrix $\bar{m}$ is positive-definite because $\bar{m}_{i}\equiv\lim_{N\to\infty}m_{ii}>0$. Therefore, $m$ is also positive-definite for sufficiently large $N$.

When the positive definiteness of $m$ is considered in Eq. (11), one knows that the sign of $\langle K_{ij}\rangle$ determines whether $\phi_{i}$ will decay or not. This is because every eigenvalue of the positive-definite matrix is positive. Thus, $\phi_{i}$ decays when $\langle K_{ij}\rangle>0$, otherwise it grows when $\langle K_{ij}\rangle<0$. Since $Q\equiv-K_{n}/K_{p}$, we get $\langle K_{ij}\rangle=K_{p}(p+(p-1)Q)$. Therefore, $\langle K_{ij}\rangle>0$ corresponds to $p>Q/(1+Q)$ and $\langle K_{ij}\rangle<0$ corresponds to $p<Q/(1+Q)$, which leads to

For $Q=1.25$, used for the numerical data shown in Figs. 1 and 2, Eq. (13) gives $p_{c}=5/9=0.55..\,$. The deviation from the numerical value $p_{c}\approx 0.6$ in the figures is supposed to be a finite $N$ effect. Note $p_{c}\approx 0.565$ of Fig. 2 for $N=800$ is closer to the prediction $p_{c}=5/9$ than the prediction in Fig. 1 for $N=200$. It is interesting that the threshold value shown in Eq. (13) is equivalent to that for the mean-field $XY$-type oscillators with random coupling disorder [19]. This means that the spatial dynamics in the system does not change the threshold value of the phase transition from the incoherent state to the fully synchronized one.

The crucial step of our argument above is the use of $\sum_{j\in C_{p}(i)}=p\sum_{j(\neq i)}$ and $\sum_{j\in C_{n}(i)}=(1-p)\sum_{j(\neq i)}$ to write Eq. (10). We here remark this is basically an annealed approximation. In fact, $K_{ij}$ is quenched by definition in the model and has never been treated in an annealed way in the numerical work. Thus the consistency between Eq. (13) and its numerical counterpart is interesting. Equation (5) gives a clue to understanding how the quenched $K_{ij}$ works as if annealed. The spatial dynamics governed by Eq. (5) is indifferent to $K_{ij}$. Based on this, one may expect the distribution of $K_{ij}$ for given $i$ has no position-dependence, i.e., the fractions for $K_{p}$ and $K_{n}$ are, respectively, $p$ and $(1-p)$ independent of position. This results in the annealed average of phase dynamics in Eq. (4). The consistency between the numerical data and the annealed average, even in the presence of quenched $K_{ij}$, means that the oscillators move to the proper places. That is, the effective annealing of $K_{ij}$ to $\langle K_{ij}\rangle$ is a characteristic property of swarmalator that can move.

We also examine that the sign of $\langle K_{ij}\rangle$ always governs the onset of the sync state in further numerical study using other $h(K_{ij})$ distributions. These observations suggest that the annealed approximation could be generally available. We remark that the effective annealing requiring proper locations of swarmalators is considered for long-term states.

In this section, we investigate the incoherent state $(R=0)$ further. In particular, we pay attention to the possibility of the nonstationary states such as the active phase wave (APW) state and the splintered phase wave (SPW) state found in Ref. [1]. Therein, the all coupling strength among the oscillators was chosen as the negative one with no disorder like $K_{ij}=K$, which corresponds to $p=0$ in the current study. We note that, in Ref. [1], various long-term states such as the async state, the APW state, the SPW state, and the static phase wave state were found depending on the strength of the negative coupling in the phase dynamics. Especially, the APW state and the SPW state are the nonstationary states where the oscillators are not static in both phase and space. It is thus interesting whether the dynamics of Eqs. (4) and (5) with quenched disorder $K_{ij}$ can induce such a nonstationary feature. Below, we again resort to the numerical data by the same scheme used for Fig. 1.

To see the possibility of the nonstationary states for $p>0$, we measure the following two quantities. One is the mean velocity $V$ defined by [1, 2]

where $v_{i}=\sqrt{{{\dot{x}}_{i}}^{2}+{{\dot{y}}_{i}}^{2}}$. The mean velocity measures the nonstatic property of the state. In other words, the finite value of $V$ means that the swarmalators move around in the plane in the long-term state. The other one is [1, 2]

where $N_{\rm{rot}}$ is the number of swarmalators circulating at least $2\pi$, and thus $U$ represents the fraction over the total swarmalators [1, 2]. Note that $U$ can play as the role of the indicator for the active phase wave state [1]. Figure 3(a) and (b) show the quantities $R$, $S$, $V$, and $U$ as a function of $p$ for $J=0.8$ and $0.2$, respectively. As $p$ increases from zero, we find that the phase transition occurs from the incoherent state $(R=0)$ to the fully synchronized state ($R=1$) at a finite value of $p_{c}$, as shown in the behavior of $R$. Interestingly, the value of $p_{c}$ for two different values of $J$ are found to be same, which means the threshold does not depend on $J$. On the other hand, the behavior of $S$, $U$, and $V$ is found to differ depending on the value of $J$ (See Fig. 3).

When $R=0$, note that the system shows a finite value of $S$ through the splintered and static phase wave. Also, the system shows the active phase wave state for small $p$ when $J$ is large. This active phase wave does not occur when $J$ is small. We find that the system shows both the stationary states and the nonstationary states depending on $p$ for a given value of $Q$ and $J$. Especially, when $J$ is large $(J=0.8)$, all five long-term states that have been reported in Ref. [1] exist depending on $p$. It is interesting to note that the nonstationary states such as active phase wave sate and splintered phase wave also exist in the presence of the positive coupling among the swarmalators in the system ($p>0$).

Above, we have seen that the result by random $K_{ij}$ is comparable with that by constant $K=K_{ij}$. The long-term state patterns, including the fully synchronized one, from both settings are similar for large $J$. It is remarkable that such analogy holds in the presence of finite faction for negative $K_{ij}$ that may bring about the frustration [21, 22] in the phase dynamics. This strongly suggests the annealed property conjectured for the onset of the sync phase at $\langle K_{ij}\rangle=0^{+}$ is also valid for non-zero $\langle K_{ij}\rangle$ cases. This motivates us to test the long-term states in the system with such quenched random couplings by mixture of different randomnesses. We are interested in whether the characters of long-term states attributed to each randomnesses will last or not.

Consider $N_{1}$ number of swarmalators in total $N$, and call them group 1 (G1). The coupling strength $K_{ij}$, therein, is assigned with $1$ for probability $p_{1}$ or $-1$ for $1-p_{1}$. Next, consider the other group, G2, of $N_{2}(=N-N_{1})$ swarmalators, where $K_{ij}=\pm 1$ is similarly assigned using probability $p_{2}$. We also consider coupling between the swarmalators in both groups as assigning $1$ to $K_{ij}$ for probability $p_{\rm m}$ or $-1$ for $1-p_{\rm m}$. Below, we use $p_{\rm m}=p_{2}$ for simplicity. When $(N_{1},p_{1})=(100,1)$ and $(N_{2},p_{2})=(100,0)$, one may expect the sync disc in G1 and the APW state in G2 if the couplings that work between the groups do not overwhelm each character. Note the average coupling in G2, $\langle K_{ij}\rangle_{2}=-1$, is in the region for the APW state in the phase diagram reported in Ref. [1]. We obtained numerical data with time unit $dt=0.01$ and consider $2\times 10^{5}$ time steps.

Figure 4(a) shows the numerical result.

Two patterns appear; one is a slightly deformed sync disc while the other is U-shaped one. We have observed the former vibrates a little (not demonstrated here) and the latter is an active phase state. The inset in Fig.4(a) is the trajectory of a swarm in the U-shape. Since it sweeps the region while changing the phase, we name the state the deformed APW (dAPW). It is interesting to note that this reminds us of the “chimera state” observed in the oscillators with identical frequency [23, 24, 25, 26]: The partial group of the swarmalators in the system exhibits the sync cluster, but the other group of the swarmalators do not join the sync cluster, showing the U-shape APW.

We vary $N_{1}$ to see a possible change in pattern. One easily see the distance between the sync disc and the dAPW increases when comparing Figs. 4(a) and 5(a). We also see the sync disc becomes more circular and static when the distance increases.

We next test whether a sync disc by $p_{1}<1$ in G1 only system ($N_{2}=0$) will still last after G2 is involved ($N_{2}>0$). The interest is a validity of the annealed approximation of $K_{ij}$ in G1. Figure 5(b) is the result for $p_{1}=0.7$, which is also composed of the sync disc and the dAPW. The sync disc is almost circular and static with enough separation from the dAPW.

We finally test whether G2 can show the SPW state while G1 remains a sync disc. As the splintered state appears for the negative coupling larger than that for APW, we use $p_{2}=0.4$, which corresponds to $\langle K_{ij}\rangle_{2}=-0.2$, the average coupling strength in G2. This value is in the region for the splintered wave in the phase diagram reported in Ref. [1]. Figure 4(b) is the result when the other settings are same as those for Fig. 4(a). There still appear sync disc and U-shaped region. This time, the swarmalators of different colors in the U-shaped region are not mixed, different from the dAPW case. The gradual change of the swarmalators’ colors in space is a characteristic of the splintered phase wave. The black solid curve near the bottom is the trajectory of a swarm, which remains near there as time goes on. As shown, it does not sweep the U region but localized in the small area while splintering. This is always the case for all swarmalators in the U-shaped region. We thus regard the U-shaped pattern as a deformed splintered phase wave (dSPW).

In the above numerical tests, the characters by each random couplings are preserved more or less in the mixed system. There apparently appear sync disc, APW, and splintered state though deformed. We here add the result is qualitatively similar for the various values of $p_{\rm m}$ (not shown here), and this does not depend on the initial condition. These observations suggest that the long-term states, demonstrated in Secs IV and V, are still robust more or less in the mixed system. Interestingly, the patterns survive the frustration taking place between the subsystems G1 and G2. A mixture with async state or of more than two groups was not tested this time. Further study on various interesting properties of the mixed deformed patterns will appear elsewhere.

We considered the population of swarmalators with random coupling strength, and explored how the coupling disorder affects the long-term states in the system. In particular, the possibility of the phase transition and the robustness of the state patterns are focused. To understand the long-term states observed in the system with quenched disorder of coupling strength, we considered the effective annealed approximation, which is mediated by the mobility of swarmalators. In the viewpoint of annealed couplings, the numerical observation in the quenched system is explained and, furthermore, such a system of the mixture of different quenched disorders is also understood.

We found that the system shows the phase transition from the incoherent state to the fully synchronized state at a certain threshold $p_{c}$, where the value of $p_{c}$ is argued in the linear stability analysis of the fully synchronized state. We also found that, in the regime of the incoherent state below the threshold, various long-term states are found to exist. Especially, the nonstationary states such as the dAPW and the dSPW besides the normal APW and SPW are discovered.

All long-term states known for the bare model in Ref. [1] are realized in the system of random coupling strength. The values of average-coupling where each state appears are similar to their counterparts in the bare model. This strongly suggests the randomly quenched coupling strengths work as if annealed. The pattern by the supposed annealed coupling is so robust that the mixture of different random couplings leads to the proper combination of each deformed patterns.

This research was supported by NRF Grant No. 2021R1A2B5B01001951 and ‘Research Base Construction Fund Support Program’ funded by Jeonbuk National University in 2021 (H.H), and by NRF Grant No. 2018R1D1A1B07049254 (H.K.L.).