Binary system modes of matrix-coupled multidimensional oscillators

The standard Kuramoto model has been instrumental in explaining synchronization and desynchronization, two emergent phenomena often observed in biological, neuronal, and physical systems. While the Kuramoto model has turned out effective with one-dimensional oscillators, real-world systems often involve high-dimensional interacting units, such as biological swarms, necessitating a model of multidimensional oscillators. However, existing high-dimensional generalizations of the Kuramoto model commonly rely on a scalar-valued coupling strength, which limits their ability to capture the full complexity of high-dimensional interactions. This work introduces a matrix, A, to couple the interconnected components of the oscillators in a d-dimensional space, leading to a matrix-coupled multidimensional Kuramoto model that approximates a prototypical swarm dynamics by its first-order Fourier harmonics. Moreover, the matrix A introduces an inter-dimensional higher-order interaction that partly accounts for the emergence of 2 d system modes in a d-dimensional population, where each dimension can either be synchronized or desynchronized, represented by a set of almost binary order parameters. The binary system modes capture characteristic swarm behaviors such as fish milling or polarized schooling. Additionally, our findings provides a theoretical analogy to cerebral activity, where the resting state and the activated state coexist unihemispherically. It also suggests a new possibility for information storage in oscillatory neural networks.


Introduction
Emergent macroscopic behavior from local or long-range interactions has been both an observation captured in biological, physical, and sociological systems [1,2], and a long-standing topic of engineering significance [3,4].Among various collective behaviors such as formation [5], consensus [6], and coverage [7], synchronization stands out as a phenomenon where two or more distinct oscillating systems adjust their rhythms or cycles to operate at the same frequency or in a specific phase relationship [8,9].A paradigmatic framework to understand the synchronization phenomena in these fields is the Kuramoto model.Originally proposed to study the phase dynamics of chemical oscillators, it exhibits a competition between the heterogeneity from the dispersed natural frequencies of the oscillators and an all-to-all sinusoidally coupling term that reconciles this heterogeneity.The Kuramoto model and its early variations not only provide phenomenological insights into collective animal behavior [10,11], circadian rhythms [12,13], and activity of neuronal networks [14] in biological systems, but are also applicable in physical and engineering contexts such as the superconducting Josephson junctions [15,16], laser arrays [17,18], and distributive clock synchronization [19].Synchronization of the model has also been examined from a Hamiltonian perspective [20,21] which associates the classical phenomenon to the emergence of quantum entanglement on the invariant manifolds.
However, due to the single dimensionality of the Kuramoto oscillators, the model in its classic form was rarely directly associated with the spatiotemporal evolution of entities moving in the three-dimensional space, such as biological and robotic swarms characterized by vector-valued elements.In recent years, high-dimensional generalizations of the Kuramoto model are being actively researched [22][23][24][25][26], including, but not limited to, the angle model on the unit sphere [23,24,27] and the swarmalator model that incorporates a distance-dependent coupling into the phase dynamics [28,29].Despite their diverse mathematical formulations, it is commonly accepted that the intensity of the interactions is modulated by a scalar-valued coupling strength.Although this approach is simple and intuitive for the extension from the one-dimensional case to the high-dimensional, it overlooks alternative possibilities for the interdependency within the dimensions of the oscillators.This oversight may lead to an oversimplified representation of the system dynamics and limit its ability to capture the full complexity of high-dimensional interactions.
In this work, we develop another high-dimensional generalization of the Kuramoto model with a matrix-valued coupling strength, where we exploit the fact that a matrix naturally represents an inter-dimensional dependence among multidimensional entities.In the proposed model, this dependence manifests as an implicit higher-order interaction that plays a crucial role in complex structures such as the brain, protein network, and social network [30][31][32][33][34]. Due to the introduction of matrix A, any component of an oscillator in a d-dimensional space will not only form pairwise interactions with components of other oscillators within the same dimension, but will also be implicitly involved in three-party interactions with components from any other dimensions.Hence, the proposed model can be viewed as a network embedded in a simplicial complex of links and triangles, connecting it to a broader context (see, for example, [35][36][37][38] and references therein).It is worthwhile to mention that the formulation with a scalar-valued coupling strength can be considered a special case of our model, where the matrix A reduces to a diagonal matrix that does not carry any higher-order interaction.
To demonstrate the applicability of our model, we have made an example out of the swarm system, with vector-valued agents and vector-valued nonlinear mappings for their interactions [39,40].Our model approximates the prototypical swarm model by the first Fourier harmonics on the d-torus.In this context, the coupling matrix A consists of Fourier coefficients to the first-order basis, thus varies accordingly to the specific choice of the vector-valued nonlinear mapping.It is expected that, by performing large-scale simulations of our model with various realizations of A, some statistical patterns may arise that align with the observations about biological swarms in nature.
Indeed, what the matrix-coupling mechanism brings to the multidimensional Kuramoto model is a unique but orderly way for synchrony and desynchrony to coexist which simulates, for example, the behaviors of the fish swarm.In this work, we report the discovery of what is referred to as the binary system mode that can be accounted as follows.Let the degree of coherence be measured on each dimension of the population.Then given all possible values of the coupling matrix A in R d×d , any dimension of the population has the potential to synchronize through a second-order phase transition when others remain coherent or incoherent.Therefore, for a d-dimensional population, the system is endowed with 2 d modes as steady states where every dimension of the population is either synchronized or desynchronized.One natural phenomenon to visualize a system mode with both a synchronized dimension and a desynchronized one is the fish schooling as a vortex when fencing off the predators [41,42].This behavior known as milling has the population remain largely static on the vertical direction, but possess varying speeds on the horizontal plane from the periphery to the center of the vortex.Another system mode representing a full synchronization can be compared to the fish with polarized schooling, where the velocities are aligned on all dimensions [43].
Upon the discovery of the aforementioned steady states, one pertinent question arises: what factors determine the emergence of one specific system mode over the others?In addressing one of the technical challenges, we reveal through stability analysis how the algebraic properties and the diagonal elements of the coupling matrix give rise to the system modes as a macroscopic, statistical pattern.The phase diagram is also derived with respect to the coupling matrix through a self-consistent analysis.Additionally, we have observed explosive dimensional phase transitions [25,35,37,44] in a numerical study demonstrated in the appendix, suggesting that the phenomenon induced by the matrix coupling mechanism may be far richer than what is presented in the current study.This reinforces the potential of the proposed model.

Motivating scenario
The standard Kuramoto model [45] describes synchronization of coupled limit-cycle oscillators, N in total, where each oscillator is represented by θ i ∈ [0, 2π) and exhibits inherent heterogeneity: where a one-dimensional oscillator θ i has a natural frequency of oscillation ω i extracted from a unimodal probability distribution.While the heterogeneity of ω i impedes synchronization, the interaction term of the model promotes it, that is modulated by a normalized, scalar-valued parameter K termed the coupling strength.With the increasing coupling strength K over a threshold K c > 0, mutual entrainment of oscillation starts to form in tiny cliques, then propagates to the entire population.In the thermodynamic limit N → ∞, this synchronization process can be properly characterized by a second-order phase transition and a bifurcation of the steady state solutions.
The model we propose in this work is an extension of (1) that adopts a matrix-valued coupling strength, i.e.
Equation ( 2) depicts a population of d-dimensional oscillators on the d-torus, denoted as Here sin(•) is a vector of sine functions, where each element represents a sine function.Each oscillator is associated with a natural frequency vector ω i = ω i1 ω i2 . . .ω id T that is allowed to vary so that ω i = ω j .The interaction between any θ i and θ j is weighed by a coupling matrix A ∈ R d×d matching the dimension of the oscillators.Additionally, we put 1/N for the equation ( 2) to be well-behaved in the thermodynamic limit N → ∞.It is then noteworthy that the proposed model reduces to the standard Kuramoto model in the one-dimensional case d = 1, where A ∈ R becomes the conventional coupling strength K.This indicates that, for there to be an actual matrix effect on the interactions, the oscillators characterized have to be at least two-dimensional.
For the remainder of this section, we motivate the model from two perspectives.Topologically, we demonstrate that the coupling matrix in (2) introduces an implicit higher-order interaction with a simple network of three oscillators, each of dimension two.Dynamically, equation ( 2) is also derived from a first harmonic approximation of a prototypical swarm model, where the coupling matrix A reflects the distribution of the nonlinear mapping on the Fourier basis.

