Volcano transition in a system of generalized Kuramoto oscillators with random frustrated interactions

In a system of heterogeneous (Abelian) Kuramoto oscillators with random or ‘frustrated’ interactions, transitions from states of incoherence to partial synchronization were observed. These so-called volcano transitions are characterized by a change in the shape of a local field distribution and were discussed in connection with an oscillator glass. In this paper, we consider a different class of oscillators, namely a system of (non-Abelian) SU(2)-Lohe oscillators that can also be defined on the 3-sphere, i.e. an oscillator is generalized to be defined as a unit vector in four-dimensional Euclidean space. We demonstrate that such higher-dimensional Kuramoto models with reciprocal and nonreciprocal random interactions represented by a low-rank matrix exhibit a volcano transition as well. We determine the critical coupling strength at which a volcano-like transition occurs, employing an Ott–Antonsen ansatz. Numerical simulations provide additional validations of our analytical findings and reveal the differences in observable collective dynamics prior to and following the transition. Furthermore, we show that a system of unit 3-vector oscillators on the 2-sphere does not possess a volcano transition.


Introduction
Numerous studies have explored emergent phenomena in complex systems across various disciplines [1].An important type of system that has received significant attention during ‡ Author to whom any correspondence should be addressed arXiv:2401.09081v1[nlin.AO] 17 Jan 2024 the last decades is coupled oscillators.Among them are large ensembles of attractively coupled oscillators that play a pivotal role in understanding transitions from incoherent states to partial synchronization [2,3].Given the similarity between the synchronization transition and ferromagnetic phase transitions [4,5], Daido introduced a model involving Kuramoto phase oscillators with heterogenous natural frequencies and random frustrated interactions, featuring both attractive and repulsive couplings between oscillators [6,7].It was conjectured that the resulting emergent phenomenon in this oscillator system might share characteristics with spin glasses [8,9], motivating to coin this peculiar state an oscillator glass.Yet, based on further examinations, the definition of an oscillator glass has been a subject of ongoing debate and dispute [10][11][12][13][14][15].
In 2018, Ottino-Löffler and Strogatz put forward a mathematically tractable model, similar to Daido's approach [10].This model of heterogeneous Kuramoto phase oscillators was formulated with a symmetric, random interaction matrix with a tunable rank.Note that the interaction matrix of the model considered by Daido had a fixed, full rank.This system featured a so-called volcano transition, which Daido suggested as indicative of an oscillator glass [6].Above the volcano transition, an incoherent state loses its stability, and a so-called volcano phase, which is characterized by partial synchronization in a system with random frustrated interactions, becomes stable and thus observable.However, the volcano phase as a partially locked state with random frustrated interactions is not characterized by the usual Kuramoto order parameter, but rather requires measuring a local degree of synchrony, i.e., a local field defined below, for each oscillator.The volcano transition involves a change in the shape of the local field distribution, shifting from being concave downwards at the origin to concave upwards, i.e., a distribution of local fields that resembles a volcano shape.In contrast to Daido's conjecture, explorations of the model with low-ranked, random interaction matrices revealed that the partial synchrony beyond the volcano transition does not exhibit characteristics of an oscillator glass.For example, the volcano phase does not exhibit an algebraic relaxation of the Kuramoto order parameter, as expected for a glassy state, but rather displays an exponential relaxation [10].Thus, the relationship between an oscillator glass and a volcano phase continues to be a topic of ongoing discussion and requires further, thorough explorations.
Realizing the significance of both nonreciprocal and random interactions for investigations related to real-world scenarios, Pazó and Gallego introduced a controllable symmetry/asymmetry parameter within the random interaction matrix to the model of Ottino-Löffler and Strogatz [11].Thus, the random, frustrated interactions are rendered nonreciprocal by this additional parameter.Employing such random interactions, a volcano transition was observed and its relation to the symmetry-controlling parameter was discussed (see Sec. 2 for details).In both of the aforementioned models [10,11], the authors derived an analytical formula to identify the critical point at which a volcano transition takes place in the thermodynamic limit.This mathematical tractability relied on the so-called Ott-Antonsen (OA) ansatz, which can be invoked due to the sinusoidal coupling of the oscillators [16][17][18][19][20]. Within the framework of the OA ansatz, the macroscopic dynamics of a system of heterogeneous Kuramoto phase oscillators can be described in terms of the so-called OA variables in the dimension-reduced manifold.Investigating the stability of incoherent states in terms of the OA variables, this approach enables the prediction of a critical point at which the volcano transition occurs.
Parallel to the studies mentioned above, the collective dynamics of ensembles of socalled generalized Kuramoto oscillators have attracted considerable attention [21][22][23][24][25][26][27][28][29][30].In particular, Lohe detailed the dynamics of generalized Kuramoto models on a compact classical Lie group and provided insights into their notable characteristics [31,32].Lohe's generalized Kuramoto oscillators can be understood as those defined on the (D − 1)sphere, with each oscillator represented as a unit vector in a D-dimensional Euclidean space.In contrast, standard Kuramoto phase oscillators are typically defined on the unit circle, with each oscillator represented as a unit vector embedded in the 2D plane (see Sec. 2 for details).Meanwhile, Barioni and Aguiar formulated the OA ansatz for higherdimensional Kuramoto oscillators [33,34].The authors proposed an expansion using higherdimensional spherical harmonics for the oscillator density function in the thermodynamic limit, mirroring Ott and Antonsen's use of a Fourier series expansion for standard Kuramoto oscillators [16][17][18][19][35][36][37].This ansatz yields the OA equation in a vector format, allowing for the description of the system's macroscopic dynamics in the dimension-reduced manifold.
In this paper, we study a higher-dimensional generalization of the aforementioned models.In Sec. 2, we first revisit the model proposed by Pazó et al., thereby introducing the random, frustrated interactions.Then, we introduce a system of heterogeneous SU(2)-Lohe oscillators, coupled through reciprocal and nonreciprocal frustrated interactions.These oscillators are parameterized by unit 4-vectors, which allows us to investigate an ensemble of heterogeneous unit 4-vector oscillators on the 3-sphere and thus to identify a volcano transition.In Sec. 3, we apply the higher-dimensional OA ansatz to this system, from which we obtain analytical predictions for the critical coupling strength as a function of the symmetry parameter.Section 4 is dedicated to validating our analytical results through numerical simulations with large ensemble sizes.Furthermore, in Sec. 5, we discuss a system of 3D unit vector oscillators for comparison with our main result of the 4D model.Finally, we summarize our results and offer an outlook for future research in Sec. 6.

