Using Phase Dynamics to Study Partial Synchrony: Three Examples

A natural extension of the PRC in the case of multiple oscillators is the collective PRC, which describes the reaction of a whole ensemble of $N$ oscillators to a perturbation. Their dynamics are measured by their average, the mean-field ${\boldsymbol{Z}=1/N\sum_{j}\boldsymbol{y}_{j}}$. For a strong enough coupling, the mean-field will also move on a limit cycle, and a perturbation of the oscillators will lead to a deviation from its stable trajectory. This means the mean-field can be seen as a complex oscillator, as demonstrated for the brain rhythm in experiments with rats velazquez2015epileptic . From the reaction of the mean-field oscillator, it is even possible to gain information about the PRC of single oscillators wilson2015determining .

Whereas the PRC $\Delta$ describes the complete reaction to a perturbation, the collective PRC can be split into two parts. The prompt PRC $\Delta_{0}$ is the immediate reaction of the mean-field and the relaxation PRC $\Delta_{r}$ the part that describes the relaxation of the oscillators after the perturbation and the resulting change in their distribution. Together these form the final PRC $\Delta_{f}$ levnajic2010phase ; hannay2015collective . Formally this can be written as

where $\tau$ is the time of the perturbation, $\varphi_{0}$ is the phase at the time of perturbation of the unperturbed system and $\bar{\varphi}$ the phase in the perturbed system. Because there does not exist a general way to describe the dynamics of the mean-field, the collective PRC and its parts have to be calculated numerically.

Numerically the collective PRC is found by measuring the change in the $k$-th period of the mean-field after the perturbation $\bar{T}_{k}$ and the unperturbed period $T$ as

After a sufficiently long time this will approach the final PRC in Eq. (5).

Consider a system of $N$ coupled Rayleigh oscillators rayleigh1945theory with the dynamics

The $\eta$ is a nonlinearity parameter, $\varepsilon$ the coupling strength and ${\dot{X}=1/N\sum_{j}\dot{x}_{j}}$. The natural frequencies $\omega_{k}$ are distributed according to a Gaussian distribution with mean $1$ and standard deviation $0.01$. For $\eta=6$ and ${N=500}$ this system has a stable limit cycle of the mean-field for ${\varepsilon\gtrapprox 0.134}$ and is in the partial synchronous regime.

The oscillators are perturbed simultaneously in the $\dot{x}$-direction with a small enough strength $P$, such that the collective PRC scales linearly with $P$ and Eq. (3) is valid. The resulting response of the mean-field is shown in Fig. 1(a). While the collective PRC is flat for small $\varepsilon$, it approaches the PRC of a single oscillator with increasing coupling strength. The approach to the PRC of a single oscillator is the expected outcome, as in the case of $\varepsilon\to\infty$, the oscillators will be fully synchronized and behave like a single unit.

An important property of the collective PRC is the needed relaxation time $t_{relax}$, as the application to real-world noisy systems becomes impossible if it is too long. When the oscillators are perturbed strongly before they fully relax, then the mean-field will not reach its limit cycle, and the phase dynamics in Eq. (3) are not applicable. To have a comparable timescale consider the synchronization time $t_{sync}$ in Fig. 1(b). It is measured as a linear approximation of the time needed to reach the stable distribution of oscillators, given by $R=\left|\boldsymbol{Z}\right|\approx const$, from a splay state, where all oscillators are distributed uniformly in the phase. The slope for the linear approximation is chosen as the value at the inflection point, i.e., the biggest $\dot{R}$. A comparison of the time scales in Fig. 1(c) shows that $t_{relax}$ (at about $40T$) is longer than $t_{sync}$. This suggests a very weak attraction to the stable distribution. The collective PRC can thus be only applied in an approximate sense.

The Kuramoto model acebron2005kuramoto ; pikovsky2015dynamics describes a weakly pairwise all-to-all coupled system of oscillators. The first order approximation in Eq. (10) is valid and the dynamics can be written as

The $\omega_{k}$ are the natural frequencies and the $\alpha_{k}$ the phase shift parameters.

