Almost global convergence to practical synchronization in the generalized Kuramoto model on networks over the n-sphere

From the flashing of fireflies to autonomous robot swarms, synchronization phenomena are ubiquitous in nature and technology. They are commonly described by the Kuramoto model that, in this paper, we generalise to networks over n-dimensional spheres. We show that, for almost all initial conditions, the sphere model converges to a set with small diameter if the model parameters satisfy a given bound. Moreover, for even n, a special case of the generalized model can achieve phase synchronization with nonidentical frequency parameters. These results contrast with the standard n = 1 Kuramoto model, which is multistable (i.e., has multiple equilibria), and converges to phase synchronization only if the frequency parameters are identical. Hence, this paper shows that the generalized network Kuramoto models for n ≥ 2 displays more coherent and predictable behavior than the standard n = 1 model, a desirable property both in flocks of animals and for robot control. The Kuramoto model describes how collective synchronization appears spontaneously in a population of rhythmic units and is typically studied on a one dimensional circle. Here, the authors generalise the Kuramoto model on higher-dimensional manifolds and show that it achieves almost global convergence to synchronization and, in even dimensions, the fully synchronized state is attainable for nonidentical frequencies, in sharp contrast with the classical one-dimensional model.

S tarting with Winfree's seminal work in the late 1960s 1 , synchronization phenomena have attracted the interests of scientists in physics and other disciplines [2][3][4] and still remain an active area of research with many open problems. One of the most successful descriptors of synchronization phenomena is the Kuramoto model, along with its many generalizations [3][4][5][6] . Part of the appeal is its wide range of applicability in explaining phenomena from the natural world and technology, including flashing fireflies 7,8 , circadian rythms 9,10 , neuronal networks 11 , power-grid networks 12 , planar vehicle coordination 13,14 , etc. The Kuramoto model addresses this diversity of topics by capturing essential aspects of weakly coupled oscillators 11 , while still being analytically and numerically tractable.
Many extensions of the Kuramoto model have been proposed. Some bridged the Kuramoto model into the field of network theory [15][16][17][18][19] , extending it from complete graphs to complex networks 3,4 . A key question in this context is how the structure of the network affects the onset of synchronization 4 . For instance, Watts and Strogatz showed that synchronization is ameliorated by the small-world property of networks 20,21 . Moreover, it was demonstrated that scale-free networks resynchronize efficiently after small perturbations 22 . Other topics of interest include characterizing the equilibria of the homogeneous network Kuramoto model (i.e., with identical natural frequencies) and their basin of attraction by theoretical and numerical means 23,24 .
Most generalizations assign a single degree of freedom to each oscillator, but there are exceptions where the oscillators evolve on nonlinear manifolds [25][26][27][28][29][30][31] . Such extensions are partly motivated by curiosity, e.g., the swarmalators 29 or theoretical results on geometric aspects of synchronization 32 , and partly to address applications like rigid-body attitude synchronization of robot swarms 33 . In particular, the generalized n-sphere model provides a unified framework for the heading orientation on S 1 , the pointing orientation on S 2 , and the full attitude of any rigid-body on SO(3) using a map from the unit quaternions on S 3 . The nsphere Kuramoto model also appears in quantum synchronization on the Bloch sphere 27,34 , in opinion consensus dynamics 35 , and in computer science applications [36][37][38] .
In this paper, we complement the previous literature by investigating the global convergence of network Kuramoto models on nonlinear manifolds. Our model extends the ndimensional sphere Kuramoto model considered by Chandra, Girvan, and Ott [39][40][41] to complex networks on ellipsoids. Like them 39 , we find qualitative differences in the systems depending on the dimension n of the sphere. Moreover, like them 41 , we address the instability of incoherent equilibria but for finite networks rather than in thermodynamic limit N → ∞. We further characterize the convergence properties. For n ≥ 2, i.e., for all dimensions except that of the standard Kuramoto model (n = 1 in our notation), we find that the systems asymptotically converges to synchronization almost globally (i.e., from almost all initial conditions), provided that the spread of the frequency matrices is sufficiently small in the norm. Moreover, for frequency matrices that are all aligned and of even dimension n, we show the asymptotical stability of a single point in the null space of the frequency matrices. These two findings contrast with the standard network Kuramoto model, which has multiple stable equilibrium sets 23,24 and can only reach phase synchronization in the homogenous case, where all-natural frequencies are identical 3 . Hence, we show that the generalized network Kuramoto model displays more coherent behavior than the standard model, an appealing property in systems such as groups of animals, and for robust and reliable performance in robot control.

