Bifurcations of clusters and collective oscillations in networks of bistable units

Networks of interacting bistable units are an important class of dynamical systems that are relevant in diverse physical contexts. While coupled units with oscillatory or chaotic dynamics have been widely considered (see, e.g., the textbook in Ref. 2), multi-component systems in which each unit has a single degree of freedom have received much less attention. Zanette³ considered globally coupled bistable elements, which interacted linearly through their common mean. The global coupling was diffusive in the sense that it vanished for uniform states. For random initial condition and a sufficiently large coupling strength, the global interaction led to a coherent motion of the entire ensemble toward one of the two stable states. The impact of additive noise to the bistable units in related network configurations was analyzed in Refs. 4–7, while Refs. 8–10 considered the collective response of bistable units that were connected through different coupling topologies. A recent study also examined the transition from a local to global coupling topology, focusing on the mathematical perspective.¹¹ 

Kouvaris and Mikhailov¹² investigated the effect of a global feedback on the dynamics of networks of diffusively coupled bistable units. The global feedback was chosen such that it altered the excitation threshold, i.e., the position of the saddle point, of the local one-dimensional reaction dynamics. In this case, localized stationary activation patterns form, the size of which can be adjusted by varying the feedback strength. The impact of a global time-delayed feedback and additive noise on an ensemble of bistable units was examined in Refs. 4 and 13, while in Ref. 14, the impact of heterogeneities on the synchronization state of globally coupled bistable electronic circuits is discussed.

Here, we consider bistable units that interact through a macroscopic observable, which is forced to take a fixed value as a global constraint. This observable can be the sum of the intrinsic state variables of the individual units or it can be the sum of a function of the individual state variables. The global constraint is enforced by allowing a parameter that controls the state of an individual unit on the equilibrium branch to adapt; thus, the global constraint can be seen as a bifurcation parameter of the coupled system. Indeed, such coupling between bistable units is realized in many physical systems: These include Li-ion batteries with phase-transition cathodes^15,16 and bistable electrochemical reactions, such as the CO oxidation on an array of Pt electrodes.^17–19 In the case of Li-ion batteries, the cathodes consist of billions of nano-particles that can be considered bistable units. Each nanoparticle can be in a Li-rich or Li-poor state depending on the chemical potential or equivalently the voltage. Thus, the battery can be seen as interacting bistable units coupled through a global constraint: When charging or discharging the battery slowly, a constraint is set to the time evolution of the total charge, while the voltage adjusts accordingly. Similarly, the CO oxidation on Pt is a prototypical electrocatalytic reaction, which exhibits bistable reaction rates and thus a bistable current–voltage characteristic. In technological applications, electrodes consist of billions of catalytically active nanoparticles on a nonreactive support that interact globally when a set current is passed through the electrode. A setup that enables the measurement of the state of an individual bistable component is an array of Pt electrodes. Such measurements reveal key dynamical properties: The bistable units activate sequentially upon a slow parameter ramp, they form clusters with most electrodes in the two stable states and at most one electrode on the unstable state of the individual electrode, and they may exhibit collective oscillations of the entire ensemble.^17–19 We recently derived a general necessary condition when a system of bistable units subject to a global constraint can become unstable in a Hopf bifurcation¹ and validated this condition in a simplified CO-electrooxidation model.

In this paper, we analyze this CO-electrooxidation model and elucidate its bifurcation properties. The behavior of the coupled units is particularly interesting when the value of the constraint is in the transitional regime between the two outer “active” and “passive” states of an uncoupled individual unit. As a globally and identically coupled system of identical units, the bifurcations are constrained by the symmetry properties of the system: The dynamical equations are $SN$-symmetric (equivariant), where $SN$ denotes the group of permutations of $N$ symbols that acts by permuting component indices; see Ref. 20 for a general introduction to dynamical systems with symmetry. The presence of symmetries gives rise to dynamically invariant subspaces, which correspond, for example, to cluster configurations where the state of some components coincides. Symmetric systems can exhibit bifurcation behavior that is nongeneric for a general dynamical system as their bifurcation behavior is constrained by the symmetry.^20,21 We will consider symmetry-breaking bifurcations where the configuration of all units being synchronized (the one-cluster configuration) loses stability. Indeed, since the state space of an individual unit is one-dimensional for the system we consider, collective oscillations cannot occur if all units are synchronized. Thus, the emergence of collective oscillations necessarily requires symmetry-breaking bifurcations away from full synchrony.