By using the complex mean-field $Z=Re^{i\theta}=1/N\sum_{j}e^{i\varphi_{j}}$ the system can be reduced to

The order parameter $R$ describes the degree of synchronization, $R=1$ denotes full synchrony and $R=0$ is an indicator for incoherence or antisymmetry.

One of the most important properties of this model is the analytical solvability. In the case of identical oscillators, the powerful Watanabe-Strogatz (WS) theory can be used watanabe1993integrability ; watanabe1994constants . It allows for the reduction of the dynamics from $N$ dimensions to just $3$ and $N-3$ constants of motion. The remaining degrees of freedom are three variables $\rho$, $\Psi$ and $\Theta$. Between them and the mean-field exists a correspondence, although they are not equal. Generally, $\rho$ and $\Theta$ are close to $R$ and $\theta$, while the final variable $\Psi$ can be seen as a measure of the system’s clustering. The correspondence between the degrees of freedom and the mean-field becomes even bigger in the thermodynamic limit when the dynamics move on the Ott-Antonsen manifold ott2008low ; ott2009long and reduced to just the mean-field variables, $R$ and $\theta$. The Ott-Antonsen manifold can only be seen as an asymptotic solution, as the transients can be very long pikovsky2011dynamics ; mirollo2012asymptotic .

When investigating multiple ($M$) groups of oscillators, the Kuramoto model can be extended to the M-Kuramoto model. The groups interact with different coupling strength and may have phase shifts,

Here $\sigma$ and $\sigma^{\prime}$ denote the different groups and $\varepsilon_{\sigma,\sigma^{\prime}}$ and $\alpha_{\sigma,\sigma^{\prime}}$ are the coupling strengths and phase shifts between groups $\sigma$ and $\sigma^{\prime}$, respectively.

Consider the system with ${\varepsilon_{\sigma,\sigma^{\prime}}=\varepsilon_{\sigma^{\prime}}}$ and ${\alpha_{\sigma,\sigma^{\prime}}=\alpha_{\sigma^{\prime}}}$ and two groups of identical oscillators, one attractive $\varphi^{a}$ and one repulsive $\varphi^{r}$

where the time was rescaled such that the attractive group has coupling strength $1$ and the repulsive ${-(1+\varepsilon)}$. Then $\varepsilon$ no longer measures the coupling strength, but the excess of repulsive coupling. These equations can be reduced further by introducing the mean-fields of the groups ${Z_{a,r}=1/N_{a,r}\sum e^{i\varphi_{a,r}}}$, as before, and a common forcing

The reduced equations are

and allow for the application of the WS theory to the system pikovsky2008partially .

In the system without natural frequencies ${\omega_{k}^{\sigma}=0}$ there exists a peculiar solitary state maistrenko2014solitary , see Fig. 2(a). When both groups have the same size ${N_{a}=N_{r}}$, the attractive group, and all the repulsive oscillators, except for one, will cluster at the same point. The remaining repulsive oscillator will be phase-shifted by $\pi$. In the case of ${\alpha_{\sigma^{\prime}}=0}$, the state exists for all values of $N$ in some region of $\varepsilon$ (which shrinks for increasing $N$), but it does not have full measure, i.e., some initial conditions may lead to a different state. For ${\alpha_{\sigma^{\prime}}\neq 0}$, the state only exists up to a critical system size and without full measure. The state’s existence is also predicted by the WS theory and is the only allowed clustered state, except for full synchrony.

Since the solitary state is a rather new discovery, there are slightly different definitions of it. The most popular one is that in a system with a natural order of the oscillators, e.g., ordered by their natural frequency, some single oscillators behave differently than the rest of the population and, more importantly, than their immediate neighbors. This is the main difference to a chimera state, where a whole subpopulation exhibits a different behavior. Kapitaniak et al. observed such a state in Ref. kapitaniak2014imperfect with metronomes coupled in a ring to their nearest neighbor and second nearest neighbor. In such a configuration, some single oscillators’ phases will differ from the rest of their synchronized neighbors. Another observation was made in Ref. hizanidis2016chimera , where superconducting quantum interference devices were placed in a one-dimensional array and coupled magnetically. In this case, the solitary oscillators exhibit a higher amplitude than the rest of the population.

