Cluster synchronization is a fundamental phenomenon in systems of coupled oscillators. Here, we investigate clustering patterns that emerge in a unidirectional ring of four delay-coupled electrochemical oscillators. A voltage parameter in the experimental setup controls the onset of oscillations via a Hopf bifurcation. For a smaller voltage, the oscillators exhibit simple, so-called primary, clustering patterns, where all phase differences between each set of coupled oscillators are identical. However, upon increasing the voltage, secondary states, where phase differences differ, are detected, in addition to the primary states. Previous work on this system saw the development of a mathematical model that explained how the existence, stability, and common frequency of the experimentally observed cluster states could be accurately controlled by the delay time of the coupling. In this study, we revisit the mathematical model of the electrochemical oscillators in order to address open questions by means of bifurcation analysis. Our analysis reveals how the stable cluster states, corresponding to experimental observations, lose their stability via an assortment of bifurcation types. The analysis further reveals complex interconnectedness between branches of different cluster types. We find that each secondary state provides a continuous transition between certain primary states. These connections are explained by studying the phase space and parameter symmetries of the respective states. Furthermore, we show that it is only for a larger value of the voltage parameter that the branches of secondary states develop intervals of stability. For a smaller voltage, all the branches of secondary states are completely unstable and are, therefore, hidden to experimentalists.

Understanding emergent behavior in collective dynamics is important in various applications, for example, neuronal firing and animal behavior. The use of controlled experimental setups provides a testbed for relevant theory on the subject. In previous work,1 it was shown that clustering behavior of four delay-coupled electrochemical oscillators in a unidirectional ring can be modeled by a network of paradigmatic Stuart–Landau oscillators with surprisingly good agreement to experimental results. In particular, depending on the choice of delay time, one can control the possible synchronization clustering patterns. In this paper, we study how the various experimentally observed clustering patterns are actually interconnected via an assortment of bifurcations, with some patterns serving as transitional states between other patterns. This additional insight into the overall structure of how the different clustering patterns are organized with respect to the control parameter provides a new perspective to the experimental observations.

Nonlinear systems of coupled oscillators form the basis of multiple areas of interdisciplinary research, from the dynamics of coupled lasers to modeling neuronal dynamics. Of particular importance and widespread interest is the phenomenon of synchronization and clustering.2–6 For example, in neuronal dynamics, different patterns of synchronization are related to normal cognitive and pathological functions of the brain.7 In power grid dynamics, analyzing synchronization and possible cluster states can provide insights into the stability of the grid.8,9

Clustering behavior is often associated with underlying symmetry properties of a system.10,11 For such systems, group theory and equivariant dynamical system theory can provide a useful tool for exploring and understanding the dynamics.12,13 For example, the symmetry properties of a network of oscillators can be used to facilitate model reductions.14 Knowledge of the symmetries can also be used to discover possible cluster synchronization patterns.9 Furthermore, the authors of Ref. 15 show that nodes in a complex network will form clusters according to their symmetry properties within the network. It is, therefore, suggested that structural connectivity of the brain could play a role in neural synchronization across distant locations.

There have also been many studies on the interplay between dynamics and symmetry in the context of networks with ring topologies, which will be the focus of this paper. For example, in Refs. 16 and 17, a type of equivariant delayed feedback control is used to target the stabilization of periodic orbits with a specified spatiotemporal pattern in rings of oscillators. In Refs. 18 and 19, rings of oscillators are used to study the gaits (i.e., walking patterns) of animals. As a result, it is suggested that transitions between different gaits can be modeled as symmetry-breaking bifurcations: for example, a horse will walk, trot, and then gallop as a result of successive bifurcations.

In more recent years, greater attention has been paid to the fact that many processes are not instantaneous, but instead possess an inherent delay. Such delays can potentially have a crucial influence on the overall dynamics of the system.20–23 In certain cases, it can turn trivial dynamics into complex dynamics.24 

In order to better understand the effects of delay and symmetry on synchronization and clustering, both experimental research and mathematical modeling have an important role to play—one can inform the other of where interesting dynamical phenomena may occur and provides clues to possible mechanisms behind the phenomena. For example, in the context of mutually coupled lasers, mathematical modeling has explained the various dynamical regimes through the existence of mutually synchronized symmetric and symmetry-broken states, which are connected through pitchfork bifurcations.25–27 Clustering has also been studied in networks of physically dissimilar mechanical and electrical oscillators,28,29 where the interplay between modeling and experiment focuses experiments on interesting parametric regions and highlights which dynamics may be most realizable experimentally.

The work presented here builds upon a previous study on the clustering patterns observed in an experimental setup of four electrochemical oscillators coupled in a unidirectional ring.1 It was shown that a delay in the coupling between each set of neighboring oscillators could control which clustering pattern would emerge. Another important parameter of the experimental setup (see  Appendix A for a brief review) is the choice of applied voltage. It is only above a critical voltage that the oscillators actually oscillate via a Hopf bifurcation. Above this critical value, two voltage regimes were studied. One lower voltage regime resulted in smooth sinusoidal oscillations. This was called the smooth regime. In the second regime, a higher voltage resulted in relaxation oscillations and was called the relaxation regime.

The clustering patterns that were observed depended on the regime. In the smooth regime, only clustering patterns with equal phase differences between neighboring oscillators (neighboring in the sense of ring network topology) were found. These are called the primary clustering patterns. A clear visualization of the clustering patterns in the smooth and relaxation regimes is crucial for appreciating the experimentally motivated modeling, as well as reconnecting the theoretical analysis of the model back to the experimental observations. Therefore, we show examples of the different primary clustering patterns in Fig. 1, as both time series in terms of current I(t) and schematic diagrams. Panels (a) and (b) show an example of an in-phase state, where all oscillators possess an equal phase. Panels (c) and (d) depict an example of the splay state, where all neighboring oscillators have a common phase difference of φ j + 1 φ j = π / 2 with their neighbors. A two-cluster state is shown in panels (e) and (f), where neighboring oscillators have a phase difference of φ j + 1 φ j = π. Similarly, panels (g) and (h) show a reverse splay state, with phase differences φ j + 1 φ j = 3 π / 2. We calculate phase via peak-to-peak linear interpolation.30 

FIG. 1.

Experimental time series and schematic diagrams of primary cluster states in the smooth regime with V 0 = 1.105V and K = 0.15: (a) and (b) in-phase with τ = 0.95 × ( 2 π ω ), (c) and (d) splay with τ = 1.25 × ( 2 π ω ), (e) and (f) two-cluster with τ = 0.50 × ( 2 π ω ), and (g) and (h) reverse splay states with τ = 0.70 × ( 2 π ω ).

FIG. 1.

Experimental time series and schematic diagrams of primary cluster states in the smooth regime with V 0 = 1.105V and K = 0.15: (a) and (b) in-phase with τ = 0.95 × ( 2 π ω ), (c) and (d) splay with τ = 1.25 × ( 2 π ω ), (e) and (f) two-cluster with τ = 0.50 × ( 2 π ω ), and (g) and (h) reverse splay states with τ = 0.70 × ( 2 π ω ).

Close modal

In the relaxation regime, the primary clustering patterns were still observed; however, additional secondary clustering patterns were observed with phase differences that were not all equal. The first example shown in Figs. 2(a) and 2(b) is a compressed two-cluster state. It resembles a two-cluster state, only that the phase difference between the two clusters is less than π. Panels (c) and (d) show an example of a compressed reverse splay state. While three of the phase differences between neighboring oscillators appear to remain equal, one of the phase differences is less than 3 π / 2. Similarly, a compressed splay state, such as the one shown in panels (e) and (f), appears to maintain three equal phase differences between neighboring oscillators, while one phase difference becomes greater than π / 2. Finally, open two-cluster states are often observed in the experiments, as depicted in panels (g) and (h). This state resembles a two-cluster state, where the oscillators within each cluster appear to have drifted apart.

