Energy landscapes and synergetic state transitions in frustrated Stuart–Landau oscillator networks: a homotopy continuation study

1 Introduction

Many natural and artificial systems spontaneously develop ordered and coherent patterns without centralized control. Synergetics, pioneered by Hermann Haken (Haken, 1977, Haken, 1983b, Haken, 1989), provides a conceptual framework for such self-organization by emphasizing how macroscopic structures emerge from the collective dynamics of many interacting components. A central message is that, near a qualitative reorganization of behavior, the high-dimensional dynamics of a complex system can often be captured by a small set of collective variables (i.e., order parameters), while the remaining degrees of freedom are effectively slaved to them. In this view, abrupt changes in global behavior are driven by shifts in control parameters that reshape the accessible macroscopic states and their stability, and fluctuations such as noise can trigger transitions among coexisting attractors.

Coupled oscillator networks constitute a prototypical arena for synergetics (Haken, 1983a): simple rhythmic units interact locally yet generate global synchronization, patterned phase-locking, and complex spatiotemporal oscillations (Strogatz, 2000; Strogatz, 2003). In some oscillator networks, collective dynamics can be organized by an energy or potential-like function, yielding a landscape-level representation of multistability and state transitions. From a synergetics perspective, the energy landscape provides a natural interpretive lens: stable macroscopic patterns correspond to minima, while unstable critical points and barriers shape the transition structure; noise can induce switching by enabling escapes across basins, and control parameters can deform the landscape, creating, destroying, or reshaping its structure. This viewpoint is especially valuable when one aims to identify and compare multiple coexisting collective states rather than focusing on a single attractor.

Recently, energy-based dynamical models and coupled oscillator networks have gained renewed attention as mechanistic models of interacting physiological systems, particularly for large-scale brain activity (Schneidman et al., 2006; Fraiman et al., 2009; Watanabe et al., 2013; Deco et al., 2017; Deco and Kringelbach, 2020; Chinichian et al., 2024; Hannay et al., 2018; Smart et al., 2023). The brain exhibits rich rhythms that may support information processing across perception, cognition, and memory (Buzsáki and Draguhn, 2004; Fries, 2005; Fries, 2015; Buzsáki and Mihály, 2023), and the coexistence of multiple metastable states provides a principled framework for discussing flexible switching among brain states (Tognoli and Kelso, 2014; Adachi and Aihara, 1997; Tsuda, 2015; Deco et al., 2019). Spin-based energy models have been influential in formalizing such metastability in terms of energy landscapes (Watanabe et al., 2013; Watanabe et al., 2014; Ezaki et al., 2017; Masuda et al., 2025). More recently, their application to health checkup data has identified multiple states and different state transitions in diabetes progression, highlighting their potential across a broader range of physiological systems (Ito et al., 2025). Meanwhile, coupled oscillator models offer an intrinsically continuous and rhythmic alternative that can express a broader repertoire of collective dynamics. These perspectives motivate studying energy landscapes in coupled oscillator networks to clarify how stable rhythmic states are organized and how noise or other possible mechanisms modulate transitions between them, as well as their potential application in networked physiology.

In this article, we investigate the energy landscape of a frustrated Stuart–Landau oscillator network and its noise-driven transitions among coexisting macroscopic states, as the frustration strength is varied across different regimes. Methodologically, we begin with an analytically tractable mode-based, phase-decoupled approximation, which suggests potentially dominant coordinates for macroscopic states; we then use homotopy continuation (Chen and Li, 2015) to track representative families of equilibria from the approximation to the full system while monitoring their energies and stability indices. In the singular cases in which the phase critical points of the initial XY-Hamiltonian-like term are non-isolated, we impose a second-order phase equilibrium condition and solve it in a constrained-solvability sense to make the continuation well posed. Extensive continuation runs reveal a partial selectivity property: an initial equilibrium with relatively high (or low) energy and stability index under the approximation tends to continue to equilibria with similarly high (or low) energy and indices in the full system. This continuation-based exploration yields a quantitative, landscape-level characterization of the typical macrostates that appear and of how their stability and energetic ordering change across frustration regimes. At an intermediate level of frustration, the energy landscape of the network exhibits an isolated global minimum and multiple 1D equi-energy troughs supported on the edges of a multigraph; for sufficiently large frustration, these troughs extend to the entire multigraph, which then constitutes the set of global minima. Based on this characterization, we demonstrate distinct noise-dependent phenomena, including noise-induced synchronization, desynchronization, and switching among macrostates; moreover, as functions of the noise level, we observe Arrhenius–Kramers-like scaling for uphill activation and diffusion-limited scaling for motion along the graph. Our results serve as a step toward an oscillator-based energy landscape analysis framework for real multivariate time series data, with potential applications in neuroscience, networked physiology, and beyond.

2 Model and methods

2.1 Stuart–Landau oscillator network: an energy-based model

We consider a network of N Stuart–Landau oscillators, with a complex-valued state s i of the ith oscillator, whose dynamics is given byfor i = 1,2 , … , N. Here, a i and b are real constants that modulate the local dynamics of a single oscillator, which undergoes a supercritical Hopf bifurcation at a i = 0 and creates a stable limit cycle s i ( t ) = a i exp ( i b t ) for a i > 0; B i represents complex-valued standard Brownian motion, and for simplicity, we assume independent, identically distributed (i.i.d.) noise across oscillators with intensity σ. The oscillators are linearly coupled via a real symmetric matrix W = [ w i j ] with zero diagonal entries because any w i i ≠ 0 should be merged with a i. The deterministic part of Equation 1 can also be expressed in polar coordinates, with s i = ρ i ⁡ exp ( i ϕ i ), asNote that we include ρ i explicitly on the left-hand side to keep the equation well defined at the phase singularity ρ i = 0, where ρ i d ϕ i / d t can be interpreted as the velocity component perpendicular to the radial axis. Moreover, Equation 2 remains invariant under the transformation ( ρ i , ϕ i ) ↦ ( − ρ i , ϕ i ± π ). In practice, we therefore allow ρ i to take negative values so that ( ρ i , ϕ i ) varies continuously across the origin. We impose the conventional restriction ρ i ∈ 0 , + ∞ and ϕ i ∈ − π , π only when a unique representation is needed.

d s i d t = s i a i + i b − | s i | 2 + ∑ j = 1 N w i j s j + σ d B i d t , (1)

d ρ i d t = ρ i ( a i − ρ i 2 ) + ∑ j = 1 N w i j ρ j ⁡ cos ( ϕ j − ϕ i ) , ρ i d ϕ i d t = b ρ i + ∑ j = 1 N w i j ρ j ⁡ sin ( ϕ j − ϕ i ) . (2)

The local dynamics in Equation 1 is a special case of the general Stuart–Landau form, in which the cubic term s i | s i | 2 may have a complex coefficient (Deco and Kringelbach, 2020; Uchiyama and Fujisaka, 2004). It is also a reduced version of the computational brain models used by Deco et al. (2017, 2019) and Deco and Kringelbach (2020), in which we use a common intrinsic angular frequency b for all oscillators instead of a distribution. A third simplification is the symmetric coupling matrix W. Under these simplifications, the deterministic dynamics in Equation 1 admits an energy-descent description in a co-rotating frame: letting ϕ i ↦ ϕ i + b t ( i = 1,2 , … , N ), the dynamics reduces to the case b = 0, where it becomes a gradient flow as follows:for i = 1,2 , … , N, with the energy function given by Equation 4:For notation in this article, we collect state variables in vector forms such as ρ = [ ρ 1 ρ 2 ⋯ ρ N ] ⊤ and ϕ = [ ϕ 1 ϕ 2 ⋯ ϕ N ] ⊤; for a scalar like H, gradients such as ∂ H / ∂ ρ and ∂ H / ∂ ϕ lie in R N, and derivatives of a vector with respect to another such as ∂ 2 H / ∂ ρ 2 and ∂ 2 H / ( ∂ ϕ ∂ ρ ) denote the corresponding Jacobian matrices (Hessian blocks for H). In addition, we write ◦ for elementwise (Hadamard) products and use the shorthand ∂ H / ( ρ ∂ ϕ ) ≔ ( 1 / ρ ) ◦ ( ∂ H / ∂ ϕ ). Then, along trajectories of Equation 3, we have:where the last equality holds if and only if ∂ H / ∂ ρ = ∂ H / ( ρ ∂ ϕ ) = 0. Thus, as indicated by Equation 5, the deterministic system in the co-rotating frame relaxes, with monotonically decreasing H, toward an equilibrium of Equation 3, or equivalently, a critical point of H ( ρ , ϕ ). Returning to the original frame, the network described by Equation 2 always converges to a frequency-synchronized solution with common angular frequency b. Accordingly, throughout the article, we focus on the gradient case b = 0 and refer to H ( ϕ , ρ ) as the energy function for the coupled oscillator network described by Equation 1, regardless of the value of b.

d ρ i d t = − ∂ H ∂ ρ i , ρ i d ϕ i d t = − ∂ H ρ i ∂ ϕ i , (3)

H ρ , ϕ = 1 4 ∑ i = 1 N ρ i 4 − 1 2 ∑ i = 1 N a i ρ i 2 − ∑ i = 1 N ∑ j < i ρ i ρ j w i j ⁡ cos ϕ i − ϕ j . (4)

d H d t = ∂ H ∂ ρ ⊤ d ρ d t + ∂ H ρ ∂ ϕ ⊤ ρ ◦ d ϕ d t = − ∂ H ∂ ρ 2 − ∂ H ρ ∂ ϕ 2 ≤ 0 , (5)

The energy function H can have a complicated critical-point structure, which depends on system parameters, most notably the coupling matrix W. Generally speaking, the equilibria of a gradient system may be isolated or non-isolated; in the latter case, they form a connected critical manifold along which the energy landscape is locally flat. Correspondingly, isolated equilibria can be asymptotically stable (local minima of H) or unstable (saddles or maxima). In contrast, non-isolated equilibria are typically semi-stable due to neutral directions tangent to the critical set. These properties are largely captured by the eigenvalue spectrum of an appropriate Hessian of H evaluated at the equilibrium. In polar coordinates, we adopt the following Hessian convention:so that the second-order Taylor expansion takes the form:which remains meaningful at the phase singularity ρ i = 0. At a critical point of H ( ρ , ϕ ), a positive definite ∇ 2 H implies an isolated minimum (a stable equilibrium), while negative eigenvalues of ∇ 2 H indicate unstable directions. Zero eigenvalues signal neutral directions and may indicate either a genuine critical manifold of non-isolated equilibria or an isolated minimum, maximum, or saddle point whose type generally must be determined by higher-order analysis. For the system given by Equation 3, the energy function in Equation 4 is invariant under the global phase shift ϕ ↦ ϕ + ε 1. Consequently, any nontrivial equilibrium with ρ ≠ 0 is non-isolated in the direction δ ϕ = 1, yielding a zero eigenvalue of ∇ 2 H. We exclude this symmetry-induced degree of freedom and assess stability considering the remaining 2 N − 1 eigenvalues; if they are all positive, the associated frequency-synchronized periodic orbit of Equation 2 is deemed a local minimum. In practice, we define the index of an equilibrium as the number of negative eigenvalues of ∇ 2 H and use it as a proxy for stability. We verify the stability of index-0 equilibria with additional zero eigenvalues by simulating perturbed trajectories. The specific expression for the Hessian in Equation 6 in a transformed frame is given by Equation S81 in the Supplementary Appendix, together with those for the whole energy homotopy, which will be described in subsequent subsections.

∇ 2 H ρ , ϕ = ∂ 2 H ∂ ρ 2 ∂ 2 H ρ ∂ ϕ ∂ ρ − Diag ∂ H ρ 2 ∂ ϕ ∂ 2 H ρ ∂ ρ ∂ ϕ − Diag ∂ H ρ 2 ∂ ϕ ∂ 2 H ρ 2 ∂ ϕ 2 + Diag ∂ H ρ ∂ ρ (6)

1 2 δ ρ ρ ◦ δ ϕ ⊤ ∇ 2 H ρ , ϕ δ ρ ρ ◦ δ ϕ ,