The solitary state has also been found numerically in simulations of coupled Lorenz oscillators shepelev2018chimera and in a ring of coupled Stuart-Landau oscillators with symmetry breaking attractive and repulsive long-range coupling sathiyadevi2019long . The observation of the existence in multiplex networks of FitzHugh-Nagumo oscillators coupled in rings with a small mismatch in the intra-layer couplings mikhaylenko2019weak , allows even for controlling strategies to tune the dynamics in, e.g., neural networks.

While the existence has been shown in many different systems, there are only a few results in investigating its emergence and stability. In a Kuramoto model with inertia, the solitary state arises from a homoclinic bifurcation and persists even in the thermodynamic limit jaros2018solitary , in contrast to it being a finite size effect in the M-Kuramoto model with attractive and repulsive interaction. Aside from differential equations, the solitary state has also been found in coupled maps. Multiplex network of non-locally coupled maps with a singular hyperbolic attractor exhibit solitary states in their transition from coherence to incoherence rybalova2017transition ; semenova2017coherenceincoherence . During the transition, more and more solitary oscillators appear, growing almost linearly with the decrease in the coupling strength. This is the result of an increase of the size of the basin of attraction of the solitary set with a decrease in the coupling, as more random initial conditions lie in this basin semenova2018mechanism . To also induce solitary states in maps with nonhyperbolic attractors, a multiplicative noise can be added to the coupling constant, thus also showing the existence of the solitary state in a noisy system rybalova2018mechanism .

These solitary state can consist of multiple solitary oscillators, but from here on, a solitary state is defined more narrowly, such that a single oscillator shows a different behavior than the rest of the population.

In Ref. teichmann2019solitary we extend the M-Kuramoto model from Ref. maistrenko2014solitary to consider non-identical groups without phase shift ${\alpha_{a}=\alpha_{r}=0}$. The oscillators in each group are identical, but there exists a difference in the natural frequencies $\omega$, i.e., in Eqs. (14) ${\omega_{a}=0}$ and ${\omega_{r}=\omega}$. In this case the solitary state changes: the cluster of the repulsive oscillator has a small phase-shift in relation to the cluster of the attractive oscillators and the phase-shift between the repulsive cluster (or the attractive cluster) and the solitary oscillator is no longer $\pi$. The solitary state also has full measure, it will always be reached, regardless of the initial condition. The region of existence of the solitary state, as well as the phase shifts between the clusters and the solitary oscillator can be calculated and fits well to numerical observations, see Fig. 2(b). The state is also not stationary, but rotates with a constant frequency

Aside from the solitary, there exist two other states in the system (Fig. 3), a fully synchronous state and a self-consistent partial synchronous state. In the fully synchronous state both clusters are fully synchronized, have a constant phase shift, and rotate with a uniform frequency. The region of existence and the stability of the state can be calculated directly from Eqs. (14). The results show that the region of stability with

is slightly smaller than the region of existence, which fits the numerical results in Fig. 3. The frequency has the same relation as the solitary state in Eq. (17). We also find numerically that the attractive group always fully synchronizes, even outside the fully synchronous state, although there is no analytical proof for this. A simple check of the stability for the attractive cluster yields the inequality

The exact dynamics of the mean-field quantities are unknown, so this cannot be solved analytically, but simulations show that ${\theta_{r}-\theta_{a}}$ is quite small and $R_{r}$ falls off very quickly with $\varepsilon$ so that this equality is always fulfilled, and the attractive group fully synchronizes.

The final state, self-consistent partial synchrony, is defined by the difference in average frequencies between the single oscillators and their mean-field clusella2016minimal . The oscillators of the repulsive group move at a frequency that is generally faster than their mean-field. The WS theory cannot be applied to the whole system at once to explain this behavior, but on each group separately pikovsky2008partially . Making the crude approximation of the thermodynamic limit as well as stationarity (numerically we find that the average frequencies of both mean-fields coincide), allows for the calculation of the mean-field variables on the Ott-Antonsen manifold. As expected, they do not fit well for small $\omega$ or $\varepsilon$, but give a surprisingly good approximation for big $\omega$ and $\varepsilon$, even for a small system of ${N_{a}=N_{r}=5}$, see Fig. 4. With the now known average values of $R_{r}$ and $\dot{\theta}_{r}$ it is possible to calculate the average frequency of the oscillators in the WS theory baibolatov2009periodically as (bars denote the time-average)