Matrix-induced higher-order interaction
Let us first examine the simplest high-dimensional case, a system of two-dimensional oscillators whose dynamics reads where {ω i1 } and {ω i2 } are the natural frequencies of oscillation on the {θ i1 } dimension and {θ i2 } dimension.
The interaction is realized through a normalized, all-to-all coupling with the matrix  where θ 11 is directly involved; while the remaining terms happen on the {θ i2 } dimension that do not concern θ 11 , which is a rare occurrence for the high-dimensionally generalized Kuramoto dynamics.For those adopting a diffusive coupling, dimension-wise dynamics typically adopt a pairwise interaction that includes the component to be updated.One possible explanation of the extra terms a12 3 [sin(θ 22 − θ 12 ) + sin(θ 32 − θ 12 )] is to view them as the result of vector additions, e.g.sin(θ 22 − θ 12 ) = sin ((θ 22 − θ 11 ) − (θ 12 − θ 11 )), indicating a joint interaction among θ 11 , θ 12 , and θ 22 .An interaction that simultaneously happens on three or more parties is called a higher-order interaction, which is topologically modeled as a simplex.It is now evident for (3) that, every nonzero off-diagonal element of A accompanies terms representing inter-dimensional higher-order interactions that each involves three parties.This is also true for the general case (2).We could now say that the proposed model is embedded in a simplicial complex structure that integrates intra-dimensional links (1-simplex) and inter-dimensional triangles (2-simplex).The higher-order interactions may serve as the microscopic foundation for the collective behaviors well observed within our model but not in other high-dimensional generalizations.

Approximation of swarm dynamics
Historically, the sine-coupling in the Kuramoto model has been put in the context of being the first harmonic approximation of a general coupling function h(•) containing higher harmonic components [46,47], that appears in The purpose of (5) was to investigate what might be fundamentally common for the mutual entrainment phenomena observed on chemical, electrical, biological, and physiological oscillators [46].This general dynamics is closely related to (1), as it was shown that some essential properties of (1) survived in ( 5), e.g.incoherence for weak coupling and bifurcation at a critical value of K, establishing (1) as a paradigmatic model for studying the synchronization phenomena in natural and artificial systems.The formulation of ( 5) is particularly meaningful for biological and robotic swarms [39,40,48] where, however, a position vector rather than a limit cycle is considered as the state for an individual.In the high-dimensional state-space, a homogeneous version of ( 5) is considered in [39,48] with the following general form: with G(•) : R d → R d being a vector-valued nonlinear function to describe the attraction and repulsion of swarming bodies x i ∈ R d .For aggregating swarms, classes of functions were examined for (6) that unanimously satisfy two rules: for the swarm to have a static center that is key to aggregation, and G(•) has a pair(s) of opposite stationary points at − → p ′ and − − → p ′ within the range of a certain distance at which the repulsion balances the attraction [39,48].The coupling function G(•) acts very much like a sine-wave on each of its dimensions, motivating us to do the following approximation with its first Fourier harmonics.Like with the standard Kuramoto model, the idea is to capture some of the intrinsic properties of (6) on a particular interval, even as the system assumes distinct formulations.Note that for now we use − → p to accentuate it being a vector in R d .It is known that [49] where f( − → n ) is the Fourier coefficient of the − → n th harmonic.To guarantee the validity of the expansion, we assume every component . Similar to the one-dimensional case, we refer to the terms with only one variable and multiplicity 1 as the first harmonics of G * r ( the entire function can be expanded as , then an approximation of G * ( − → p ) with its zeroth and first harmonics3 yields Since G(•) and G * (•) are real functions, there is in fact i G1 ∈ R d×d .By inserting (7) into (6), we demonstrate that model (6) with a nonlinear interaction function G(•) can be approximated by the following dynamics, where the matrix Ã = i G1 of Fourier coefficients couples the sine-valued vectors.Although due to the expansion into Fourier series, it is more appropriate to consider each dimension of x i modulo 2π.We therefore adopt the heterogeneous Kuramoto framework with rescaled time and random offsets.This allows us to propose the following matrix-coupled, multidimensional phase oscillator model which is the other way to motivate model (2).We note that the multistable solutions to be introduced in the following sections can also be observed in a homogeneous setting, whose stability is readily inferred from the linear stability analysis of the current model.For that reason, our focus will be on the heterogeneous population for the remainder of this work.The proposed dynamics is a direct high-dimensional generalization of the standard Kuramoto model (1).The resemblance between (1) and ( 9) lies not only in their mathematical formulation, but also in their role as the first harmonic approximation of a widely applicable N-body dynamics.Therefore, we expect the matrix-coupled multidimensional Kuramoto model to have real-world instances beyond the swarming scenario presented here.
→ n , one can also consider including the multivariate terms, e.g.those with bases exp (i (s j p j + s k p k )) , s j = ±1, s k = ±1, j, k ∈ d.However, any basis containing more than one variable will correspond to a coefficient matrix that is not square, introducing significant complexity to the dynamics.Therefore we deem the current level of approximation suitable for developing numerical and theoretical analysis.

Dimensional order parameters
In this section, we introduce the dimensional order parameters to derive a mean-field characterization of the model in arbitrary dimensions.The complex order parameter for each dimension is defined as e iθ jr = σ r e iΨr , r ∈ {1, 2, . . ., d} , (10) where the magnitude satisfies 0 ⩽ σ r ⩽ 1 and Ψ r denotes the average phase of θ jr projected onto the complex plane.When σ r → 1, it indicates a high degree of synchrony on the rth dimension of the oscillators, meanwhile σ r → 0 indicates an incoherent distribution of {θ jr } where their complex projections cancel out.
In our example of the two dimensional oscillators (3), the equation of motion can then be adapted to which demonstrates that the instantaneous frequencies are, in fact, directly controlled by the {θ i1 } mean-field and {θ i2 } mean-field.The mean-field actions are then modulated through the matrix elements of the same row.
In general, it is also possible to derive a compact form of (2) with the order parameters for d dimensions, i.e.
Here 12) is essential in determining the fixed point of an individual and the distribution of the population.We see from (12) that it is not always possible for θ i to be fixed on every dimension, where the natural frequency ω ir is dominating over the mean-field influence.When θi = 0 is indeed solvable, it is ideal to have an invertible A to obtain a unique, closed-form expression of the fixed point θ * i ; we refer to such oscillators as being fully synchronized (by the mean-field).For oscillators that are not fixed on all of their dimensions, it is expected that they will be relatively in motion to the entrained population, for which we refer to them as the drifting oscillators.

Numerical results
Our main objective is to investigate the effect of varying coupling matrices on the order parameter of each dimension, thus to identify the qualitatively distinct states of the system induced by this coupling mechanism.Now that the parameter space is augmented from the conventional R to R d×d , there are obviously numerous ways that each element of the coupling matrix can be adjusted.Our approach involves scaling an invertible, real symmetric initial matrix by an incremental real number, commencing at zero.This enables us to preserve, or at least trace, the algebraic characteristics of the initial matrix while observing the complete transition the system undergoes, from a state of no matrix coupling effect to one where such coupling is particularly strong.We also make the assumption that the natural frequency ω i of an oscillator does not distinguish between different dimensions, i.e. ω i = ωi 1 d , whereas ωi for i ∈ {1, . . ., N} follow a unimodal probability distribution.