2.2 Mode-based approximation: A coarse description

In this subsection, we introduce a mode-based, phase-decoupled approximation that provides a coarse description of the energy landscape. The goal is to capture the main features of the stationary statistics of the system as described by Equation 1 and to obtain a convenient starting point for the subsequent, more accurate landscape exploration. For this purpose, we transform the state variable into an orthogonal eigenbasis E = [ e 1 e 2 ⋯ e N ] obtained from the eigendecomposition of the linear part of Equation 1,where a = [ a 1 a 2 ⋯ a N ] ⊤ and d = [ d 1 d 2 ⋯ d N ] ⊤ collects the corresponding eigenvalues. Let z = E ⊤ s denote the projection of the original state variable s onto the eigenbasis and introduce polar coordinates ( r , θ ) again for the new state variable z, by z = r ◦ exp ( i θ ). Then, the energy function H can be rewritten in ( r , θ ) following from the coordinate transformations aswhere superscript † denotes the conjugate transpose. By further definingwhere ( ⋅ ) ◦ 2 denotes elementwise squaring and δ i j is the Kronecker delta, the energy can be written compactly asBecause an orthogonal change of coordinates preserves the gradient-flow structure, the deterministic dynamics in ( r , θ ) can be obtained analogously to Equation 3 by differentiating H ( r , θ ) with respect to r and θ. We interpret the eigenbasis projections z = [ z 1 z 2 ⋯ z N ] ⊤ as oscillation modes; accordingly, we refer to the gradient dynamics induced by Equation 9 as mode dynamics as it describes the evolution of the mode amplitudes and phases. Finally, when W + Diag ( d ) has repeated eigenvalues, different choices of eigenmodes in E yield different representations of the mode dynamics. We may have the following strategies to handle this ambiguity: (i) choose a basis that distributes oscillation power as evenly as possible across oscillators, for example, by Riemannian minimization of ∑ i , j e i j 4 over the orthogonal manifold; (ii) apply random orthogonal rotations within the degenerate eigenspaces to circumvent certain singularities (as used later); or (iii) enforce the continuity of E with respect to some control parameter, implemented numerically by a small perturbation to E to lift the degeneracy.

W + Diag a = E Diag d E ⊤ ,

H r , θ = 1 4 ∑ i = 1 N ρ i 4 − 1 2 ∑ i = 1 N a i ρ i 2 − 1 2 ∑ i = 1 N ∑ j = 1 N w i j ρ i ρ j ⁡ exp i ϕ i − ϕ j = s = ρ ◦ exp i ϕ 1 4 ∑ i = 1 N s i † s i 2 − 1 2 s † Diag a s − 1 2 s † W s = s = E z 1 4 ∑ i = 1 N ∑ j 1 = 1 N ∑ j 2 = 1 N e i j 1 e i j 2 z j 1 † z j 2 2 − 1 2 z † Diag d z = z = r ◦ exp i θ 1 4 ∑ i = 1 N ∑ j = 1 N e i j 2 r j 2 + ∑ j 1 = 1 N ∑ j 2 ≠ j 1 e i j 1 e i j 2 r j 1 r j 2 ⁡ cos θ j 1 − θ j 2 2 − 1 2 r ⊤ Diag d r , (7)

B = E ◦ E = E ◦ 2 , W k = w i j k = e k i e k j , C θ = c i j θ = cos θ i − θ j − δ i j , (8)

H r , θ = 1 4 ∑ k = 1 N B r ◦ 2 k + r ⊤ W k ◦ C θ r 2 − 1 2 d ⊤ r ◦ 2 . (9)

Now, assume that the couplings are sufficiently “random” so that, for each k, the matrix W ( k ) has balanced positive/negative entries, and the phase-dependent cross term, r ⊤ ( W ( k ) ◦ C ( θ ) ) r, is negligible compared with ( B r ◦ 2 ) k. This leads us to neglect the interaction among modes via the phase and approximate the energy bywhich further yields the approximate mode dynamics as an amplitude system,with the constraints r ≥ 0. Here, we let ≥ and ≫ denote elementwise “no less than” and “strictly greater than,” respectively, while x > y means x ≥ y with x ≠ y. If B ⊤ B is nonsingular (hence positive definite), Equation 11 in both the deterministic and stochastic settings admits the following properties:

H approx r = 1 4 B r ◦ 2 ⊤ B r ◦ 2 − 1 2 d ⊤ r ◦ 2 , (10)

d r d t = − ∂ H approx r ∂ r = − B ⊤ B r ◦ 2 − d ◦ r (11)

• Equilibria via active sets. Letting ψ : = r ◦ 2, the critical points of H approx ( r ) in Equation 10 satisfy B ⊤ B ψ − d ◦ ψ = 0; that is, each component either vanishes ( ψ i = 0 ) or satisfies the linear equation ( B ⊤ B ψ ) i = d i. The equilibria of system Equation 11 can thus be described by choosing a support, or “active set” of indices with ψ i > 0 and then solving a principal sub-system. This is expressed by where the superscripts (1) and (2) above denote components in two complementary subspaces and ( B ⊤ B ) ( 1,1 ) denotes the corresponding principal submatrix. Because every principal submatrix of a positive definite matrix is positive definite, there are at most 2 N candidate solutions (one per choice of the active set), and only those with ψ ( 1 ) ≫ 0 are feasible equilibria.
• Uniqueness of the stable equilibrium and mode sparsity. In terms of ψ, the approximate energy is the convex quadratic function given by Equation 13: which is strongly convex on the convex set 0 , ∞ N. Thus, it has a unique minimizer ψ ∗ ∈ 0 , ∞ N, corresponding to a unique stable equilibrium r ∗ = ψ ∗ ◦ 1 2 of Equation 11; all other feasible critical points are saddles. If ( B ⊤ B ) − 1 d ≫ 0, then the minimizer ψ ∗ = ( B ⊤ B ) − 1 d (the full-support solution to Equation 12); otherwise, ψ ∗ coincides with the feasible solution to Equation 12 that is nearest to ( B ⊤ B ) − 1 d in the Mahalanobis distance induced by B ⊤ B, lying on the boundary of 0 , + ∞ N. This suggests that, under Equation 11, some mode amplitudes tend to vanish to minimize the energy, so the stable state is typically supported on only a few oscillation modes. From the perspective of synergetics, these active modes act as order parameters that can organize highly coherent behavior such as synchronization.
• Steady state under noise as a truncated Gaussian in ψ . As the i.i.d. Gaussian noise in Equation 1 is invariant under orthogonal transformations, the stationary probability density in z admits a Gibbs form P approx ( z ) ( z ) ∝ exp − 2 H approx ( z ) ( z ) / σ 2. Under the phase-decoupled approximation H approx ( z ) ( z ) = H approx ( r ) = H approx ( ψ ) ( ψ ) with z ◦ z † = r ◦ 2 = ψ, integrating out θ and invoking ∏ i r i d r i = 2 − N ∏ i d ψ i easily lead to P approx ( ψ ) ( ψ ) ∝ exp ( − 2 H approx ( ψ ) ( ψ ) / σ 2 ), that is, on ψ ∈ 0 , + ∞ N, where Z ( B , d , σ ) is the normalization factor or partition function. Thus ψ follows a truncated multivariate normal distribution on the positive orthant given in Equation 15, with the previously mentioned ψ ∗ attaining the maximum density. We expect this distribution to concentrate near ψ ∗ as σ → 0, revealing clear dominant modes. For sufficiently large σ, fluctuations may overwhelm the deterministic dynamics and degrade the organized behavior.

Although the approximate mode dynamics (Equation 11) admits the characterization above, the associated partition function Z ( B , d , σ ) for the truncated normal distribution shown by Equation 14 generally has no closed-form expression. Closed-form evaluation of Z ( B , d , σ ) and marginal density functions is possible only for special choices of the quadratic form. In our previous work (Li et al., 2022), an even simpler mean-field approximation has been introduced by taking all entries of B = [ e i j 2 ] to be approximately 1 / N, which yields the following: for ψ ∈ 0 , + ∞ N, where O is the all-one matrix. The amplitude dynamics induced by Equation 16 for r (i.e., replacing ψ by r ◦ 2 and then differentiating with respect to r) reduces toWhen d 1 > d 2 > ⋯ > d N + > 0 ≥ d N + + 1 ≥ ⋯ ≥ d N, this system has exactly N + nontrivial equilibria, each supported on a single mode (all other components vanish). When d 1 > 0, all equilibria are saddles except r ∗ = [ N d 1 0 ⋯ 0 ] ⊤, which is the unique stable equilibrium of the system given by Equation 18 and corresponds to the global minimizer of H approx ( ψ ) ( ψ ) in Equation 16; this reflects an extreme form of mode sparsity, where only one order parameter survives, and corresponds to synchronization on the pattern encoded by the leading eigenvector e 1. Accordingly, analytic examination of marginal densities of P approx ( ψ ) ( ψ ) in Equation 17 shows that a ring-like distribution can arise only in the leading mode coordinate z 1. Furthermore, this ring-like distribution collapses to a unimodal distribution when the leading eigenvalue d 1 is below a threshold d t h ( d 2 , d 3 , … , d N ; N , σ ), which provides an approximate criterion for a stochastic analogue of oscillation death (see Section 3.8 and Li et al. (2022) for details).

H approx ψ ψ = 1 4 N ψ ⊤ O ψ − 1 2 d ⊤ ψ , (16)

P approx ψ ψ = 1 Z d , N , σ exp − 1 2 N σ 2 ψ ⊤ O ψ + 1 σ 2 d ⊤ ψ , (17)

d r d t = − r ⊤ r N 1 − d ◦ r . (18)

2.3 Homotopy continuation: exploring the full landscape

The approximate analysis above provides a coarse picture of the Stuart–Landau oscillator network, highlighting the statistical dominance of a single organized oscillatory macrostate and its degradation under strong noise. However, the full energy function, Equation 4 or Equation 9, and its associated gradient dynamics may exhibit a richer critical-point structure, including multiple local minima and saddles that can mediate synergetic transitions. Resolving the full landscape typically requires numerical tools, and we use homotopy continuation in this work. Consistent with the preceding formulation, we work with H ( r , θ ) in Equation 9 and introduce the parametric homotopy as Equation 19 where the homotopy parameter τ ∈ R continuously deforms the landscape from the approximation ( τ = 0 ) to the full model ( τ = 1 ). Starting from a critical point ( r 0 , θ 0 ) of H ( r , θ ; 0 ), we attempt to trace the corresponding family of critical points as τ varies; that is, we compute a continuation curve ( r ( τ ) , θ ( τ ) , τ ) in the augmented state space, which satisfies the equilibrium conditions An important point here is that the equilibrium condition for the phase θ would normally be ∂ H / ( r ∂ θ ) = 0. However, the approximate energy landscape H ( r , θ ; 0 ) = H approx ( r ) is independent of θ so that ∂ H / ( r ∂ θ ) = 0 holds trivially at τ = 0 and imposes no constraint on θ 0. This violates the nonsingularity condition required by the implicit function theorem for local continuation. We therefore impose the phase equilibrium condition up to the first order in τ using the rescaled form of Equation 20b, which remains equivalent to the original condition for τ ≠ 0 while imposing a critical constraint at τ = 0. For continuation, we differentiate the equilibrium conditions (Equations 20a, b) with respect to τ and set the total derivatives to 0. This yields the tangency conditions where we write H in place of H ( r , θ ; τ ) for brevity. Equations 21a,b form an implicit ordinary differential equation for ( r ( τ ) , θ ( τ ) ) along the solution manifold; it is also referred to as the homotopy system. Numerically integrating the homotopy system traces the desired continuation curves, allowing us to transport critical points from τ = 0 to τ = 1 and thereby compare the approximate and full energy landscapes.

H r , θ ; τ = 1 4 ∑ k = 1 N B r ◦ 2 k + τ r ⊤ W k ◦ C θ r 2 − 1 2 d ⊤ r ◦ 2 , (19)