FIG. 2.

Experimental time series and schematic diagrams of secondary cluster states in the relaxation regime with V 0 = 1.2V and K = 0.10: (a) and (b) compressed two-cluster with τ = 0.65 × ( 2 π ω ), (c) and (d) compressed reverse splay with τ = 0.81 × ( 2 π ω ), (e) and (f) compressed splay with τ = 1.06 × ( 2 π ω ), and (g) and (h) open two-cluster states with τ = 0.51 × ( 2 π ω ).

FIG. 2.

Experimental time series and schematic diagrams of secondary cluster states in the relaxation regime with V 0 = 1.2V and K = 0.10: (a) and (b) compressed two-cluster with τ = 0.65 × ( 2 π ω ), (c) and (d) compressed reverse splay with τ = 0.81 × ( 2 π ω ), (e) and (f) compressed splay with τ = 1.06 × ( 2 π ω ), and (g) and (h) open two-cluster states with τ = 0.51 × ( 2 π ω ).

Close modal

A well-calibrated model of Stuart–Landau oscillators was found to match the transitions between different clustering patterns, as well as their common frequency, in the smooth regime with very good accuracy. Furthermore, an extension to the classic Stuart–Landau oscillators with experimentally derived coupling functions was found to agree remarkably well with the experimental results in the relaxation regime. The results presented in Ref. 1 effectively took the form of bifurcation diagrams in terms of the coupling delay parameter, with the branches of various clustering patterns in the experiment being accurately matched by stable branches of corresponding solutions to the model. However, because the focus of that study was on the accurate reproduction of experimental results, the analysis of the model was limited to only the stable behavior. While successful, this approach left important open questions. It was still unclear how these states lose their stability. Perhaps more intriguing is the question: Where do the secondary states in the relaxation regime come from, when there is no evidence of their existence in the smooth regime?

In this paper, we conduct a systematic bifurcation analysis of the model for the unidirectional ring of delay-coupled electrochemical oscillators that reveals the roles of the unstable solutions and various bifurcation types in the organization of clustering patterns with respect to the coupling delay. More specifically, we consider each type of secondary cluster state observed in the experiments and identify how they relate to other cluster states. A useful tool in this context is the investigation of the symmetry groups of the various cluster states. We find that each secondary cluster state can be interpreted as a transitional cluster state, en route between certain primary states via symmetry-breaking bifurcations. Primary states are associated with symmetry groups of order 4, while secondary states have symmetry groups of order 2 or 1. Finally, we find that although the general interconnected structure of clustering patterns is preserved in the smooth regime, all secondary cluster states become unstable and can, therefore, not be directly observed in experiments.

In Sec. II, we introduce the model of delay-coupled oscillators in both a general and experiment-specific framework, where the electrochemical oscillators are modeled by Stuart–Landau oscillators. We also investigate the phase space and parameter symmetries of the model. The results of systematic analysis of the transitions between cluster states in the relaxation regime are presented in Sec. III. Finally, we discuss the results, how they relate to dynamics in the smooth regime, and further open questions in Sec. IV.

In this section, we will present the model behind the numerical analysis. First, we will review the general framework of coupled Stuart–Landau oscillators in the context of cluster synchronization. In Sec. II B, we generalize the coupling scheme in the model for cases of non-sinusoidal coupling. Then, in Sec. II C, we describe clustering solutions from the model and discuss their symmetries.

Let us consider N supercritical Hopf normal forms, also known as Stuart–Landau oscillators, that are unidirectionally coupled in a ring configuration with transmission delays,
z ˙ j ( t ) = [ λ + i ω ( 1 + i γ ) | z j ( t ) 2 | ] z j ( t ) + K z ( j + 1 ) mod N ( t τ ) ,
(1)
where z j ( t ) = r j ( t ) e i φ j ( t ) C denotes the jth complex dynamical variable, j = 1 , , N. λ, ω 0, and γ are real-valued constants, where γ couples the frequency to the oscillation amplitude. The parameters K R and τ denote the coupling strength and delay, respectively. In this study, we will focus on the formation of synchronized behavior as the delay is varied. In polar coordinates, that is, using radius r j and phase φ j variables, Eq. (1) can be rewritten as follows:
r ˙ j ( t ) = [ λ r j ( t ) 2 ] r j ( t ) + K r j + 1 ( t τ ) cos [ φ j + 1 ( t τ ) φ j ( t ) ] ,
(2a)
φ ˙ j ( t ) = ω γ r j ( t ) 2 + K r j + 1 ( t τ ) r j ( t ) sin [ φ j + 1 ( t τ ) φ j ( t ) ] ,
(2b)
where all indices have to be taken modulo N throughout this paper.
As shown in Ref. 1, this simple model reproduces the correct clustering patterns and transitions between different patterns observed in the experiments for the smooth regime very well. However, in contrast to the smooth regime, the relaxation regime is not near the Hopf bifurcation that initiates oscillations. The coupling between oscillators cannot, therefore, be assumed to be sinusoidal, and we require a more general representation of the coupling. We rewrite the model using general interaction functions H i, i { r , φ },
r ˙ j ( t ) = [ λ r j ( t ) 2 ] r j ( t ) + K r j + 1 ( t τ ) H r [ φ j + 1 ( t τ ) φ j ( t ) ] ,
(3a)
φ ˙ j ( t ) = ω γ r j ( t ) 2 + K r j + 1 ( t τ ) r j ( t ) H φ [ φ j + 1 ( t τ ) φ j ( t ) ] .
(3b)
In this paper, we use the radial and angular interaction functions as determined in Ref. 1 (cf. Sec. V therein). In short, the experimentally obtained functions H r and H φ can each be approximated by a fifth-order Fourier series, which yields the following system of equations:
r ˙ j ( t ) = [ λ r j ( t ) 2 ] r j ( t ) + K r j + 1 ( t τ ) ( n = 0 5 a n , r cos { n [ φ j + 1 ( t τ ) φ j ( t ) ] } + b n , r sin { n [ φ j + 1 ( t τ ) φ j ( t ) ] } ) ,
(4a)
φ ˙ j ( t ) = ω j γ r j ( t ) 2 + K r j + 1 ( t τ ) r j ( t ) ( n = 0 5 a n , φ cos { n [ φ j + 1 ( t τ ) φ j ( t ) ] } + b n , φ sin { n [ φ j + 1 ( t τ ) φ j ( t ) ] } ) ,
(4b)
where all indices have to be taken modulo N. The Fourier coefficients a n , r, b n , r, a n , φ, and b n , φ that best match experimentally observed interaction functions can be found in  Appendix B. The values of the other parameters in the model are chosen to match the experiments, as discussed in Ref. 1. Note that γ = 0 throughout, so that the frequency does not depend on the radius.

We investigate the occurrence of cluster states that are characterized by a common collective frequency Ω m. The cluster states come in two flavors. Primary states have a collective amplitude and equal phase lags Δ φ m = 2 π m / N between neighboring oscillators; i.e., r j r 0 , m and φ j ( t ) = Ω m t + j Δ φ m. Secondary states differ in amplitudes and phase lags but still oscillate with a common frequency.

The integer m = 0 , , N 1 labels the specific primary cluster state. Following Ref. 13, these states may be characterized as discrete rotating waves when m is coprime to N, discrete standing waves when m = 0, and a discrete alternating wave when m = N / 2. In our case of N = 4, m = 0 corresponds to the in-phase state [cf. Fig. 1(a)], m = 1 to the splay state [cf. Fig. 1(b)], m = 2 to the two-cluster state [cf. Fig. 1(c)], and m = 3 to the reverse splay state [cf. Fig. 1(d)].