In our bifurcation analysis, we concentrate specifically on the structure imposed by the invariant subspaces due to symmetry. Here the one-, two-, and three-cluster subspaces are crucially important as they contain all equilibria of the system. We give a detailed bifurcation analysis of the steady states in these subspaces and highlight bifurcations that give rise to and stabilize limit cycles as they induce collective oscillations observed in experiments. Specifically, the paper is organized as follows. In Sec. II, we summarize the model equations. In Sec. III, we consider the one-dimensional, fully synchronized dynamics where the states of all units take the same value and form one cluster. In Sec. IV, we discuss symmetry-breaking bifurcations away from full synchrony to two-cluster equilibria and bifurcations of the two-cluster equilibria; these can give rise to collective oscillations (that are transversely unstable for most cluster sizes). In Sec. V, we outline bifurcations of three-cluster equilibria and indicate that the transversely unstable collective oscillations within the two-cluster subspace can be stabilized in a transcritical bifurcation. We conclude with some remarks in Sec. VI.

We consider the following simplified model for CO-electrooxidation on Pt microelectrodes. Its derivation from an established, more detailed model version can be found in the supplement of Ref. 1. The intrinsic state $xk$ of unit $k$ evolves according to

where $y$ imposes a global constraint and $a$ and $b$ are positive constants; throughout this paper, we fix $a=0.05$⁠, $b=0.01$ as in Ref. 1. In terms of the CO dynamics, $xk$ is the CO-covered fraction of the surface and $u$ the electrode potential. The macroscopic observable $ytot$ is the total CO oxidation current, and $y$ denotes the mean current per electrode (or nanoparticle). In our setting, $y$ is the bifurcation parameter, while $u$ adjusts such that the system can attain the preset mean value. Note that this type of coupling corresponds to a control that fixes the value of a function of the mean state of the ensemble and allows the “natural control parameter”—in our case $u$—to change. This situation is different from the situation considered in Ref. 3, where the main bifurcation parameter of the bistable unit is also the bifurcation parameter of the globally coupled ensemble and thus set to a constant value.

First, consider a single uncoupled unit, $N=1$⁠, with state $x1=x$⁠; the corresponding equilibrium dynamics $x˙=0$ are shown in Fig. 1. When $u$ is treated as a parameter, $x$ is bistable in $u$ with equilibrium branch as in Fig. 1(a); there are two saddle-node (SN) bifurcations that bound the region of bistability. For ease of notation, we subdivide the branch of equilibria into three segments specified by the points where these SN bifurcations occur and the branch turns over: The active segment (“a”), where the current is high, the middle (“m”) segment, and the passive segment (“p”), where the current is low (see Fig. 1). We denote the value of $u$ corresponding to the transitions between the passive and middle segment with $up→m$ and the transition between the middle and active segment with $um→a$⁠. The transitions between segments play an explicit role in the bifurcation analysis as we will discuss below. For the equilibrium values of $x$⁠, there is a single-valued, but non-monotonic relationship between $u$ and $y$⁠, depicted in Fig. 1(b). By contrast, the value of $x$ at equilibrium depends monotonically on $y$ [Fig. 1(c)].

The dynamical equations (1) can be rearranged to the explicit form

with the electrode potential

which is determined by the intrinsic states $xk$ and the bifurcation parameter $y$⁠. If we define the two functions

the dynamical equations take the form

Configurations where the states of the units form three clusters are particularly interesting because all equilibria of Eq. (5) necessarily lie in a three-cluster subspace. To see this, consider the fixed points $xk∗$ of Eq. (1) with $u=u∗$⁠, which satisfy

Since this cubic polynomial has at most three distinct roots, the coefficients of any equilibrium of (5) take at most three distinct values and thus form a cluster configuration.

As a globally coupled network of