Governing Equations with Random Frustrated Interactions
For the exploration of a volcano transition, one can consider the system of heterogeneous Kuramoto phase oscillators introduced in Ref. [11].Each oscillator where N is the number of oscillators, J ∈ R + is the (global) coupling strength, and P i (t) denotes complex local forcing fields (see below).The natural frequency (ω i ∈ R) of oscillator i is selected from any unimodal distribution, which we take to be centered at zero, exploiting the rotational symmetry of Eq. ( 1).More precisely, we employ a Gaussian distribution where ∆ is the standard deviation, which is fixed at ∆ = 1 throughout this paper.Pazó et al. observed a volcano transition with an interaction matrix (M i j ) = M ∈ R N×N which characterizes quenched random disorder in terms of interactions between oscillators.The interaction matrix was obtained as follows [11]: Each oscillator is assigned two timeindependent random vectors where L ≥ 2 is the dimension of the interaction vectors.Their components are selected randomly and independently as ±1 with probability 0.5, and the matrix of random frustrated interactions is defined as where η ∈ [−1, 1] is a bifurcation parameter controlling the weight between S and A in the interaction matrix M. Here, the symmetric matrix S is defined by inner products S i j := ⟨u i , u j ⟩ − ⟨v i , v j ⟩ that characterizes a reciprocal interaction whereas the asymmetric matrix A given by A i j := ⟨u i , v j ⟩ − ⟨v i , u j ⟩.It determines an anti-reciprocal interaction between oscillators i and j.Both satisfy S ii = A ii = 0 for i = 1, ..., N (no self-interaction).It is worth noting that in Ref. [10] the authors only investigated symmetric, reciprocal interactions between oscillators, corresponding to η = 1 in Eq. (3).However, for η < 1, the interaction is rendered nonreciprocal.Additionally, the interaction matrix encompasses both attractive (M i j > 0) and repulsive (M i j < 0) couplings.All the detailed properties of the above interaction matrix are summarized in Refs.[10,11].Throughout this paper, we maintain a constant value of L = 2.This choice adheres to the condition L ≪ log 2 N.This also ensures that, for the sake of analytical tractability, the interaction matrix (M i j ) is intentionally configured as a low-ranked quenched disorder, rather than being fully ranked, i.e., rank(M) = 2L except for η = 0 [11].
In the above setup, the detection of a volcano transition can be accomplished by assessing radial distributions of complex local fields, which are defined by Note that the global Kuramoto order parameter is defined by which cannot be used to observe a volcano transition since it remains close to zero for all coupling strength J.This explains a prominent difference between well-known partial synchronization in a system of globally coupled oscillators and a volcano phase of oscillators with random frustrated interactions.The latter can be detected by measuring a local degree of coherence on the connectivity topology whereas the former is described by the global order parameter.At the transition, the radial distribution h(r i ) of the local fields becomes concave up at the origin (the most probable value of r i is greater than zero) from being concave down (the most probable value of r i is zero) [6,10,11].Furthermore, employing the OA ansatz [16][17][18][19][20], an analytical prediction of the critical coupling strength J c was derived in Refs.[10,11], using the linear stability analysis of incoherent states in the thermodynamic limit: which is a function of the parameter η.Above J c , the incoherent state is destabilized.Hence, it was shown that a volcano transition can occur in a system of heterogeneous Kuramoto oscillators with nonreciprocal as well as reciprocal interactions [11].Moreover, Eq. ( 6) demonstrates that the likelihood of observing the volcano transition decreases as the symmetry parameter η becomes smaller, and eventually goes to zero as η → 0 + .As outlined above, the previous explorations of the volcano transition primarily focused on systems of heterogenous Kuramoto phase oscillators, employing Eq. ( 1) or similar systems [6,10,11].In this paper, we consider a system of so-called generalized Kuramoto oscillators proposed by Lohe (for details of the below description, see Ref. [31]) that consists of non-Abelian Kuramoto oscillators where each oscillator Here, † indicates the Hermitian conjugate.Note that in a broader context, one has the flexibility to consider various classical compact Lie groups for this generalization, e.g., U i ∈ U(n) or SO(n) [31].However, our focus in this paper is exclusively on the subgroup SU(2).In Eq. ( 7), one can choose the local dynamics to be dictated by a time-independent Hermitian matrix with zero trace, denoted as H i ∈ su 2 (Lie algebra).Then, the right-hand side of Eq. ( 7) retains Hermitian properties.This implies that the dynamics of the oscillators are restricted to evolve on the compact manifold of SU(2) as long as the initial condition U i (0) ∈ SU(2) holds.As a consequence, U i U i † is a constant of motion and U i (t) ∈ SU(2) for all t > 0. In essence, the local dynamics is nothing but a finite-dimensional Schrödinger's equation, as it can be expressed by i Ui = H i U i .It is noteworthy that the Abelian Kuramoto model in Eq. ( 1) can be revisited using U i = e −iφ i ∈ U(1) and setting . Furthermore, as pointed out in Ref. [31], the equation exhibits a chiral SU(2) × SU(2) covariance.This can be seen when considering the transformations U i → SU i T and H i → SH i S † , with S, T denoting constant unitary matrices.In such cases, SU i T also satisifies Eq. ( 7) if U i is a solution.
Given our exclusive focus on the SU(2) model, we can parameterize U i ∈ SU(2) using a real 4-vector with unit magnitude on the 3-sphere since the group SU( 2) is diffeomorphic to the 3-sphere [38].This parameterization reads for i ∈ [N] where σ k for k = 1, 2, 3 are Pauli matrices and I n ∈ R n×n is the identity matrix [31].Substituting Eq. ( 8) into Eq.( 7), one can obtain the governing equations of the unit 4-vectors on S 3 that read ẋi Also, Eq. ( 9) guarantees that if ∥x i (0)∥ = 1 holds, then ∥x i (t)∥ = 1 for t > 0, namely, ∥x∥ := ⟨x, x⟩ is a constant of motion.Moreover, if we restrict the system to 2-vector oscillators (x i = (cos φ i , sin φ i ) ⊺ ∈ S 1 ), then we regain the usual (2D) Kuramoto model as in Eq. ( 1).In Eq. ( 9), the local dynamics is governed by a real, anti-symmetric natural frequency matrix for each oscillator, i.e., where ⊺ denotes the transpose of a matrix.These natural frequency matrices, as described in Eq. ( 19) of Ref. [31], do not represent the most general form of elements of so 4 .Instead, they span su 2 , corresponding to an SU(2) subgroup of the SO(4) group.This correspondence arises from the local isomorphism SO(4) ∼ = SU(2) × SU(2) and the presence of the chiral covariance [31].Nevertheless, for the sake of generality, we will utilize the most general form of elements from so 4 to assign each oscillator an anti-symmetric natural frequency matrix.Thus, a natural frequency matrix of each oscillator for our system is given by for i ∈ [N] and satisfies Ω ⊺ i = −Ω i .We also assume that the natural frequency matrices are distributed according to a distribution G(Ω).As in Ref. [39], we here choose the natural frequency matrix distribution G(Ω) as follows: each component ω a i for each a = 1, ..., 6 is selected randomly and independently from a Gaussian distribution in Eq. ( 2) with ∆ = 1.This method allows us to exploit rotational invariance and so the anti-symmetric matrices are set to have zero mean.This fact will be useful later in Sec.3.2.Lastly, the random frustrated interactions M i j in Eq. ( 9) are defined in Eq. (3).Also in this system, the global Kuramoto order parameter, which is here defined as a center of mass of oscillators on S 3 , namely, remains close to a zero vector for J ≥ 0 and so it is not indicative of a volcano transition.Therefore, we define a local field vector as for i ∈ [N], similarly to Eq. ( 4).Given its four-dimensional distribution, we will measure its partial information by examining distributions on individual planes in R 4 such as h(x This approach provides an effective means to observe a volcano transition in our system.See numerical results in Sec. 4 for details.

Theoretical Result: Critical Coupling Strength
In this section, we derive an analytical expression for the critical coupling strength J c in the non-Abelian Kuramoto model ( 9) at which the volcano transition occurs, i.e., at which the incoherent state loses its stability.

Higher-dimensional Ott-Antonsen Ansatz
We start by applying the so-called higher-dimensional OA ansatz, recently introduced in Refs.[33,34], to our system (9) with random frustrated interactions.Although we are interested in 4D oscillators, we here consider any even-dimensional oscillators defined on the (D − 1)-sphere for D ≥ 2, i.e., x i ∈ S D−1 , in order to provide a general description.The odd-dimensional cases are considered in Sec. 5.
The governing equations ( 9) of the non-Abelian Kuramoto oscillators describe a finitesized ensemble.However, to apply the OA ansatz, we need to consider the thermodynamic limit N → ∞.The state function of this system is the density of oscillators f (r, Ω, w;t), specified by the continuous natural frequencies Ω and the 2 2L interaction vectors w := (u ⊺ , v ⊺ ) ⊺ ∈ R 2L which are discrete parameters for a continuous unit vector oscillator r = r(θ 1 , • • • , θ D ) ∈ S D−1 .Since the oscillator density function is defined on the sphere S D−1 , it can be expanded by higher-dimensional spherical harmonics [34] such that similar as the oscillator density function for the Abelian Kuramoto oscillators can be expanded in Fourier series [16][17][18][19][20]. Furthermore, in the thermodynamic limit, the continuous local field vector is defined by where M(w, w ′ ) is the (w, w ′ )-th matrix element of M ∈ R 2 2L ×2 2L [11].Due to the conservation of the number of oscillators, the oscillator density function is governed by the continuity equation that reads [33,34] where ∇ S f denotes the component of the gradient of f along the surface S D−1 of the unit ball [39].
The OA ansatz for the usual 2D Kuramoto model assumes that the n-th Fourier coefficient is given by the n-th power of the first coefficient.This first coefficient constitutes the OA variable and describes the system's macroscopic dynamics [16,18,19].Likewise, Barioni et al. suggested the higher-dimensional Ott-Antonsen ansatz constructed similarly to the 2D case.Following their ansatz, we assume a manifold where the coefficients in Eq. ( 13) are given by where the vector and ρ ∈ S D−1 [34].The symbol with an overbar indicates the complex conjugate.Plugging Eq. ( 16) into Eq.( 13), the oscillator density is written as a higher-dimensional normalized Poisson kernel such that where S D := 2π D/2 Γ(D/2) as pointed out in Refs.[30,34].Hence, in the higher-dimensional OA manifold, the local field vector defined in Eq. ( 14) reads for each interaction vector w since the normalized Poisson kernel leads to where ρ = ρ ∥ρ∥ is a normalized OA vector and r = η ρ + n with ⟨ρ, n⟩ = 0 [21,30].Finally, substituting Eq. ( 18) and Eq. ( 17) into the continuity equation ( 15), one can derive the higherdimensional OA equation [34]: Note that Lohe obtained a similar dimension-reduced governing equation for the identical higher-dimensional oscillators using the higher-dimensional Watanabe-Strogatz transformation [32].

The Critical Coupling Strength
In this subsection, we obtain the critical coupling strength exploiting the higher-dimensional OA ansatz in Sec.3.1.The volcano transition occurs at the parameter value at which the incoherent state loses its stability.Therefore, we determine the linear stability of the incoherent state, as also done in Refs.[10,11].Let us consider an infinitesimal perturbation vector β with ∥β∥ ≪ 1 off the incoherent state ρ 0 = 0. Substituting ρ(Ω, w,t) = ρ 0 + β(Ω, w)e λt into Eq.( 20), we obtain the following linearized equation: for each interaction vector w.As described in Refs.[34,39], there exists a real orthogonal matrix O that transforms the anti-symmetric natural frequency matrix Ω into a blockdiagonalized form [40].Under the transformation Ω → O ⊺ ΩO, the natural frequency matrix reads where all the off-diagonal blocks are zero matrices and each diagonal block is given by Applying the transformation β → Oβ in Eq. ( 21), we obtain In Eq. ( 23), one can choose ) is a joint distribution of the frequencies {ω i } and F(O) denotes a uniform distribution of orthogonal matrices (see Appendix C of Ref. [39]).Considering each subspace with W i and using a complex representation, we rewrite the perturbation vector as where b i ∈ C 1 for i = 1, ..., D/2.Then, the linearized equation ( 23) reads for each w We define now bn This allows us to algebraically rearrange Eq. ( 25) and integrate it over the joint distribution of the frequencies on both sides.We then arrive at Note that at the critical point, the real part of the eigenvalue λ goes to zero [11].To exploit this, let us consider the first component b1 (w).We obtain b1 (w) = J where the value of g(0) is given by [39] g(0 for D ≥ 2 if we consider the zero mean distribution of natural frequency matrices.Since all the other components obey the same equation ( 26), we obtain a matrix form of an eigenvalue equation that determines the stability of the incoherent state: where I is the identity matrix.The nontrivial solutions of Eq. ( 29) ( b ̸ = 0) are given by eigenvectors of M corresponding to nonzero eigenvalues [11].The nonzero eigenvalues of M are given by √ η2 2L and − √ η2 2L with multiplicity L [10,11].Thus, the volcano transition can only be detected provided that the system satisfies η > 0. Finally, we provide the analytical form of the critical coupling strength: which is a function of the symmetry parameter η.Comparing Eq. ( 30) to Eq. ( 6), we observe that both have the same form except for g(0) ̸ = g(0) due to the higher-dimensional nature of our vector oscillators.Moreover, for D = 2, we regain the critical coupling strength for the Abelian Kuramoto model in Eq. ( 6).

Numerical Results: Volcano Transition for 4D Kuramoto Oscillators
In Sec. 3, we derived the critical coupling strength (30) as a function of the symmetry parameter η in the thermodynamic limit using the higher-dimensional OA ansatz.In this section, we validate our analytical findings through numerical simulations in a system of unit 4-vector oscillators as given in Eq. ( 9), considering a large value of N that satisfies L ≪ log 2 N to have a low-rank interaction matrix.In Fig. 1 (a), a phase diagram in the η −J plane is illustrated together with the analytically determined critical coupling strength J c (η). Akin to the Abelian Kuramoto model [11], we observe the divergence of the volcano phase boundary as η → 0 + , a region where no volcano transition is expected to occur.Additionally, the phase diagram shows that a volcano transition can be observed not only with symmetric and reciprocal interactions (η = 1) but also with asymmetric and nonreciprocal interactions (0 < η < 1).In Fig. 1 (b-c), we present normalized distributions for the components of local field vectors, focusing on fully symmetric interactions (η = 1).We numerically assess 2D distributions for the third (P 3 i ) and fourth (P 4 i ) components of the local field vectors (12), both below and above the critical coupling point.Notice that the distributions for the first (P 1 i ) and second (P 2 i ) components display the same characteristics (not shown here).As depicted in Fig. 1 (b), the distribution of local field vectors for J < J c exhibits a bell-shaped unimodal pattern centered at the origin.This observation confirms that, for J < J c , the most probable value of the local fields is indeed a zero vector.On the contrary, as shown in Fig. 1 (c), the zero vector is not the most likely local field vector for J > J c ; instead, the distribution exhibits its peak at a nonzero vector.This distinctive characteristic, where the behavior differs below and above the critical point, is a hallmark of volcano transitions [6,10,11].Next, we shift our focus to nonreciprocal interactions characterized by 0 < η < 1.To explore this, we also conducted numerical simulations, however, employing nonreciprocal random interactions (η = 0.8).As illustrated in Fig. 1 (d-e), similarly, the distribution of local fields follows a bell-shaped unimodal feature for J < J c whereas it possesses its peaks at nonzero local fields for J > J c , which also indicates a detection of a volcano transition for nonreciprocal interactions.
In Fig. 2, cross sections of distributions of local fields are depicted with various coupling strengths for η = 1 (reciprocal), η = 0.8 and η = 0.4 (nonreciprocal).For both reciprocal and nonreciprocal interactions, the distributions are bell-shaped and centered at the origin below the critical coupling strength, which statistically demonstrates that the incoherent state is stable.As the coupling J increases up to J c , this bell-shaped distribution becomes ever broader while the peak at the origin becomes lower.Above the critical point J > J c the incoherent state is destabilized and the peak of distributions appears at a non-zero point.Finally, we examine the positions of the maximum as a function of the coupling strength in order to illustrate where a volcano transition takes place.
In Fig. 3, the position of the maximum ∥P * ∥ is depicted as a function of the coupling strength J for three different symmetry parameters, within our numerical capacity.The  12) for reciprocal interactions with η = 1.0 and (b) J = 1.8, (c) J = 2.8 and for nonreciprocal interactions with η = 0.8 and (d) J = 2.0, (e) J = 2.6.Note that J c = 2.12769 for η = 1.0 and J c = 2.37883 for η = 0.8.We perform numerical integration of Eq. ( 9) using a fourth-order Runge-Kutta method with a time step of dt = 0.05 and a system size of N = 500.After discarding transient steps, we collect the data between t = 1500 and t = 2000 for each ensemble, and then compute the average over 100 samples, each having distinct sets of natural frequencies and random interaction matrices.maximum ∥P * ∥ stays near zero for J < J c where the incoherent state is stable.However, it bifurcates off the zero value close to the value of J c analytically determined in Eq. (30).Note that beyond the critical point, we observe certain fluctuations as the coupling strength increases.This observation might arise from either the limited number of ensembles or the limited number of oscillators.
As discussed in Sec. 3, the volcano transition is characterized by the loss of stability of the incoherent and the emergence of partial synchrony.To illustrate this phenomenon, we examine the time-averaged velocity of each component of the oscillators, defined as A blue histogram: J = 2.8 > J c ; an orange histogram: J = 1.8 < J c .All figures are obtained from reciprocal interactions (η = 1).However, similar results for nonreciprocal interactions with η < 1 (results not shown here) can be found.We perform numerical integration of Eq. ( 9) using a fourth-order Runge-Kutta method with a time step of dt = 0.05 and a system size of N = 500.We average the instantaneous velocities of each component for ∆t = 1000 up to T max = 3000.
In Fig. 4, we present density distributions of time-averaged velocities for each component, both below and above the transition point J c .Below the transition (orange), incoherent states are evident, showing no substantial locking behavior and exhibiting a relatively broad distribution around zero.Conversely, above the volcano transition (blue), a significant portion of the oscillators is locked, resulting in a density distribution of time-averaged velocities that is highly concentrated around zero.This observation corroborates the fact that the incoherent state indeed becomes destabilized, giving rise to a partially locked state for J > J c .It is noteworthy that a similar outcome was observed in a system of Abelian Kuramoto oscillators in Ref. [10].
For a more detailed analysis of the properties of observable dynamical states, we explore the degree of correlation between each oscillator x i and the local field vector P j .As previously noted, the quenched random interaction matrix M encompasses both attractive (M i j > 0) and repulsive (M i j < 0) couplings between oscillators.Consequently, in a partially locked state for J > J c , oscillator x i tends to align with the local field vector P j when it is coupled attractively with oscillator x j , i.e., M i j > 0. Conversely, in the case of repulsive coupling, the oscillator x i tends to be in anti-alignment with the local field vector P j [10].In this analysis, we explore the cosine similarity between an oscillator x i and the local field vector P j for all i, j ∈ [N] defined by Note that when oscillator i is aligned with the local field j, then S c (x i , P j ) = 1, indicating that the angle between them is 0. In contrast, S c (x i , P j ) = −1 implies that the angle between  31)) are plotted against the interaction matrix element M i j for (a) J = 1.8 < J c and (b) J = 2.8 > J c .All figures are obtained for reciprocal interactions (η = 1).However, similar results for nonreciprocal interactions with η < 1 (results not shown here) can be found.We perform numerical integration of Eq. ( 9) using a fourth-order Runge-Kutta method with a time step of dt = 0.05 and a system size of N = 4000.them is ±π, i.e., they are in anti-alignment.In Fig. 5, we display the densities of S c (x i , P j ) as a function of the interaction matrix element M i j for both (a) J < J c and (b) J > J c .Below the volcano transition, we observe that the local field vectors are not correlated to the oscillators that characterize the incoherent state.This can be noticed from the distribution of alignments, which shows vertical stripes [10], as depicted in Fig. 5 (a).However, as shown in Fig. 5 (b), above the critical point, the density of S c (x i , P j ) versus M i j exhibits peaks at approximately S c (x i , P j ) = 1 for the attractive coupling M i j > 0 and at around S c (x i , P j ) = −1 for the repulsive coupling M i j < 0. This observation demonstrates how the oscillators are aligned with respect to the local fields in the emergence of a partially locked state (a volcano phase) [10].
In previous studies [10,[12][13][14][15], it was conjectured that an oscillator glass might be characterized by nonexponential or algebraic relaxation of the global Kuramoto order parameter Γ(t) defined in Eq. (5).To investigate the decay of the global Kuramoto order parameter in our system, we tracked ∥Γ(t)∥ as given in Eq. ( 11) over time, commencing from an initially synchronized state.In Fig. 5 (c), the decay of the norm of the global Kuramoto order parameter is shown as a function of time with different coupling strengths.For both cases of J < J c and J > J c , we observed only an exponential relaxation of the norm of the center of mass vector as time passes.This outcome is in line with the findings in the Abelian Kuramoto model.In Ref. [10], the authors showed that a volcano phase exhibits an exponential relaxation of the global order parameter for a low-rank interaction matrix, i.e., L ≪ log 2 N. Here, we arrive at the same conclusion: The presence of volcano phases in a system of non-Abelian Kuramoto oscillators (9) does not necessarily indicate the existence of a glassy state for L ≪ log 2 N.

Absence of Volcano Transition for 3D Kuramoto Oscillators
So far, we have explored a volcano transition that occurs in a system of heterogeneous unit 4-vector oscillators with reciprocal and nonreciprocal random interactions.Various studies have highlighted distinctions between even and odd dimensions in the collective dynamics of systems comprising higher-dimensional oscillators.For example, globally coupled heterogeneous oscillators show different characteristics in terms of the onset of partial synchronization.For odd dimensions, a partially locked state discontinuously occurs, no matter how small the global coupling strength is [28,34,39].Furthermore, when we consider odd dimensions, higher-dimensional, identical oscillators with a phase-lag rotation matrix are known to possess additional complete synchronization as fixed points at the north and south poles, regardless of the phase-lag parameter [30].This particular distinction arises due to a specific direction of a rotation in an odd dimension.An even-dimensional rotation is an isoclinic rotation consisting of planes of rotations.On the contrary, for odd dimensions, there is a so-called rotational axis perpendicular to the planes of rotation.In many cases, oscillators are inclined to this axis such that one can have complete synchronization at the poles [30], or a partially locked state along the natural frequency vector [28,31,34,39], regardless of a given parameter.
In previous studies on volcano transitions [6,7,10,11], a system of Abelian Kuramoto oscillators was investigated with random frustrated interactions.In the context of the current paper, such Abelian Kuramoto oscillators are 2-dimensional vector oscillators on the unit circle.Thus, one can expect a volcano transition from an incoherent state to a volcano phase, as for 4D oscillators in the previous sections.In this section, we demonstrate that no volcano transition occurs for odd dimensions since oscillators are partially locked along the average natural frequency vector, irrespective of the global parameter J > 0. We study a system of unit 3-vectors with the same random interactions.Note that equation (7) for SU(2) corresponds to the unit 4-vector model on S 3 .To consider a system of unit 3-vector oscillators on S 2 , Lohe restricted the system to follow where U i ∈ SU(2) for i ∈ [N].In this case, the chiral covariance is given by S = T [31].Then, oscillators in Eq. ( 32) can be parameterized by the same method as given in Eq. ( 8), however, setting x 4 i ≡ 0 identically.Thus, one obtains the following system of unit 3-vectors [31,33]  for all i ∈ [N]: where the natural frequency vector reads ω i = (ω 1 i , ω 2 i , ω 3 i ), each component of which is selected randomly and independently from an unimodal distribution [31].The OA ansatz in Sec.3.1 can be applied to this system with D = 3, which leads to the following evolution equations of the OA vector and its norm [33]: We now look for an equilibrium solution of Eq. (34).The second equation indicates that an equilibrium solution should be either ∥ρ∥ = 1 or ρ ⊥ P. When we consider the coupling strength J as a very small parameter, then d dt ρ ∼ ω × ρ as J → 0 + , which thus implies that an equilibrium should satisfy ρ ∼ ±ω.Employing small perturbations ρ = ±ω + δ ρ and ∥ρ∥ = 1 + δ ρ, one can show that the solution ρ ∼ ±ω is stable for ⟨ ρ, P⟩ > 0 since the perturbation is governed by d dt δ ρ = − J 2 ⟨ ρ, P⟩δ ρ [33].Thus, we draw a conclusion that there is no volcano transition for D = 3: The locked state, i.e., either ρ ∼ ω or ρ ∼ −ω, becomes already stable for J → 0 + , regardless of the value of η.
This conjecture is verified by numerically solving Eq. ( 33) for large N.In Fig. 6, distributions of components of local fields are depicted for very small coupling strengths: (a) J = 0.01 and (b) J = 0.05.For both cases, we observe a volcano-shaped distribution of local fields.The most probable value of the local fields does not occur at the origin but rather the peak of the distribution is detected at a non-zero vector.This numerical observation is consistent with the above analytical conjecture that no volcano transition can be observed in a system of unit 3-vector oscillators on the 2-sphere.

Summary and Outlook
In this paper, we generalized a system of Kuramoto phase oscillators to higher-dimensional Kuramoto vector oscillators with both reciprocal and nonreciprocal frustrated interactions proposed in Ref. [11].This generalization also predicted a volcano transition in which the incoherent state becomes unstable and a partially locked state emerges.We derived an analytical expression of the critical coupling strength at which the volcano transition occurs using the so-called higher-dimensional OA ansatz in the thermodynamic limit [34].With numerical simulations of large-sized ensembles of oscillators, we validated our analytical results for 4D oscillators.The obtained collective behaviors exhibited distinct dynamical characteristics below and above the volcano transition, respectively.A significant portion of oscillators above the volcano transition becomes entrained, and their vectors align or antialign with local field vectors, depending on whether the interactions between oscillators are attractive or repulsive.In contrast, the incoherent states below the threshold point do not exhibit these characteristics.Furthermore, we proposed the absence of a volcano transition in a system of 3D Kuramoto oscillators.This conjecture was supported by our observation of a stable partially locked state as the coupling strength approached zero from the positive side.We also confirmed this conjecture through numerical simulations.Note that another study on higher-dimensional phase oscillators with random frustrated interactions was recently reported, where the authors considered an equilibrium in the mean-field equation for vector spin models on random networks with high connectivity, featuring an arbitrary degree distribution and links with random weights [41].
In fact, our analytical outcome (30) is not confined to 4D Kuramoto oscillators alone.Therefore, one has the flexibility to numerically substantiate this conclusion using an ensemble of oscillators with any even dimensionality, such as D = 6, 8, • • • .This could provide further insights into how the dimensionality of the oscillators influences observable collective dynamics.Similarly, the conjecture of the absence of a volcano transition can be extended to encompass investigations for any odd-dimensional oscillators in future research.Moreover, in this paper, we focused solely on low-rank interaction matrices for the sake of analytical convenience.Nevertheless, just as in Refs.[6,10,11], one can employ full-ranked interaction matrices to investigate a volcano transition and its relation with glassy oscillator states, which still remains controversial.
Recently, generalized Kuramoto oscillators have received considerable attention, with a specific emphasis on the dynamical and spectral properties of their collective behaviors.On the one hand, there has been a trend toward exploring higher-dimensional generalizations to investigate phenomena like synchronization or chimera states [21][22][23][24][25][26][27][28][29][30].On the other hand, there has been a growing interest in the complexification of oscillator dynamics [42][43][44][45].While these two generalizations possess distinct characteristics, it is worth considering future investigations into any generalization of Kuramoto oscillators to explore collective dynamics within a variety of interaction setups.This could encompass scenarios such as systems with other random frustrated interactions than those discussed here, oscillators arranged on a ring geometry, or within graph topologies, as these could offer avenues for a closer approximation to real-world systems.

Figure 1 .
Figure 1.(a) Phase diagram obtained from the analytic Eq. (30).The dashed line indicates the critical coupling strength as a function of the symmetry parameter η.For comparison, we display critical coupling strengths that are numerically obtained (black solid dots).(b-e) Normalized distributions of components of local fields (12) for reciprocal interactions with η = 1.0 and (b) J = 1.8, (c) J = 2.8 and for nonreciprocal interactions with η = 0.8 and (d)J = 2.0, (e) J = 2.6.Note that J c = 2.12769 for η = 1.0 and J c = 2.37883 for η = 0.8.We perform numerical integration of Eq. (9) using a fourth-order Runge-Kutta method with a time step of dt = 0.05 and a system size of N = 500.After discarding transient steps, we collect the data between t = 1500 and t = 2000 for each ensemble, and then compute the average over 100 samples, each having distinct sets of natural frequencies and random interaction matrices.

Figure 2 .Figure 3 .Figure 4 .
Figure 2. Cross-sections of normalized distributions of local fields (12) for both reciprocal (the first row) and nonreciprocal (the second and third row) random interactions.Note that J c = 2.12769 for η = 1.0,J c = 2.37883 for η = 0.8, and J c = 3.36418 for η = 0.4.Orangecolored distributions are obtained for J < J c while blue-colored distributions are for J > J c .The numerical methods used are the same as those given in the caption of Fig. 1.

Figure 5 .
Figure 5. Densities of S c (x i , P j ) (see Eq. (31)) are plotted against the interaction matrix element M i j for (a) J = 1.8 < J c and (b) J = 2.8 > J c .All figures are obtained for reciprocal interactions (η = 1).However, similar results for nonreciprocal interactions with η < 1 (results not shown here) can be found.We perform numerical integration of Eq. (9) using a fourth-order Runge-Kutta method with a time step of dt = 0.05 and a system size of N = 4000.After discarding transient steps, we collect the data between t = 1000 and t = 2000.(c) Time evolution of the norm of the global Kuramoto order parameter(11) for η = 1.0 in a log-log scale.Gray dashed guideline indicates an algebraic relaxation.All the considered cases exhibit an exponential relaxation: J = 1.0 < J c (blue), J = 6.0 > J c (green), and J = 15.0 > J c (red).
Figure 5. Densities of S c (x i , P j ) (see Eq. (31)) are plotted against the interaction matrix element M i j for (a) J = 1.8 < J c and (b) J = 2.8 > J c .All figures are obtained for reciprocal interactions (η = 1).However, similar results for nonreciprocal interactions with η < 1 (results not shown here) can be found.We perform numerical integration of Eq. (9) using a fourth-order Runge-Kutta method with a time step of dt = 0.05 and a system size of N = 4000.After discarding transient steps, we collect the data between t = 1000 and t = 2000.(c) Time evolution of the norm of the global Kuramoto order parameter(11) for η = 1.0 in a log-log scale.Gray dashed guideline indicates an algebraic relaxation.All the considered cases exhibit an exponential relaxation: J = 1.0 < J c (blue), J = 6.0 > J c (green), and J = 15.0 > J c (red).

Figure 6 .
Figure 6.Distributions of components of the local fields (12) for (a) J = 0.01 and (b) J = 0.05.The numerical methods used are the same as those given in the caption of Fig. 1.