Let us briefly discuss the symmetries present in the dynamical system (4) for the case N = 4. As we will see throughout our analysis, it is helpful to relate the different clustering patterns observed in the experiments to their symmetry groups. This will aid us in our goal of understanding how they are related to each other and how the system transitions from one state to another.

There are two phase space symmetries, namely, the symmetry under a global shift of all phase variables by a common phase φ s, as well as the symmetry under cyclic permutations of the coordinate indices. More precisely, if we are given a particular solution ( r ( t ) , φ φ ( t ) ) T, with r ( t ) = ( r 1 ( t ) , , r 4 ( t ) ) T and φ φ ( t ) = ( φ 1 ( t ) , , φ 4 ( t ) ) T, then the functions defined via
( r ^ ( φ s , j ) ( t ) φ φ ^ ( φ s , j ) ( t ) ) = ( G j r ( t ) G j φ φ ( t ) + φ s b )
(5)
are also solutions of (4) for all φ s [ 0 , 2 π ) and j { 0 , , 3 }. Here, G is the unidirectional coupling matrix defined in Eq. (A4), which cyclically permutes coordinate indices j j 1 and b = ( 1 , 1 , 1 , 1 ) T. A pair of the form ( ϕ s , j ), therefore, performs a symmetry operation on the set of all solutions of (4). The set of all symmetry operations ( φ s , j ) forms a group with neutral element ( 0 , 0 ) and composition rule,
( φ s , 1 , j 1 ) ( φ s , 2 , j 2 ) = ( φ s , 1 + φ s , 2 mod 2 π , j 1 + j 2 mod 4 ) .
(6)
Let us denote this group by H. We note that H is isomorphic to the direct sum U ( 1 ) Z 4, where U ( 1 ) = { e i ϕ : ϕ [ 0 , 2 π ) } is the unitary group in one dimension and Z 4 = Z / 4 Z = { 0 , 1 , 2 , 3 } is the cyclic group with four elements. In mathematical terms, we can also express the above observations by saying that the dynamical system (4) is H-equivariant.31 
Now that we have identified the phase space symmetry group of our system, we can use it to classify the different solutions of the system based on their symmetry. We are particularly interested in cluster states, where all oscillators share a common frequency. For example, the splay state (i.e., the primary state with m = 1 introduced in Sec. II A) is found to be invariant under the action of the element ( π / 2 , 1 ) H since φ ^ j = φ j + 1 + π / 2 = φ j. More generally, a primary state of index m is invariant under a subgroup
H 4 , m = { ( m j π / 2 , j ) : j = 0 , , 3 } H .
(7)
Note that the order of all H 4 , m is four, and they are isomorphic to Z 4. The secondary states are less symmetric, and the three options for invariant subgroups are
H 2 , 0 = { ( 0 , 0 ) , ( 0 , 2 ) } ; H 2 , 1 = { ( 0 , 0 ) , ( π , 2 ) } ; H 1 , 0 = { ( 0 , 0 ) } .
(8)
H 2 , 0 and H 2 , 1 are isomorphic to Z 2, and H 1 , 0 is the trivial group with only the neutral element. The groups of the clustering patterns observed in the experiments are summarized in Table I. Note that the pattern corresponding to H 2 , 1, which was not observed in the experiments, will be discussed in Sec. III D.
TABLE I.

Symmetry groups of cluster states.

Pattern Group Representation
In-phase  Z 4  H4,0 
Splay  Z 4  H4,1 
Two-cluster  Z 4  H4,2 
Reverse splay  Z 4  H4,3 
Compressed two-cluster  Z 2  H2,0 
Compressed splay  Z 1  H1,0 
Compressed reserve splay  Z 1  H1,0 
Open two-cluster  Z 1  H1,0 
Pattern Group Representation
In-phase  Z 4  H4,0 
Splay  Z 4  H4,1 
Two-cluster  Z 4  H4,2 
Reverse splay  Z 4  H4,3 
Compressed two-cluster  Z 2  H2,0 
Compressed splay  Z 1  H1,0 
Compressed reserve splay  Z 1  H1,0 
Open two-cluster  Z 1  H1,0 
In addition to the phase space symmetry, the dynamical system (4) also possesses an important parameter symmetry for cluster states. For a particular τ, let us assume that we are given a solution of (4) of the form
( r ( t ) φ φ ( t ) ) = ( r 0 φ φ 0 + Ω t b ) ,
(9)
with r 0 , φ φ 0 R 4 and Ω R. With c = ( 0 , π / 2 , π , 3 π / 2 ) T, it then follows that
( r ^ ( t ) φ φ ^ ( t ) ) = ( r 0 φ φ 0 + c + Ω t b )
(10)
is a solution of (4) for τ ^ = τ + π / ( 2 Ω ). This can be seen by noting that
φ ^ j + 1 ( t τ ^ ) φ ^ j ( t )
(11)
= φ 0 , j + 1 + ( j + 1 ) π 2 ( φ 0 , j + j π 2 ) Ω ( τ + π 2 Ω )
(12)
= φ j + 1 ( t τ ) φ j ( t ) .
(13)

It is important to note that because of the presence of the vector c in the transformation, the symmetry group of the cluster state changes. For example, an in-phase state at delay τ, which is invariant under H 4 , 0, will become a splay state at τ ^ with symmetry H 4 , 1. More generally, H 4 , m becomes H 4 , m + 1, and the groups H 2 , 0 and H 2 , 1 get exchanged.

Repeated application of the above transformation yields a parameter symmetry group, which is isomorphic to Z. Explicitly, for every j Z, the transformation
( r ^ 0 φ φ ^ 0 τ ^ Ω ^ ) = ( r 0 φ φ 0 + j c τ + j π 2 Ω Ω )
(14)
describes a parameter transformation between cluster states in the ( τ , Ω ) plane. As our original system is only defined for positive τ, we will restrict ourselves to that case. We also note that stability is not necessarily conserved under this parameter symmetry. This is in contrast to phase space symmetry, where stability information is preserved.

In the following, we explain the nature of the experimentally observed, secondary cluster states by exploring the transitions between different primary cluster states. We demonstrate these transitions by means of bifurcation diagrams for the model given by Eq. (4) with N = 4. We will demonstrate that already, this small number gives rise to rich dynamical scenarios.

The results in this paper are obtained numerically with the continuation software DDE-Biftool.32,33 Throughout this section, we fix the system parameters as λ = 2.89, ω = 2.43, and K = 0.189, as in the previous work on the same experimental setup.1 In order to facilitate a thorough analysis, we reduce the model by replacing the phase variables with phase-difference variables (see  Appendix C for details). In particular, phase differences with respect to oscillator 4 are considered so that oscillator 4 becomes the reduced system’s rotating frame of reference. This means that periodic solutions are represented by relative equilibria; for example, a two-cluster periodic solution is represented by a relative equilibrium, where all phase differences are π. Therefore, in the bifurcation diagrams below, bifurcations of period solutions are represented by their counterpart bifurcations of relative equilibria; e.g., torus bifurcations are represented by Hopf bifurcations.

Figure 3 shows the bifurcation diagram of primary states in the ( τ , Ω )-plane; that is, we vary the coupling delay and depict the collective frequency of the synchronized cluster state. Red plus signs and blue crosses denote Hopf and pitchfork bifurcations, respectively, and stability/instability of the states are indicated by solid/dotted line styles. For better visualization, we use a color code to highlight the phase difference between oscillators 1 and 4; that is, φ 1 φ 4 [ 0 , 2 π ]. Therefore, the blue, green, yellow, and red curves correspond to in-phase, splay, two-cluster, and reverse splay states, respectively. The different states are connected via the parameter symmetry (14), and this explains the striking similarity between the various branches.

FIG. 3.