Suppose that the units form $M$ clusters such that the $k$th cluster contains $Nk$ units with identical state $ξk$⁠. Evidently, we should have $N=∑j=1MNj$⁠. Define the relative cluster size $nk=NkN$ and write $C(n1,…,nM)$ for the cluster subspace where the first $N1$ oscillators form cluster $1$⁠, the next $N2$ oscillators cluster $2$⁠, and so on. The effective dynamics on $C(n1,…,nM)$—and by symmetry any other cluster configuration with these cluster sizes—are $M$-dimensional, and the state $ξk$ of cluster $k∈{1,…,M}$ evolves according to

While $nk$ only takes finitely many values for finite networks, in the limit of infinitely many units, $N→∞$⁠, we can see $nk$ as a continuous parameter.

The simplest cluster configuration $C(1)$ is a single cluster; that is, all units are synchronized with $xk=ξ1=ξ$⁠. With $P,Q$ as above, the synchronized dynamics of (5) are given by

Equation (8) has only one single equilibrium point

for any $y∈(0,1)$⁠, which is identical to the curve of an individual element depicted in Fig. 1. Restricted to the synchronized subspace, the system is monostable for all parameter values; i.e., for each $y$⁠, there is exactly one single equilibrium. (If $u$ is given instead of $y$⁠, however, there can be three equilibria.) Linearizing Eq. (8) at the equilibrium $ξ∗$ yields

and thus, $ξ∗$ is always stable within $C(1)$⁠.

The stability of $ξ∗$ as an equilibrium of the full system (5) depends on the $N−1$ eigenvalues transverse to $C(1)$⁠. Due to the symmetry, all transversal eigenvalues are the same; they are all negative on the “a” and “p” segments and positive on the “m” segment. This implies that the equilibrium is stable on the “a” and “p” segments and unstable on the “m” segment. The change of stability of the equilibrium $ξ∗$ happens in $SN$-equivariant transcritical bifurcations at $up→m$ and $um→a$ where all eigenvalues transverse to $C(1)$ pass through zero simultaneously. Note that at these transcritical bifurcations, an uncoupled individual bistable unit undergoes a SN bifurcation if $u$ is considered a bifurcation parameter [cf. Eq. (1a)]; that is, the bifurcation happens where all synchronized units jointly transition from the passive to middle segment of the branch or from the middle to active segment of the branch.

Now, we consider two-cluster configurations $C(n1,n2)$ with relative sizes $n1$ and $n2$⁠, respectively. The two-dimensional dynamics on $C(n1,n2)$ are given by

The intersection of all possible two-cluster subspaces contains the fully synchronized subspace ${ξ1=ξ2}$⁠. Note that for $n1$ and $n2$ fixed, the system has a parameter symmetry

which corresponds to exchanging the clusters. Since $n1+n2=1$⁠, the parameter $n1$ fully determines the system, and the parameter symmetry can also be written as

Thus, without loss of generality, we may assume $n1≥n2$⁠. If the clusters are of equal size, $n1=n2$⁠, the parameter symmetry yields a system symmetry $(ξ1,ξ2)↦(ξ2,ξ1)$ on the cluster subspace.

The branch of fully symmetric equilibria in $C(1)$ interacts with branches of two-cluster solutions in $C(n1,n2)$ in transcritical bifurcations. Note that we say that two (distinct) branches of equilibria or limit cycles interact if they are involved in the same bifurcation; for example, two branches of equilibria interact in a simple transcritical bifurcation. All branches and their bifurcations described throughout the paper were calculated numerically using AUTO-07P.²² In the following, we will first focus on the specific relative cluster sizes $n1=0.8$⁠, $n2=0.2$ to illustrate the bifurcation behavior. Figure 2 shows the equilibria branches in the $(u,y)$-plane as the bifurcation parameter $y$ is varied. The black line shows the fully symmetric equilibrium $ξ==(ξ∗,ξ∗)$—with $ξ∗$ as in Sec. III—within $C(1)⊂C(n1,n2)$⁠; this branch is the same as the one shown in Fig. 1(b).