Observing the binary system modes
Under the above assumptions, the key finding of this research is demonstrated in figure 2 from an experiment we performed on N = 5000 oscillators.The heterogeneous population is mediated respectively by the four randomly generated matrices A 1 , A 2 , A 3 , A 4 from a uniform initial distribution.In each case of A = kA u , u ∈ {1, 2, 3, 4}, we adiabatically increase k from 0 to 3.2 in 22 steps and calculate the average of σ 1,2 in a span of time, when the order parameters are considered to have reached their equilibrium.We note that this scaling is fundamentally different from the strengthening K in other generalized Kuramoto models, as the matrix elements are inherently d(d + 1)/2 independent variables undergoing changes in the parameter space.As a result, figure 2 reports that with A 1 and A 4 , both the {θ i1 } dimension and the {θ i2 } dimension go through qualitatively similar transitions either from being largely incoherent, due to the limit size effect, to a complete frequency synchronization (σ 1,2 → 1), or to a complete desynchronization (σ 1,2 → 0).Meanwhile with A 2 and A 3 , the transitions on {θ i1 } dimension and {θ i2 } dimension go to opposite directions, i.e. σ 1 → 1, σ 2 → 0 or σ 1 → 0, σ 2 → 1 as the elements of A are simultaneously and sufficiently increased.Take θ i1 as the toroidal angle and θ i2 as the poloidal angle, we can project these combinations of σ 1,2 onto the 2-torus  for visualization, i.e. each oscillator has the cartesian coordinates 3 demonstrates how the system has set into four qualitatively distinct modes of distribution and motion, each to a configuration of A u .
Motivated by what is observed from the two-dimensional population, we then apply the same method on an ensemble of three dimensional Kuramoto oscillators, coupled through matrices A u ∈ R 3×3 , producing a total of eight combinations of σ 1,2,3 as depicted in figure 4. At this stage, what we have shown with two dimensional and three dimensional population is that, the matrix coupling enables the multidimensional Kuramoto oscillators to separate the transition to coherence/incoherence on each dimension, which we call an independent dimensional phase transition.

Binary system modes as multistable solutions
Is any of the above system mode the only solution we can observe from a fixed coupling matrix?Among the majority of the generalized Kuramoto models, it is not uncommon to have systems that are multistable at a particular coupling strength or within a range of the coupling strength; therefore the basins of attraction become decisive in the emergent equilibria [50,51].In figure 4, our experiment on N = 100 matrix-coupled oscillators indicates that, with a fixed value of the matrix A = kA u , the system modes may indeed be limitedly attracting to the 10 4 uniformly distributed initial values with different k-values.The criteria that determine which system mode the data correspond to, even when σ r may be significantly between (0, 1) due to small k, is detailed in the next subsection.
Notably, in our experiments with A 1 and A 4 , there is no display of multistability in contrast to those with A 2 and A 3 , where all system modes except mode 11 are identified.Yet interestingly, when we gather the 2668 samples of θ(0) from the A 2 experiment that converge to mode 01 at k = 20, and start to progressively increase k from zero with these initial values, the dimensional phase transitions induced are all towards the equilibrium mode 10 as those in figure 3(b).In fact, all 10 4 initial values have yielded the above result with ) in a combinatorial sense, the corresponding system modes σ3σ2σ1 also runs from 000 (panel (a)) to 111 (panel (h)).Right: line graph illustrates the percentage of the 10 4 uniformly distributed initial values that converge to the same system mode at kAu to that in figure 2, for u ∈ {1, 2, 3, 4}, k ∈ {1, 1.5, 2, 20, 40}.Bar chart displays the average values of σ1, σ2 at each k that are considered to have chosen mode 11 with A1 (two bars on the left), and those considered to have chosen mode 10 with A2 (two bars on the right).The complete data can be found in the appendix.
the adiabatically increased k ∈ [0, 2.3].The same occurred to our examination of the 160 samples that converge to mode 10 with kA 3 , k = 20, where the system eventually converges to mode 01 with incremental k, as all 10 4 initial values do.This suggests that the system mode to be observed has to do with our practice of the experiment which, in this work, is to simultaneously scale all the elements of the matrix from zero.The scaling process gives rise to a unique, definitive outcome out of the many possible stable solutions.Subsequently in the article, our theory will reveal the exact multistable solutions to expect from a coupling matrix, and the reason why the system mode we observe emerges from the scaling experiment.

Simulation specifications
The four initial matrices adopted to generate figure 3  .The differential equation ( 2) was numerically integrated using a 5th order Runge-Kutta formula, with relative error tolerance 10 −5 and step length ∆t = 0.01.The natural frequency ω i = ωi 1 2 of the population follows the standard Lorentzian distribution g(ω) = 1/[π(ω 2 + 1)].At each value of k, we integrated until the order parameter σ 1 (σ 2 ) had settled around an equilibrium, then took the average of σ(t) for the last 1500 integration units.The experiment performed on the three dimensional population has adopted the same convention.
In our test of the 10 4 initial values, we determine the system mode a particular θ(0) is attracted to based on the numerical result in figure 3 and the following experimental fact: whenever a dimension of the population is tending to an incoherent solution from an almost (but not exactly) uniform initial distribution, be it with N = 1000 or with N = 5000, its order parameter σ r significantly overcomes the limit size effect which causes σ r = ϵ > 0 at k = 0, and becomes almost vanishing for k > k critical .For N = 100, there is approximately ϵ = 0.1.Consequently, we set up the criterion that if σ 1 σ 2 is comparable to the data in figure 3 we consider θ(0) to be attracted to mode 10.The criteria for mode 01 and mode 00 are similarly defined, with Note that all the initial value tests are performed on the same population, i.e. the natural frequencies {ω i } are controlled.

Theoretical results
In this section, we present an overview of the theoretical results addressing several problems arising from the numerical experiments.These include clarifying the mechanism governing the system mode and deriving the steady state solutions of the dimensional order parameters.We then revisit the standard Kuramoto model to show that some essential theories developed on it can be fully recovered in our framework.Meanwhile we find that the opposite is not true: system modes with both synchronized and desynchronized dimensions, which are intricately related to the matrix-induced higher-order interaction, will not be observed with a scalar-valued coupling adopted by the standard K model.

Linear stability analysis
The key problem at hand is what determines the system mode when the matrix coupling is sufficiently strong.Our findings suggest that it is the combined effect of the principle submatrices and the diagonal elements of A: with the increasing matrix-coupling effect, dimensions with the larger diagonals are earlier in exhibiting phase transitions; any subsequent dimension joins in this order only if the principle submatrix of A considering that dimension remains positive definite.
To delve deeper into this explanation, we outline several simple rules derived from the stability analysis (refer to appendices B and D), then show how together, they associate a coupling matrix to a system mode at sufficiently large k.

Linear stability of the system modes
Under the premise that A is real symmetric and invertible, the linear stability analysis reveals that: (1) The system mode that corresponds to a full coherence (σ is linearly stable at the end of the transition A = kA initial if and only if the initial matrix is positive definite.(2) If the system admits mode at the end of the transition A = kA initial , then the principle submatrix consisting dimensions {i τ +1 , i τ +2 , . . ., i d } of A is positive definite.
(3) If there exists a diagonal element a rr < 0, then σ r → 0, the rth dimension of the population will remain incoherent for the entire transition.(4) If the initial matrix is negative definite, then the incoherent solution is stable on every dimension, and we observe σ ) Given a d dimensional population, the 2 d system modes exist under the matrix coupling mechanism.
Here (1) and ( 2) are obtained through a linearization of the system of finite population around the fully coherent solution, followed by our analysis on the distribution of the Jacobian eigenvalues; (4) is an inference from (3) given that A is real symmetric, both of which stem from our analysis on the linear stability of the fully incoherent solution in the thermodynamic limit.To prove statement (5), consider statements (1) and (3) and an arbitrary combination of σ 1 σ 2 • • • σ d , where S = {i 1 , i 2 , . . ., i τ } ⊂ {1, 2, . . ., d} denotes dimensions that are incoherent, and S = {i τ +1 , i τ +2 , . . ., i d }, coherent.If we construct A initial in such a way that a rr < 0 for r ∈ S, while the reduced system that leaves out these dimensions, expressed by equation (23), is characterized by an A ′ that is rid of the corresponding rows and columns to r ∈ S. Then by ensuring A ′ to be positive definite, we know for r ∈ S there is coherence due to (1).