Bifurcation diagram of primary states in the ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. The curve color indicates the phase difference between oscillators 1 and 4. Red plus signs/blue crosses denote Hopf/pitchfork bifurcations. Parameters are λ = 2.89, ω = 2.43, and K = 0.189.

FIG. 3.

Bifurcation diagram of primary states in the ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. The curve color indicates the phase difference between oscillators 1 and 4. Red plus signs/blue crosses denote Hopf/pitchfork bifurcations. Parameters are λ = 2.89, ω = 2.43, and K = 0.189.

Close modal

One can see how the different cluster states change in their collective frequency and how the stable branches that are experimentally accessible (see Ref. 1 for full branches) are connected by unstable branches. Stable branches connect the maximum and minimum of the collective frequency curves and also appear at intermediate frequencies around the intrinsic frequency ω. The transitions between stability and instability occur via pairs of Hopf and pitchfork bifurcations in very close proximity to each other. We will use Fig. 3 as an overall reference plot and explore how secondary states emerge in that bifurcation diagram. For this purpose, we will display the curves of the respective secondary states on top of the relevant primary states and adjust the range of the τ-axis accordingly.

As a brief guide, we consider the secondary states: compressed two-cluster state (Sec. III A), compressed reverse splay state (Sec. III B), compressed splay state (Sec. III C), and open two-cluster state (Sec. III D). Aiming for an intuitive understanding, we also provide examples of solutions along the branches in the phase plane, which depict steps along the transition between different primary states.

Let us start by considering the compressed two-cluster state, shown schematically in Fig. 2(b). This is a configuration, where oscillators group in alternating pairs {1,3} and {2,4}. The compression arises from a phase lag between the two pairs that is different from π.

Figure 4 illustrates how compressed two-cluster states emerge as connecting branches between the in-phase (blue) and two-cluster solutions (yellow). The bifurcation diagram highlights the recursive nature of these connecting branches: every time the in-phase solutions encounter a pitchfork bifurcation, a branch of compressed two-cluster solutions is born, which terminates at a pitchfork bifurcation of two-cluster solutions. This diagram also highlights how the phase differences along the primary branches are fixed. For example, the phase difference φ 1 φ 4 is either fixed at 0 or π. On the other hand, along the secondary branches, the phase differences vary, providing the first clue to their role as transitional states. Here, the phase difference φ 1 φ 4 varies continuously between 0 and π or π and 2 π.

FIG. 4.

Bifurcation diagram showing transitions between the two-cluster and in-phase states in the ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Other parameters as in Fig. 3.

FIG. 4.

Bifurcation diagram showing transitions between the two-cluster and in-phase states in the ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Other parameters as in Fig. 3.

Close modal

Figure 5 provides a closer look at the transition between in-phase and two-cluster states. Panel (a) is a zoomed-in version of the bifurcation diagram in the ( τ , Ω )-plane shown in Fig. 4. Here, we see clearly that the branch of compressed two-cluster states undergoes pitchfork bifurcations, creating intervals of stable solutions, which can be observed experimentally. Indeed, the stable interval of compressed two-cluster states with approximately τ [ 0.63 , 0.76 ] was shown experimentally in Ref. 1 [see Fig. 8(a) therein]. However, it is not immediately clear from panel (a) that this branch emerges from pitchfork bifurcations, and that there are in fact two branches of compressed two-cluster states that overlap each other in the ( τ , Ω )-plane. Therefore, in panel (b), we show the same diagram in the ( τ , r 4 )-plane. This allows us to distinguish the different states via their respective radius. Now, we see the two branches of compressed two-cluster states (with Z 2-symmetry) that connect the in-phase and two-cluster branches (with Z 4-symmetry). Considering the color scheme, we see that one transition involves the phase difference φ 1 φ 4 increasing from π to 2 π, while the other involves the same phase difference decreasing from π to 0.

FIG. 5.

Bifurcation diagram showing a transition between the two-cluster and in-phase state in (a) the ( τ , Ω )-plane and (b) the ( τ , r 4 )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Open and closed circles denote the solutions shown in the phase plane in panels (c)–(e) and (f)–(h), respectively. The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

FIG. 5.

Bifurcation diagram showing a transition between the two-cluster and in-phase state in (a) the ( τ , Ω )-plane and (b) the ( τ , r 4 )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Open and closed circles denote the solutions shown in the phase plane in panels (c)–(e) and (f)–(h), respectively. The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

Close modal

Panels (c)–(e) of Fig. 5 are example solutions, corresponding to the open circles in the bifurcation diagrams, and panels (f)–(h) are the symmetrically related solutions, corresponding to the closed circles. Panels (c) and (f) both appear similar to a two-cluster state, except that now, the φ j + 1 φ j = φ j φ j 1 symmetry is broken as a result of the pitchfork bifurcation. Panels (d) and (g) are examples where the two clusters of oscillator pairs have moved closer to each other and are now stable. This example is typical of solutions observed experimentally and corresponds to the schematic diagram shown in Fig. 2(b); note the quantitative agreement of the delay times with the experimental observation. Finally, in panels (e) and (h), the two clusters are so compressed that they are almost in-phase.

In the context of symmetry, the compressed two-cluster states are invariant under the group H 2 , 0. This is a subgroup of both H 4 , 0 and H 4 , 2, which are the symmetry groups of the in-phase state and the two-cluster state, respectively. It, therefore, follows that the compressed two-cluster state can connect to those two states through pitchfork bifurcations.

Figure 6(a) shows a branch of compressed reverse splay states, together with branches of reverse splay and in-phase states. The example solutions in panels (b)–(d) depict a gradual transition from a reverse splay to an in-phase configuration. The stable solution shown in panel (c) is typical of experimental observations and corresponds to the schematic diagram shown in Fig. 2(d); again, note the quantitative agreement of τ with the experiment.

FIG. 6.

Bifurcation diagram showing a transition between the reverse splay and in-phase state in the (a) ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Circles denote the solutions shown in the phase plane in panels (b)–(d). The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

FIG. 6.

Bifurcation diagram showing a transition between the reverse splay and in-phase state in the (a) ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Circles denote the solutions shown in the phase plane in panels (b)–(d). The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

Close modal

In this case, however, the branch of secondary cluster states is not connected to the primary states branches by a simple pair of pitchfork bifurcations. This is clear from the fact that the primary branches are Z 4, while this branch of secondary cluster states possesses only Z 1 symmetry. Figure 7 provides further details on how the branch of compressed reverse splay states is related to the branch of reverse splay states. Panels (a) and (b) show the reverse splay branch losing stability at a pitchfork bifurcation. At this initial pitchfork bifurcation, the φ j + 1 φ j = φ j φ j 1 symmetry is lost, and two branches of solutions with Z 2 symmetry are born. Then, very shortly afterward, further pitchfork bifurcations take place, at which the φ j + 2 φ j = φ j φ j 2 symmetry is lost, resulting in four branches with Z 1 symmetry. This is in agreement with the equivariant branching lemma,13 which characterizes the isotropy subgroups of bifurcating solution branches. Only with both sets of symmetries gone, can all oscillators approach each other to form the compressed reverse splay state. As a result of the sets of pitchfork bifurcations, the compressed reverse splay branch shown in Fig. 6(a) is actually four symmetrically related curves overlapping in the ( τ , Ω )-plane. Figures 7(c)7(f) show an example of four symmetrically related solutions in the phase plane, each with the same parameters and collective frequency Ω as Fig. 6(c).

FIG. 7.

Pitchfork bifurcations, denoted by blue crosses, of the reverse splay state in (a) the ( τ , Ω )-plane and (b) the ( τ , r 4 )-plane. Solid/dotted curves represent stable/unstable solutions. Panels (c)–(f) are solutions resulting from the pitchfork bifurcations for the same parameters and collective frequency Ω as Fig. 6(c). Other parameters as in Fig. 3.