Results and discussion
Kuramoto models. Kuramoto first proposed a network of N harmonic oscillators on a plane, each with a natural frequency ω i 5 . An uncoupled oscillator satisfies € x i ¼ ω 2 i x i . Set x i ¼ cos θ i for unit amplitudes and phase angle θ i , then _ θ i ¼ ω i . The standard Kuramoto model of a complex network of coupled oscillators is given by where i = 1, …, N and k ij is the coupling strength between oscillator i and j. The coupling parameters k ij are weakly positive and symmetric, k ij = k ji ≥ 0. The network is recovered from k ij as a graph G ¼ ðV; EÞ where the nodes are the oscillators and the edges are node pairs with strictly positive coupling k ij . The Kuramoto model is usually expressed in polar coordinates as in Eq. (1) but it can also be formulated in Cartesian coordinates x i ¼ ½cos θ i sin θ i T 2 R nþ1 , where n = 1 is the dimension of the circle, which yields the dynamics where the matrix Ω 2 R nþ1 nþ1 is skew-symmetric. The Kuramoto model in Eq. (2) can easily be generalized. By changing the dimensions n to any strictly positive integer, we obtain x i 2 R nþ1 with ∥x i ∥ 2 = 1 and Ω i as a skew-symmetric (n + 1) × (n + 1) matrix. This yields the n-dimensional network Kuramoto model. For further perspective on the generalizations of Eq. (2), we consider the network Kuramoto model with homogeneous frequencies on the n-sphere. It is derived from Eq. (2) with Ω i = Ω j = Ω for all fi; jg 2 E. In this special case, like for the standard network Kuramoto model with homogenous frequencies, there is no loss of generality in setting Ω i = 0. In the following, we refer to this system as the homogeneous model. Note that the agent dynamics still differ with respect to k ij , which compose the effect of the network. The standard homogeneous (and heterogeneous) Kuramoto model can be interpreted as a gradient descent flow 42,43 . The homogeneous Kuramoto model on S n can also be interpreted as a gradient descent flow on the sphere surface: where v is an orthogonal projection on the tangent space T x i S n of the n-sphere at x i , and V is a potential function that measures the distance to the synchronized state.
Eq. (3) allows us to consider the natural frequency terms Ω i in Eq. (2) as perturbations of a homogeneous model obtained when all frequencies are equal. The advantage of such an approach is that we can use the large arsenal of tools from the optimization that are available for studying gradient descent flows. This idea leads to our main result, which is valid for frequencies Ω i that have a small spread around their mean frequency 1 N ∑ N i¼1 Ω i . This approach has also allowed us to generalize Eq. (3) to a large number of spaces including the Stiefel manifold 44 and general Riemannian manifolds 45 .
We remark on a third generalization of the Kuramoto model which turns out to be equivalent to the n-sphere model. It is easily verified that for any y i 2 E n ¼ fz 2 R nþ1 j z > Az ¼ 1g, the ellipsoid E n is invariant under see ref. 46 . The change of variables  28,47 . The global analysis of dispersed equilibria, as carried out in this paper, is challenging since it cannot rely on a reduction of a subset of the n-sphere to R n . Our main results establish that any dispersed equilibrium point of Eq. (2) is unstable in the sense of Lyapunov. Practical synchronization yields convergence to the desired set, i.e., synchronization 48 . Since all other equilibria are unstable, over time almost all trajectories turn to practically synchronized configurations. Most of the trajectories which originate near dispersed equilibria diverge from them. The only exceptions are the equilibria themselves and their (small) basins of attraction. The set N of dispersed equilibria and their basins of attraction have measure zero on ðS n Þ N . This implies that an initial condition that is uniformly distributed on ðS n Þ N has zero probability of being in N .
In order for our main results to hold, we require that the spread of frequencies Ω i is small in the following sense: where K ¼ min fi;jg2E k ij and n ≥ 2. The dependence of the bound in Eq. (5) on ðx i Þ N i¼1 can be removed. To do this, for dispersed configurations ðx i Þ N i¼1 , consider the left-hand side of Eq. (5): it can be upper bounded by a constant, obtained by replacing x i ⋅ x j in the right-hand side by cos π=N. The value x i Á x j ¼ cos π=N is found through optimization as explained in Supplementary Note 4. It holds for a path graph with N nodes that are connected as an arc. This arc minimizes the bound in Eq. (5) for dispersed configurations, while the total angle between agents satisfies ∑ N i¼1 acos x i Á x j ¼ π. Note that the condition Eq. (5) holds for the homogeneous model Eq. (3) except for phase synchronization 33 .
Interpretation. The result in Eq. (5) is surprising in light of what we know about the standard network Kuramoto model Eq. (1). Recall that Eq. (1) is multistable and has many other stable equilibrium sets besides synchronization 23 . These include the set  of q-twisted configurations of the homogeneous model Eq. (3), given by where the network is given by a cycle graph On the 2-sphere, we can define q-twisted configurations as equilibria of the homogeneous network Kuramoto model: where R is a rotation by 2qπ/N around a vector perpendicular to {x 1 , …, x N }, and addition of subindices is modulo N. In contrast to the q-twisted equilibria of the standard model 23 , our main result implies that T 2 q is an unstable equilibrium manifold of limited interest.
Since q-twisted configurations defined by T 2 q are unstable equilibria of Eq. (2), one might conjecture that, on the 2-sphere, there exist other configurations with analogous properties. The qtwisted configurations have the geometric appearance of regular polygons so regular polyhedra may be the equivalent of q-twisted configurations on the 2-sphere. Candidates for such configurations are the Platonic solids (the convex regular polyhedra in R 3 ) and Archimedian solids, as well as polyhedra corresponding to molecules of the fullerene family 49 . These configurations can be inscribed in a 2-sphere, such that their vertices lie on the sphere surface, see Fig. 2. Moreover, their edges have a certain symmetry that makes the configurations into equilibria of the homogeneous model Eq. (3), i.e., the sum of agent positions ∑ fi;jg2E k ij x j projected on the tangent space T x i S n is zero. By continuity of solutions 50 , there exist perturbed versions of these configurations that are equilibria of the heterogeneous model. However, our main result implies that they are unstable provided (Eq. 5) is satisfied.
To understand why there is a difference between the cases of n = 1 and n ≥ 2, it helps to consider the homogeneous model Eq. (3). Since the flow evolves along the negative gradient of V, it follows that V is decreasing with time. Also note that for k ij ≠ 0, by definition of V, Let us consider a 1-twisted configuration on the circle or 2sphere, i.e., holds for all future times t ∈ [0, ∞). In particular, for the standard Kuramoto model to go from a configuration in T 1 q to a synchronized state, there must be a time τ and a pair of oscillators that satisfy |θ i (τ) − θ j (τ)| = π. This is not possible for large N, as it would contradict Eq. (8).
lies on a great circle, which we without loss of generality assume to be the equator. We can perturb the agents towards the north pole without increasing V. To see this, use spherical coordinates By decreasing the s parameter from π/2 at the equator to 0 at the north pole we continuously deform the great circle into a point while staying on the sphere, see Fig. 3. The value of the potential function strictly decreases: V ¼ 4N sin s sin 2 π=N ! 0 as s → 0. For a more intuitive understanding of this result, let again the graph be the cycle given by Eq. (7). We may consider the agents as beads on a string, where each bead is an agent and the string is composed of all pairwise geodesic curves (great circles) that connect any two neighboring agents. This idea is valid in the much more general setting of gradient descent flows like Eq. (3) on Riemannian manifolds M 45 . These flows are generalizations of the homogeneous Kuramoto model with zero natural frequencies. The fact that V is decreasing can again be used to bound ∥x i − x i+1 ∥. For sufficiently small agent distances, the pairwise geodesics are unique and change continuously. This implies that the string, i.e., the closed curve composed of all pairwise geodesics, changes continuously. If M is multiply connected, then by definition of this property, there exists at least one closed curve which cannot be continuously deformed to a point. Hence we only need to choose our initial conditions so that the string approximates this curve and the topology of the manifold will prevent synchronization.
The manifold property of being simply connected versus multiply connected is one thing that distinguishes the circle from all other spheres (we do not consider the 0-sphere S 0 ¼ fÀ1; 1g). The circle is multiply connected since it cannot be continuously deformed to a point. The spheres are simply connected since any closed curve can be continuously deformed to a point as we did with the equator in Eq. (9) and Fig. 3. In this sense, synchronization on the circle is more similar to synchronization on SO(n), which is multiply connected, than to synchronization on S n . Note that the two manifolds are diffeomorphic S 1 ffi SOð2Þ through the linear map Phase synchronization. We have a second result that applies to spheres S n of even dimension n. Define convergence to a point, Fig. 3 Convergence to synchronization. Evolution of a configuration with 9 agents connected by a cycle graph on the sphere S 2 . The configuration is continuously deformed from the equator into a point on the north pole (i.e., synchronization) by letting the latitude parameter s in the map given by (Eq. 9) decrease from s = π/2 to s = 0. An equivalent version of this continuous deformation is not possible for the standard Kuramoto model since the circle is a simply connected manifold.
i.e., x i = x j for all {i, j} in E, as phase synchronization. We show that phase synchronization is possible on S n in the case of nonidentical natural frequencies. Like our first result Eq. (5), this finding distinguishes the network Kuramoto model on S 1 and S 2 . Phase synchronization on S 1 requires identical natural frequencies, ω i = ω j for all {i, j} in E. A key difference between this result and Eq. (5) is that it does not constrain the magnitude of the frequency norms in any way. However, like the findings of 39,51 it depends on the parity of the dimension n of the sphere. Let kerΩ i denote the null space of the natural frequency matrix Ω i . If there exists a nonzero vector v such that then for any initial condition which satisfies it holds that x i → v as t → ∞. Moreover, v is an asymptotically stable equilibrium. Note that if the frequency matrices satisfy the synchronization condition given by our result, Eq. (5), then the convergence to phase synchronization is almost global. Simulations also indicate that condition (Eq. 11) is redundant. To understand why it's needed for the proof, see Supplementary Note 5. To interpret the synchronization condition, consider the 2-sphere. If e.g., we chose v = e 3 , it means Initialize all agents in the northern hemisphere, then all agents converge to phase synchronization in the north pole. Note that a random Ω i ∈ so(n + 1) has full rank if n is odd, whereby v = 0 and (Eq. 11) cannot hold. For even n however, the dimension of kerΩ i is always at least 1.
It should be noted that condition (Eq. 10) refers to a special case that would not appear if the frequency matrices are selected uniformly at random. Rather, it suggests that all frequency matrices are in a sense somewhat similar. This may occur in a real-world system where the drift term in Eq. (2) is due to an exterior cause. For example, if the agents are fish they may be caught in a stream, or if they are drones the external influence could be wind. Such situations can also be modeled by our assumption (Eq. 5) which assumes the frequencies have a small spread around a large mean value. Another scenario is a single agent with a nonzero frequency matrix, perhaps due to a malfunction in an otherwise homogeneous robotic system. From a theoretical perspective, condition (Eq. 10) is interesting as a counterexample to what one might assume based on the standard Kuramoto model since it shows a counter-intuitive property of the high-dimensional model.
To understand the mechanics behind the result (Eq. 11), it is helpful to consider a two-agent case on S 2 , where Ω 1 yields a rotation around the vector e 3 . Consider three cases: (i) If the agents are in sync, then agent 2 will break away from agent 1 due to the drift term Ω 1 x 1 . The only exception to this is if x 1 belongs to kerΩ 1 , i.e., x 1 = e 3 ; (ii) If the agents are on different latitudes θ 1 ≠ θ 2 in the spherical coordinates Eq. (9), then the latitude θ i of the agent i with the lowest latitude will increase. In sum, the agent with the lowest latitude θ i moves closer to e 3 ; (iii) If the agents have the same latitude θ 1 = θ 2 , then the drift term Ω 1 x 1 in the dynamics of agent 1 will cause it to rotate at this latitude. However, both agents will also move closer to the north pole. This is due to the coupling terms being aligned with the tangent space of the great circle connecting the two agents (see Fig. 4). The great circle connecting two points of the same latitude has a segment at a higher latitude than the common latitude of the two points themselves.