The locations where the fully symmetric equilibrium solution undergoes transcritical bifurcations are indicated by yellow dots. At these bifurcation points, the symmetric branch interacts with the branch of two-cluster equilibria (dark blue line). It is useful to subdivide the branch of two-cluster equilibria into segments: We write “aa” if both clusters are on the active segment, “am” if the first cluster is on the active and the second one on the middle segment, etc. This leads to the labeling shown in Fig. 2, and the boundaries of the segments are given by the local maxima and minima corresponding to $u∈{up→m,um→a}$⁠. This labeling yields some intuition of what happens at the transcritical bifurcations: At the first transcritical bifurcations, the two-cluster equilibrium transitions from “pm” to “mp” and the clusters switch roles as the majority cluster (⁠$n1=0.8$⁠) switches from “p” to “m” and the minority cluster (⁠$n2=0.2$⁠) from “m” to “p.” Similarly, at the second transcritical bifurcation, there is a transition from “ma” to “am.” Following the two-cluster equilibrium branch further, it folds over twice in SN bifurcations (green dots) and undergoes a Hopf bifurcation at $y≈0.8$ (pink dot). This gives rise to a branch of limit cycles (light blue) that is stable within $C(n1,n2)$ and ends in a homoclinic bifurcation close to the transcritical bifurcation at $y≈0.9$⁠.

There is actually a series of bifurcations close to the second transcritical bifurcation as the parameter $y$ is varied as shown in Fig. 3. When the Hopf bifurcation occurs at $y≈0.8$⁠, there are three equilibria in $C(n1,n2)$⁠, the symmetric equilibrium $ξ=$ in $C(1)⊂C(n1,n2)$ as well as nonsymmetric equilibria $ξ+=(ξ1+,ξ2+)$ with $ξ1+<ξ2+$ and $ξ−=(ξ1−,ξ2−)$ with $ξ2−<ξ1−$⁠; the former gives rise to the stable limit cycle. At the transcritical bifurcation, $ξ−$ interacts with $ξ=$ and satisfies $ξ1−<ξ2−$ after the bifurcation point. This equilibrium then collides with the limit cycle in a homoclinic bifurcation at $y≈0.91$⁠. Finally, the two equilibria $ξ−$ and $ξ+$ meet in the fold where the nonsymmetric branch in Fig. 2 folds over itself. Note that the transcritical bifurcation must happen before the homoclinic bifurcation to allow for this bifurcation scenario: It cannot occur on the “ma” segment before the transcritical bifurcation but has to occur on the “am” segment (but before $ξ−$ and $ξ+$ merge in the SN bifurcation).

Considering transverse stability beyond $C(n1,n2)$⁠, we realize that the branch of two-cluster equilibria undergoes transcritical bifurcations whenever $u∈{up→m,um→a}$⁠. Hence, there are six in total, the two indicated by yellow dots in Fig. 2, where all bifurcating branches lie within $C(n1,n2)$⁠, and four additional ones where the two-cluster equilibria interact with three-cluster equilibria; we will discuss these in more detail in Sec. V.

To elucidate the dynamics and bifurcations for other relative cluster sizes, we continued the bifurcations in two parameters, $y$ and $n1$⁠. The resulting two-parameter bifurcation lines are shown in Fig. 4. Note that the symmetry of the bifurcation diagram reflects the parameter symmetry (12). The transcritical bifurcation lines are horizontal since the bifurcation condition does not depend on the relative cluster sizes: Consider equilibria on $C(n1,n2)$⁠, that is, solutions of Eq. (6) with two distinct real roots $ξ1∗≠ξ2∗$ (a nonsymmetric equilibrium on $C(n1,n2)$⁠). Now suppose that these roots coincide at a bifurcation for $u=u∗$⁠, defined in (3), which means that (6) has a triple root $ξ1∗=ξ2∗=ξ∗$ (an equilibrium on $C(1)⊂C(n1,n2)$⁠). Then, the associated value of $y$ is

independent of $n1$⁠. Furthermore, for any given relative cluster sizes $(n1,n2)$⁠, there is a branch of equilibria in $C(n1,n2)$⁠. These branches arise in saddle-node bifurcations and interact with the one-cluster equilibrium in the transcritical bifurcation. The more asymmetric the clusters are, the larger the distance of the saddle-node bifurcation to the transcritical bifurcation in the bifurcation parameter $y$ is. The Hopf and homoclinic bifurcations of the two-cluster states described above extend to a large range of cluster sizes from the most asymmetric distribution $n1=1$ (or, equivalently, $n1=0$⁠) to $n1≈0.4$ (respectively $n1≈0.6$⁠). They meet in a codimension-two Bogdanov–Takens bifurcation point at $(y,n1)≈(0.907,0.34)$ [$(y,n1)≈(0.907,0.66)$⁠, respectively]. For increasing (respectively decreasing) cluster size $n1$⁠, the Hopf bifurcation curve approaches the saddle-node bifurcation curve again. In this region of parameter space, the limit cycle solution grows extremely quickly as in a Canard phenomenon²³ before it is destroyed in a branch of homoclinic bifurcations.