∂ H r , θ ; τ ∂ r = 0 , (20a)

∂ H r , θ ; τ τ r ∂ θ = 0 . (20b)

d d τ ∂ H r , θ ; τ ∂ r = ∂ 2 H ∂ r 2 d r d τ + ∂ 2 H ∂ θ ∂ r d θ d τ + ∂ 2 H ∂ τ ∂ r = 0 , (21a)

d d τ ∂ H r , θ ; τ τ r ∂ θ = ∂ ∂ r ∂ H τ r ∂ θ d r d τ + ∂ ∂ θ ∂ H τ r ∂ θ d θ d τ + ∂ ∂ τ ∂ H τ r ∂ θ = 0 , (21b)

Next, we briefly introduce the implementation of continuation. To handle possible fold bifurcations on the continuation curve where ∥ d r / d τ ∥ or ∥ d θ / d τ ∥ becomes infinite, we adopt an arc-length l-parametrization, instead of the above τ-parametrization, by multiplying Equation 21a and Equation 21b with d τ / d l. The homotopy system is then expressed aswhere ( r ′ , θ ′ , τ ′ ) = ( d r / d l , d θ / d l , d τ / d l ), or more briefly,where A = [ A i j ] ( i , j = 1,2 ) and b = [ b 1 ⊤ b 2 ⊤ ] ⊤ denote the corresponding blocks of derivatives in the augmented coefficient matrix in Equation 22. Moreover, this 2 N × ( 2 N + 1 ) system only determines the tangential direction of the continuation curve, where we can apply the normalizationto specify a unit tangential vector, with its orientation chosen by Equations 23, 24 are solved and then integrated using the conventional approach of Euler’s predictor and generalized Newton’s corrector, with the algorithm details given in Section 2 of Chen and Li (2015). Generally, both orientations in Equation 25 should be attempted when continuing to τ = 1. Here, however, for each initial configuration ( r 0 , θ 0 , 0 ) with τ ′ ≠ 0, we retain only the orientation for which τ ′ > 0 because any continuation curve that starts with τ ′ < 0 and reaches τ = 1 must bend back and cross τ = 0 at another “initial” configuration ( r 0 , θ 0 , 0 ). In addition, note that the continuation curves in this study should come in pairs due to the reflection symmetry ( r , θ ) → ( r , − θ ), which entails accompanying sign reversals of A 12, A 21, b 2, and θ ′ for Equations 23–25 to remain satisfied.

∂ 2 H ∂ r 2 ∂ 2 H ∂ θ ∂ r ∂ 2 H ∂ τ ∂ r ∂ ∂ r ∂ H τ r ∂ θ ∂ ∂ θ ∂ H τ r ∂ θ ∂ ∂ τ ∂ H τ r ∂ θ r ′ θ ′ τ ′ = 0 , (22)

A 11 r , θ ; τ A 12 r , θ ; τ b 1 r , θ ; τ A 21 r , θ ; τ A 22 r , θ ; τ b 2 r , θ ; τ r ′ θ ′ τ ′ = 0 , (23)

r ′ 2 + θ ′ 2 + | τ ′ | 2 = 1 (24)

sgn det A 11 A 12 b 1 A 21 A 22 b 2 r ′ ⊤ θ ′ ⊤ τ ′ = + 1 or − 1 . (25)

After introducing the homotopy continuation method, we examine the specific form of our energy homotopy and its derivatives. For convenience, we expand H ( r , θ ; τ ) in the homotopy parameter τ as Equation 26: where H 1 ( r , θ ) in Equation 27b is written compactly using V ( r ), which is defined asThis expansion clarifies appropriate initialization as τ → 0. The initial amplitudes r 0 at τ = 0 are chosen as the critical points of H approx ( r ), which have already been prepared in Section 2.2. The phase initialization is then determined at the first order by H 1 ( r , θ ) because the phase equilibrium condition (Equation 20b) yields Equation 29:Note that, the function H 1 ( r 0 , θ ) in θ, for a fixed r 0, has exactly the form of an XY Hamiltonian, which is extensively studied in statistical physics. We therefore obtain θ 0 according to Equation 30:In this study, we do not seek to enumerate all critical points; instead, for each chosen r 0 of interest, we sample random phases θ, apply the generalized Newton corrector within the θ-subspace to obtain possibly multiple candidate θ 0, and then perform homotopy continuation up to τ = 1.

H r , θ ; τ = H 0 r + τ H 1 r , θ + τ 2 H 2 r , θ , (26)

H 0 r = 1 4 B r ◦ 2 ⊤ B r ◦ 2 − 1 2 d ⊤ r ◦ 2 = H approx r , (27a)

H 1 r , θ = 1 2 ∑ k = 1 N B r ◦ 2 k r ⊤ W k ◦ C θ r = 1 2 r ⊤ V r ◦ C θ r , (27b)

H 2 r , θ = 1 4 ∑ k = 1 N r ⊤ W k ◦ C θ r 2 , (27c)

V r = ∑ k = 1 N B r ◦ 2 k W k = E ⊤ Diag B r ◦ 2 E . (28)

∂ H r , θ ; τ τ r ∂ θ = ∂ τ r ∂ θ τ H 1 r , θ + τ 2 H 2 r , θ = ∂ H 1 r , θ r ∂ θ + O τ . (29)

∂ H X Y θ r 0 ∂ θ = 0 , with H X Y θ : = 1 2 r 0 ⊤ V r 0 ◦ C θ r 0 . (30)

Finally, we provide full expressions for the homotopy system (Equation 23) for completeness. Although the separation of functions in r and θ in Equations 27a–c facilitates our physical understanding of the energy homotopy, the coefficient matrix blocks in the homotopy system have far more complex expressions. Nevertheless, we find that by introducing the quantities defined in Equation 31, along with previously defined C ( θ ) in Equation 8 and V ( r ) in Equation 28, H ( r , θ ; τ ) can be rewritten simply in a “quadratic” form as Equation 32 (see also Equation S68 in the Supplementary Appendix):Accordingly, the coefficient matrix blocks in Equation 23 are given by where we have In Equations 33a,b, 34a–c, the explicit dependence of the relevant quantities on r and θ has been suppressed for compactness of notation. We refer the reader to the Supplementary Appendix for detailed derivations of these expressions, together with the necessary preliminaries provided therein.

R = Diag r , S θ = s i j θ = sin θ i − θ j , G r , θ = E ◦ E R C θ , F r , θ = E ◦ E R S θ , (31)

H r , θ ; τ = 1 4 r ⊤ B R + τ G ⊤ B R + τ G r − 1 2 r ⊤ Diag d r . (32)

A r , θ , τ = A 0 r , θ + τ A 1 r , θ + τ 2 A 2 r , θ , (33a)

b r , θ , τ = b 0 r , θ + τ b 1 r , θ = G ⊤ B R + B R ⊤ G F ⊤ G r + τ 2 G ⊤ G r 0 r , (33b)

A 0 = Diag B ⊤ B r ◦ 2 − d + 2 R B ⊤ B R 0 − V ◦ S + 2 F ⊤ B R V ◦ C R − Diag V ◦ C r , (34a)

A 1 = V ◦ C + D i a g B ⊤ G r + 2 G ⊤ B R + 2 B R ⊤ G V ◦ S R + 2 B R ⊤ F R − Diag V ◦ S r − S ◦ E ⊤ Diag G r E + 2 F ⊤ G C ◦ E ⊤ Diag G r E R + 2 F ⊤ F R − Diag G ⊤ G r , (34b)

A 2 = C ◦ E ⊤ Diag G r E + 2 G ⊤ G S ◦ E ⊤ Diag G r E R + 2 G ⊤ F R + Diag F ⊤ G r 0 0 . (34c)

2.4 Additional singularity at initialization: multiple branches

In the previous subsection, a first-order phase equilibrium condition is introduced for the specific form of homotopy (Equation 19) so that the initial phases can be determined from the critical points of an XY Hamiltonian and numerically computed using, for example, the Newton–Raphson method. Following from the implicit function theorem, isolated initial phase configurations can be obtained when the corresponding Jacobian block, given by Equation 35,has rank N − 1, after excluding the degree of freedom associated with the global phase shift. However, this full-rank condition does not necessarily hold, especially in the presence of additional symmetries. When A 22 ( 0 ) ( r 0 , θ 0 ) is singular, the first-order phase equilibrium condition alone is insufficient to determine isolated phase equilibria because flat directions may remain along which the energy gradient varies only at order O ( τ 2 ). In such cases, the O ( τ 2 ) term in homotopy Equation 26 must be taken into account to identify the phase configurations that persist as isolated continuation branches. This motivates the second-order phase equilibrium condition introduced below. In this subsection, we consider the singular cases where ∂ H X Y / ( r 0 ∂ θ ) ≡ 0 on a manifold of dimension greater than one; this implies that the matrix A 22 ( 0 ) ( r 0 , θ ) has a locally constant rank strictly less than N − 1 along this manifold. Here, we assume that the Jacobian block for r at ( r , τ ) = ( r 0 , 0 ), that is, A 11 ( 0 ) ( r 0 ), is nonsingular.

A 22 0 r 0 , θ 0 = ∂ ∂ θ ∂ H X Y θ r 0 ∂ θ | θ = θ 0 = V ◦ C R − Diag V ◦ C r | r = r 0 θ = θ 0 , (35)

2.4.1 Strong singularity

We first consider the extreme case in which the XY Hamiltonian, H X Y ( θ ) = 1 2 r 0 ⊤ V ( r 0 ) ◦ C ( θ ) r 0, is independent of θ. Equivalently, the phase equilibrium condition satisfies ∂ H X Y ( θ ) / ( r 0 ∂ θ ) = − V ( r 0 ) ◦ S ( θ ) r 0 ≡ 0 for all θ, so the Jacobian block A 22 ( 0 ) ( r 0 , θ ) vanishes identically. In this strongly singular case, the initial phase θ 0 cannot be determined by the preceding equilibrium condition, and the continuation with Equation 23 cannot be conducted in the usual way. To de-trivialize continuation at τ = 0, we further rescale the phase equilibrium condition (Equation 20b) and impose the second-order form ∂ H ( r , θ ; τ ) / ( τ 2 r ∂ θ ) = 0. Taking the limit τ → 0 yieldswhere we have used ∂ H 1 ( r 0 , θ ) / ( r 0 ∂ θ ) | θ = θ 0 = 0 for the second equality. Equation 36 shows that the phase equilibrium condition at τ = 0 in the singular case coincides with the phase tangency condition (Equation 21b). In the strongly singular setting, this leads towhere θ and d r / d τ are unknowns whose solutions yield the initial phase θ 0 and the initial tangent ( d r / d τ ) | τ = 0. This system is linear in d r / d τ, while θ enters through the coefficients in a nonlinear fashion. Because A 11 ( 0 ) ( r 0 ) is invertible, we can eliminate d r / d τ and obtain the residual (suppressing the explicit dependence on r 0 for brevity):After fixing the global phase-shift degree of freedom, the nonlinear system f s ( θ ) = 0 comprises N − 1 equations in N − 1 unknowns and can be solved for possibly multiple solutions θ 0, using iterative methods such as Newton’s method and the Gauss–Newton method.

lim τ → 0 ∂ H τ 2 r ∂ θ = lim τ → 0 ∂ H 1 τ r ∂ θ + ∂ H 2 r ∂ θ = d d τ ∂ H 1 r ∂ θ + ∂ H 2 r ∂ θ | τ = 0 = ∂ ∂ r ∂ H 1 r ∂ θ d r d τ + ∂ ∂ θ ∂ H 1 r ∂ θ d θ d τ + ∂ ∂ τ τ ∂ H 2 r ∂ θ | τ = 0 = ∂ ∂ r ∂ H τ r ∂ θ d r d τ + ∂ ∂ θ ∂ H τ r ∂ θ d θ d τ + ∂ ∂ τ ∂ H τ r ∂ θ | τ = 0 = 0 , (36)

A 11 0 r 0 A 21 0 r 0 , θ d r d τ + b 1 0 r 0 , θ b 2 0 r 0 , θ = 0 , (37)