Simulations
Practical synchronization. We are interested in how geometric and topological properties of the manifold affect the convergence of Kuramoto models. To this end, we study a large family of manifolds that include the n-spheres, the generalizations of S 1 ffi SOð2Þ to the special orthogonal group SO(n), and other manifolds in-between. This family is the Stiefel manifolds St(p, n) defined by where p ≤ n. The Stiefel manifolds are rectangular orthogonal matrices; the manifold St(p, n) is composed of the first p columns of all orthogonal matrices Q ∈ O(n). The following diffeomorphisms show the versatility of the Stiefel manifold framework Sð1; n þ 1Þ ffi S n , St(n − 1, n) ≅ SO(n), and St(n, n) ≅ O(n). The Kuramoto model on Stiefel manifolds can be obtained as a generalization of the gradient flow formulation for the n-sphere model, Eq. (3). Along those lines, the following model is introduced in 44 : Eq. (13) can be expressed as a coupled oscillator system by the addition of frequency terms, where Ω ∈ so(n), Ξ ∈ so(p). See refs. 44,52 for more details. For the homogeneous model, we are interested in the measure μðBðSÞÞ of the basin of attraction BðSÞ of the synchronization manifold S. To find an equivalent of μðBðSÞÞ for practical synchronization we introduce an order parameter The factor 1/p ensures that R = ∥R∥ F = 1 when all matrices are equal, S i = S j for all fi; jg 2 E. By taking the average value of R at a large time T over many simulations we get an index of convergence to synchronization for Eq. (14) that is similar to what μ BðSÞ ð Þis for Eq. (13). Note that R is a generalization of the order parameter of the standard Kuramoto model on S 1 3 .
The details of numerical integration of Eq. (13) are given in the "Methods" section. Results of the simulation are displayed in Table 1. Note that almost global convergence to synchronization on St(p, n) is observed for pairs (p, n) that satisfy p ≤ n − 2. That almost global convergence holds for pairs that satisfy p ≤ 2n/3 − 1 has been proved analytically 44 . However, we conjecture that it also holds for p ≤ n − 2, which corresponds to all the simply connected Stiefel manifolds 53 . As we argued in the Interpretation section based on ref. 45 , almost global convergence for all connected graphs is impossible on a manifold that is multiply connected. Failures to reach consensus occur when p ∈ {n − 1, n}, i.e., for the special orthogonal group and the orthogonal group. For the case of O(n), the probability that all agents belong to the same connected component (det Q ¼ 1 or det Q ¼ À1) is 2 −4 ≈ 0.06, which explains the numbers on the diagonal where p = n.
Many results that guarantee synchronization on nonlinear space are limited to geodesically convex sets, e.g., hemispheres of the n-sphere 28,47,48,54 . This does not scale well for large N since the probability that all N agents belong to a hemisphere decreases exponentially with N. However, even though a theoretical guarantee does not scale well, actual system performance may still do. This is the case on S n for n 2 Nnf1g for the homogeneous model since S is almost globally stable 33 .
It is of interest to investigate how μðBÞ scales with N when the consensus manifold is not almost globally stable. Figure 5 reveals that on Stð1; 2Þ ffi S 1 , μðBÞ decreases steadily with N (see also ref. 23 ). On St(2, 3) ≅ SO(3), St(3, 4) ≅ SO(4), and St(4, 5) ≅ SO(5) the value of μðBÞ settles on 0.5. Note that there is not any immediate reason to expect that lim N!1 μðBÞ even exists. For each N = 5k, where k 2 N, there is a new q-twisted equilibrium set that is stable on St(n − 1, n). However, there does not appear to be any large difference between μðBÞ when N = 5k for k 2 N compared to N − 1 except for the case of N = 5. The size of St(p, n) N grows with N. The average time required to satisfy our simulation termination criteria increases with N. This introduces a dependence on N in the bias towards underestimating μðBÞ, but the effect is small as we have investigated by increasing the termination time T.
We compare the size of the synchronization basin, μ BðSÞ ð Þ with the average order parameter from M simulations, 1=M ∑ M i¼1 R, see Fig. 5. The way the μ and R curves depend on the number of agents N are strikingly similar. Inspect condition Eq. (5) for n = 2, a network given by the cycle graph C N , and a constant RHS obtained by setting hx i ; x j i ¼ cos π=N. This yields which suggest that Ω i should scale as OðK=N 4 Þ. However, the use of hx i ; x j i ¼ cos π=N in the condition is conservative since it is derived for a worst-case scenario. In simulations, we find that scaling Ω i as 1/N suffices to approximately replicate the μ BðSÞ ð Þ curves, see Fig. 5.
Phase synchronization. We use simulations to study the order parameter R given by Eq. (15) as a function of the coupling parameters K ¼ min fi;jg2E k ij under the assumptions of the synchronization condition in Eq. (10). We draw Ω 1 such that ðΩ 1 Þ st 2 Nð1; 1Þ if s < t and ðΩ 1 Þ st 2 NðÀ1; 1Þ if s > t. Then we set Ω i = ω i Ω 1 where ω i ∈ N(0, 1) for i = 2, …, N. Generically, this results in skew-symmetric matrices of maximal rank. However, that means rank Ω i = n for even n and rank Ω i = n + 1 for odd n. In particular, this implies the existence of a v satisfying Eq. (10) for even n but not for odd n. We find that it is not necessary to limit the initial conditions to a hemisphere as we required for our analytic result; the simulated system is observed to converge to phase synchronization from (almost) all initial conditions if a single joint nullspace vector exists. We use the complete graph G ¼ K N since it allows for a comparison of our findings with similar results in refs. 39,51 . For this reason, we also use the normalization factor 1/N in front of the coupling terms in Eq. (2). In particular, we find that the phase transition from disorder to synchronization is explosive and discontinuous for even n but not for odd n where it is continuous, see Fig. 6. We expected this difference based on Eq. (10) since the full rank of generic skewsymmetric matrices with even dimension n + 1 (the dimension of the matrix, not the sphere) prevents condition Eq. (10) from applying. We also find that practical synchronization holds when small perturbations are added to each of the frequency matrices (note that the condition of Eq. (5) does not need to be satisfied).