Onset of dimensional phase transition
To derive an analytical expression of the bifurcation point A critical = k critical A initial , we have adopted the method in [52] and set up a multi-variable Fourier formulation, which yields similar characteristic equations on the eigenmodes to that of the classic model, except the coupling strength that determines the stability of incoherence is now replaced by the diagonal elements of the matrix.The characteristic equations then predict that, should there be a phase transition towards synchronization on σ r , it would happen when the diagonal element a rr is scaled across a critical value of 2 π g(0) .Considering the diagonal elements are not necessarily identical, the bifurcation of each σ r understandably occurs at a distinct stage of the scaling, hence the differences in k critical for σ 1,2 in figure 3, and for σ 1,2,3 in figure 4. As k is adiabatically increased from zero, the dimension with the larger diagonal is earlier in reaching the bifurcation point.
Combining the above results, we have an important inference that predicts the system mode emerging at the end of the scaling experiment.Suppose the scaled coupling matrix A initial has unrepeated diagonal elements a i1 < a i2 < . . .< a i d with a i d > 0, and denote ÃI its principle submatrix defined by dimensions in set I. The d dimensions then have different critical values for bifurcation As k is adiabatically increased to k i d , dimension i d must be the first to bifurcate to the synchronized branch and thus initializes I = {i d }.Then at each critical value k ir larger than the previous, bifurcation on dimension i r will occur if and only if ÃI∪{ir} is positive definite, which determines whether to update I ← I ∪ {i r } or not.The procedure continues until {i 1 , i 2 , . . ., i d } is exhausted, and one obtains the system mode that takes 1 for i r ∈ I, and 0 otherwise, that will emerge at a sufficiently large value of k.For the above A initial , I also bears the description that, given I = m, ÃI has the largest diagonal sum among all of A initial 's mth order principle submatrices that are positive definite.In the insets, we compared the critical value of k from the data and that predicted by the self-consistent analysis and the stability analysis (marked with cyan diamond).The steady state solutions are obtained by solving equations ( 22) and (34) using Matlab, fsolve.For equation (34), the initial evaluation of the equations is set as σ1 = σ2 = 0.3 for k ⩽ 0.55, while it approaches the experiment data for k > 0.55.

Self-consistent analysis
The second question we are concerned with is the analytical description to the synchronization branch that may appear on any dimension with a second-order phase transition.Here, we employ a self-consistent analysis, wherein oscillators with fixed points on dimension r are subject to distinct velocity fields, and therefore need to be treated separately regarding their contributions to σ r (see appendix C for further details).
A comparison of the simulation result to the calculated synchronization branch is presented in figure 5. Notice the close fit between the analytically derived coherent σ r and the simulation data on the N = 5000 population.We have also marked up predictions of k critical from the steady state solution and the stability analysis in the insets.The theoretical predictions show the best agreement with the data when only one of the two dimensions is synchronizing, reducing the system to the classic Kuramoto model.When more than one dimensions is undergoing phase transition, both predictions of k critical slightly deviate from the data and also from each other.

Reiterate the standard Kuramoto model
Previously, as we motivate the matrix-coupled multidimensional model, there has been some interesting parallels between equation (2) and the standard Kuramoto model, particularly evident in the mathematical formulation and a first harmonic approximation derivation.This may be one of the first hints that some substantial parts of our theoretical results can be traced back to the theory developed on the Kuramoto model.As each oscillator has a one-dimensional phase θ i ∈ [0, 2π), the bistability demonstrated by the Kuramoto model corresponds to the two system modes where either the single dimension of the oscillators admits synchronization (σ 1 → 1), or remains desynchronized (σ 1 → 0).We know that the Kuramoto model (1) has no stable synchronized solution for K < 0, and the phase transition does not happen until the coupling strength goes across the threshold K = K c > 0. This is compatible with our discovery on the multidimensional case that the initial matrix A initial has to be positive definite for the full-synchronization solution to be stable, so that we observe and that a negative definite definite matrix, and K < 0, negative definite.In our interpretation, K is also the only diagonal element corresponding to the single dimension of all the oscillators, there is no surprise that it shares the same analytical expression at the critical point as the diagonals during the dimensional phase transition, i.e.K c = (a 11 ) critical = 2 π g(0) .What is fundamentally different with the multidimensional case d ⩾ 2 due to a matrix-valued coupling strength can be illustrated through the following example.We may consider a population of multidimensional phase oscillators θ i ∈ [0, 2π) d coupled with the conventional scalar K, namely This setup is equivalent to the proposed model ( 2) with a diagonal matrix A = diag{K, K, . . ., K} which is real-symmetric and, provided K = 0, invertible.A positive K then corresponds to a positive definite coupling matrix A, and a negative K renders A negative definite.According to our theory, the system (13) exhibits only the fully synchronized mode (11  2) that can only be induced by a matrix coupling explicitly addressing the inter-dimensional interactions.Recall from the two-dimensional example formulated by equation ( 3) that every inter-dimensional interaction is an implicit higher-order interaction among three variables, we safely infer that the matrix-induced higher-order interaction is the microscopic root for the coexistence of synchronization and desynchronization in the system modes.