In the full system (5), all these bifurcations happen simultaneously in a large number of invariant subspaces. On the one hand, for fixed $(n1,n2)$⁠, the bifurcations happen simultaneously in the symmetric copies of $C(n1,n2)$⁠. On the other hand, recall that $(n1,n2)$ parameterize invariant subspaces of different cluster sizes. The transcritical bifurcations for all the different $(n1,n2)$ happen at the same value of $y$ given by (13). Thus, the symmetric equilibrium undergoes a transcritical bifurcation simultaneously in all these invariant subspaces as one would expect in our symmetric setting. In the case of symmetric clusters, $n1=n2=0.5$⁠, these degenerate to a symmetric pitchfork bifurcation, in line with the reflection symmetry $(ξ1,ξ2)↦(ξ2,ξ1)$ on $C(n1,n2)$⁠; cf. Fig. 4.

Now, consider the three-cluster subspace $C(n1,n2,n3)$ with three clusters of relative sizes $n1+n2+n3=1$ whose states $ξ1,ξ2,ξ3$ evolve according to

The invariant three-cluster space $C(n1,n2,n3)$ contains the invariant two-cluster spaces $C(n1+n2,n3)$ as ${ξ1=ξ2}$⁠, $C(n1,n2+n3)$ as ${ξ2=ξ3}$⁠, and $C(n1+n3,n2)$ as ${ξ3=ξ1}$⁠, on which the dynamics are given by (11). As above, any permutation of the three clusters yields a parameter symmetry. If $n1=n2$⁠, $n2=n3$⁠, or $n3=n1$⁠, then (14) has a symmetry with respect to the transpositions that swap the corresponding clusters. Thus, without loss of generality, we may assume $n1≥n2≥n3$⁠. If $n1=n2=n3=13$ (i.e., all clusters are of equal size), then Eq. (14) is equivariant with respect to the full permutation of three elements.

Since any equilibrium is constrained to a three-cluster configuration due to condition (6), any branch of steady states is contained in a cluster subspace $C(n1,n2,n3)$⁠. Figure 5 shows the equilibrium branches in $C(n1,n2,n3)$ for relative cluster sizes $n1=0.8$⁠, $n2=0.16$⁠, and $n3=0.04$⁠. We have $C(1)⊂C(n1,n2,n3)$⁠; therefore, the branch of symmetric equilibria in Fig. 1 appears (black line). Note that $n2+n3=0.2$⁠, and thus, $C(0.8,0.2)⊂C(0.8,0.16,0.04)$ for ${ξ2=ξ3}$ and the equilibrium branches shown in Fig. 2 reappear in Fig. 5 (blue line). In addition, $C(0.8,0.16,0.04)$ also contains the two-cluster subspaces ${ξ1=ξ3}$ and ${ξ1=ξ2}$⁠; these correspond to $C(0.84,0.16)$ and $C(0.96,0.04)$⁠, respectively. The steady state bifurcations can be read in Fig. 4, and the corresponding branches of equilibria are depicted in Fig. 5 as purple and teal lines, respectively. Finally, there is one branch of three-cluster equilibria. Taken together, for $u∈[up→m,um→a]$⁠, there are $33=27$ equilibria because each $xj$ can take one of three values while $x˙j=0$⁠, $u˙=0$⁠).

