Transitional cluster dynamics in a model for delay-coupled chemical oscillators

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.⁷ 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.¹⁴ Knowledge of the symmetries can also be used to discover possible cluster synchronization patterns.⁹ 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.²⁴ 

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.¹ 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.³⁰ 

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.

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.

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)= Ω mt+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.

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.

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.¹ 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.

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π$⁠.

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$⁠.

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.

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,¹³ 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).

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.

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.

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.

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)].

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 $τ$⁠.

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.

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.

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.³⁴ 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/Hg₂SO₄/K₂SO₄ (sat) reference electrode with a 3M H₂SO₄ 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.³⁵ The startup or shutdown of one oscillator does not alter the behavior of the others, and oscillator dynamics are independent when uncoupled.

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$⁠.

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.

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.