Discussion
In this paper, we introduced a new high-dimensional generalization of the Kuramoto model, which gives rise to a novel macroscopic phenomenon where any combination of coherence/incoherence of the dimensions of the population is possible.Usage of the term 'coupling matrix' has been pervasive in previous literatures in network science [53][54][55], or even sometimes referring to the connectivity matrix in brain theory [56], but caution is due where it stands for an adjacency matrix as the cases above that just arranges the scalar-valued pairwise coupling strengths into a matricial form.An explicit matrix effect on the diffusive coupling between high-dimensional Kuramoto oscillators was rarely studied save for a few exceptions [26,57], not to say that in the representation of equation ( 2), to the best knowledge of the authors.
The key finding of this work is the multistable solutions known as the binary system modes, and their correspondence to the coupling matrix's algebraic properties.We have also identified the matrix-induced higher-order interactions as the microscopic root for the intermediate system modes where synchronization and desynchronization coexist by dimension.In deriving these results, we have several theoretical contributions: (1) for a d-dimensional population, we proved the existence of 2 d system modes under matrix coupling where each dimension of the population can be fully coherent or fully incoherent; (2) necessary and/or sufficient conditions are established between the 'positiveness' of the coupling matrix, or the lack thereof, and the system modes observed at the end of the scaling experiment; (3) steady state solutions of the order parameters are derived for arbitrarily many dimensions that tend towards synchronization through a second-order phase transition, along with a complete phase diagram for the two dimensional population; (4) analytical estimation of the bifurcation point on each dimension is obtained through a multi-variable Fourier analysis.
In the numerical study elaborated in the appendix, we discover an explosive synchronization/desynchronization across the dimensions of the population, when we try to show that the system modes are interchangeable through matrix manipulation.This is somewhat exceptional since it was believed that with Kuramoto-based interaction models, fast switching between the states must happen when there are either degree-correlated natural frequencies [58], a global feedback from the order parameter [25,44], or higher-order effects with explicit three-way or even four-way interactions [35][36][37].Meanwhile, our model admitted explosive synchronization/desynchronization without the above mechanisms, which may provide fresh insight into the understanding of such phenomena that have been related to epileptic seizures, bistable perception and other behaviors of the brain [59][60][61].Our experiment also suggests a simple way of controlling the macroscopic state to switch between pairs of system modes which may turn out effective with further investigations in applications.
The main limitation of this work is that, the experiments and the analysis are performed under the premise that the coupling matrix is real symmetric.Although the self-consistent analysis and the multi-variable Fourier analysis do not touch on this condition and remain applicable even when the coupling matrix is asymmetric, the linear stability condition of the fully coherent solution is indeed based on this property, so do all the subsequent conclusions on the relationship between the system mode and the coupling matrix.This limitation stems from the block structure of the Jacobian matrix ( 14) that complicates the analysis of its spectrum distribution.Existing results often require strong mathematical condition on the blocks or the producted matrices, while many of them do not extend to ensure the stability of the overall matrix.Another limitation to consider is the assumption made about the natural frequency vector ω i = ωi 1 d .Such a form does not affect the linear stability analysis with the Jacobian matrix, but it does play a role in the calculation of the coherent branch and the multi-variable Fourier analysis.Should the components of ω i also follow some random distribution, we would anticipate a distinct and considerably more complicated procedure to calculate the dimensional order parameters.
The proposed model is distinct from other network dynamics by the unique macroscopic phenomena it displays.For the majority of the variations of the Kuramoto model, including the multi-layer [25], the higher-order [35,37], and the high-dimensional [23,24,44] ones, great effort was dedicated to the synchronization phenomenon and the nuance in the phase transition toward it, which is also a main focus of this work.But in terms of desynchronization that is indeed discussed on a high-dimensional population, a direct comparison to figure 3 is difficult to draw when [24,62] defined the population on the unit sphere.If we loosely apply the idea to a population on the 3D unit sphere that one angular variable is synchronized while the other is not, in analogy to mode 10 or mode 01 in this work, then as in [24], all of these incoherent states, together with an infinite continuum of distributions whose centroid is the origin, collapse into a single representation of a zero-magnitude order parameter.According to the precessing equation d in this case and the distribution G(W i ), the aforementioned system modes are not steady state distributions that can be observed in their system.Another idea reminiscent of our discovery is the isolated desynchronization in cluster synchronization studied on a rather general form of the complex network [63][64][65], where a population of high-dimensional oscillators is divided into a number of subsets that each evolves along a unique, synchronized trajectory.On occasion, as the mentioned works proved, synchronized solutions may lose stability for some of these subsets while others carry on, which can be detected by identifying the symmetries in the network.Although this description bears certain similarity to our phenomenon, a rigorous comparison reveals the difference considering the population studied in this work is also high-dimensional.Thus by their definition, no cluster is formed for any of the system modes with one of the dimensional order parameters being zero, as it indicates desynchronization for the entire population.
Apart from potentially contributing to the understanding of biological swarms as mentioned in the Introduction and the Motivating Scenario, further applications of our results may be found, as we believe, in the field of neuroscience where a desynchronized state in part of an ensemble is particularly meaningful.A notable example is the unihemispheric slow-wave sleep observed in birds and aquatic mammals including the dolphin [66,67], where the two hemispheres of the brain alternate between the resting state (synchronized state) and the activated state (desynchronized state), the two behaviors existing independently at the same time.The analogy provided by our numerical experiment is that, since the synchronization and desynchronization of the two dimensions are separated, their combinations suggest four qualitatively distinct modes of the system that are possible to switch between one and another through manipulating the matrix elements, which may facilitate our understanding of the mentioned phenomenon.On the other hand, in terms of the information storage and processing of oscillatory neural networks (ONN) [68,69] that was introduced based on the Kuramoto model [70], our finding suggests a new way to reconcile the phase-relation logic of the ONN with the boolean logic adopted by the traditional von Neumann machine, where the former encodes information with the in-phase (logic '0') or anti-phase (logic '1') distributions of the oscillators.Here we have shown that through a set of statistical measures, the dimensions as digits naturally carry boolean information.

Appendix B. Linear stability of the coherent solution
Let us first present the theoretical ground on statement (1).Before, we made the observation that when the weight matrix A is real symmetric, a full synchronization emerges at a positively distributed spectrum of A. This has been observed in oscillators of dimension two as well as dimension three.By linearizing the whole system around the coherent solution (σ r → 1, r = 1, . . ., d), we will show that this correspondence of λ(A) to the order parameters is readily generalized into arbitrary dimension d.To see this, consider the Nd dimensional dynamic system (2).To perform the linear stability analysis for the entire system, we would like to evaluate the Jacobian matrix with respect to the state variable θ in a closed form.This is realized by identifying with i = j ∈ {1, . . ., N}, r = s ∈ {1, . . ., d}, and a rr , a rs are the elements of the coupling matrix A. Then by sequencing θ dimension-wise as [θ 11 , θ 21 , . . ., θ N1 , . . ., θ 1d , θ 2d , . . ., θ Nd ], the Jacobian matrix can be written as where L r , r ∈ {1, . . ., d} are the graph Laplacians of networks G r for the rth dimension of the population.Because of ( 2), the networks G r = (V r , E r , W r ) are (i) complete, (ii) undirected, and (iii) weighed by cos(θ jr − θ ir ) for any edge {ir, jr} ∈ E r .Formalizing this with the incidence matrix H for the complete graph, there is When θ jr − θ ir < π 2 , G r , r ∈ {1, . . ., d}, are weighed by strictly positive scalars that render the Laplacians L r positive semi-definite, with nullity(L r ) = 1. Let we denote the eigenvalues of A ∈ R d×d as λ 1 , λ 2 , . . ., λ d , and that of Q as υ 1 , υ 2 , . . ., υ Nd , while the Jacobian has eigenvalues µ 1 , µ 2 , . . ., µ Nd .We will first look into some of the properties of the eigenvalues of P, Q respectively, before studying their combined effect on the product of P, Q, which leads to the distribution of µ 1 , µ 2 , . . ., µ Nd and an explanation of the independent dimensional transition to synchronization.
We first notice that for P = A ⊗ I N , its eigenvalues are λ 1 , λ it means that the condition that θ * i exists is where • ∞ denotes the infinity norm for vectors.Note that Σ −1 exists because we are going to evaluate the stability of the coherent solution.In our experiment, the natural frequency vector ω i is set to be ωi 1 d , and the coupling matrix A is scaled from an initial matrix A u by a factor k. Therefore ( 15) is also One sees that with finite N, as long as k is large enough, the fixed point θ * i always exists.Moreover, for oscillators θ i and θ j , Again, for large enough k there will be |θ * ir − θ * jr | < π 2 for arbitrary i, j ∈ {1, . . ., N}, i = j, and for arbitrary dimension r ∈ {1, . . ., d}.Therefore, the graph Laplacians L r , r ∈ {1, . . ., d} are positive semi-definite, so is the block diagonal matrix Q.We now have a simple bound for υ 1 , υ 2 , . . ., υ Nd that is υ ⩾ 0.
For the Jacobian matrix J(θ) = − 1 N PQ whose eigenvalues are denoted as µ 1 , µ 2 , . . ., µ Nd , we mention that these eigenvalues are real when P, Q satisfy the above conditions.This is easy to prove considering PQ and QP have the same non-zero eigenvalues.Suppose u is an eigenvector of QP with respect to eigenvalue Then for µ ′ i = 0, since P is Hermitian, u, Pu is real, we have that µ ′ i , i = 1, . . ., dN − d are real, so are µ 1 , µ 2 , . . ., µ Nd .To further specify the distribution of the Jacobian eigenvalues, it is necessary to introduce the following lemma from [71, theorem 1].