f s θ : = b 2 0 θ − A 21 0 θ A 11 0 − 1 b 1 0 θ . (38)

2.4.2 Weak and general singularity

We now consider the general singular scenario in which the Jacobian block A 22 ( 0 ) ( r 0 , θ ) is rank-deficient, with locally constant rank at most N − 2 on a manifold in the θ-subspace. This occurs when the condition expressed by Equation 30 admits additional degrees of freedom; for instance, two subgroups of oscillation modes may have internally fixed phase patterns while the relative phase difference between the groups can be arbitrary, yielding A 22 ( 0 ) ( r 0 , θ ) with a rank of N − 2. In this setting, the same second-order rescaling argument as Equation 36 applies, leading to the systemat τ = 0, where θ, x, and y are unknowns whose solutions provide θ 0, ( d r / d τ ) | τ = 0, and ( d θ / d τ ) | τ = 0, respectively. Because A 11 ( 0 ) ( r 0 ) is invertible, we can eliminate x and write Equation 39 aswhere the dependence on r 0 is suppressed for brevity. For a given θ, the above equation poses a linear solvability problem in y: letting y ∗ ( θ ) denote any least-squares solution of y for a given θ, we can define the residual as Equation 41, so that the solvability of Equation 40 for y is equivalent to the equation f ( θ ) = 0. Because θ 0 must also lie on the critical manifold for phase equilibrium, we can finally formulate the initialization as the following coupled system in θ: Note that, for a fixed θ, the solvability of Equation 40 for y generally imposes N − 1 − rank A 22 ( 0 ) independent compatibility conditions (after excluding global shifts). In contrast, the local manifold given by Equation 42b itself has dimension N − 1 − rank A 22 ( 0 ). Thus, Equations 42a,b supplies precisely N − 1 constraints to determine isolated admissible phases θ 0 in generic situations. The strongly singular case is recovered as the special case where A 22 ( 0 ) ( θ ) ≡ 0 so that Equation 42b becomes the full θ-subspace and f ( θ ) reduces to f s ( θ ) in Equation 38. Equations 42a,b can be solved iteratively, for example, using a constrained Gauss–Newton-type method, in which each iteration computes a trial step δ θ and performs a line search on it. This trial step is given by the Karush–Kuhn–Tucker (KKT) system in Equation 43:where K ( θ ) is a Jacobian proxy that satisfies ∂ f ( θ ) 2 / ∂ θ = 2 K ( θ ) ⊤ f ( θ ), and μ is an auxiliary vector returned. Further details about this method and explicit expressions for the relevant quantities are provided in the Supplementary Appendix.

A 11 0 r 0 0 A 21 0 r 0 , θ A 22 0 r 0 , θ x y + b 1 0 r 0 , θ b 2 0 r 0 , θ = 0 (39)

A 22 0 θ y + b 2 0 θ − A 21 0 θ A 11 0 − 1 b 1 0 θ = 0 , (40)

f θ : = A 22 0 θ y ∗ θ + b 2 0 θ − A 21 0 θ A 11 0 − 1 b 1 0 θ , (41)

f ( θ ) = 0 , (42a)

∂ H X Y ( θ ) r 0 ∂ θ = 0 . (42b)

K θ ⊤ K θ A 22 0 θ ⊤ A 22 0 θ 0 δ θ μ = − K θ ⊤ f θ ∂ H X Y θ / r 0 ∂ θ , (43)

2.4.3 Mode reduction

When the initial amplitude r 0 has a support (an active set) S with the corresponding inactive set S c = { 1,2 , ⋯ , N } \ S, it may happen that the modes indexed by a subset of S c become structurally degenerate: the corresponding columns of both A 22 ( 0 ) ( θ ) and K ( θ ) vanish identically on a local critical manifold. Hence, the associated phase components do not affect either condition in Equations 42a,b. This results in a singular KKT system and a non-isolated solution set, but this situation can be handled simply by disregarding these degenerate modes and restricting the computation to the remaining ones. In practice, such degenerate modes typically remain silent along the continuation curve, allowing initialization and tracing to be performed as usual in the reduced system. We refer to this as a mode reduction, in which mode sparsity persists beyond the phase-decoupled approximation and into the entire homotopy system.

Because r 0, θ 0, and ( d r / d τ ) τ = 0 are now fixed, our remaining concern would be the determination of d θ / d τ at τ = 0. To this end, we differentiate the corresponding equilibrium condition to consider the phase tangency condition at the second order. Taking the limit τ → 0 and using ∂ H 1 / ( r ∂ θ ) = 0 at ( r , θ ) = ( r 0 , θ 0 ), we obtain Equation 44:or more simply, the following equation for d θ / d τ:Here, the first term is known, whereas the third term depends implicitly on ( d θ / d τ ) τ = 0 and is algebraically cumbersome. In practice, we evaluate this term numerically (e.g., by a central difference approximation); detailed analysis of Equation 45 is beyond the scope of this article. Because tangents such as ( d θ / d τ ) τ = 0 are used only in the Euler predictor at the start of continuation, moderate errors are typically corrected by the subsequent generalized Newton corrector, which projects each predicted point back onto a continuation curve near ( r 0 , θ 0 , 0 ). Finally, because θ 0 may not be unique, multiple admissible initial phases may give rise to multiple continuation branches emanating from the same critical amplitude r 0 of the approximate landscape H 0 ( r ), some of which can be tracked to the full system at τ = 1.

lim τ → 0 d d τ ∂ H τ 2 r ∂ θ = lim τ → 0 ∂ ∂ r ∂ H τ 2 r ∂ θ d r d τ + ∂ ∂ θ ∂ H τ 2 r ∂ θ d θ d τ + ∂ ∂ τ ∂ H 1 τ r ∂ θ + ∂ H 2 r ∂ θ = lim τ → 0 1 τ ∂ ∂ r ∂ H 1 r ∂ θ + ∂ ∂ r ∂ H 2 r ∂ θ d r d τ + 1 τ ∂ ∂ θ ∂ H 1 r ∂ θ + ∂ ∂ θ ∂ H 2 r ∂ θ d θ d τ + ∂ ∂ τ ∂ H 1 τ r ∂ θ = ∂ ∂ r ∂ H 2 r ∂ θ d r d τ + ∂ ∂ θ ∂ H 2 r ∂ θ d θ d τ | τ = 0 + lim τ → 0 d d τ ∂ H 1 τ r ∂ θ = ∂ ∂ r ∂ H 2 r ∂ θ d r d τ + ∂ ∂ θ ∂ H 2 r ∂ θ d θ d τ + d 2 d τ 2 ∂ H 1 r ∂ θ | τ = 0 = 0 , (44)

A 21 1 r 0 , θ 0 d r d τ | τ = 0 + A 22 1 r 0 , θ 0 d θ d τ + d 2 d τ 2 ∂ H 1 r ∂ θ | τ = 0 = 0 . (45)

3 Numerical results

3.1 Basic setting for a frustrated network

Here, we consider a network of N = M 2 coupled oscillators arranged on a 2D square lattice with the periodic boundary condition in both directions (i.e., a toroidal grid). Each oscillator interacts with its four nearest neighbors; this local coupling is encoded by a symmetric, hollow, 4-regular adjacency matrix L ∈ R N × N. Moreover, to incorporate uniform global frustration, which represents long-range mean-field inhibition or antiferromagnetic interaction, we augment the local coupling with an all-to-all interaction term I − O, where O is the all-one matrix, and I is the identity. The total coupling matrix is given in Equation 46: where κ , κ f ∈ R control the overall coupling strength and relative strength of frustration, respectively. The interplay between L and I − O introduces competition between local synchronization and global desynchronization, yielding a rich structure that shapes the collective dynamics of the network. Throughout this article, we fix the overall coupling strength at κ = 0.1 and set the intrinsic oscillation magnitude a i = 1 for all i and the common intrinsic frequency b = π / 50 (period 100).

W = κ 1 4 L + κ f N I − O , (46)

The competing influences of local excitatory and global inhibitory couplings are reflected by the eigenvalue spectrum of the total coupling W. The lattice adjacency matrix L is diagonalizable in the Fourier basis and has well-known eigenvalues given by Equation 47:for j , l ∈ { 0,1,2 , … , M − 1 }; the largest eigenvalue is λ 0,0 ( L ) = 4, with corresponding eigenvector 1 (the uniform mode). The frustration term has eigenvalues { λ ( I − O ) } = { 1 − N , 1,1 , … , 1 }, where the single eigenvalue 1 − N also corresponds to an eigenvector 1. Because L and I − O commute, the eigenvalue spectrum of the total coupling is given in such a structured way as Equation 48:In other words, as κ f increases, the ground mode is penalized by the global frustration, leading to a decreasing eigenvalue, while all the other eigenvalues are shifted upward uniformly. We consider N = 4 2 in this article, with a schematic of the lattice structure shown in Figure 1a and the change in eigenvalues with κ f shown in Figure 1b.

λ j , l L = 2 ⁡ cos 2 π j M + 2 ⁡ cos 2 π l M (47)

λ j , l κ − 1 W = 1 + κ f 1 N − 1 , for j , l = 0,0 , 1 2 cos 2 π j M + cos 2 π l M + κ f N , otherwise . (48)

For N = 4 2, it is easily verified that a common orthonormal basis that simultaneously diagonalizes L and I − O (hence W, for any κ f) can be chosen as the normalized Hadamard matrix, whose entries all have an absolute value 1 / 4. The Hadamard basis, denoted by E H, corresponds to modes for which the oscillation power is evenly distributed across oscillators. It can be obtained by Riemannian minimization of ∑ i , j e i j 4 over the orthogonal manifold. Figure 1c showcases the 2D patterns of these oscillation modes, corresponding to reshaped column vectors of E H, which is shown in Figure 1d. Finally, a general orthonormal basis E, obtained by random rotations of E H within the eigenspaces corresponding to repeated eigenvalues, is shown in Figure 1e; this representative orthonormal basis E, and other randomly sampled ones, will be used to isolate the initial amplitudes r 0, as discussed later.

3.2 Approximate equilibria (H 0 critical points)

Next, we investigate the coupled oscillator lattice starting from the equilibria r 0 to Equation 11, which are used to initialize the continuation of critical points of the homotopy (Equation 19) at τ = 0. With the Hadamard basis E H, both matrices B H ( = E H ◦ E H ) and B H ⊤ B H reduce to N − 1 O. This coincides with the naïve mean-field model in Li et al. (2022), where Equation 18 admits exactly N + nontrivial equilibria when W + Diag ( a ) has exactly N + distinct positive eigenvalues, and Equation 17 describes the stationary distribution even when the eigenvalues are repeated. However, the drawback of this singular structure in our setting is that the active-set system ( B H ⊤ B H ) ( 1,1 ) r ( 1 ) ◦ 2 = d ( 1 ) in Equation 12 typically has no isolated solutions for critical points r 0 of H 0 ( r ), except in the eigenspace of a single eigenvalue. This renders the Jacobian block A 11 ( 0 ) ( r 0 ) in Equation 39 singular, violating the fundamental requirement for continuation. To avoid such a degeneracy, we construct H 0 ( r ) using the randomly generated orthonormal basis E shown in Figure 1e; with probability one, B ⊤ B then has 32,003 nonsingular principal submatrices, while H ( r , θ ) in Equation 9 remains invariant.

We calculated the solutions r 0 to Equation 12 and plotted the corresponding energy values H 0 ( r 0 ) in Figure 2, for κ f = 0.1, 0.4, and 0.7. Each solution is labeled with an integer in { 0,1,2 , … , 2 16 − 1 }, whose base-2 representation specifies the active set (i.e., the indices of nonzero components of r 0), with the rightmost bit corresponding to the leading vector component. For example, a solution of the form [ 0 r 2 r 3 r 4 r 5 0 ⋯ 0 ] ⊤ with r 2 , r 3 , r 4 , r 5 > 0 is labeled by ( 11110 ) 2, that is, 30. Under the current parameter setting, a random basis E typically yields roughly ( 3.3 ± 0.6 ) × 1 0 3 feasible solutions r 0; in our calculation, the numbers of solutions are about 4.3 × 1 0 3. Although varying κ f changes the range of eigenvalues in d, Figure 2 shows that the overall energy patterns of the approximate equilibria remain broadly similar. Among solutions with labels in 2 n , 2 n + 1 ( 1 ≤ n ≤ 14 ), solution #2 n typically has a particularly high H 0. Meanwhile, solutions that are close to “all-one” (up to scaling) within the selected eigenspace tend to achieve lower H 0, especially the small-numbered ones with nonzero leading modes (e.g., solutions #1 and #30).