These equilibria now (ex)change their stability properties as they bifurcate. In addition to the transcritical bifurcations of the fully synchronized equilibrium, the two-cluster equilibria also undergo transcritical bifurcations in which three-cluster equilibria are involved. More specifically, each of the three two-cluster subspaces contained in $C(0.8,0.16,0.04)$⁠, that is, $ξ1=ξ2$⁠, $ξ2=ξ3$⁠, $ξ1=ξ3$⁠, contains a transcritical bifurcation for each of the two double root of the equilibrium equation (6). This yields six transcritical bifurcations in addition to the transcritical bifurcations of the fully synchronized equilibrium that already appeared in Fig. 2. These bifurcations link to the physical interpretation as in Sec. IV: At each transcritical bifurcation, the value of $u$ is identical with one of the values at which the individual unit undergoes a saddle-node bifurcation when $u$ is considered a bifurcation parameter. Indeed, it can be shown that when the transcritical bifurcation involves an equilibrium with the larger of the two $u$-values, i.e., $u=up→m$⁠, the components of the three-cluster state that are on the middle and passive branch segment switch the segment, while the component on the active segment remains unchanged. For example, if $ξjx$⁠, $j∈{1,2,3}$ denotes the $j$th component of a three cluster that is on the segment $x∈{p,m,a}$⁠, then an equilibrium $ξ∗=(ξ1a,ξ2m,ξ3p)$ on the “amp” segment transitions to $ξ∗=(ξ1a,ξ2p,ξ3m)$ on the “apm” segment. Correspondingly, at the lower value of $u$⁠, the components with values on the middle and active segment switch segment, while the one on the passive segment remains there, such as $ξ∗=(ξ1a,ξ2m,ξ3p)$ transitions to $ξ∗=(ξ1m,ξ2a,ξ3p)$⁠.

We can evaluate the values at which the transcritical bifurcations occur explicitly as in Sec. IV; in general, they depend on the cluster sizes $n1$⁠, $n2$⁠, $n3$⁠. Consider equilibria on $C(n1,n2,n3)$⁠, that is, solutions of Eq. (6) with three distinct real roots $ξ1∗≠ξ2∗≠ξ3∗≠ξ1∗$ (a nonsymmetric equilibrium on $C(n1,n2,n3)$⁠). Now, suppose two of these roots coincide at a bifurcation for $u=u∗$⁠, defined in (3) so that Eq. (6) has one real double root $ξd$ and one real single root $ξs$⁠. Without loss of generality, we assume that $ξ1∗=ξs$⁠, $ξ2∗=ξ3∗=ξd$ on $C(n1,n2+n3)⊂C(n1,n2,n3)$⁠. Then, the associated value of $y$ is

which is linear in $n1$⁠. The cluster at $ξs$ is what distinguishes this from the previous lower-dimensional case, where this “bystander” was not present. By symmetry, the condition for the transcritical bifurcation in the other two-cluster subspaces depends linearly on $n2$ and $n3$⁠, respectively. Thus, the transcritical bifurcations are given by lines as the cluster sizes are varied, shown in Fig. 6. In this figure, we see that eight transcritical bifurcations exist for every $n1$⁠. The value $n1=0.8$ corresponds to Fig. 5. We see the two transcritical bifurcations in $C(n1,n2+n3)$ following Eq. (15). Furthermore, $n2$ and $n3$ are chosen to be linearly dependent on $n1$ (see the caption of Fig. 6); therefore, the corresponding bifurcations lie on straight lines, too.

For the full system with $N$ identical units, this means that for each $ℓ∈{2,…,N}$ (i.e., the combined size of cluster 2 and 3, which fall together at the bifurcation), there are $(Nℓ)∑k=1ℓ−112(ℓk)$ transcritical bifurcations at $y=u∗[Q(ξs)−Q(ξd)](N−ℓ)/N+u∗Q(ξd)$⁠. These bifurcation points lie in $(Nℓ)$ different subspaces, depending on which $xk$ are at $ξs$ and $ξd$⁠. If you consider bifurcations as identical when they can be switched between by index permutation, then there is only one such bifurcation for each combination of $ℓ∈{2,…,N}$ (i.e., combined size of cluster 2 and 3) and $k∈{1,…,ℓ−1}$ (i.e., size of cluster 2). Furthermore, many of these transcritical bifurcations lie on the same points in the parameter space and in the phase space such that in total, there are actually just $N−1$ multi-branch bifurcation-points at every $u∗$⁠, i.e., one for each possible size of cluster 1 (the cluster at $ξs$⁠) and one for each $ℓ∈{2,…,N}$⁠. The number of branches is $2+2∑k=1ℓ−112(ℓk)$⁠. These bifurcations are invariant under intra-cluster index-permutations and are very similar to the equivariant bifurcations we saw in lower dimensions. Here, however, they can occur at various different values of $y$⁠. In fact, for $N→∞$⁠, they occur at each of the two $u∗$-values in the entire $y$-intervals from the transcritical bifurcations almost to the respective fully symmetric state.