FIG. 7.

Pitchfork bifurcations, denoted by blue crosses, of the reverse splay state in (a) the ( τ , Ω )-plane and (b) the ( τ , r 4 )-plane. Solid/dotted curves represent stable/unstable solutions. Panels (c)–(f) are solutions resulting from the pitchfork bifurcations for the same parameters and collective frequency Ω as Fig. 6(c). Other parameters as in Fig. 3.

Close modal

The connection between the compressed reverse splay branch and the in-phase branch in Fig. 6(a) also requires a more detailed inspection. Figure 8(a) provides a zoomed-in view of the bifurcation diagram. Here, the vertical axis represents the mean collective frequency of the four oscillators, Ω ¯, because the bifurcation diagram now involves periodic solutions. At the Hopf bifurcation, denoted by the red plus sign, the in-phase solutions lose stability and a branch of unstable periodic solutions emerge. To emphasize, these are periodic solutions of the phase differences of the oscillators. In the original non-rotating frame of reference used in Eq. (4), these solutions correspond to solutions on a torus. The periodic solutions, born at the Hopf bifurcation, fulfill the Z 4 symmetry condition φ j + 1 ( t ) = φ j ( t + T / 4 ), where T is the period. This condition persists as they approach the four saddles on the four symmetrically related Z 1 branches of compressed reverse splay states, where a heteroclinic bifurcation takes place, denoted by the green diamond. A periodic solution approaching the heteroclinic bifurcation is shown in panel (b) as a projection onto the ( φ 1 φ 4 , φ 2 φ 1 )-plane, where green circles are the saddles. Panels (c)–(f) show zoomed-in plots of the saddles in the phase plane, providing a clear illustration of the symmetry. Finally, we see in panel (g) that the time series appears to plateau at the values of φ 1 φ 4 that correspond to the four saddles. This shows that, as expected with increasing period toward a heteroclinic bifurcation, the periodic solution begins to spend more time in the vicinity of the saddles.

FIG. 8.

(a) Bifurcation diagram showing the connection between the in-phase branch and secondary cluster states involving periodic solutions in the ( τ , Ω ¯ )-plane, where Ω ¯ is the mean collective frequency. Solid/dotted curves represent stable/unstable solutions. The plus sign, triangle, and diamond denote a Hopf, fold, and heteroclinic bifurcation, respectively. (b) A periodic solution close to the heteroclinic bifurcation projected onto the ( φ 2 φ 4 , φ 1 φ 4 )-plane, where green circles represent saddles. Panels (c)–(f) show close-ups of the four saddles shown in panel (b) in the ( x , y ) phase plane. (g) Time series of φ 1 φ 4 corresponding to the solution shown in panel (b). Dotted lines correspond to the four saddles. Other parameters as in Fig. 3.

FIG. 8.

(a) Bifurcation diagram showing the connection between the in-phase branch and secondary cluster states involving periodic solutions in the ( τ , Ω ¯ )-plane, where Ω ¯ is the mean collective frequency. Solid/dotted curves represent stable/unstable solutions. The plus sign, triangle, and diamond denote a Hopf, fold, and heteroclinic bifurcation, respectively. (b) A periodic solution close to the heteroclinic bifurcation projected onto the ( φ 2 φ 4 , φ 1 φ 4 )-plane, where green circles represent saddles. Panels (c)–(f) show close-ups of the four saddles shown in panel (b) in the ( x , y ) phase plane. (g) Time series of φ 1 φ 4 corresponding to the solution shown in panel (b). Dotted lines correspond to the four saddles. Other parameters as in Fig. 3.

Close modal

The compressed splay state arises in an analogous fashion to the compressed reverse splay state, except that here the transition is between the in-phase and splay primary branches. Figure 9(a) shows the corresponding bifurcation diagram in the ( τ , Ω )-plane, with a branch of compressed splay states connecting the two primary state branches. In this case, the connection to the in-phase branch is via multiple pitchfork bifurcations, i.e., the mechanism demonstrated in Fig. 7, while the connection to the splay branch is via a branch of periodic solutions, i.e., the mechanism demonstrated in Fig. 8.

FIG. 9.

Bifurcation diagram showing a transition between the in-phase and splay state in the (a) ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Circles denote the solutions shown in the phase plane in panels (b)–(d). The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

FIG. 9.

Bifurcation diagram showing a transition between the in-phase and splay state in the (a) ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Circles denote the solutions shown in the phase plane in panels (b)–(d). The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

Close modal

Panels (b)–(d) show example solutions along the compressed splay branch, indicated by the circles in panel (a). Again, the stable solution, shown in panel (c), agrees with the schematic diagram of experimentally observed compressed splay states shown in Fig. 2(f) for the same value of τ. The striking similarity between Figs. 6(a) and 9(a) is again explained through the connection of the branches via the parameter symmetry of (14). However, we stress that this connection only holds for cluster states and not limit cycles. It also does not preserve stability of states. This explains that observed minor difference between Figs. 6(a) and 9(a), which indicates that here, the transition branch encounters a pair of Hopf bifurcations.

FIG. 10.

Bifurcation diagram showing a transition between the reverse splay and splay state in (a) the ( τ , Ω )-plane and (b) the ( τ , r 4 )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Open and closed circles denote the solutions shown in the phase plane in panels (c)–(e) and (f)–(h), respectively. The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

FIG. 10.

Bifurcation diagram showing a transition between the reverse splay and splay state in (a) the ( τ , Ω )-plane and (b) the ( τ , r 4 )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Open and closed circles denote the solutions shown in the phase plane in panels (c)–(e) and (f)–(h), respectively. The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

Close modal

In the model, the simplest way for open two-cluster states to emerge is along branches connecting the splay and reverse splay branches. Figure 10 shows two (symmetrically related; Z 2) open two-cluster branches in (a) the ( τ , Ω )-plane and (b) the ( τ , r 4 )-plane that connects to the ( Z 4) primary branches at pitchfork bifurcations. Note that the curves of solutions take on the exact same forms as the curves in Fig. 5, although different cluster types are involved. This is due to the parameter symmetry (14). While in Fig. 5 the connecting branch has H 2 , 0 symmetry, the branch in Fig. 10 has H 2 , 1 symmetry. Therefore, throughout the transition from reverse splay to splay state in panels (c)–(e) and (f)–(h) of Fig. 10, we see how oscillators 1 and 3 maintain a phase difference of π, as do oscillators 2 and 4. However, the phase differences between neighboring oscillators shift, allowing pairs of oscillators to approach each other. In panels (d) and (g), we see an example of pairs of neighboring oscillators with a zero phase difference, resembling a case of a “closed” two-cluster state. Note that this state is different from the primary two-cluster state, where the clusters are given by pairs of oscillators { 1 , 3 } and { 2 , 4 } and all radii are equal.

Although the secondary branch between the reverse splay and splay branches in Fig. 10 provides a mechanism for solutions that could be described as open two-cluster, they do not correspond to the open two-cluster states observed in the electrochemical oscillator experiments. Comparison with the schematic diagram in Fig. 2(h) reveals that the phase difference of π between oscillators 1 and 3, as well as 2 and 4 is not present in the open two-cluster states observed experimentally.

Figure 11 shows the transition branch between the two-cluster and reverse splay branches that provides the mechanism for the open two-cluster states, as they appear in the experiments. The details of the connections between this open two-cluster branch and the primary branches are analogous to Figs. 6 and 9. Generally, the phase differences between the paired oscillators in the two-cluster state increase [see Fig. 11(b)] until they approach the reverse splay state [see Fig. 11(e)]. Along the way, the clustering pattern appears as in the experiments [cf. panel (c) with Fig. 2(h)].

FIG. 11.

Bifurcation diagram showing a transition between the two-cluster and reverse splay state in the ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Open circles denote the solutions shown in the phase plane in panels (b)–(e). The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