Varying the frustration strength κ f modifies the details of the energy pattern in Figure 2. A notable change is the relative ordering of solutions #1 and #30: solution #1 has the lowest H 0 when κ f < 0.5, while the energy gap between them decreases as κ f increases. Once κ f > 0.5, solution #30 attains a lower H 0 than solution #1, and a series of odd-labeled solutions emerge (e.g., #3, #5, … , #29, #93, … , #1551). These solutions are accompaniments to solution #1, taking the form [ r 1 ε 2 ε 3 ⋯ ε 11 0 ⋯ 0 ] ⊤ with 0 ≤ ε 2 , ε 3 , … , ε 11 ≪ r 1. This trend aligns with the intuition that the preferred configuration of the coupled oscillators shifts from full to partial synchronization, as global frustration increases.

3.3 Homotopy continuation instances

Now, we perform homotopy continuation from selected solutions r 0 and examine how the equilibria of H 0 ( r ) (i.e., H ( r , θ ; 0 )) continue to those of H ( r , θ ; 1 ). We focus mainly on solutions #1 and #30, which attain the lowest and second-lowest values of H 0, and also consider solutions #28, #24, and #32 as representative higher-energy examples. For κ f = 0.7, we additionally investigate the small-odd-numbered solutions accompanying solution #1, for example, solution #25, which turns out to be weakly singular.

3.3.1 Solution #1

This is a highly degenerate case exhibiting strong singularity and mode reduction. At τ = 0, for any θ, all blocks in the coefficient matrix in Equation 22 vanish except A 11 ( 0 ) ( r 0 ). In particular, condition Equation 37 implies d r / d τ | τ = 0 = 0, while θ 0 remains arbitrary. Moreover, the above-mentioned vanishing blocks remain zero for finite τ. Consequently, the continuation reduces to d r / d τ ≡ 0, and the homotopy system effectively collapses to the first mode, with equilibrium #1 persisting for all τ ∈ [ 0,1 ]. Numerical calculation shows that this equilibrium maintains an index of 0 (i.e., no negative eigenvalues) throughout τ ∈ [ 0,1 ] when κ f = 0.1 and 0.4; for κ f = 0.7, the index is initially 8, decreases to 4 at a small τ, and then remains 4 up to τ = 1. The energy values are H = − 4.758, − 4.516, and − 4.280 for κ f = 0.1, 0.4, and 0.7, respectively.

3.3.2 Solution #30

Similar to solution #1, solution #30 at τ = 0 remains invariant for any choice of orthonormal basis E. It is also a strongly singular case in the sense that A 22 ( 0 ) ( r 0 , θ ) ≡ 0 for any θ. On the other hand, the choice of E affects the admissible initial conditions that satisfy Equations 37, 45. With the basis E in Figure 1e and κ f = 0.4, we could find 50 initial phase configurations: 8 with collinear phases and 21 non-collinear configurations together with their sign-flipped counterparts (arising from the reflection symmetry). Here, we fix θ 2 = 0 to remove the global phase-shift degree of freedom. In all successful runs, the homotopy system reduces to at most eight modes, with a potential support { 2,3,4,5,12,13,14,15 }; the continuation curves that reach τ = 1 can be grouped into the following five types:

• The homotopy system is reduced to four modes ( r 2 , r 3 , r 4 , r 5 ≠ 0 ); the coefficient matrix along the continued curve has rank 19 (A 22 ( 0 ) has rank 3) in a neighborhood of τ = 0;
• An 8-mode reduction, with non-collinear initial phases and vanishing r 12, r 13, r 14, and r 15 at τ = 1;
• An 8-mode reduction, with collinear initial phases and vanishing r 12, r 13, r 14, and r 15 at τ = 1;
• An 8-mode reduction, with non-collinear initial phases and non-vanishing r 12, r 13, r 14, and r 15 at τ = 1;
• An 8-mode reduction, with collinear initial phases and non-vanishing r 12, r 13, r 14, and r 15 at τ = 1.

Here, by “collinear,” we mean that the relevant components of θ are { 0 , π }-locked. For cases (ii)–(v), the coefficient matrix along the continued curves has rank 23 (A 22 ( 0 ) has rank 7) in a neighborhood of τ = 0.

Figure 3 summarizes typical continuation trajectories from solution #30 with κ f = 0.4. For each previously mentioned type (i)–(v), we plot the mode amplitudes r, phases θ, the equilibrium’s (negative) eigenvalues, and the energies H ( r , θ ; τ ) as functions of the homotopy parameter τ ∈ [ 0,1 ]. The quantity H 0 ( r ) is also evaluated along the same curve and plotted together with H ( r , θ ; τ ) to illustrate the magnitude of the θ-dependent cross terms in Equation 19 disregarded under the phase-decoupled approximation. The continuation in (i) leads to an equilibrium with no negative eigenvalues at τ = 1, which can be confirmed as an energy local minimum, while index-1 saddles are found in cases (ii) and (iii). These three equilibria share the same energy value H = − 4.431 (which equals the initial energy H 0), and have one, three, and three additional zero eigenvalues, respectively, implying a degeneracy that will be discussed later. Meanwhile, index-3 and index-5 saddles with higher energies are obtained in (iv) and (v). To aid interpretation, we transform back to the Hadamard basis E H and visualize the continuation curves in mode-wise polar coordinates in Figure 4a. We then further reconstruct the corresponding oscillator configurations in the original coordinates (see Figure 4b). As a result, the equilibria obtained at τ = 1 admit particularly simple representations in the Hadamard basis: case (i) involves two modes, while cases (ii) and (iii) each reduce to a single mode. Specifically, (ii) and (iii) correspond to the fifth and second Hadamard patterns in Figure 1c, while the local minimum in (i) is a quadrature superposition of the third and fourth Hadamard patterns, forming four locally synchronized clusters as shown in panel (i) of Figures 4a,b. We note that these results are representative, and other similar minima and saddles also exist owing to symmetry.

Similar to the case of solution #1, the frustration strength κ f has a profound influence on the continuation from solution #30. Calculations show that, although the continuation curves r ( τ ) and θ ( τ ) differ only slightly from those in Figure 3, significant changes appear in the eigenspectra of the equilibria when κ f increases from 0.1 to 0.7, as shown in Figures 5a–c. In all cases, increasing κ f lifts certain eigenvalues toward positive values; consequently, the number of negative eigenvalues tends to decrease, resulting in increased stability (in the sense of fewer unstable directions in the state space). For example, the index-1 equilibrium at τ = 1 in Figure 5a (case (i)) with κ f = 0.1 becomes index-0 in Figure 5b with κ f = 0.4. Similar changes occur in (ii) and (iii) for stronger frustration κ f = 0.7 in Figure 5c, when the index-5 equilibrium in (v) becomes index-4.

3.3.3 Solution #25

This solution appears for sufficiently strong frustration such as κ f = 0.7, where A 22 ( 0 ) has a rank of 13 (and the coefficient matrix has rank 29), corresponding to a weak singularity. Two additional degrees of freedom emerge for satisfying the equilibrium condition to the first order: θ 16 is unconstrained, and θ 1 , θ 2 , … , θ 15 fall into two internally collinear groups with a free inter-group phase offset. We could find eight initial θ 0 satisfying the equilibrium condition to the second order, among which two pairs of quadrature-locked θ 0 extend to τ = 1. In these cases, the homotopy system exhibits no mode reduction, and the equilibria are initially index-4 saddles that become index-2 saddles at τ = 1. A typical continuation result, including the final oscillator configuration, is illustrated in column (i) of Figures 6a,b.

3.3.4 Solutions #28, #24, and #32

These solutions all correspond to singular initialization where A 22 ( 0 ) has a rank of 7 (the coefficient matrix has rank 23). All additional degrees of freedom reside in the phases of single modes. The homotopy system reduces to nonsingular 8-mode systems, with support { 2,3,4,5,12,13,14,15 } for solutions #28 and #24, or { 1,6,7,8,9,10,11,16 } for solution #32. For solution #28, we can find three collinear θ 0 that satisfy the second-order equilibrium condition and extend to τ = 1. The equilibrium has index 16 initially and continues to index 10–16 at τ = 1. For solution #24, two such collinear θ 0 are found; the equilibrium has index 24 initially and continues to index 20–23 at τ = 1. For solution #32, only one θ 0 is found; starting with index 30 and H 0 = − 1.567, the equilibrium still has a high index of 29 and a high energy H = − 0.494 when reaching τ = 1. Varying κ f within [0.1,0.7] makes no fundamental difference. Typical continuation results, including the corresponding oscillator configurations, are shown for κ f = 0.7 in Figures 6a,b as cases (ii), (iii), and (iv).

3.4 Partial selectivity of the continuation

From the preceding results, we have observed a possible selectivity in the energy and index of the equilibria reached by continuation: starting from an approximate solution r 0 at τ = 0 with relatively high (or low) energy, the equilibria obtained at τ = 1 also tend to have relatively high (or low) energy and indices. To examine this effect, we performed homotopy continuation from a collection of approximate solutions randomly sampled from the full set. Specifically, for κ f = 0.4, we resampled about 430 approximate solutions (i.e., 10 % of the total) using importance sampling to obtain an approximately uniform distribution of the initial energy H 0. Among these, about 360 equilibria r 0 with about 1.2 × 1 0 4 initial phases θ 0 (among a total of 1.3 × 1 0 4) produce about 2.6 × 1 0 3 equilibria at τ = 1. The results are summarized in the hexagonal binning plots shown in Figure 7.

Figure 7a shows that, in the target system ( τ = 1 ), the index and energy of an equilibrium are positively associated with a roughly linear trend; both global and local index-0 minima are also evident from the plot. Figures 7b,c further indicate that both the target energy H ( r , θ ; τ ) τ = 1 (or H ( τ = 1 ) for brevity) and the target index increase with the initial energy H 0 ( r 0 ) (or, H ( τ = 0 )). Finally, Figure 7d shows the comparison of the target and initial indices: although the relationship exhibits substantial scatter, a clear positive correlation remains. While Figure 7 is based on continuation at κ f = 0.4 with E chosen as in Figure 1e, qualitatively similar behavior is observed for other random bases and other values of κ f of interest. Together, these observations suggest partial selectivity of this eigenspace-based homotopy continuation procedure, in the sense that it preferentially finds equilibria with higher (or lower) energy and, correspondingly, higher (or lower) index, depending on the initialization energy. Such selectivity may greatly facilitate the exploration of low-index critical points, including minima and index-1 saddles, which usually represent macrostates and transition mediators.

The observed partial selectivity is closely related to how the approximation is constructed. The phase-decoupled approximation is formulated in the eigenspace coordinates of the linear part of the dynamics and is not an arbitrary surrogate. In this representation, the effect of the linear coupling term, which is central to shaping the collective energy landscape, is already partially incorporated, albeit at the cost of a modified higher-order structure. Therefore, the approximate and full landscapes remain structurally close in the sense relevant to continuation, and the homotopy acts as a structured deformation rather than a complete reorganization, along which some properties, such as the energy and index, undergo only limited changes. Consequently, equilibria with a relatively low or high energy and index in the approximation tend to continue to equilibria with similar properties in the full system.

3.5 Coverage of low-index equilibria