Lemma 1. Let A be a positive definite or semidefinite Hermitian matrix of order n with the smallest eigenvalue m and the largest eigenvalue M. Let B be a Hermitian matrix of order n with eigenvalues
Then AB has real eigenvalues Λ 1 ⩽ . . .⩽ Λ n , and it holds for suitable factors Θ v , which lie between m and M, that Since both P and Q are square, PQ actually has exactly the same eigenvalues as QP, for which Q is positive semi-definite, P is Hermitian and invertible.Let the eigenvalues of PQ and QP be sequenced as Nd } where each λ ′ can be either positive or negative, lemma 1 suggests that for each µ ′ v , there exist a λ ′ v of the same order and a Θ v that has 0 ⩽ Θ v ⩽ υ max , such that Also, since λ ′ = 0, there are d times that Θ v takes the value zero.We then conclude that Θ v = 0 only happens when (i) |λ ′ v | is the smallest, if all those in λ ′ have the same sign, or when (ii) λ ′ v is the smallest positive eigenvalue or the largest negative eigenvalue, if λ ′ has both positive and negative elements.Since we are dealing with a large population, it is reasonable to assume that N d.Consider (i) where 1 and Θ N+1 = 0, then µ ′ N+1 = 0. Since there exists j ∈ {1, 2, . . ., N} where Θ j = 0, λ j = λ 1 > 0, there is µ ′ N+1 < µ ′ j which is a contradiction.The same applies when λ ′ are all negative and when they can be both.Equation ( 18) then captures the full spectrum of the Jacobian at the coherent equilibrium.
For the nonzero eigenvalues of QP and PQ, since Θ v > 0, the sign of µ ′ v is completely determined by the sign of λ ′ v , which is just the eigenvalue of A. We can then sum up the above discussion into the following theorem.
Theorem 1.Consider the Jacobian matrix in (14) where A is invertible, real-symmetric, and blkdiag{L 1 , . . ., L d } is positive semidefinite.Then the nonzero eigenvalues of the Jacobian are of the same sign if and only if the eigenvalues of A are of the same sign.
Under the same premise, combined with equation (18), it is then derived that the fully coherent solution is linearly stable if and only if A is positive definite.In contrast, should the coupling matrix had any negative eigenvalue, the system would have at least one dimension that could not synchronize.The experiment goes further to exhibit that this partly incoherent state is actually a fully incoherent solution for some of the dimensions, but its opposite for the rest.
It is then readily inferred that if the system exhibits at the end of the transition A = kA initial , then the principle submatrix consisting dimensions {i τ +1 , i τ +2 , . . ., i d } of A must be positive definite.Because the solution 0 0 is essentially about the fully coherent dimensions {i τ +1 , i τ +2 , . . ., i d } that decouple from the fully incoherent ones {i 1 , i 2 , . . ., i τ }, we can set σ i1 , σ i2 , • • • , σ iτ as zero and analyze the reduced system from (12) where θ i is also rid of the {i 1 , i 2 , . . ., i τ } dimensions.If the principle submatrix is not positive definite, the fully coherent solution on the reduced system will in turn be unstable, then mode 0 0 is unstable on the original system.

Appendix C. The coherent branch and the onset of synchronization
Despite the ultimate system modes we can now relate to the properties of the coupling matrix, it remains a challenge to analytically derive their phase diagrams where one or more dimensions may tend to synchronization.
In one of our experiments that is not shown here, it was observed that with A = kA initial , the transition in dimension {θ i1 } or {θ i2 } experiences several minor jumps at lower values of k for even populations as large as N = 1000, before gaining coherence consistently with each increased k.Also, due to the limit size effect, one needs to tune k slightly below zero to see the onset of synchronization.Yet for N = 5000, the transitions of σ 1,2 are smooth enough to show prospect of fitting analytical results.These transitions are monotone, have a clear bifurcation at a critical k above zero, and are in all comparable to that of the second-order.For this reason, we seek to tackle the σ(A) solution in the continuum limit N → ∞, where the order parameters are in turn r ∈ {1, . . ., d}.The distribution function n r (θ r ) can be broken into two parts depending on if its represented population are fixed on the rth dimension; given that the natural frequencies are identical on every dimension, we write To specify the distribution of the oscillators on each dimension, we now employ a self-consistent analysis where it is assumed that the order parameters ρ r = σ r e iΨr , r ∈ {1, . . ., d} are fixed, therefore the drifting oscillators -oscillator without a fixed point on all of its dimensions -must form a stationary distribution in the space [0, 2π) d .We now have a density function f(θ 1 , . . ., θ d , ω) = f(θ, ω) that is independent of time.Also assume that Ψ r coincides with the symmetry center of g(ω), this gives ρ r = σ r , and From here, it is essential to separate the treatment of each system mode where there may be σ r = 0, because it affects our way of obtaining the solution of θr = 0 from ( 21) and our interpretation of the solution.Now we have divided the population into two categories, the fully synchronized oscillators and the drifting ones, then system modes with σ r = 0 automatically excludes the existence of the former.But the experiments producing results like σ 1 = 0, σ 2 = 1 validate that, the oscillators being drifting does not mean they have zero contribution to any of the order parameters, i.e. they are impossible to form certain degree of coherence on any dimension.Moreover, as we shall see, even for system modes as σ 1 = 1, σ 2 = 1, the drifting oscillators demand further division to yield an exact theoretic prediction.We will first give the calculation of the order parameters on the two dimensional population, then discuss our treatment on arbitrary dimensions.When one of the order parameters remains zero due to the coupling matrix and the subsequent instability of the coherent branch, the other dimension in fact degenerates to the classic Kuramoto model, e.g. when σ 2 = 0, (11) turns into θi1 = ω i1 − a 11 σ 1 sin (θ i1 − Ψ 1 ) , except the classic coupling strength K is now substituted by the matrix element a 11 .There is then naturally as the expression of the coherent branch of σ 1 [52].As we scale up A 2 with k > 0, a 11 is strengthened in the manner of the classic K, thus it is not hard to expect that the transition the order parameter goes through conforms with the classic case.This simple example also suggests that for population with higher dimensions, we need to leave out the dimensions with σ r = 0 for a given system mode, and analyze the reduced model which, compared with (21), involves an A ′ and a Σ ′ of dimension d ′ < d, and A ′ is obtained by taking out the corresponding rows and columns of A to σ r = 0.The analysis of the reduced system can then build much on the following calculation of mode 11.
During the transitions of σ 1 , σ 2 from 0 to 1, there are maximally three types of oscillators with different contributions to the order parameter -the fully synchronized, the orbiting, and the fully incoherent oscillators.For the first category, consider the model ( 21) and solve for θr = 0 for r ∈ {1, 2} simultaneously, we would have on each dimension a synchronization domain, defined as which is always symmetric about the origin despite the value of r.But for the oscillators to have fixed points on both dimensions, the natural frequency must satisfy In fact, given that θr = ω − N j =1 a rj σ j sin(θ j ) where σ j = 0, θ can either be stationary on all of its dimensions or none of its dimensions, and any oscillator with ω ∈ R\D is considered a drifting oscillator.For ω ∈ D, we can now derive the fixed point on each dimension, Since for the entrained population, put ( 20), (26), and ( 27) in (19), there is where γ = arcsin( ), sin ∆ r = ω σr A −1 1 2 r .For population with ω ∈ R\D, the natural frequency exceeds at least one of the synchronization domains D r , and under all circumstances, the oscillators are relatively in motion to the mean field.However, it does not imply that there is no solution for θr = 0.Here we propose another domain which we refer to as the original domain on the rth dimension, defined as The original domain depicts the most general condition where θr = 0 has a solution.Then, by definition, there is We also denote The contribution of the drifting population to σ r distinguishes between D o \D and R\D o , where we call the former orbiting oscillators, because when projected onto the unit sphere, they are constantly drifting on a circle along the z-axis that rotates azimuthally for a small angle after each period (when N is relatively small, as N = 1000, and the θ 1 mean field is not stationary).We expect this circle, or the orbit, to be slightly oscillatory about a fixed plane when N → ∞.Oscillators with ω ∈ R\D o constitute the fully incoherent population, we first treat this simpler case where the continuity equation (40) needs to be evaluated with ∂f/∂t = 0, this requires where • 2 is the Euclidean norm, and C(ω, A, Σ) is the normalization constant that guarantees ˆπ −π ˆπ −π f (θ, ω) dθ 1 dθ 2 = 1 for every ω ∈ R\D o .Because each D o r is symmetric about the origin, so is R\D o , the constant should be invariant against ω → −ω.We can then write n fi r (θ r ) and plug it into (19) so that the contribution of the fully incoherent oscillators is Without loss of generality, for r = 1, substituting (31) becomes where sin(θ ′ ) = − sin(θ ′ 1 ) sin(θ 2 ) T .Notice how under the 2-norm of the denominator, each transformation θ 2 → −θ 2 or ω → −ω in the first term is a change of the sign before sin(θ 2 ) or ω, and this arrangement can be found in the second term.With C(ω, A, Σ)g(ω) being invariant to ω → −ω, the two terms in (32) are equal, and the contribution of the fully incoherent oscillators to σ r is actually zero.Since to deal with the orbiting oscillators, the relationship between different domains needs to be specific, we use our experiment with initial matrix A 1 as an example, where the data shows , where we each sampled two oscillators and tracked their velocities over time, as displayed in figure 6.We did this because ostensibly, when ω ∈ D o 2 \D, both θ1 = 0 and θ2 = 0 are viable, but then it is undecided which of the two dimensions should be fixed for such oscillators.Experiment shows in figure 6 that in this case, there is θ1 = 0 instead of θ2 = 0.This is sensible because ω is well within the margin D o 1 , and also, if there are more than two dimensions to the oscillators, it is not reasonable to randomly pick one from θ1 = 0, . . ., θd = 0 to be unsolvable.For the rest ω ∈ D o \D o 2 , there is obviously θ1 = 0, θ2 = 0. Based on the above analysis, where θ * 1 = arcsin( 1 a11σ1 (ω − a 12 σ 2 sin θ 2 )) and C(ω, A, Σ) is defined correspondingly.Though θ * 1 is not stationary, it guarantees that θ1 is close to zero (as is observed from the experiment), and it allows θ 1 to oscillate about arcsin( ω a11σ1 ), for which θ 1 is stationary when averaged over time.Therefore the contribution of the orbiting oscillators to σ 1 is where 1 \D o 2 exhibit similar behaviors, except their desynchronized θi2 never really coincide with that of the entrained population.θ281 and θ958, however, experience a bottleneck of the θi2 velocity when they pass through the synchronized population as if they are themselves entrained, but will eventually move on.and the last line of (33) stands because the imaginary part Combining our analysis of ( 28), ( 31) and (33), the order parameters for mode 11 are derived as In fact, for the d dimensional population, the calculation of the synchronization branches follows similar procedures as illustrated in figure 7. We mentioned that the dimensions with σ r = 0 must be left out for a given system mode, and the reduced system ( 23 25) and ( 29) for the d-dimensional case are plotted on the |ω| axis, with an assumed relative order.In calculating σ r ′ , the population with ω ∈ D o r ′ (in the bracket) is fixed on dimension r ′ , and therefore contributing to σ r ′ despite that most of these oscillators are unfixed on some other dimensions thus overall drifting.Each segment on the axis within D o r ′ corresponds to oscillators with a unique arrangement of their dimensions being fixed or unfixed, therefore corresponds to its own velocity field and density distribution.Integral (37) is then considered for all the relevant frequency segments to obtain σ r ′ .Note that f 0 here represents equation ( 27) in the d-dimensional case.
where F r is obtained by solving Ār ′ d is the principle sub matrix of A containing the r ′ th to the dth dimension.Since A is positive definite, all of its principle minors are positive, which means Ār ′ d should be invertible, and therefore, the expression in ( 35) is uniquely obtained.All fixed points are essentially the linear combination of ω, sin θ 1 , sin θ 2 , . . ., sin θ r ′ −1 .We then write the density function for where the velocity θr is also a linear combination of ω, sin θ 1 , sin θ 2 , . . ., sin θ r ′ −1 .Then the contribution of 1/2 = 0, for the same reason as (32).This means σ o r should only take into account frequency intervals that allow the oscillators to be fixed on dimension r.And for r > r ′ − 1, the contribution on the above interval is Summing (37) over r ′ ⩽ r, we arrive at the expression of σ r in terms of the orbiting oscillators σ o r .Combining σ o r with σ fi r = 0 and we obtain the full equation of the synchronization branch of σ r with arbitrary dimensional population.