As seen above for the two-cluster subspaces, also, three-cluster equilibria might undergo a Hopf bifurcation. In Fig. 6, the locations of the Hopf bifurcations of the three-cluster and the different two-cluster equilibria are depicted in the $(y,n1)$-plane. The labels of the different curves indicate the segments of the three components of the equilibria from which the limit cycles bifurcated. All of them contain at least one cluster on the middle segment. As shown in Ref. 1, in the full system, any steady state with more than one element on the middle branch is unstable. Hence, even if the Hopf bifurcations are supercritical and the limit cycle branches off equilibria that are stable in the cluster subspaces, they are unstable in the full system. The limit cycles that emerge in a Hopf bifurcation will, therefore, not be observable in systems containing a large number of individual elements. However, as we will demonstrate in Subsection V B, the three-cluster limit cycles might interact with a two-cluster limit cycle in a transcritical bifurcation, thereby stabilizing the limit cycle in the two-cluster subspace.

Consider the phase portraits in the three-cluster subspace $C(0.8,0.16,0.04)$ depicted in Fig. 7. At $y=0.78625$⁠, there are two fixed points and two limit cycles, which were both born at slightly smaller values of $y$ in supercritical Hopf bifurcations. One limit cycle and one fixed point are three-cluster limit sets, while the other two limit sets lie in the invariant $C(n1,n2+n3)$ subspace; i.e., $ξ2=ξ3$⁠. The two-cluster limit cycle is stable within the two-cluster subspace but unstable in the direction perpendicular to it, while the three-cluster limit cycle is stable within the invariant three-cluster subspace. As $y$ is increased to $y=0.79$⁠, we see that the limit cycles have increased in size and have moved closer together. At $y≈0.79625$⁠, there is a transcritical bifurcation of the limit cycles; i.e., they coincide on the two-cluster subspace $C(n1,n2+n3)$⁠. Then, at $y=0.8$⁠, the three-cluster limit cycle is on the other side of the $C(n1,n2+n3)$ plane. In the transcritical bifurcation, the limit cycles have exchanged their stabilities perpendicular to the two-cluster subspace. Hence, the two-cluster limit cycle is now stable within the three-cluster subspace.

Analogous to the equivariant transcritical bifurcation of equilibria, in the full system, the transcritical bifurcation of limit cycles is an equivariant transcritical bifurcation. As above, the cluster sizes are not true bifurcation parameters but rather enumerate different invariant subspaces in the same system. Thus, the transcritical bifurcation takes places simultaneously in all the different subspaces, which are related to each other by index permutation, and in the corresponding subspaces, which correspond to different cluster sizes $n2$ and $n3$ at fixed $n2+n3=1−n1$⁠. In other words, in the full system with $N$ identical components and sufficiently large $N$ at $y=0.79625$⁠, multiple limit cycles from different three-cluster subspaces meet in the same two-cluster limit cycle, each changing the stability of a different Floquet exponent of the two-cluster limit cycle. Since all transversal eigenvalues (i.e., intra-cluster eigenvalues) of the Jacobian matrix are degenerate due to symmetry,¹ it actually affects all of them in the same way. As a consequence, in the full system, the equivariant transcritical bifurcation of two-cluster limit cycles stabilizes the latter and macroscopic two-cluster oscillations should be observable in arbitrarily large systems, i.e., even in the thermodynamic limit. We verified this conclusion numerically for a system with $N=20$ units and the other parameters as in Fig. 7.

In the CO-electrooxidation model we analyzed in this paper, the fully synchronized one-cluster equilibrium loses stability in an $SN$-equivariant transcritical bifurcation where $N−1$ real eigenvalues change sign simultaneously. Such bifurcations can arise, for example, as a perturbation of a pitchfork bifurcation²¹ where the invariant subspace persists. Indeed, the pitchfork bifurcation for identical cluster sizes (⁠$n1=n2$⁠) perturbs to a pair of saddle-node and transcritical bifurcation for nonidentical cluster sizes. While the symmetric transcritical bifurcation happens simultaneously for all cluster sizes, the parameter values of the saddle-node bifurcation are distinct for different relative cluster sizes as seen in Fig. 4: The most asymmetric two-cluster state is generated first, then increasingly symmetric ones, and finally the symmetric $(12,12)$-cluster in a pitchfork bifurcation. The symmetry-breaking bifurcations from two- to three-cluster equilibria occur in a similar fashion, just that there is an additional cluster as an “bystander.” More generally, if the phase-space geometry allows, $SN$-equivariant transcritical bifurcation can also have a global flavor.²⁴ 