FIG. 11.

Bifurcation diagram showing a transition between the two-cluster and reverse splay state in the ( τ , Ω )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, and triangles denote Hopf, pitchfork, and fold bifurcations, respectively. Open circles denote the solutions shown in the phase plane in panels (b)–(e). The colored circular arcs highlight the phase differences between oscillators 1 and 4. Other parameters as in Fig. 3.

Close modal

Comparing Figs. 11(c) and 11(e), it becomes clear that along the branch, oscillators 2 and 3 must pass each other. Panel (d) shows an example of a stable solution where oscillators 2 and 3 are very close together. Now, if we return to the experimental output and consider the open two-cluster state just before it disappears as τ is being increased, we indeed find a cluster state (overlooked in Ref. 1, now shown in Fig. 12) that matches Fig. 11(d) and exists for the same value of τ.

FIG. 12.

Experimental time series (a) and a schematic diagram (b) of a secondary cluster state in the relaxation regime: open two-cluster state with τ = 0.62 × ( 2 π ω ).

FIG. 12.

Experimental time series (a) and a schematic diagram (b) of a secondary cluster state in the relaxation regime: open two-cluster state with τ = 0.62 × ( 2 π ω ).

Close modal

We have analyzed a mathematical model of delay-coupled oscillators to gain a deeper understanding of experimental results of a system of delayed-coupled electrochemical oscillators. Our analysis, conducted by numerical continuation, details the bifurcations involved in the appearance and disappearance of branches of various cluster types in experiments. Furthermore, we elucidate the role of the secondary cluster states as transitions between certain primary states. While primary states have symmetry groups of order 4, secondary cluster states have symmetry groups of smaller order, and the connection between them is mediated through pitchfork bifurcations.

We explicitly demonstrate where each of the four types of experimentally observed secondary clusters belongs along the transitional routes between certain primary states. As emphasized earlier, in each case of a secondary state, what we show in the above figures are only the parts of the curves needed to explain the relevant transition. In fact, the curves are more complicated and intertwined. For example, Fig. 13 shows the complete branch with Z 1 symmetry to which the compressed splay states in Fig. 9 belong. Interestingly, it reveals that this single secondary branch actually provides continuous routes between all primary states via a plethora of fold bifurcations since each extremum of the branch connects to one of the primary branches (not shown for clarity) by the same mechanisms demonstrated in Figs. 7 and 8.

FIG. 13.

The entire branch of secondary cluster states shown only partially in Fig. 9. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs and triangles denote Hopf and fold bifurcations, respectively. Other parameters as in Fig. 3.

FIG. 13.

The entire branch of secondary cluster states shown only partially in Fig. 9. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs and triangles denote Hopf and fold bifurcations, respectively. Other parameters as in Fig. 3.

Close modal

Our results have focused on the clustering patterns, as observed in the relaxation oscillation regime. But how do these results relate to the smooth oscillation regime where the oscillators first begin to oscillate upon increasing the potential V 0? Do the secondary states also exist in this regime? Why were the secondary states not observed there?

Figure 14 is the bifurcation diagram for the smooth oscillation regime in the ( Ω ¯ , τ )-plane. It reveals that, according to the model, the primary cluster states, observed in the experiments, do possess the same bifurcations as in the relaxation regime. Furthermore, similar to in our results above, we find branches of secondary cluster states that connect the different primary state branches via pitchfork bifurcations and small branches of periodic solutions. However, this interesting behavior of secondary states is always unstable. This provides an explanation for the lack of secondary states being observed in the experiments. In the smooth regime, the secondary states are unstable and cannot be observed directly. It is only when the potential V 0 in the experimental setup is increased further, moving the system away from the smooth regime, that the bifurcation structure of the secondary state branches becomes more complex and creates pockets of stable solutions.

FIG. 14.

Bifurcation diagram of smooth oscillators in the ( Ω ¯ , τ )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, triangles, and diamonds denote Hopf, pitchfork, fold, and heteroclinic bifurcations, respectively. Parameters: λ = 1.1025, ω = 3.4228, and K = 0.3.

FIG. 14.

Bifurcation diagram of smooth oscillators in the ( Ω ¯ , τ )-plane. Solid/dotted curves represent stable/unstable solutions. Curve color indicates the phase difference between oscillators 1 and 4. Plus signs, crosses, triangles, and diamonds denote Hopf, pitchfork, fold, and heteroclinic bifurcations, respectively. Parameters: λ = 1.1025, ω = 3.4228, and K = 0.3.

Close modal

When we adapt the original model (2) for general coupling forms in model (3), we tailor the relevant terms to mimic the coupling observed in the experiments in the relaxation regime, as outlined in Sec. II B. This turns out to be sufficient for achieving a very good agreement between the model and experiments. For example, in this paper, we match examples of various cluster solutions from the model to experimental observations for the same value of τ. Furthermore, the model reproduces the bifurcation structure inferred from the experiments with good quantitative agreement regarding oscillation frequency and τ (see Figs. 6 and 8 of Ref. 1). Therefore, in contrast to the modeling of the coupling, the modeling of the local dynamics of the oscillators, which is not adapted in model (3), appears to be less important for the study of the clustering behavior of the network.

Further work could include bridging the gaps between the two snapshots we have of the smooth regime ( V 0 = 1.105) and the relaxation regime ( V 0 = 1.2). The key to this would be studying how the interaction functions evolve between the regimes. Another obvious avenue for further work would be to investigate the influence of the coupling strength K, which would undoubtedly be crucial for synchronization behavior.

The observation of the consistently close pairing of Hopf and pitchfork bifurcations along the curves of primary states could warrant further study. Similar observations are made in other systems.26,27 How generic is this observation in relation to cluster synchronization? Is it even a dynamical requirement that a Hopf bifurcation occur in close proximity to the pitchfork bifurcations along the primary branches?

While we have found very good agreement between the secondary clustering patterns in the model and the experimental observations, there is one pattern in the model that was not observed experimentally, namely, the open two-cluster state, as shown in Figs. 10(d) and 10(g). Of course, this could be due to the model not capturing some aspect of the network of coupled oscillators that exclude this particular state from observations. Alternatively, one possible explanation, which would require a closer study, is that those solutions simply have relatively small basins of attraction compared to other co-existing solutions. For example, the range of τ, for which the secondary branch in Fig. 10 is stable, lies within the range of τ for which the in-phase branch is stable. Perhaps, when the oscillators transition between cluster states, some secondary branches are overlooked because of the size of their basin of attraction at the time of the transition.

Finally, we have remarked throughout the paper that the unstable solutions cannot be directly observed in the experiments. One possible approach to indirectly observe these states experimentally would be to apply Pyragas control in order to stabilize the unstable solutions.34 This would require implementing a delayed self-feedback to each of the four oscillators and setting this delay time equal to the period of the unstable cluster state to be stabilized. Assuming the model to be accurate also for the unstable states, this period would be equal to 2 π / Ω. One challenge might be finding appropriate control strengths to successfully stabilize the unstable solutions; nonetheless, it would be fascinating to explore how well the model captures not only the stable, but also unstable clustering patterns of the electrochemical oscillations.

This work was dedicated in memory of John L. Hudson.

This project was supported by the University College Cork in the framework of the SEFS New Connections Grant Award scheme. PH acknowledges further support by Deutsche Forschungsgemeinschaft (DFG) under Project ID 434434223—SFB 1461. We are also grateful to two anonymous referees for valuable feedback, which helped us improve the manuscript.

The authors have no conflicts to disclose.