Appendix D. Linear stability of the incoherent solution
Here we are interested in the stability of the state where the oscillators are incoherent on all dimensions, i.e. σ 1 = σ 2 = 0. We aim to derive an explicit estimation of the critical value of the coupling matrix from which the system starts the transition towards coherence on each dimension, under the presumption that such a transition indeed occurs.Assume that N → ∞, we consider the population as a flow on a two-dimensional plane with 2π-periodic bounds.For our purpose, the flow has a perturbed density function where ϵ is small, and f(θ 1 , θ 2 , ω, t)dθ 1 dθ 2 denotes the fraction of oscillators with natural frequency ω on both dimensions, that are found in the infinitesimal area Following the method of [52], given that the oscillators are conserved, we set up the continuity equation where the velocity field is that of equation (11).Inserting ( 11) and ( 38) into (40) gives To perform Fourier analysis, we expand η(θ 1 , θ 2 , ω) into Fourier series Note that for real-valued functions, c m,n (ω, t) = c * −m,−n (ω, t).As we evaluate the order parameter σ 1 e iΨ1 for N → ∞, it is noted that only term c −1,0 (ω, t)e −iθ1 + c.c. contributes, c.c. representing the complex conjugate of the previous term.This similarly applies to where harmonics c 0,−1 (ω, t)e −iθ2 + c.c. alone contribute to the order parameter.It is further derived that and We separate the terms that governs the evolution of η(θ 1 , θ 2 , ω, t) in ( 41) from the others, where η stands for the rest of the harmonics.Equation (39) implies that c 0,0 = 0, therefore η 0 can be considered the fundamental mode of η where the disturbances on the θ 1 dimension and the θ 2 dimension are decoupled.Inserting ( 43)-( 47) into ( 41) and linearizing result in where c(ω) = ´+∞ −∞ c(ω)g(ω)dω.Solving (48) with respect to e iθ1 and e iθ2 , we derive and To determine the stability of solution (38), set c 1,0 (ω, t) = b 1,0 (ω)e κ1t , c 0,1 (ω, t) = b 0,1 (ω)e κ2t , it is obvious that ( 49) and ( 50) have the same discrete spectrum as that of the classic model derived in [72], which means κ 1 , κ 2 are derived from the characteristic Equations In this work we have adopted the standard Lorentzian distribution with g(ω) = 1 π(ω 2 +1) , this leaves the critical value of a 11 (a 22 , resp.) at which the effect of η(θ 1 , θ 2 , ω, t) on the θ 1 (θ 2 , resp.) dimension renders (38) unstable to be Similarly, one will find that only the first harmonics contribute to η in the sense of (40) Recall that in the numerical experiment, we scale all the elements of A by a real number k and increase it slightly at each step.The coupling matrix A is required to be real symmetric, but the diagonal elements are not necessarily identical.This means that whichever is positive and larger will reach its critical value first, and lead to the phase transition to synchronization on its corresponding dimension.Equations ( 51) and ( 52) also suggest that if a diagonal element a rr is negative, the characteristic equation on κ r will be unsolvable for a discrete spectrum, leaving only the continuous spectrum on the imaginary line [52], which means σ r = 0 remains neutrally stable against the disturbance from η on the rth dimension for the entire time.As a result, when a rr < 0, the rth dimension of the population will not see a transition toward synchronization with an increasing positive k.