Methods
Proof sketch. A full derivation of the result in Section "Practical synchronization" is given in Supplementary Note 2-4. This section provides a brief sketch of the main ideas. For the sake of simplicity, we focus on the stronger assumption 2K Nðn þ 1Þ ðn À 1 À cos π=NÞ Á ð1 À cos π=NÞ: Note that we do not require ∑ N j¼1 Ω j ¼ 0 in Eq. (5) for Eq. (16) to hold, rather we just let each Ω i be small in the norm.
Our approach is based on a matrix perturbation result 55 that relates the spectrum of three matrices. Let A(x) and B(x) be the linearization matrices of (Eq. 2) and its particularization for Probability measure, μ BðSÞ ð Þ2 ½0; 1, of the basin of attraction BðSÞ of the sync manifold S on the Stiefel manifold St(p, n) with parameters p, n for a cycle network with 5 agents given by the graph C 5 defined by Eq. (7). The calculation of μðBÞ is done by Monte Carlo integration using M = 10 4 samples of the uniform distribution on St(p, n) for each pair (p, n), see "Methods" section for details. We set the threshold value ε = 0.01 in Eq. (20) and the final time T = 200. Rows in the table fix p, columns fix n. An empty cell indicates that there is no Stiefel manifold for that pair (p, n). Bold font indicates that (p, n) satisfies 2 3 n À 1<p p À 2 whereby our previous result 44 does not apply. Ω = diag(Ω 1 , …, Ω N ), using skew-symmetry of Ω and symmetry of B(x). The result on matrix perturbations 55 yields jRe αðxÞ À βðxÞj ≤ k AðxÞ À BðxÞk 2 where Re αðxÞ and β(x) are the eigenvalues with the largest real part of A(x) and B(x) respectively. A lower bound for β(x) is known 33 . We use techniques from optimization to derive the following inequality which holds at any dispersed equilibrium x, Suppose that the assumption (Eq. 16) holds, then where we used the result (Eq. 17). From jRe αðxÞ À βðxÞj < βðxÞ we get Re αðxÞ > 0. It follows that the equilibrium x is exponentially unstable.
Simulations. We perform simulations to investigate how the basin of attraction BðSÞ ¼ fðS i Þ N i¼1 2 Stðp; nÞ N j lim t!1 Φðt; ðS i Þ N i¼1 Þ 2 Sg: For the heterogeneous model, we have a practical synchronization manifold D of (14) defined by The existence of such a manifold for sufficiently small Ω i , Ξ i follows by a perturbation theory argument 48 . The probability measure μðBðSÞÞ is the fraction of initial conditions for which the system in Eq. (13) converges to the synchronization manifold S. The measure μðBðSÞÞ can be calculated by Monte Carlo integration: À! a:s: μðBÞ as M ! 1; where 1 : St(p, n) N → {0, 1} is the indicator function and ðS k i Þ N i¼1 for each k ∈ {1, …, M} is a sample drawn from the uniform distribution on St(p, n) N . A uniform sample S is found by drawing X 2 R n p such that each element of X is independent and identically normally distributed N ð0; 1Þ and forming S ¼ XðX > XÞ for a threshold value ε ∈ (0, ∞) at a fixed time T. The tensor ΦðT; ðS k;0 Þ N k¼1 Þ is obtained by integrating Eq. (13) using the function ode45 in MATLAB. If Eq. (20) is satisfied at T, we count this as a case of convergence to S. There are potential  15); each point corresponds to 10 4 simulations. The graph is the cycle graph C N given by (7). We set the coupling parameters between neighboring agents i and j to k ij = 1. The (s, t)-elements of the frequency matrix ðΩ i Þ st 2 N ð1; 1=NÞ 2 if s < t and ðΩ i Þ st 2 N ðÀ1; 1=N 2 Þ if s > t, where N ðÁ; ÁÞ denotes a normal distribution. Fig. 6 The transition from incoherence to synchronization depends on the parity (odd/even) of the dimension n of the n-sphere S n . The order parameter R is given by Eq. (15). The uniform coupling parameters k ij = K between agent i and j range from negative, where R = 0 indicates incoherence, to positive, where R = 1 indicates synchronization. We take the average value of R from M = 10 3 system evolutions per point. We use N = 100 agents. To find the discontinuous phase transition for small K values we have chosen a large final time T = 2000.
issues with long simulation times causing numerical errors to accumulate so that either the system leaves the Stiefel manifold or is perturbed from B c into B. We have dealt with such issues. In particular, there might be a bias towards underestimating the value of μðBÞ due to T < ∞, but the effect is small for large T.

Data availability
The data used in this paper is available from the GitHub account of the corresponding author: https://github.com/johanmarkdahl/CommPhys2021. The data is also available on request to the corresponding author.