To better understand the multistability of the system, we usually aim to identify as many low-index equilibria as possible. However, in the presence of repeated eigenvalues of W, initialization based on a single randomly rotated E may reveal only a subset of equilibria in the full system, as illustrated by the foregoing results. Moreover, the minima and index-1 saddles shown in Figure 3 possess additional zero eigenvalues and identical energy, indicating degeneracy and possibly non-isolation. We, therefore, perform a coverage analysis of the low-index equilibria obtained by continuation from low-energy approximate solutions across repeated continuation trials with different random bases E (up to roughly 200 choices of E and 8 × 1 0 3 trials). To identify equilibria while properly accounting for the invariance under global phase shifts, we define the distance between two equilibria z 1 = r 1 ◦ exp ( i θ 1 ) and z 2 = r 2 ◦ exp ( i θ 2 ) bywhere the minimum is attained at γ = γ ∗ = − arg ( z 1 † z 2 ). In other words, D ( z 1 , z 2 ) is the minimum Euclidean distance in state space between the two deterministic oscillatory solutions represented by z 1 and z 2. It is also invariant under the orthonormal transformation z = E ⊤ s. We adopt a numerical threshold value D t h = 1 × 1 0 − 3 to determine whether a newly obtained equilibrium should be recorded as distinct. Because there is no prior knowledge of which equilibrium will be obtained in each trial, we treat the process as a finite-population discovery model under sampling with replacement. In each continuation trial, a particular equilibrium type is drawn with probability p. The count q of distinct equilibria discovered then increases with the number of draws (trials) m and eventually saturates. The expected value of q ( m ) is given bywhere n denotes the total number of distinct equilibria of interest. Here, the number of draws m may be defined either as the number of continuation trials that successfully reach τ = 1 or as the total number of continuation trials including those that terminate without yielding an equilibrium, with corresponding interpretation of p. We adopt the former convention in the following analysis. For convenience, we use subscripts on D and q to label the starting solution for continuation, and we use superscripts to denote the indices of the target equilibria in the discussion below.

D z 1 , z 2 : = min γ ∈ R z 1 − exp i γ z 2 = z 1 2 + z 2 2 − 2 | z 1 † z 2 | = r 1 2 + r 2 2 − 2 r 1 ◦ r 2 ⊤ ⁡ exp i θ 1 − θ 2 , (49)

E q m = n 1 − 1 − p n − 1 m ≈ n 1 − exp − m p n − 1 , (50)

At the frustration strength κ f = 0.4, we focus on continuation from approximate solution #30, which persists under any randomly rotated E and consistently has the lowest initial energy except for the non-branching starting solution #1. Figure 8a plots the distinct counts q 30 0 ( m ) and q 30 1 ( m ) of index-0 and index-1 equilibria obtained after m continuation trials, shown as mean values with bootstrapped sample standard deviations. Both q 30 0 ( m ) and q 30 1 ( m ) increase with little tendency of saturation, and the resultant target equilibria all have two or four zero eigenvalues (i.e., one or three eigenvalues modulo the global phase shift); fitting Equation 50 is either unsuccessful or yields quite large n, suggesting non-isolated equilibria. Figure 8b shows the histogram densities of the distances D 30 0,0 among index-0 equilibria, D 30 1,1 among index-1 equilibria, and D 30 0,1 between index-0 and index-1 equilibria, all of which exhibit continuous distributions. The noticeable density for small values of D 30 0,1 suggests adjacency between the index-0 and index-1 sets, while the multimodal distribution of D 30 0,0 and D 30 1,1 suggests isolated substructures within them. For comparison, Figures 8c,d show a saturating sequence q 30 2 ( m ) and discrete values of D 30 2,2 for index-2 equilibria continued from solution #30, all with no additional zero eigenvalues, suggesting that they are isolated. The empirical q 30 2 ( m ) profiles are well described by the discovery model with n = 96 and p = 0.106; for a typical sequence, we also computed the posterior probability that n > q 30 2 ( m ) under this model and found that it decreases to about 0.02. This suggests that it is unlikely that more than these 96 index-2 equilibria remain undiscovered.

Coverage analysis is also conducted for low-index equilibria at κ f = 0.7, as summarized in Figures 8e,f. In this case, the energy ordering of approximate solutions #1 and #30 reverses, as shown previously in Figure 2. Continuation from the minimum-energy solution #30 again produces 96 index-2 equilibria, together with a set of index-0 equilibria whose distance histogram resembles the overall one in Figure 8b. This suggests that the stabilization shown in Figure 5 has occurred for all previously identified index-1 equilibria. Meanwhile, we also investigate continuation from the three approximate solutions with the next-lowest energies after #30 under various random choices of E (with labels ranging from ∼ 10 to ∼ 20000, collectively denoted by “OL”). As shown in Figure 8e, we find 24 isolated index-1 equilibria and an additional 80 isolated index-2 equilibria, for which the posterior probabilities that more remain undiscovered decrease to nearly 0 and 0.012, respectively. These index-1 saddles, all with H = − 4.328, may lie on the ridges separating multiple subsets of the equi-energy index-0 minima at H = − 4.447. Our observations in this subsection also clarify that the discovery behavior can be meaningfully characterized only when the target equilibria are isolated modulo the global phase shift, as indicated by the absence of additional zero eigenvalues. We note that, in such cases, Equation 50 can be fitted once a clear saturation trend appears, even before full saturation is reached, to estimate the total number n of distinct equilibria.

3.6 Energy landscapes: embedding-based visualization

After we have identified the low-index equilibria, we next investigate their geometric organization using manifold embedding (nonlinear dimensionality reduction), using the pre-computed distances defined in Equation 49. In principle, manifold embedding seeks a low-dimensional Euclidean representation that preserves these distances as well as possible. Commonly used methods include t-distributed stochastic neighbor embedding (tSNE), multidimensional scaling (MDS), and Isomap. Figures 9a–c show the embeddings of the low-index equilibria at κ f = 0.4, obtained using two-dimensional tSNE, Isomap, and MDS, respectively. We find that tSNE, which is well suited to capturing local structures, reveals four six-spoked, asterisk-shaped sets of index-1 equilibria and 12 segments of connected index-0 equilibria, although the global arrangement of these sets appears irregular. In contrast, Isomap reveals a cleaner global structure, namely, the skeleton of a 3-simplex, or a complete graph on four vertices ( K 4 ), where the index-1 equilibria occupy regions near its vertices and the index-0 equilibria lie in each of its six edges. The geometry of the full equi-energy set of index-0 and index-1 equilibria is most clearly recovered by MDS, which yields a complete multigraph with four vertices and 12 edges; this structure is captured even more clearly by three-dimensional MDS, as shown in Figures 9d,e.

These results are also consistent with the number of additional zero eigenvalues found previously in Figure 5, namely, one flat direction for the index-0 equilibrium and three for the index-1 equilibria. We also confirm that the vertices represent the Hadamard modes or basis vectors in E H corresponding to the second largest eigenvalue of L (with multiplicity 4), and a point on an edge represents the superposition of two Hadamard modes with certain amplitudes and phase offset Δ θ; this has been suggested by previous results shown in Figure 4. At κ f = 0.7, three-dimensional MDS reveals a similar geometric organization, in which the entire graph is stabilized, and newly appearing isolated index-1 equilibria emerge around it, as shown in Figures 9f,g. We therefore infer that, at intermediate frustration, the energy landscape mainly consists of an isolated global minimum (i.e., the full synchronization) and 12 flat 1D troughs (i.e., proper superposition of two Hadamard modes with sufficiently large | Δ θ |), separated by unstable breaches. At sufficiently large frustration, these breaches heal, and the troughs merge into a connected set of global minima.

3.7 Simulated dynamics: synergetics and transitions

3.7.1 Spatiotemporal dynamics under noise

The results above show that varying the frustration strength κ f profoundly deforms the energy landscape of the coupled oscillator network. Under sufficiently weak frustration, the fully synchronized state is the unique energy minimum. As κ f increases, the synchronized state gradually loses stability, while some low-index, low-energy saddles on a complete multigraph progressively stabilize. At an intermediate frustration level (e.g., κ f = 0.4), saddles in the middle portions of the edges become index-0 and coexist with the synchronized state as additional local minima. For stronger frustration (e.g., κ f = 0.7), all saddles on the graph turn into index-0 equilibria, whereas the synchronized state turns into a saddle. In contrast, the high-index saddles undergo no qualitative change over this range, and we do not consider larger values of κ f in this article. Such an energy-based structure supports rich spatiotemporal dynamics, including emergent (synergetic) rhythms and noise-driven state transitions.

To illustrate this, we numerically integrated Equation 1 using the Euler–Maruyama scheme with time step δ t = 0.01. Figures 10, 11 show representative spatiotemporal patterns, together with the corresponding time evolution of the energy and the projections (or their amplitudes) onto the Hadamard basis vectors, to aid visualization of transitions. In Figures 10a–d, we apply transient noise over the interval t ∈ [ 100,150 ]. For κ f = 0.4, the network initialized at the index-0 equilibrium (i) in Figure 3 remains near its initial state under a small perturbation (σ = 0.05; Figure 10a); meanwhile, the amplitudes of the two constituent Hadamard modes and their phase offset undergo slight changes before and after the perturbation. This is consistent with a local minimum in a trough, where small transient noise induces a slight displacement along the flat direction. By contrast, when initialized at the index-1 saddle (ii) in Figure 3, which is in a breach, the system lingers there for some time but readily relaxes to synchronization after a small perturbation (Figure 10b). A substantially larger perturbation is required to induce synchronization from the index-0 local minimum (σ = 0.3; Figure 10c), which drives the trajectory out of the current energy trough. For κ f = 0.7, synchronization itself is a saddle, and even small noise induces desynchronization (σ = 0.05; Figure 10d), in contrast to the noise-induced synchronization in Figure 10b. In Figure 11, stationary i.i.d. Gaussian noise is applied over the entire time range. For κ f = 0.4 and relatively large noise σ = 0.3, the spatiotemporal dynamics appears to be consistent with a Gibbs (Boltzmann) distribution: the global energy minimum (synchronization) is visited transiently from time to time, interspersed with visits to other patterns. For κ f = 0.7 and intermediate noise ( σ = 0.12 ), the network exhibits persistent wandering on the graph of troughs: typically, two Hadamard modes attain appreciable amplitudes at the same time, indicating a mixed state on an edge between two vertices, while the energy is approximately maintained within [ − 4.4 , − 4.3 ].

3.7.2 Stationary transition kinetics

In this subsection, we conduct a simple investigation of the system’s transition mechanisms through the transition kinetics of its stationary dynamics, focusing in particular on their scaling with the noise level. We adopt a simple relative-distance-based partition of the state space, with a hysteresis rule to reduce flickering. Specifically, let D i ( t ) denote the distance from the network state to the index-0 equilibrium set S i; when D i ( t 1 ) < α min j ≠ i D j ( t 1 ), the system is classified as entering macrostate i at time t 1. It is then labeled as exiting macrostate i at the first time t 2 > t 1 such that D i ( t 2 ) > β min j ≠ i D j ( t 2 ). The sojourn time T s, defined as the time spent in a macrostate from entry to exit, and the first-passage time T p from macrostate i to j, defined as the time from entry into macrostate i to the first subsequent entry into macrostate j, can then be analyzed statistically. We numerically integrated Equation 1 with time step δ t = 0.1 over a duration up to Δ t = 2 × 1 0 5 (after discarding a transient length of 5 × 1 0 3), and obtained the following results.