Appendix E. Mode switching
The dimensional phase transitions in the main text are induced by adiabatically scaling an initial matrix A initial with k ⩾ 0.Here we show another method to modulate the elements of the coupling matrix.The purpose is to answer, with a case study, if the 'system modes' exhibited as our main finding are possible to switch between one another purely through manipulating the coupling matrix.
Recall the proposed model and the system modes encoded into binary digits 00, 01, 10, 11 encouraged by the combination of σ 1 σ 2 .We now find representative matrices M 1 , M 2 , M 3 , M 4 for these modes, which means A = M i in (54) along with a proper initial value settles the system into σ 1 σ 2 = 00/01/10/11.By presenting the following method on adjusting the representative matrix, we demonstrate that the mode of the system also switches correspondingly, and introduce some of the many open problems to be studied of this model.
For the experiment, we realize the switchings on N = 100 oscillators whose natural frequencies ω i1 = ω i2 = ω i are drawn from the standard Lorentzian distribution.Every switching process indicated by (a)∼(l) in figure 8 starts with the numerical integration of system (54) from θ i (0) = 0 and A = M i , before its transition to A = M j which is broken down into 51 steps, i, j ∈ {1, 2, 3, 4} and can be specified by so that M 0 = M i = a i 11 a i .The idea is to keep the change of matricial elements gradual with the incremental matrix Λ, and therefore an effective way of compensating the differences between M i and M j is to make a linear interpolation, with Λ = (M j − M i )/50.At each step, an average of the steady state order parameter σ 1 (σ 2 ) over time is then evaluated as a data point in figure 9. However just as there are innumerous ways to select the four representative matrices, it is beyond the scope of this work to study the switching diagram for all possible configurations.Actually, we will show one of the cases where the interpolation scheme fails to accomplish all twelve switchings.The purpose is to suggest the fluidity of the 'modes' we established as qualitatively different states of the system, the interchangeability granted by the few extra degrees of freedom from the matricial coupling.Three examples are demonstrated in figure 9 for their switching processes under (55), with both asymmetric and symmetric matrices.Specifically, the representative matrices for (i) are M 1 = −13.13    We observe that the assigned representative matrices fulfills all twelve switchings for case (i) in figure 9. However for (ii.d) which is from 10 to 01, the switching does not happen with direct interpolation between M 3 = 16 30 30 −4 and M 2 = 16 30 30 20 , and the symmetry is broken between (ii.d) and (ii.c); actually, one needs to extend the increasing of a 22 till around 31 to induce the switching.We take a different route in (iii.d) by first evenly decreasing a 11 while increasing a 22 of M 3 to 0 30 30 8 in 26 steps, then increasing evenly both a 11 and a 22 in the remaining steps to M 2 which has allowed the desired switching from 10 to 01. (i) and (iii) verify that, as the representative matrix M i morphs into M j in the twelve scenarios, there are cases where the system completes a switching from M i 's corresponding mode to that of M j 's with either continuous transitions or explosive transitions to synchronization/desynchronization.As we pointed out in Discussion, previous discoveries of the explosive phase transitions are accompanied by modifications of the Kuramoto model such as degree-correlated natural frequencies, a global feedback from the order parameter, or higher-order effects whose mathematical formulation contains explicit three-way or higher nonlinear terms.In contrast, the model we studied does not explicitly involve these traits.A preliminary explanation is that the matrix coupling mechanism has indeed incorporated higher-order interactions which give rise to the explosive dimensional phase transitions, but the inter-dimensional three-way interaction appears as a two-way interaction term in the model as we commented about figure 1. Future investigations can then focus on three problems, namely necessary and sufficient stability conditions of the system modes, an adjustment scheme of the coupling matrix A that guarantees the switching of all pairs of modes, and the conditions that guarantee the first-order phase transitions and the second-order ones.

Figure 1 .
Figure 1.The matrix coupling on a simple network.Here we illustrate how the three oscillators and their two dimensions are interacting from the perspective of variable θ11.Gray, solid lines indicate pairwise interactions between components of the same dimension, which take place on the 1-simplex; inter-dimensional interactions that θ11 is involved in are indicated with triangles regarded as the 2-simplex.The first-row elements of A weigh the influences from these interactions on θ11.

Figure 3 .
Figure3.Visualization of the system modes on a two-dimensional torus.Take θ i1 as the toroidal angle and θ i2 as the poloidal angle, the population of oscillators (black circles) are projected onto the surface of a 2-torus with major radius Ra = 10 and minor radius R b = 6.4.We illustrate with 100 oscillators from the original ensemble N = 5000 whose initial conditions were uniformly distributed.The four panels correspond to the four modes the system eventually settled in with the initial matrices A1, A2, A3, A4 at k = 10.Panel (a) displays mode 11 where almost all oscillators are fully synchronized (in frequency) for both elements θ i1 and θ i2 ; the oscillators remain relatively static and travel at a constant, but minuscule speed of the mean field.Panel (b) shows, at mode 10, the oscillators are only synchronized on the θ i1 element and stay disarranged on a ring that moves slowly on the toroidal direction.Panel (c) corresponds to mode 01 where the population is only synchronized on the θ i2 element and is distributed on a ring that shrinks or stretches depending on the common poloidal angle.Panel (d) demonstrates how the population is incoherent on both dimensions at mode 00, where oscillators drift across the surface and do not form any groups or show any particular pattern.

Figure 4 .
Figure 4. Transition of σ1,2,3 for the three-dimensional oscillators and the initial value experiment.Left: panels (a)∼(h) are simulation results of equation (3) on N = 1000 three-dimensional oscillators with coupling matrices

Figure 5 .
Figure5.The steady state solutions and simulation data of dimensional phase transitions.The four panels demonstrate the order parameter magnitudes σ1,2 undergoing phase transitions under mode 11, mode 10, and mode 01.In the insets, we compared the critical value of k from the data and that predicted by the self-consistent analysis and the stability analysis (marked with cyan diamond).The steady state solutions are obtained by solving equations (22) and (34) using Matlab, fsolve.For equation(34), the initial evaluation of the equations is set as σ1 = σ2 = 0.3 for k ⩽ 0.55, while it approaches the experiment data for k > 0.55.

Figure 7 .
Figure 7. Schematic illustration of the calculation of the dimensional order parameters.In this figure, the domains defined by equations (25) and (29) for the d-dimensional case are plotted on the |ω| axis, with an assumed relative order.In calculating σ r ′ , the population with ω ∈ D o r ′ (in the bracket) is fixed on dimension r ′ , and therefore contributing to σ r ′ despite that most of these oscillators are unfixed on some other dimensions thus overall drifting.Each segment on the axis within D o r ′ corresponds to oscillators with a unique arrangement of their dimensions being fixed or unfixed, therefore corresponds to its own velocity field and density distribution.Integral(37) is then considered for all the relevant frequency segments to obtain σ r ′ .Note that f 0 here represents equation(27) in the d-dimensional case.

Figure 8 .
Figure 8. Switching diagram of σ1σ2.The diagram illustrates the combination of σ1 and σ2 and their interchanging processes indexed from (a) to (l).

Figure 9 .
Figure 9. Mode switching based on direct interpolation of matricial elements.The order parameters σ1,2 are plotted against the interpolation steps s. (a)∼(l) correspond to the switching processes marked in figure 8.
A = a 11 a 12 a 21 a 22 ∈ R 2×2 , that acts on the sum of vectors [sin(θ j1− θ i1 ) sin(θ j2 − θ i2 )] T .For clarity, we consider a simple connected network of N = 3 as demonstrated in figure1.As we break down the expression of θ11 to see what information θ 11 uses to update its state, there is θ11 = ω 11 + a 11 3 [sin (θ 21 − θ 11 ) + sin (θ 31 − θ 11 )] + a 12 3 [sin (θ 22 − θ 12 ) + sin (θ 32 − θ 12 )] , (4) which involves variables from the same dimension with θ 11 that are θ 21 , θ 31 , and those from the other dimension, θ 12 , θ 22 , θ 32 .We observe that θ 11 forms an explicit pairwise interaction with the former, namely a 11 3 [sin (θ 21 − θ 11 ) + sin (θ 31 − θ 11 )] if a multivariate function F : R d → C meets F( • • • 1) d and the fully desynchronized mode (00 • • • 0) d out of the 2 dpossibilities.This implies the absence of co-existence of synchronization and desynchronization across different dimensions in terms of system modes.The rationale behind this is that with the scalar K, equation (13) in fact depicts layers of one-dimensional Kuramoto oscillators θ ir , r ∈ d evolving independently of each other, where there is no interaction or coupling across the dimensions of the population.It can then be concluded that, compared with the standard Kuramoto model, the intermediate system modes (those except (11 • • • 1) d and (00 • • • 0) d ) are what is truly new about equation (

Table 2 .
Mean values of the [σ1, σ2] that are considered to have converged to the system mode on the left column.

Table 3 .
Standard deviation of the [σ1, σ2] that are considered to have converged to the system mode on the left column.