Equations that describe globally coupled identical one-dimensional units also directly relate to more general contexts, such as coupled oscillators. If a general $SN$-equivariant system undergoes a $SN$-equivariant transcritical bifurcation, one can describe the dynamics on the corresponding center manifold through a suitable normal form. This approach has been used to describe the dynamics of identical Stuart–Landau oscillators at a symmetry-breaking bifurcation.^25–27 Truncating to appropriate order, one obtains a cubic $SN$-equivariant system of coupled one-dimensional variables $xk$^20,26–30 that describe the dynamics on the center manifold. In the limit of two large clusters and one cluster of vanishing size, the dynamics is characterized by a web of heteroclinic orbits between the two-cluster equilibria.²⁷ These cubic vector fields show the same bifurcation behavior close to the symmetry-breaking bifurcations of the system analyzed here, including the transcritical and saddle-node bifurcations, as well as the secondary bifurcations where the two-cluster solutions lose stability; cf. Refs. 25 and 26. A system of coupled one-dimensional units with a cubic vector field, however, is unable to capture some of the more intricate secondary bifurcation behavior, such as Hopf and homoclinic bifurcations, described here. These require the inclusion of higher-order terms, such as the coupling of (1) of our CO-electrooxidation model that is beyond cubic order, or a higher-order approximation on the center manifold.

The dynamical features described in this paper arise in a mathematical model that captures the experimental control of a physical system. From this perspective, the bifurcations discussed are common to many real-world systems. Yet, we made one simplification, which will never be strictly fulfilled in a physical system, namely, that all bistable units are identical. While the permutational symmetry is broken if the units are slightly heterogeneous, hyperbolic equilibria and limit cycles—which are the main focus of our analysis—will persist. The exact behavior of the branches at the (symmetric) bifurcation points will of course differ as these bifurcations break up into a series of generic bifurcation points; how exactly depends on the heterogeneity. However, the qualitative behavior that is relevant for real-world systems is preserved: In all ensemble steady states, each of the bistable units will be on one of its three branch segments. We can still define a cluster as a group of units that are on the same segment, and a unit transits between the segments at its turning points of $u$⁠, i.e., $up→m$ or $um→a$⁠. However, for the individual units, the turning points do not necessarily occur at identical values of $u$ as the branches deform. We verified that Hopf and homoclinic bifurcations still occur in a model of heterogeneous bistable units. Consequently, our bifurcation analysis facilitates a qualitative understanding of experimental observations also for nearly symmetric systems and the conditions for the occurrence of oscillations. Consider, e.g., Li-ion batteries. Here, a global constraint, corresponding to the constant current in our model, is imposed on the total charge. The hysteresis in the cell voltage $U$ over the total charge, which is measured when charging and discharging Li-ion batteries with phase-transition cathodes even at infinitely slow rates [14], could be seen as the “shadow” of the transcritical bifurcations at the turning points of the charge vs voltage curve of the individual nano-sized storage particles in the symmetric system.

Another physical system that shows similar symmetry-breaking bifurcations, including a Hopf bifurcation, are certain electronic circuits, which, in addition, can be considered superlattice models.³¹ While in that case the governing equations of the individual elements were two-dimensional, our studied revealed that stable, macroscopic collective oscillations may also exist in ensembles of coupled one-dimensional bistable units. This might lead to a reinterpretation of the origin of oscillations observed in many-particle systems, where macroscopic oscillations have been interpreted as the result of synchronization of oscillating units.

The authors thank A. Bonnefont for input concerning heterogeneous ensembles. C.B. acknowledges support from the Institute for Advanced Study at the Technical University of Munich through a Hans Fischer fellowship. This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through “e-conversion” Cluster of Excellence (Grant No. EXC 2089/1-390776260).

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

The authors have no conflicts to disclose.