For κ f = 0.4, we adopt a coarse-grained description in which the union of all energy troughs is treated as a single macrostate (ET), together with the macrostate corresponding to the global minimum (GM), using α = 0.7 and β = 1.0. Figures 12a,b show the sojourn time in GM, T s GM, the first-passage time from GM to ET, T p GM → ET, and the ratio of sojourn times T s GM / T s ET as functions of 1 / σ 2, plotted on log–log and semi-log scales, respectively. For relatively small σ, both T s GM / T s ET and T p GM → ET can be well described by Equation 51:as predicted by simple effective two-state models with effective energy differences ( Δ H ) eff ( 1 ) and ( Δ H ) eff ( 2 ). In practice, these effective energy differences are estimated by fitting a straight line to ln ( T s GM / T s ET ) versus 1 / σ 2 or ln ⁡ T p GM → ET versus 1 / σ 2 over the low-noise regime, with ( Δ H ) eff estimated as one-half of the fitted slope. In the present case, we obtained ( Δ H ) eff ( 1 ) = ( 5.45 ± 0.17 ) × 1 0 − 2 over 1 / σ 2 ∈ [ 36,100 ] and ( Δ H ) eff ( 2 ) = ( 5.49 ± 0.05 ) × 1 0 − 2 over 1 / σ 2 ∈ [ 9,81 ]. These results are consistent with an Arrhenius–Kramers law, indicating thermally activated escape. The effective energy differences are smaller than the actual energy gap of 0.0846 between the global minimum and the equi-energy graph, likely because of the coarse partition of the state space to define the two macrostates. In contrast, for sufficiently large σ, T s GM and T p GM → ET appear to follow power-law relations given in Equation 52:with C 1 = 1.01 ± 0.03 and C 2 = 0.703 ± 0.027 over 1 / σ 2 ∈ [ 1,6 ]. This behavior, with an exponent close to 1, suggests diffusion-dominated motion under strong noise. In contrast, the smaller exponent C 2 may reflect the increasing difficulty of reaching ET as the accessible state space and the effective transition region expand for large σ. Moreover, the sojourn time T s ET shown in Figure 12c exhibits diffusion-like behavior over a broad range 1 / σ 2 ∈ [ 1,49 ]. A fit of the form ( 1 / σ 2 ) C 3 gives C 3 = 1.08 ± 0.04, indicating diffusion-limited escape from ET. Meanwhile, the decrease of T s ET for 1 / σ 2 > 49 is possibly attributable to the imperfect state-space partition, whereby thermally driven excursions outside the troughs can be counted as stays in, or escapes from, ET.

T s GM σ T s ET σ ∝ exp 2 Δ H eff 1 σ 2 , T p GM → ET σ ∝ exp 2 Δ H eff 2 σ 2 , (51)

T s GM σ ∝ 1 σ 2 C 1 , T p GM → ET σ ∝ 1 σ 2 C 2 , (52)

For the larger frustration κ f = 0.7, at which the entire multigraph forms a connected trough, no natural distinction between metastable states remains. Nevertheless, we can still examine the residence within prescribed regions, or quasi-states. To do so, we naïvely partition the state space to form 16 quasi-states, associated with the four vertices and 12 edge midpoints, using the common parameters α = 0.95 and β = 1.0. The sojourn time T s VT near a vertex and the first-passage time from a vertex to one of its nearest edge midpoints T p VT → MD are also shown in Figure 12c. Both quantities exhibit nearly linear behavior in the log–log plot over almost the full investigated range 1 / σ 2 ∈ [ 1,100 ], consistent with the previously mentioned scaling relation with fitted exponents C 4 = 0.966 ± 0.008 and C 5 = 0.979 ± 0.008, respectively. This supports the interpretation of diffusion motion on the graph. Overall, despite possible inaccuracies introduced by the naïve partition of the state space, these results suggest an Arrhenius–Kramers-like scaling for uphill escape at low noise and a diffusion-dominated scaling for motion at large noise or on the equi-energy graph; these mechanisms are also consistent with the previously obtained energy landscapes.

3.8 Phase-decoupled approximation revisited

Finally, we return to the phase-decoupled approximation used to initiate the continuation and revisit its performance in describing the statistical properties of the oscillator network dynamics. Specifically, we examine the modeled properties of z H and compare them with those computed from the projections z H = E H ⊤ s with simulated s. Here, we fix κ f = 0.4 and focus on the following aspects:

• Simultaneous diagonalization. The phase-decoupled approximation implicitly assumes simultaneous diagonalization of the coupling matrix W and the covariance matrix Cov ( Re ( s ) ) or Cov ( Im ( s ) ) because the amplitude dynamics (Equation 11) does not favor any particular phases and thus predicts uncorrelated components of z H. Figure 13a shows the ranges of the absolute Pearson correlation coefficients for pairs of components of Re ( s ) and Re ( z H ) = E H ⊤ Re ( s ) obtained from the simulation. Even in the deterministic limit, approximately 95 % of these correlation coefficients after projection onto E H have absolute values below 0.04, supporting this assumption.
• Variance of Re ( z H ) , that is, 1 2 E [ ψ H ] ( ψ H = r H ◦ r H ) . The approximate energy model (Equation 13) induces the exponential-family distribution P approx ( ψ ) ( ψ ) ∝ exp ( − 2 H approx ( ψ ) ( ψ ) / σ 2 ) for ψ. This remains valid even when B ⊤ B is singular, provided that the distribution exists, that is, it is integrable. For the network and basis E H considered here, this gives rise to the degenerate distribution Equation 17, which can be sampled numerically using a Metropolis–Hastings sampler with a ridge-perturbed random-scan Gibbs proposal, or treated analytically as in Li et al. (2022). Figure 13b shows the comparison of E [ ψ i H ] (i = 1,2,6,12,16 for different eigenvalues) calculated from simulated z H = E H ⊤ s with the corresponding values obtained from the approximation (either analytically or by sampling). The phase-decoupled approximation captures the qualitative changes in the variances but exhibits a substantial lag in the growth of the dominant mode as the noise level decreases, similar to the results in Li et al. (2022).
• Covariance Cov ( Re ( s ) ) . Under the phase-decoupled approximation, the covariance matrix of the oscillators is given by Cov ( Re ( s ) ) = 1 2 E H Diag E [ ψ H ] E H ⊤. We compare this with the covariance matrix obtained from simulation, using the dimension-normalized affine-invariant Riemannian metric (AIRM) distance for positive definite matrices M 1 , M 2 ∈ R n × n, defined as Equation 53:where λ i are the generalized eigenvalues defined by M 2 v = λ M 1 v, and with det ( ⋅ ) denoting the determinant. In other words, Equations 54a, b decompose the log-spectrum underlying the AIRM distance into its mean and standard deviation, which represent the overall scale change and the anisotropic deformation, respectively. Figure 13c shows d scale and d shape between the simulated covariance Cov ( Re ( s ) ) and that obtained from the approximation, as functions of 1 / σ. We observe notable scale mismatches in both the large- and small-noise regimes. The error accompanying the slow growth of the dominant mode primarily stems from the scale, whereas the shape error remains comparatively small.
• Probability integral transform (PIT) diagnostics for truncated normality. Although the probability distribution for ψ given by Equation 14 may be degenerate (when B ⊤ B is singular), the conditional distribution of each ψ i given the remaining components is still a truncated normal shown by Equation 55: For each sampled ψ i from the simulation, we then compute u i : = F i ( ψ i | ψ j , j ≠ i ), where F i denotes the cumulative distribution function of the truncated normal above. If the approximation were exact, the resulting u i would follow the uniform distribution on [ 0,1 ] for all i = 1,2 , … , N. Figure 13d shows the histogram of all u i values computed for ψ i H at σ = 0.25. The distribution deviates from uniformity with a bias toward smaller values and a mean of 0.33. Across different i and a range of σ, the mean value of u i lies roughly in the interval [ 0.25 , 0.45 ], indicating that the conditional truncated normality implied by the phase-decoupled approximation is not exactly followed.
• Threshold σ t h for stochastic oscillation death. It is natural to expect that, as the noise level σ increases, the deterministic synchronized oscillation is eventually overwhelmed by noise and becomes invisible. In other words, the distribution of the dominant mode z 1 H changes from ring-like to unimodal. Under the phase-decoupled approximation, the marginal distribution of z 1 H predicts that the largest eigenvalue d 1 of W + Diag ( a ) should satisfy d 1 > d t h ( d 2 , d 3 , … , d N ; N , σ ) for the distribution to remain ring-like; for fixed d, this yields a threshold value σ t h for σ. Figure 13e shows the comparison of the theoretical threshold 1 / σ t h = 5.15 with the result of Hartigan’s dip test for unimodality, which yields a p-value below 0.05 when 1 / σ ≳ 4.9 (thereby rejecting unimodality), in close agreement with the approximation. Meanwhile, the projections on other Hadamard basis vectors, such as z 2 H and z 6 H, always have a unimodal distribution.

Overall, our numerical investigations in the present and previous studies suggest empirically that the phase-decoupled approximation performs well quantitatively when couplings and parameter heterogeneity are weak, and the noise level is at least intermediate, namely, when Equation 56 is satisfied:where ρ ( ⋅ ) denotes the spectral radius and a ̄ > 0 and d ̄ > 0 denote the mean values of the components of a and d, respectively (a ̄ = d ̄ because W has zero diagonal), and Equation 57 holds:These results show that, although the phase-decoupled mode approximation is not exact, it yields a probability distribution that captures the statistical properties of the full system at a rough semi-quantitative level. We therefore view this approximation as a surrogate system that is energetically and statistically close to the full system, making it a plausible basis for constructing the homotopy used in the continuation.

ρ W + Diag a − a ̄ I a ̄ = d − d ̄ 1 ∞ d ̄ ≲ 0.05 , (56)

σ a ̄ ≳ ‖ d − d ̄ 1 ‖ ∞ d ̄ , or σ ≳ d ̄ ‖ d − d ̄ 1 ‖ ∞ . (57)

4 Conclusion and discussions

4.1 Energy-based coupled oscillator network and its homotopy continuation

The Stuart–Landau oscillator network with homogeneous frequency and symmetric coupling provides a simple yet expressive gradient system for modeling collective dynamics from the viewpoint of its energy function. In this work, we revisited this classical setting and analyzed its energy landscape using a combination of a mode-based, phase-decoupled approximation and homotopy continuation:

• We introduced the parametric homotopy Equation 19 (and its equivalent form Equation 32) in the eigenbasis of the coupling matrix, which continuously deforms a tractable mean-field mode-amplitude system at τ = 0 into the full system at τ = 1.
• We developed a systematic approach to initialize continuation at τ = 0 from a singularity, arising from non-isolated equilibria in the phase subspace (after excluding the global phase-shift degree of freedom), by incorporating a second-order phase equilibrium condition.
• We numerically tracked how the approximate equilibria split into multiple equilibria of the full system, and how their energy values and indices change with the frustration strength κ f in a 2D toroidal lattice with local coupling and global frustration.
• By launching extensive continuations, we revealed a partial selectivity of our proposed homotopy: approximate solutions with relatively high (or low) energy and indices at τ = 0 continue to target equilibria at τ = 1 that also have relatively high (or low) energy and index.
• We identified the metastable energy landscape of the 2D lattice, at an appropriate frustration strength, as consisting of an isolated global minimum and a family of 1D energy troughs arranged along the middle portions of the edges of a multigraph, with transitions involving both the Arrhenius–Kramers and diffusion-dominated mechanisms.

In our investigation, varying the frustration strength κ f reshapes the energy landscape and hence the set of equilibria. Under weak frustration, the fully synchronized state is the unique minimum, whereas increasing κ f stabilizes some low-index saddles into additional local minima and destabilizes the fully synchronized state into a saddle. Consequently, multiple sets of isolated and non-isolated low-energy minima, including synchronized and asynchronous states, may coexist over an intermediate range of κ f. In the spirit of synergetics, these minima can be represented by ordered macroscopic patterns described by a few collective modes, while the index characterizes the stability and instability of the corresponding order parameters. Noise then induces switching between ordered states, as reflected in the changes in the dominant Hadamard modes in our simulations. These results imply that the high-dimensional coupled oscillator dynamics is effectively governed by a small set of slow collective variables, in line with the so-called slaving principle.

At the same time, several theoretical and algorithmic questions concerning the homotopy method remain open as follows:

• Origin of selectivity. The proposed homotopy method uses a phase-free amplitude system as the starting system, while the phase initialization relies on the conventional Newton–Raphson solution of the leading XY Hamiltonian (in a generic case) or the constrained-solvability problem (in singular cases). Compared with a purely local Newton–Raphson search in the full ( r , θ ) space, the homotopy method exhibits a certain selectivity, which facilitates the search for low-index equilibria. Although we intuitively attribute this selectivity to the structural closeness between the approximate and full systems, achieved by incorporating the coupling effects into the eigenspace coordinates, this selectivity is at present largely empirical, and its theoretical origin and performance remain to be clarified more rigorously.
• Completeness of homotopy. In homotopy continuation, a homotopy is said to be complete if, after tracking all start solutions, every isolated solution of the target system is reached by at least one continuation path. We have not yet established whether the proposed homotopy is complete in this sense, namely, whether every isolated equilibrium of the target system, or at least every low-index one, is reachable from some equilibrium of the surrogate approximation. Moreover, when W + Diag ( a ) has repeated eigenvalues, and B ⊤ B is singular, the equilibria found by continuation depend on the choice of eigenbasis E. Although randomization of E may exhaust the set of reachable equilibria, it remains unknown whether some equilibria are still inaccessible under this procedure. These questions call for further theoretical study of the proposed homotopy.
• Computational complexity and scalability. A systematic analysis of the computational complexity and scalability of the present continuation framework remains critical for the broader applicability of the method. In particular, it would be important to clarify how the computational cost, numerical stability, and equilibrium coverage vary with lattice size, and whether changes in system size introduce additional failure modes in branch tracking or correction. It would also be of interest to investigate how these aspects depend on the symmetry structure of the system as symmetry may affect both the organization of equilibria and the behavior of the continuation procedure.

All of these issues constitute important directions for future research.

4.2 Energy landscape analysis based on coupled oscillator models

The idea of organizing dynamics in terms of an energy landscape has a long-standing history. In statistical physics and chemistry, the potential-energy landscape framework for viscous liquids and glasses, originating with Goldstein’s barrier picture (Goldstein, 1969) and developed into the “inherent-structure” formalism (Stillinger and Weber, 1982), interprets slow relaxation as thermally activated motion among local minima separated by energy barriers (Stillinger, 1995; Debenedetti and Stillinger, 2001). In neural computation, the Amari–Hopfield neural network (Amari, 1972; Hopfield, 1982) is endowed with an energy (Lyapunov) function, whose local minima are attractors encoding stored patterns, and basins of attraction implement content-addressable recall (Hopfield, 1982). Recently, its extensions with higher-order interactions have been applied to modern machine learning tasks (Demircigil et al., 2017; Krotov and Hopfield, 2018; Ramsauer et al., 2021). In developmental and stem-cell biology, Waddington’s epigenetic landscape (Waddington, 1957) is now interpreted in dynamical-systems terms, with cell fates as attractors and fate decisions as transitions between basins, in the state space of gene regulatory networks (Kauffman, 1993; Bhattacharya et al., 2011). This picture has been made quantitative, considering quasi-potential and flux descriptions of non-equilibrium gene regulatory dynamics (Wang et al., 2011; Zhou et al., 2012) and birth-death-based decompositions into steady-state and pluripotency landscapes (Shi et al., 2022), along with data-driven methods that infer attractors and branching differentiation trajectories directly from high-dimensional single-cell transcriptomic data (Guo and Zheng, 2017; Shi et al., 2019; Zhang et al., 2019). More broadly, there is a growing trend to model disease as abnormal energy-minimum attractors and develop early-warning indicators for impending shifts from healthy to pathological states (Huang et al., 2009; Fard and Ragan, 2017; Uthamacumaran, 2021; Chen et al., 2012; Aihara et al., 2022).

In contemporary neuroscience, a prominent energy-based framework for modeling brain-wide neural dynamics, especially intrinsic or spontaneous activity, is the energy landscape analysis (Watanabe et al., 2014; Ezaki et al., 2017) implemented on the basis of a kinetic Ising model of spins (Schneidman et al., 2006; Fraiman et al., 2009; Watanabe et al., 2013). Such an Ising model can be viewed as an Amari–Hopfield network with stochastic updates (Sompolinsky, 1988), where the occurrence of binary activity patterns follows a Gibbs (Boltzmann) distribution over the energy. Energy minima thus correspond to frequently visited brain states, and thermally driven barrier crossings generate transitions among them. This approach has been applied to binarized fMRI and other population-level data to identify repertoires of brain states and their transition structure and reveal possible altered state accessibility in psychiatric and neurological disorders, including autism (Watanabe and Rees, 2017), Alzheimer’s disease (Xing et al., 2024), aphasia (Watanabe et al., 2025), and others. It should be noted that energy landscape analysis is also applicable beyond neurodynamics, for example, to the modeling of disease (Fard and Ragan, 2017). Recently, an energy landscape analysis of health checkup data clarified multiple pathways to diabetes development in obese and non-obese subjects (Ito et al., 2025).

Our key motivation for revisiting the energy-based formulation of the Stuart–Landau oscillator network is to use it as a basis for a continuous energy landscape analysis, in contrast to the popular binary approach. This perspective becomes natural when oscillatory behavior is of primary interest, and it also rests on a substantial neuroscientific background. Computational whole-brain models based on coupled oscillators, including Stuart–Landau oscillators, naturally describe collective rhythms (full and partial synchronization), multistability, metastability, and criticality. They have been successfully applied to account for brain dynamics (Deco et al., 2017; Cabral et al., 2017; Deco and Kringelbach, 2020), quantify brain flexibility (Chinichian et al., 2024), and predict disease-related alterations (Demirtaş et al., 2017; Idesis et al., 2024). We speculate that, compared with Ising-based landscapes, oscillator-based energy landscapes offer several potential advantages:

• Continuous oscillatory states. Oscillator networks model oscillations in continuous state variables and are suited for faithful modeling of brain signals such as fMRI and electroencephalography (EEG), as well as other oscillatory time series. In particular, they allow stable system states to be described as limit-cycle attractors, rather than static binary patterns.
• Differential structure. With a continuous Euclidean state space, equilibria and their stability can be characterized by standard linear spectral methods, and pathways can be treated as curves in this space. This avoids combinatorial searches over the { ± 1 } N hypercube and the need to analyze high-dimensional transition operators or Boolean Jacobians over F 2.
• Explicit noise scale. Due to a lack of local dynamics of the node in the spin network, Ising-based approaches typically recover landscapes and couplings only relative to an overall noise level (temperature). In oscillator networks, the noise level is encoded in the diffusive wandering of trajectories and can be estimated explicitly, allowing the effects of couplings and temperature to be distinguished.

With these potential advantages, a practical workflow for energy landscape analysis based on coupled oscillator networks may parallel the now-familiar Ising-based framework (Masuda et al., 2025) while extending it to continuous oscillatory dynamics. In general, it may include the following steps:

• Multivariate time series from the nodes of a physiological network are preprocessed by standard procedures such as normalization, detrending, and filtering, depending on the recording modality and the frequency range of interest. When the number of channels is large, clustering, parcellation, or other dimension-reduction methods may be introduced to obtain a tractable set of effective variables.
• The processed data are then fitted with the energy-based Stuart–Landau oscillator model (Equation 1), for example, by Monte Carlo maximum likelihood estimation. In contrast to the Ising setting, this formulation allows not only the effective coupling matrix W but also the local parameters a and the noise intensity σ to be inferred from data.
• Using the inferred W and a, one constructs the corresponding energy landscape and investigates its low-energy equilibrium structures, particularly local minima, which represent macrostates, and low-index saddles, which may mediate transitions between macrostates. The homotopy method developed in the present study is intended to contribute to this stage of the analysis.
• The identified states may then be further analyzed in terms of their mutual distances (Equation 49), low-dimensional embeddings, and transition statistics, thereby providing a coarse-grained description of the system’s dynamical organization in a form suitable for comparison across conditions, subjects, and physiological states.

Therefore, the present study forms a preparatory step toward an oscillator-based energy landscape framework. By invoking more systematic methods for finding low-index equilibria and connecting paths (E and Zhou, 2011; Quapp and Bofill, 2014; Yin et al., 2020), a more comprehensive analysis of states and transition paths in the energy landscape of Stuart–Landau oscillator networks will become accessible, paving the way for applications to real-world systems and empirical data.

4.3 State transitions in energy-based models and beyond

Although a landscape with deterministic energy-descent dynamics conveniently models stable states in complex systems, it usually serves only as a diagnostic reference once we seek to describe transitions among these states because the ideal energy-descent structure is then broken. In addition to thermal activation, which yields a steady-state distribution over states, many mechanisms lead to state changes even at the deterministic level. Recently, there has been growing interest in modeling transitions among metastable states by a sequence of saddles connected by heteroclinic orbits in explicitly non-gradient systems (Afraimovich et al., 2004a; Afraimovich et al., 2004b; Rabinovich et al., 2008; Meyer-Ortmanns, 2023). Such a setting, in a more complicated form, leads to chaotic itinerancy (Tsuda, 2001; Kaneko and Tsuda, 2003; Tsuda, 2015; Inoue et al., 2020) and transitive chaotic synchronization (Kanamaru et al., 2023). Here, we can distinguish several mechanisms that induce transitions among states defined by an underlying landscape, including but not limited to the following:

• Direct perturbations of the state. Transitions can be caused by perturbations that act directly in the state space. One extreme is transient forcing by external inputs or pulses that push the system into another basin; the other extreme is persistent stochastic forcing that drives ongoing escapes, as in kinetic Ising models and stochastic gradient systems.
• Non-gradient flux-driven wandering. Here, the deterministic dynamics in the original state space is Markovian but no longer gradient: the underlying landscape still defines states, while an additional flux component introduces circulatory flow (Wang et al., 2011). Such a situation may arise from parameter inhomogeneities, asymmetric couplings (e.g., in sequential associative memory (Amari, 1972; Sompolinsky and Kanter, 1986)), or other mild non-gradient perturbations. In many models of this type, original minima become metastable waypoints that are organized by heteroclinic orbits, as mentioned above.
• Delay-induced transitions. When the vector field depends on time-delayed state variables, the deterministic dynamics is non-Markovian and, except in special cases, admits no finite-dimensional gradient–flux decomposition. Nonetheless, delay can break monotone descent in the underlying energy landscape and induce transitions between corresponding states (Kleinfeld, 1986; Janson and Marsden, 2021).
• Fast or slow landscape deformation, or tipping. In this case, the energy landscape itself changes in time with a drifting control parameter, externally prescribed or internally generated. If the parameter slowly crosses a bifurcation point, a minimum disappears or loses stability, and the system jumps to another state (bifurcation-induced tipping, B-tipping). If the parameter stays in a regime where the minimum exists but changes too fast for the system to track it, the trajectory can lose the branch and fall into a different basin (rate-induced tipping, R-tipping) (Kuehn, 2011; Ashwin et al., 2017). In particular, B-tipping is frequently adopted in modeling cell differentiation processes and disease progression (Bhattacharya et al., 2011; Aihara et al., 2022).
• Strongly non-gradient chaotic itinerancy and transitive chaotic synchronization. Here, the deterministic dynamics is strongly modified so that the original minima merely correspond to preferred regions (attractor ruins) in the state space, embedded in a more intricate heteroclinic skeleton than in the second case above. The resulting dynamics exhibits chaotic itinerancy and transitive chaotic synchronization: stationary epochs near one quasi-stable attractor ruin are separated by irregular transitions to others. Representative examples include the chaotic associative memory model (Adachi and Aihara, 1997; Kitajima et al., 2003), in which chaotic neurons replace spins in the Amari–Hopfield network and coupled Stuart–Landau oscillators of a more general form (Uchiyama and Fujisaka, 2004), where its equivalent Ginzburg–Landau map shows strongly nonrelaxational dynamics.

The above mechanisms can act separately or in combination in shaping state transitions. Although these mechanisms have been well studied, our interest lies in modeling and predicting real-world systems while explicitly accounting for them. Developing and generalizing oscillator-based energy landscape analysis to incorporate them may lead to more realistic and mechanistic interpretations of complex systems and remains an important direction for future work.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2026.1778380/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Eckehard Schöll, Technical University of Berlin, Germany

Reviewed by: Mohamed Ridha Znaidi, University of Southern California, United States

Bo-Wei Qin, Fudan University, China,