which shows clearly the expected difference between the average frequency of the mean-field and the single oscillators.

Even in a simple model of phase coupled oscillators in the weak coupling limit, it is possible to find interesting dynamical states. This can be extended further by considering more complex coupling schemes than simple all-to-all coupling. For even richer dynamics, higher-order coupling terms need to be considered.

Extending the first order phase approximation for coupled oscillators in Eq. (10) to higher orders needs knowledge of the phase dynamics of the uncoupled units. One of the oscillators with a known analytical phase is the Stuart-Landau oscillator pikovsky2003synchronization ; landau1944problem . Its dimensionless dynamics in a coupled system are

where $A$ is a complex amplitude, $\omega$ is the frequency and $\gamma$ is the non-isochronicity parameter and determines $\Phi(A)$. Using ${A=Re^{i\theta}}$ leads to

From there it follows for the phase of the uncoupled system ${\varphi=\theta-\gamma\ln(R)}$ and

In Ref. matheny2019exotic , eight nanomechanical systems (NEMS) with dynamical equations resembling the Stuart-Landau oscillator were coupled in a circle. The oscillators were connected with their two neighbors, and the coupling term consisted of the average of them. The resulting phase reduction up to the second order in the coupling strength contained terms not present as physical links in the coupling scheme, such as ${\sin(\varphi_{k+2}-\varphi_{k})}$, ${\sin(\varphi_{k+2}-2\varphi_{k+1}+\varphi_{k})}$ and similar terms in the opposite direction. These non-structural terms only appear in the second-order approximation and would not be recovered in the typical first-order phase reduction. The observed system exhibits many dynamical states, such as traveling waves or weak chimeras, which would not necessarily be expected in the first-order phase approximation. Numerical simulations of the phase model show good agreement with the experimental results, pointing to the importance of the non-structural coupling terms in complex synchronization behavior.

Similarly, in Ref. leon2019phase a system of coupled Stuart-Landau oscillators was investigated numerically. The oscillators were coupled all-to-all to their mean-field. In the phase reduction up to the second-order in $\varepsilon$ appeared terms of the form ${\sin(\varphi_{m}+\varphi_{n}-2\varphi_{k})}$, where $m$, $n$ and $k$ are all possible combinations, coupling three oscillators. They determined that these terms are necessary to explain the system’s full dynamics but are again not present in the first-order phase dynamics, which only contains the existing pairwise connections.

In both Refs leon2019phase ; matheny2019exotic with coupled Stuart-Landau-like oscillators, the first-order approximation of the phase dynamics yields a Kuramoto-like model. This allows the novel terms to be seen as an extension of the Kuramoto model in Eq. (11) to higher coupling strength.

Based on the observations for the NEMS we investigated the higher order phase reduction for a system of three coupled nonidentical Stuart-Landau oscillators in Ref. gengel2020high . In difference to Ref. leon2019phase one pairwise connection is missing. The coupling scheme is a line, as can be seen in Fig. 5 and is given by

where the $c_{j,k}$ are of ${\mathcal{O}(1)}$ and the $\beta_{j,k}$ are phase shifts between the two oscillators. To extend the phase reduction beyond the first order in $\varepsilon$ we use a perturbation method, where the $R$ and $\dot{\varphi}$ are functions of the phases and expand them as a power series

In the dimensionless form the limit cycle of the uncoupled oscillator has the amplitude ${R_{0}=1}$, as can be easily checked in Eq. (22). Inserting these assumptions in Eqs. (23), the dynamics of $r^{(1)}$, $r^{(2)}$, $\psi^{(1)}$, $\psi^{(2)}$ and so on can be found by gathering the powers of $\varepsilon$. The full calculation is rather long, so please see Ref. gengel2020high for a detailed explanation. Here it will suffice to say that to find the phase dynamics in the second order in $\varepsilon$, a partial differential equation has to be solved using a Fourier series to represent the coupling function $G_{k}$. This representation is motivated by the fact that $G_{k}$ has to be $2\pi$-periodic for all the phases. The resulting reduction for the first oscillator up to the second order in $\varepsilon$ is