Andrew Keane: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Alannah Neff: Data curation (equal); Formal analysis (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Karen Blaha: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Andreas Amann: Methodology (equal); Formal analysis (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Philipp Hövel: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The motivation for the theoretical study presented in this paper is to provide a thorough understanding of the transitions between different clustering patterns during experiments with electrochemical oscillators. The observed clustering patterns are introduced in Ref. 1, as well as numerous examples of transitions that occur as the control parameter; i.e., the coupling delay, is varied. Some additional (quantitatively different) runs are shown in this paper in Figs. 1 and 2, and a (qualitatively different) run is shown in Fig. 12. Therefore, for the convenience of the reader, we provide a brief summary of the experiment here. For further details, we refer the reader to Ref. 1.

We conduct experiments performed in an electrochemical cell consisting of four 1-mm-diameter Ni working electrodes (99.98% pure), a Pt mesh counter electrode, and an Hg/Hg2SO4/K2SO4 (sat) reference electrode with a 3M H2SO4 electrolyte. We electrically couple the four electrodes in a unidirectional ring and measure the electrochemical dissolution currents I j , j = 1 , 2 , 3 , 4 , via zero resistance amperemeters (ZRAs). We select four oscillators with similar uncoupled frequencies from an array of 64 oscillators. We can vary the character of the oscillators with our choice of applied voltage. Lower voltages closer to the Hopf bifurcation are nearly harmonic; higher voltages produce higher harmonic oscillations that appear less sinusoidal. These two dynamical behaviors are the above-mentioned smooth and relaxation oscillations, respectively. In these experiments, we observe smooth and relaxation oscillators for an applied voltage of V 0 = 1.105 and 1.2 V, respectively.

Negligible intrinsic electrical interactions exist between uncoupled oscillators.35 The startup or shutdown of one oscillator does not alter the behavior of the others, and oscillator dynamics are independent when uncoupled.

We introduce coupling of the form
V j ( t ) = V 0 + δ V j ( t ) ,
(A1)
where V j , j = 1 , 2 , 3 , 4, denotes the voltage between working and reference electrodes and δ V j is the change in the circuit potential of the jth element due to feedback. The feedback voltages are
δ V j ( t ) = K n = 1 4 g j n [ V n ( t τ ) R P I ^ n ( t τ ) ] ,
(A2)
where R P = 650 Ω is a resistance, K denotes the overall coupling gain, and τ is the coupling time delay; we impose time delayed coupling via the real-time data acquisition system with a multichannel potentiostat. I ^ j is the normalized current measured by the ZRAs such that
I ^ j ( t ) = A ¯ A j ( I j ( t ) I ¯ j ) ,
(A3)
where A j and I ¯ j are the amplitude and mean current of oscillator j and A ¯ is the mean amplitude of the population, n = 1 4 I n max / 4.
The coupling structure is determined by g j n, which belongs to the adjacency matrix G. We apply unidirectional coupling in a four-member ring with
G = ( 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 ) ,
(A4)
which is implemented via the multichannel potentiostat.

The Fourier coefficients of the radial and angular interaction functions, which are used in Eq. (4), are obtained in Ref. 1 and are summarized in Table II for the reader’s convenience. They are scaled to normalize the angular interaction function max | H φ | = 1.

TABLE II.

Fourier coefficients of Eq. (4).

Radial interaction function
a 0 , r = 0.455 79   
a 1 , r = 0.979 48  b 1 , r = 1.823 54 
a 2 , r = 0.361 10  b 2 , r = 0.079 63 
a 3 , r = 0.297 24  b 3 , r = 0.548 54 
a 4 , r = 0.058 46  b 4 , r = 0.090 98 
a 5 , r = 0.115 58  b 5 , r = 0.092 51 
Angular interaction function   
a0,φ = 0   
a 1 , φ = 0.006 10  b 1 , φ = 0.316 22 
a 2 , φ = 0.358 11  b 2 , φ = 0.290 20 
a 3 , φ = 0.253 41  b3,φ = −0.055 85 
a 4 , φ = 0.135 41  b 4 , φ = 0.007 99 
a 5 , φ = 0.071 83  b 5 , φ = 0.004 25 
Radial interaction function
a 0 , r = 0.455 79   
a 1 , r = 0.979 48  b 1 , r = 1.823 54 
a 2 , r = 0.361 10  b 2 , r = 0.079 63 
a 3 , r = 0.297 24  b 3 , r = 0.548 54 
a 4 , r = 0.058 46  b 4 , r = 0.090 98 
a 5 , r = 0.115 58  b 5 , r = 0.092 51 
Angular interaction function   
a0,φ = 0   
a 1 , φ = 0.006 10  b 1 , φ = 0.316 22 
a 2 , φ = 0.358 11  b 2 , φ = 0.290 20 
a 3 , φ = 0.253 41  b3,φ = −0.055 85 
a 4 , φ = 0.135 41  b 4 , φ = 0.007 99 
a 5 , φ = 0.071 83  b 5 , φ = 0.004 25 
In order to simplify the numerical analysis, we rewrite the model (3) with N = 4 in terms of phase differences relative to oscillator 4: φ ^ j = φ j φ 4. We approximate the terms in the interaction functions of Eq. (3) as φ j + 1 ( t τ ) φ j ( t ) = φ ^ j + 1 ( t τ ) + φ 4 ( t τ ) φ ^ j ( t ) φ 4 ( t ) φ ^ j + 1 ( t τ ) φ ^ j ( t ) φ ˙ 4 ( t ) τ and rewrite the model as
r ˙ j ( t ) = [ λ r j ( t ) 2 ] r j ( t ) + K r j + 1 ( t τ ) H r [ φ ^ j + 1 ( t τ ) φ ^ j ( t ) φ ˙ 4 ( t ) τ ] ,
(C1a)
φ ^ ˙ j ( t ) = ω γ r j ( t ) 2 + K r j + 1 ( t τ ) r j ( t ) H φ [ φ ^ j + 1 ( t τ ) φ ^ j ( t ) φ ˙ 4 ( t ) τ ] φ ˙ 4 ( t ) , j = 1 , 2 , 3.
(C1b)

Equation (C1b) now represents the phases of oscillators 1, 2, and 3 within the rotating frame of oscillator 4. Therefore, the system only requires seven dimensions [i.e., j = 1 , , 4 for Eq. (C1a) and j = 1 , 2 , 3 for Eq. (C1b)]. Note that when j = 3, the term φ ^ j + 1 ( t τ ) in Eq. (C1b) becomes zero.

Of course, there is also the variable φ ˙ 4 ( t ), which must be solved for every time Eq. (C1) is evaluated. This is done by solving the nonlinear implicit equation
0 = ω γ r 4 ( t ) 2 + K r 1 ( t τ ) r 4 ( t ) H φ [ φ ^ 1 ( t τ ) φ ˙ 4 ( t ) τ ] φ ˙ 4 ( t ) .
(C2a)
Since DDE-Biftool runs in Matlab, we use the function fzero to solve Eq. (C2).

This reduction raises the question of whether the above approximation φ 4 ( t τ ) φ 4 ( t ) φ ˙ 4 ( t ) τ has a significant effect on the stability of the solutions. Therefore, we calculated a sample of the above results without the reduction and confirm no noticeable difference between the two sets of results. As an additional check, we also confirm the stability of various solutions by means of numerical simulation.

1.
K.
Blaha
,
J.
Lehnert
,
A.
Keane
,
T.
Dahms
,
P.
Hövel
,
E.
Schöll
, and
J. L.
Hudson
, “
Clustering in delay-coupled smooth and relaxational chemical oscillators
,”
Phys. Rev. E
88
,
062915
(
2013
).
2.
J.
Juang
and
Y.-H.
Liang
, “
Cluster synchronization in networks of neurons with chemical synapses
,”
Chaos
24
,
013110
(
2014
).
3.
M.
Lodi
,
F.
Della Rossa
,
F.
Sorrentino
, and
M.
Storace
, “
Analyzing synchronized clusters in neuron networks
,”
Sci. Rep.
10
,
16336
(
2020
).
4.
P. R.
Protachevicz
,
M.
Hansen
,
K. C.
Iarosz
,
I. L.
Caldas
,
A. M.
Batista
, and
J.
Kurths
, “
Emergence of neuronal synchronisation in coupled areas
,”
Front. Hum. Neurosci.
15
,
663408
(
2021
).
5.
M. C.
Soriano
,
J.
García-Ojalvo
,
C. R.
Mirasso
, and
I.
Fischer
, “
Complex photonics: Dynamics and applications of delay-coupled semiconductors lasers
,”
Rev. Mod. Phys.
85
,
421
(
2013
).
6.
Y.
Han
,
S.
Xiang
, and
L.
Zhang
, “
Cluster synchronization in mutually-coupled semiconductor laser networks with different topologies
,”
Opt. Commun.
445
,
262
267
(
2019
).
7.
A.
Schnitzler
and
J.
Gross
, “
Normal and pathological oscillatory communication in the brain
,”
Nat. Rev. Neurosci.
6
,
285
296
(
2005
).
8.
A. E.
Motter
,
S. A.
Myers
,
M.
Anghel
, and
T.
Nishikawa
, “
Spontaneous synchrony in power-grid networks
,”
Nat. Phys.
9
,
191
197
(
2013
).
9.
L. M.
Pecora
,
F.
Sorrentino
,
A. M.
Hagerstrom
,
T. E.
Murphy
, and
R.
Roy
, “
Cluster synchronization and isolated desynchronization in complex networks with symmetries
,”
Nat. Commun.
5
,
4079
(
2014
).
10.
B. D.
MacArthur
,
R. J.
Sánchez-García
, and
J. W.
Anderson
, “
Symmetry in complex networks
,”
Discrete Appl. Math.
156
,
3525
3531
(
2008
).
11.
P. S.
Skardal
, “
Symmetry and symmetry breaking in coupled oscillator communities
,”
Eur. Phys. J. B
92
,
1
11
(
2019
).
12.
P.
Chossat
and
R.
Lauterbach
,
Methods in Equivariant Bifurcations and Dynamical Systems
(
World Scientific Publishing Company
,
2000
), Vol. 15.
13.
M.
Golubitsky
,
I.
Stewart
, and
D. G.
Schaeffer
,
Singularities and Groups in Bifurcation Theory: Volume II
(
Springer Science & Business Media
,
2012
), Vol. 69.
14.
B.
Pietras
and
A.
Daffertshofer
, “
Network dynamics of coupled oscillators and phase reduction techniques
,”
Phys. Rep.
819
,
1
105
(
2019
).
15.
V.
Nicosia
,
M.
Valencia
,
M.
Chavez
,
A.
Díaz-Guilera
, and
V.
Latora
, “
Remote synchronization reveals network symmetries and functional modules
,”
Phys. Rev. Lett.
110
,
174102
(
2013
).
16.
I.
Schneider
, “
Delayed feedback control of three diffusively coupled Stuart–Landau oscillators: A case study in equivariant Hopf bifurcation
,”
Philos. Trans. R. Soc. A
371
,
20120472
(
2013
).
17.
I.
Schneider
and
M.
Bosewitz
, “
Eliminating restrictions of time-delayed feedback control using equivariance
,”
Discrete Contin. Dyn. Syst. A
36
,
451
467
(
2016
).
18.
J.
Collins
and
I.
Stewart
, “
Hexapodal gaits and coupled nonlinear oscillator models
,”
Biol. Cybern.
68
,
287
298
(
1993
).
19.
J. J.
Collins
and
I. N.
Stewart
, “
Coupled nonlinear oscillators and the symmetries of animal gaits
,”
J. Nonlinear Sci.
3
,
349
392
(
1993
).
20.
F. M.
Atay
,
Complex Time-Delay Systems: Theory and Applications
(
Springer
,
2010
).
21.
A.
Zakharova
,
I.
Schneider
,
Y.
Kyrychko
,
K.
Blyuss
,
A.
Koseska
,
B.
Fiedler
, and
E.
Schöll
, “
Time delay control of symmetry-breaking primary and secondary oscillation death
,”
Europhys. Lett.
104
,
50004
(
2013
).
22.
T.
Erneux
,
J.
Javaloyes
,
M.
Wolfrum
, and
S.
Yanchuk
, “Introduction to focus issue: Time-delay dynamics,”
Chaos
27
,
114201
(
2017
).
23.
A.
Otto
,
W.
Just
, and
G.
Radons
, “
Nonlinear dynamics of delay systems: An overview
,”
Philos. Trans. R. Soc. A
377
,
20180389
(
2019
).
24.
R. C.
Calleja
,
A.
Humphries
, and
B.
Krauskopf
, “
Resonance phenomena in a scalar delay differential equation with two state-dependent delays
,”
SIAM J. Appl. Dyn. Syst.
16
,
1474
1513
(
2017
).
25.
S.
Yanchuk
,
K. R.
Schneider
, and
L.
Recke
, “
Dynamics of two mutually coupled semiconductor lasers: Instantaneous coupling limit
,”
Phys. Rev. E
69
,
056221
(
2004
).
26.
H.
Erzgräber
,
B.
Krauskopf
, and
D.
Lenstra
, “
Compound laser modes of mutually delay-coupled lasers
,”
SIAM J. Appl. Dyn. Syst.
5
,
30
65
(
2006
).
27.
E.
Clerkin
,
S.
O’Brien
, and
A.
Amann
, “
Multistabilities and symmetry-broken one-color and two-color states in closely coupled single-mode lasers
,”
Phys. Rev. E
89
,
032919
(
2014
).
28.
F.
Della Rossa
,
L.
Pecora
,
K.
Blaha
,
A.
Shirin
,
I.
Klickstein
, and
F.
Sorrentino
, “
Symmetries and cluster synchronization in multilayer networks
,”
Nat. Commun.
11
,
3179
(
2020
).
29.
K.
Blaha
,
R. J.
Burrus
,
J. L.
Orozco-Mora
,
E.
Ruiz-Beltrán
,
A. B.
Siddique
,
V.
Hatamipour
, and
F.
Sorrentino
, “
Symmetry effects on naturally arising chimera states in mechanical oscillator networks
,”
Chaos
26
,
116307
(
2016
).
30.
C. G.
Rusin
,
H.
Kori
,
I. Z.
Kiss
, and
J. L.
Hudson
, “
Synchronization engineering: Tuning the phase relationship between dissimilar oscillators using nonlinear feedback
,”
Philos. Trans. R. Soc. A
368
,
2189
2204
(
2010
).
31.
Z.
Balanov
,
W.
Krawcewicz
, and
H.
Steinlein
,
Applied Equivariant Degree
(
American Institute of Mathematical Sciences
,
Springfield, MO
,
2006
), Vol. 1.
32.
K.
Engelborghs
,
T.
Luzyanina
, and
G.
Samaey
, “DDE-BIFTOOL: A MATLAB package for bifurcation analysis of delay differential equations,” TW Report, Vol. 305, 2000.
33.
J.
Sieber
,
K.
Engelborghs
,
T.
Luzyanina
,
G.
Samaey
, and
D.
Roose
, “DDE-BIFTOOL manual—Bifurcation analysis of delay differential equations,” arXiv:1406.7144 (2014).
34.
K.
Pyragas
, “
Continuous control of chaos by self-controlling feedback
,”
Phys. Lett. A
170
,
421
428
(
1992
).
35.
Y.
Zhai
,
I. Z.
Kiss
, and
J. L.
Hudson
, “
Control of complex dynamics with time-delayed feedback in populations of chemical oscillators: Desynchronization and clustering
,”
Ind. Eng. Chem. Res.
47
,
3502
3514
(
2008
).