with $a_{k;\boldsymbol{l}}^{(j)}$ and $b_{k;\boldsymbol{l}}^{(j)}$ denoting the coefficients for cosine and sine terms respectively and $k$ being the index of the oscillator. A vector of integers $\boldsymbol{l}$ contains the coefficients before the phases and $j$ the power in $\varepsilon$.

Using the perturbation method up to the second order in $\varepsilon$ yields again non-structural terms in Eq. (26), e.g. terms of the form ${\sin(\varphi_{3}-\varphi_{1})}$. Aside from these, there exist additional terms coupling all three oscillators ${\sin(2\varphi_{2}-\varphi_{1}-\varphi_{3})}$ and second harmonics ${\sin(2\varphi_{2}-2\varphi_{1})}$. The coefficients of the second order also contain the frequency difference between the oscillators, whereas the first order terms do not.

A Fourier Ansatz is used to measure the coupling terms numerically and verify the analytical results. The phase dynamics are ${2\pi}$-periodic, so they can be written as a multidimensional Fourier series

where ${\boldsymbol{\varphi}=(\varphi_{1},\varphi_{2},\ldots,\varphi_{N})}$ is a vector of all the phases, $\boldsymbol{l}$ an N-dimensional vector of integers and ${\boldsymbol{\varphi}\cdot\boldsymbol{l}=\sum_{j}\varphi_{j}l_{j}}$ the scalar product. The $a_{k,\boldsymbol{l}}$ and $b_{k,\boldsymbol{l}}$ are Fourier coefficients. This form resembles the analytical results in Eq. (26) and finding the relevant coupling terms is reduced to fitting the Fourier coefficients $a_{k;\boldsymbol{l}}$ and $b_{k;\boldsymbol{l}}$ up to some maximum value with ${\left|l_{j}\right|\leq m}$. Because the coupling function is real, ${a_{k;-\boldsymbol{l}}=-a_{k;\boldsymbol{l}}}$ and the number of coefficients to consider halves. Still, the number of terms that need to be fitted scales like $m^{3}$, which increases the necessary number of data points for a good fit very rapidly with $m$ and leads to the curse of dimensionality.

In the case of partial synchrony, the Fourier coefficients can be fitted onto one long time series, although some initial transient has to be integrated over before the dynamics reach the torus spanned by the phases. For a good fit, the trajectory should be long enough to cover the whole torus. However, this is not the case in the synchronous regime; instead, the dynamics will settle on a single synchronized trajectory. The integration has to be stopped before reaching this synchronous trajectory and restarted with different initial conditions until the torus is sufficiently filled.

Fitting the coefficients for different coupling strengths $\varepsilon$ yields the power series

and a similar series for $b_{k;\boldsymbol{l}}$. The coupling function is then reconstructed by comparing these series and Eq. (27). The results of this method fit very well to the analytical prediction, as can be seen in Fig. 6.

The numerical method to reconstruct the phase can also be used in cases where the phase is unknown, but then the phases and their derivatives have to be calculated numerically. For finding the phase, an autonomous copy of each oscillator is integrated. It evolves for a number of its autonomous periods $T$ until it reaches its limit cycle. The phase after the relaxation is then also the phase of the point before the relaxation. To find the derivative, one observes the infinitesimal time step $dt$ of the perturbed system and sets it in relation to a different time step $\overline{dt}$ on the limit cycle. This time difference $\overline{dt}$ is determined by the motion on the limit cycle and the derivative $\boldsymbol{f}$ in Eq. (1). Using the relation of $dt$ and $\overline{dt}$ allows the calculation of the phase derivative, even if the oscillator is perturbed far from its limit cycle. For a full explanation and the resulting equation, see Ref. gengel2020high .