Finding equilibria of the finite size Kuramoto model amounts to solving a nonlinear system of equations, which is an important yet challenging problem. We translate this into an algebraic geometry problem and use numerical methods to find all of the equilibria for various choices of coupling constants K, natural frequencies, and on different graphs. We note that for even modest sizes (N ∼ 10–20), the number of equilibria is already more than 100 000. We analyze the stability of each computed equilibrium as well as the configuration of angles. Our exploration of the equilibrium landscape leads to unexpected and possibly surprising results including non-monotonicity in the number of equilibria, a predictable pattern in the indices of equilibria, counter-examples to conjectures, multi-stable equilibrium landscapes, scenarios with only unstable equilibria, and multiple distinct extrema in the stable equilibrium distribution as a function of the number of cycles in the graph.

The Kuramoto model is a fascinating model proposed in 1975 to study synchronization phenomena1 that have gained attention due to its applicability from various scientific communities, including biology, chemistry, physics, and electrical engineering. This model has been used to study various phenomena including neural networks, chemical oscillators, Josephson junctions and laser arrays, power grids, particle coordination, spin glass models, and rhythmic applause.2–4 The objective of this article is to introduce an algebraic geometric interpretation of the problem of computing equilibria for the Kuramoto model. Our interpretation requires solving a system of polynomial equations, and we discuss techniques from numerical algebraic geometry that are well suited for these problems. Computing equilibria at various points in parameter space and analyzing their stability can provide new insight in synchronization phenomena.

The Kuramoto model is defined as a system of autonomous ordinary differential equations as

dθidt=ωiKNj=1Nai,jsin(θiθj)fori=1,,N,
(1)

where K is the coupling strength, N is the number of oscillators, Ω = (ω1,…,ωN) is the vector of intrinsic natural frequencies, and ai,j = aj,i ∈ {0, 1} is the (i, j) entry of the adjacency matrix of the coupling graph which we assume is undirected. The natural frequencies ωi indicate how the system oscillates in the absence of any dissipation or exogenous forces. Without loss of generality, we assume i=1Nωi=0. In particular, if ω¯:=i=1Nωi/N, then the transformation to a rotating frame, θiθiω¯t leads to the rescaled natural frequencies ω̃i=ωiω¯ satisfying i=1Nω̃i=0.

The equilibrium conditions are dθidt=0 for all i. This system of equations has an O(2) freedom, i.e., for any α ∈ (−π, π], the equations are invariant under replacing all θi with θi + α. This rotational symmetry leads to the continua of equilibria. To remove this O(2) freedom resulting in finitely many equilibria, we fix one of the angles, say, θN = 0, and remove the equation dθNdt=0 from the system. The remaining system consists of N − 1 nonlinear equations in N − 1 angles. Since i=1Ndθidt=i=1Nωi=0, no information is lost by removing one equation from the system. See Appendix for an equivalent formulation removing the O(2) freedom.

Provided the coupling strength K is strong enough, the oscillators will synchronize as t. In this setup, a critical coupling Kc(N) exists at which the number of stable equilibria switches from 0 to a nonzero value. In the special case of N and long-range (all-to-all) coupling ai,j = 1, one may analytically compute Kc(N). However, for the finite size Kuramoto model, such an analysis may turn out to be very difficult. In particular, finding all equilibria, analyzing stability, and finding Kc(N) are known to be prohibitively difficult for a finite but large oscillator population.

Since we are working with undirected coupling graphs, i.e., ai,j = aj,i for all i, j, we point out that the equilibria of system (1) can also be viewed as the stationary points of the potential energy landscape drawn by the mean-field XY model with an exogenous perturbation term

V(θ)=K2Ni,j=1Nai,j(1cos(θiθj))i=1Nωiθi,
(2)

whose gradient reproduces the right-hand side of Eq. (1). Hence, in the following, we use the words equilibria and stationary points interchangeably.

All stationary points of the finite N mean-field XY model (i.e., the Kuramoto model with homogeneous frequencies) were identified in Ref. 5. Building on Ref. 5, all stationary points of the one-dimensional nearest-neighbour XY model (i.e., the Kuramoto model with local coupling) for any given N with either periodic or anti-periodic boundary conditions have been found.6–9 

Using these solutions, a class of stationary points of the two-dimensional nearest-neighbour XY model10 and the XY model with long-range interactions11 were built and analyzed (see also Refs. 12 and 13). In Refs. 6, 7, 12, and 14, all of the stationary points for small lattices were found using algebraic geometry methods. Bounds on the number of equilibria15 as well as some counterintuitive examples to plausible conjectures16 have been reported for the same model in the domain of power systems.

In Ref. 5 for the finite N mean-field XY model, and in Refs. 6 and 10 for the nearest neighbour XY models, it was shown that there were exponentially many isolated stationary solutions as N increases. Moreover, even after breaking the global O(2) symmetry, there were continuous solutions at the maximum value of the energy which were independently observed and termed as incoherent manifolds in Ref. 17.

In Refs. 18–20, necessary and sufficient conditions are given for fixed points to exist for the finite size Kuramoto model for complete and bipartite graphs, and explicit upper and lower bounds of Kc(N) for these systems were also computed, followed by providing an algorithm to compute Kc(N). For the complete graph, similar results were presented in Refs. 21 and 22, where it was additionally shown that there is exactly one single stable equilibrium for K > Kc(N). In Ref. 23, an analogous result was shown for acyclic graphs, short cycles, and complete graphs as well as combinations thereof. In the case of homogeneous natural frequencies networks with sufficiently high nodal degrees, the only stable fixed point is known to be the phase-synchronized solution.24 In Ref. 25, all of the stable synchronized states were classified for the one-dimensional Kuramoto model on a ring graph with random natural frequencies in addition to computing a lower bound on Kc(N) (see also Refs. 26, 44, and 45).

We initiate an approach to study the Kuramoto model by using an algebraic geometry interpretation of the equilibria and studying synchronization. As mentioned above, upon assuming, without loss of generality, that ωN=i=1N1ωi, fixing θN = 0, and removing the Nth equation, we have

0=ωiKNj=1Nai,jsin(θiθj)fori=1,,N1.
(3)

We use the identity sin(xy)=sinxcosysinycosx and substitute si=sinθi and ci=cosθi to transform the N − 1 equations into polynomials that are coupled by the Pythagorean identity si2+ci2=1. This results in a system of 2(N − 1) polynomials in 2(N − 1) variables

0=ωiKNj=1Nai,j(sicjsjci)0=si2+ci21.
(4)

Given a fixed integer N, our first goal is to find all real solutions of polynomial system (4) for a given choice of Ω, K, and adjacency matrix A = [ai,j]. We accomplish this using a technique from numerical algebraic geometry via parameter homotopies, which we briefly review below. Once the system is solved, we determine which solutions (if any) are stable steady-state solutions by analyzing the eigenvalues of the Jacobian at each real solution.

For a fixed N, we can interpret (4) as the system FN(s, c; K, Ω, A) = 0 with variables s and c and parameters K, Ω, and A. This leads to using parameter homotopies27 for solving, which is a two-step process. First, in the ab initio phase, one computes all solutions for a sufficiently random set of parameters. Then, in the parameter homotopy phase, one solves for given parameters by deforming from the parameters selected in the ab initio phase. The ab initio phase is performed once, while the parameter homotopy phase is performed for each set of parameters. The following is a short introduction with more details provided in Ref. 28, Chap. 7 and Ref. 29, Chap. 6.

This process is a generalization of other standard homotopy methods, such as a total degree homotopy or multi-homogeneous homotopy frequently used, e.g., Refs. 28, 29, and 32–40. The method has also recently used in solving the power-flow systems.30,31 For example, consider solving a system of m polynomial equations in m variables, f(x)=(f1(x),,fm(x))=0 where x=(x1,,xm). For a total degree homotopy, the space of interest is the space of all polynomial systems g(x) of m polynomials in m variables where di:=degfi=deggi. Since a random element in this space has the Bézout number of roots, namely, i=1mdi, the system G=(G1,,Gm)=0, where Gi(x)=xidi1 is sufficiently random in this space.

In this case, the ab initio phase for solving G = 0 is trivial. Then, in the parameter homotopy phase, one can consider the homotopy H(x,t)=(1t)f(x)+γtG(x) for a random γ. Now, the solutions of H(x, t) = 0 at t = 1 are known and we want to compute the solutions of H(x, 0) = 0. Using a numerical predictor-corrector method implemented in Bertini,41 we track each of the solution paths from t = 1 to t = 0. The number γ ensures that such path tracking will obtain all of the isolated complex solutions of f(x) = 0, from which the real and nonreal solutions can be identified. For nonsingular solutions, this identification can be done certifiably.42 

Now, returning to our case FN, suppose, for simplicity, that we want to study the behavior of the solutions of FN = 0 as a function of K, that is, we also fix Ω and A. In the ab initio phase, we pick a random K=K0 and compute the solutions S0 of FN(s,c;K0,Ω,A)=0 using a total degree homotopy described above via Bertini.

In the parameter homotopy phase, we consider solving for various choices of parameters K, e.g., K = 2–100. For each choice of K, we simply track the solution paths of FN(s,c;K·(1t)+t·K0,Ω,A)=0 starting at t = 1 with solutions S0. Each of these computations is relatively inexpensive, thereby making it practical to compute the solutions at hundreds of different choices for K.

The resulting solutions are sorted based on whether they are real or not, and their stability is investigated as follows. From the sine si and cosine ci values, we compute the corresponding angle θi. The eigenvalues of the Jacobian of (1) are used to analyze stability. The index is the number of negative eigenvalues so that real solutions with index zero correspond to stable steady states.

We start with the most prominent and well studied case corresponding to the undirected complete graph, namely, ai,j = 1. We select Ω = (ω1ωN) to be N equidistant numbers, namely, ωi=1+(2i1)/N. Figures 1 and 2 summarize the number of real-valued solutions to the polynomial system (4) for N = 3–18 at various values of K. In general, we see that the number of distinct real solutions tends to grow as K increases. However, there are some surprising exceptions to this expected monotonicity property, e.g., the N = 9 case has 328 real solutions at K = 60 and 326 real solutions at K = 70. Our results suggest that eventually the number of real solutions stabilizes to a constant for each N. For N = 18, we computed 140 356 real solutions when K = 50; computing these without our polynomial formulation would be very difficult.

FIG. 1.

Number of equilibria for the case of equidistant natural frequencies on the complete graph at 1 ≤ K ≤ 2.

FIG. 1.

Number of equilibria for the case of equidistant natural frequencies on the complete graph at 1 ≤ K ≤ 2.

Close modal
FIG. 2.

Number of equilibria for the case of equidistant natural frequencies on the complete graph at 2 ≤ K ≤ 100.

FIG. 2.

Number of equilibria for the case of equidistant natural frequencies on the complete graph at 2 ≤ K ≤ 100.

Close modal

Once the real solutions are found, we turn our attention to determining which ones are stable. By analyzing the Jacobian, we find that there is exactly a single stable steady state solution for N = 3–18 at each K value tested, with the exception of small values of K which result in no real solutions. This result is consistent with theoretic findings for the complete graph.19–22 Next, for each N, we use the results of our stability analysis to determine numerical upper and lower bounds on Kc(N) as follows. We acquire an upper bound on Kc(N) by determining the smallest tested value of K for which a real steady state is found. Similarly, we find a lower bound on Kc(N) by determining the largest tested value of K for which no real steady states are found. These results are presented in Figure 3. We verify that these bounds for Kc(N) are consistent with the known explicit bounds from Ref. 4, Corollary 6.7. For N ∈ {3, 4}, our computed lower bound is less than the explicit lower bound due to the coarse resolution of tested numeric values for Kc(N).

FIG. 3.

Bounds for Kc(N) in the case of equidistant natural frequencies on the complete graph. Known explicit bounds from Refs. 4, Corollary 6.7 are shown for comparison.

FIG. 3.

Bounds for Kc(N) in the case of equidistant natural frequencies on the complete graph. Known explicit bounds from Refs. 4, Corollary 6.7 are shown for comparison.

Close modal

In Ref. 16, it is conjectured that if there is a stable equilibrium, then there is a stable equilibrium satisfying |θiθj|<π/2 for all neighbors {i, j}. We checked this conjecture and we quickly found counterexamples. For N = 3 and K = 1.15, there is one stable equilibria at (θ1,θ2,θ3)(1.6820,0.8410,0). Since θ1 and θ3 are coupled with |θ1θ3|>π/2, we reject the conjecture.

For larger N, the unique stable equilibria can have angles even more spread out, such as (θ1,,θ15)(2.0981,1.8867,,0.2114,0) when N = 15 and K = 1.3. These interesting numerical observations prompt us to analytically reject the conjecture in Ref. 16; in Ref. 4, it is shown that the worst-case equilibrium configuration for K = Kc and n > 3 corresponds to a tripolar distribution of all angles spanning a half-circle. From Refs. 21 and 22, we know that this is the only stable equilibrium. These two arguments suffice to reject the conjecture in Ref. 16.

As mentioned earlier, our stability analysis is based on computing the index corresponding to each real solution. Figure 4 shows a histogram of these values for the case K = 100. We observe that when j is small enough relative to N, the number of real solutions with index j is exactly (Nj). In particular, for 3 ≤ N ≤ 4 this behavior occurs for j ≤ 1 and, for 5 ≤ N ≤ 18, this behavior occurs for j 25(N1) . Based on these results, we conjecture that this phenomenon occurs in general: for equidistant natural frequencies, if 0 ≤ jN and K ≫ 0, then the number of equilibria with index j is expected to be (Nj). This conjecture can be used when N may be too large to compute all solutions. For example, we expect 75 287 520 real equilibria of index 5 for N = 100.

FIG. 4.

Number of equilibria with given index for equidistant natural frequencies and a complete graph with K = 100.

FIG. 4.

Number of equilibria with given index for equidistant natural frequencies and a complete graph with K = 100.

Close modal

Although our investigation initially focused on the complete graph being the most well studied case, we also performed a preliminary analysis of the coupling arrangement defined by an undirected cyclic graph. Cyclic graphs are known for having multi-stable equilibria,4 entirely unstable equilibrium landscapes16 which is distinct from acyclic or complete graphs.23,43,45

For N = 10, we used the same equidistant natural frequencies mentioned earlier with Figure 5 showing the number of equilibria and the number of stable equilibria at each integer K = 0–100.

FIG. 5.

Number of equilibria and number of stable equilibria for a cyclic coupling arrangement when N = 10.

FIG. 5.

Number of equilibria and number of stable equilibria for a cyclic coupling arrangement when N = 10.

Close modal

For some values of K, the system possesses only unstable equilibria. In particular, when K is 13, 14, or 15, there are 76, 164, and 260 equilibria, respectively, all of which are unstable. Since the critical points of the potential (2) correspond to power flow in a transmission network, a very interesting technological implication of these results is that there are power demands which can be met only with unstable equilibria. We find one stable equilibria at K = 16 and conclude 15 ≤ Kc(10) ≤ 16. Figure 5 also indicates multistability for some values of K, with at most three stable equilibria for each investigated sample. These results are in contrast to the complete graph case, in which exactly one stable equilibrium exists whenever the system has real-valued solutions.

Next, we take a closer look at the geometric configuration of angles occurring at the stable equilibria shown in Figure 5. For relatively small K such as K = 16, 17, and 18, the configuration of angles at stable equilibria shows no discernible structure. For K = 35 there are two stable equilibria, with one exhibiting phase sync, i.e., angles clustered around 0, and one exhibiting a splay state, i.e., angles approximately uniformly distributed on [0, 2π). For the splay state, the θi decrease on [0, 2π).

For each K = 36–100, there are three stable equilibria, with each case consisting of one phase sync and two splay states. In one of these splay states θ1,…,θ9 are arranged in increasing order on [0, 2π), and in the other θ9,…,θ1 are arranged in decreasing order on [0, 2π). In Figure 6, we depict the three stable equilibria at K = 100. As K increases, the steady state phase sync gradually becomes more tightly clustered, with the angle range decreasing from ∼0.9432 for K = 35 to ∼0.3256 for K = 100. This is in accordance with the asymptotic result that exact phase sync is a critical point of the Kuramoto potential (2) as K.4 

FIG. 6.

Configuration of θ1,…,θ10 for the three stable equilibria occurring at K = 100 for a cyclic graph with N = 10.

FIG. 6.

Configuration of θ1,…,θ10 for the three stable equilibria occurring at K = 100 for a cyclic graph with N = 10.

Close modal

We now turn our attention to random coupling arrangements. In these computations, we choose each ai,j as 0 or 1 according to a pre-determined probability P while setting aj,i = ai,j. In other words, [ai,j] is the adjacency matrix of a symmetric Erdös-Rényi random graph with coupling probability P, where we restrict our attention to undirected, connected graphs. For these numerical experiments, we use the same equidistant natural frequencies discussed earlier.

First, we fix N = 8 and K = 100 and investigate the number of stable equilibria for sparse random graphs. For each c = 0, 1, …, 20, we construct 100 random graphs having exactly c cycles. We achieve this by generating adjacency matrices of symmetric Erdős-Rényi random graphs until we have 100 instances of connected graphs with the desired number of cycles. Depending on the graph, we find 1, 2, or 3 stable equilibria. The averages of these results are shown in Figure 7 showing that the sparsest (i.e., acyclic) as well as sufficiently dense graphs have exactly one stable equilibrium confirming the analytic results for these two extremal topologies.4,21–24 Since cycles are known to display multi-stability (see Figure 5), it is to be expected that there is more than one stable equilibrium for small cycle numbers. Quite surprisingly, multiple distinct extrema can be observed in the stable equilibrium distribution in Figure 7.

FIG. 7.

Number of stable equilibria for random graphs according to number of cycles in the case of N = 8 and K = 100.

FIG. 7.

Number of stable equilibria for random graphs according to number of cycles in the case of N = 8 and K = 100.

Close modal

Next, we fix N = 10 and investigate how Kc(N) depends on the density of the graph. For each coupling probability P = 0.25, 0.375, …, 0.875, we generated 100 random graphs. For each graph, we computed equilibria and determine stability at the following values of K: for P = 0.25 we use K = 1–20; for P = 0.375 we use K = 1–15; for P = 0.5 we use K = 1–10; and for each P = 0.625, 0.75, and 0.875, we use K = 1–5. In all cases, we find a value of K such that at least one stable equilibrium occurs at K while no stable equilibria occur at K – 1, thereby allowing us to estimate bounds on Kc(10). Figure 8 shows these results sorted according to the number of cycles in the graph.

FIG. 8.

Average bounds on Kc(N) for random graphs according to number of cycles in the case of N = 10.

FIG. 8.

Average bounds on Kc(N) for random graphs according to number of cycles in the case of N = 10.

Close modal

For undirected and connected graphs, a theoretic lower and a conjectured theoretic upper bound are4,23

N·maxi{ |ωi|degi }Kc(N)N· BTLω ,

where degi=jai,j is the degree of node i, B is the oriented incidence matrix, and L is the network Laplacian matrix. These bounds are shown in Figure 8 for comparison. We verified that our 600 individual results as well as the averages shown are consistent with these theoretic bounds and validate their accuracy. Figure 8 indicates that the numerical upper bound is tighter than the theoretic bound on average, while the numerical lower bound tends to be weaker. Since the tightness of our computed bounds depends on the resolution of K values tested, one could check more refined K values to obtain tighter numerical bounds for particular graphs of interest.

We have interpreted the problem of finding all of the equilibria of the Kuramoto model as an algebraic geometry problem by translating the stationary equations into a system of polynomial equations. The coupling constant, natural frequencies, and entries of the adjacency matrix of the corresponding graph can all be viewed as parameters of the system. With this unification, we show that polynomial homotopy continuation combined with parameter homotopies cannot only find all isolated equilibria at each parameter point, but can also scan over a large set of parameters without having to re-solve the system at each instance.

For complete graphs, we find unexpected non-monotonicity in the number of equilibria as K increases. Our stability analysis leads us to a conjecture regarding the expected number of equilibria of certain index. Our investigation of angle configuration at steady states leads us to reject a conjecture from Ref. 16. For cyclic graphs, we note certain coupling constants K result in only unstable equilibria, and we show multi-stable scenarios with a multi-modal stable equilibrium distribution depending on the number of cycles. Our experiments on random graphs shed insight into the Kuramoto model in a broader sense and show the importance of the coupling probability. All of our results are consistent with the available theoretic bounds which are quite loose at times.

Perhaps more importantly for complex systems areas, with the algebraic geometry interpretation and computational methods, we have developed a reliable and direct way of bounding the critical coupling Kc(N) for synchronization in the finite size Kuramoto model that does not rely on any specific frequency distribution or graph topology. We are in the process of writing a more specific algorithm to compute Kc(N) up to an arbitrary precision based on the algebraic geometry interpretation and parameter homotopy computations.

By using this algebraic geometric interpretation, one can also study situations where there is an infinite number of critical points, e.g., the homogeneous case when ωi = 0 for all i. The resulting polynomial system will have positive-dimensional solution components which can be computed via numerical algebraic geometry, e.g., see Ref. 29, Chap. 8. Analysis of such components using an algebraic geometric viewpoint is beyond the current scope of this paper but could be addressed in the future.

Additionally, we aim to extend our investigations to directed graphs, negative weights, cosine coupling, and statistical properties of random graph models, as well as verify conjectures for variations of the Kuramoto model as studied in Ref. 43.

D.M. would like to thank Carlo Laing and Steven Strogatz for their feedback at the initial stages of this work. N.S.D., J.D.H., and D.M. were supported by DARPA Young Faculty Award and Sloan Research Fellowship with N.S.D. and J.D.H. additionally supported by NSF DMS-1262428. F.D. was supported by in part by ETH Zürich startup funds.

In Sec. I, we remarked that the equilibrium conditions dθidt=0 form a system of N equations with O(2) freedom. That is, for any α ∈ (−π,π], the equilibrium conditions are invariant under the rotation θiθi + α for all i. We addressed this O(2) freedom by fixing θN = 0 and removing the Nth equation from the system. We claimed no information was lost in this simplification since iωi=0. In this Appendix, we follow different steps to remove the O(2) freedom and show that both formulations are indeed equivalent.

For any undirected coupling graph A = [ai,j], (1) implies that idθidt=iωi. In particular, equilibrium, i.e., dθidt=0 for i = 1,…,N, can occur only if iωi=0.

We proceed to solve dθidt=0 for i = 1,…,N by considering ϕi,j=θiθj for each i, j = 1,…,N. Since the system has an O(2) freedom, it is natural to work with differences of phases rather than the phases themselves.

With this setup, we have the system

ϕi,j=θiθjforeachi,j=1,,N0=ωiKNj=1Nai,jsinϕi,jfori=1,,N.

As written, both the linear and trigonometric conditions include redundant information. For example, ϕi,i=0,sinϕi,i=0, and ϕi,j=ϕj,i. Thus, we can immediately restrict to variables ϕi,j with 1 ≤ i < j ≤ N. Upon observing that ϕi,j=ϕi,Nϕj,N, we are left with N − 1 variables: ϕi,N for i = 1,…,N − 1.

For i = 1,…,N − 1, this setup yields the following:

0=ωiKNj=1Nai,jsinϕi,j=ωiKN((j=1N1ai,jsin(ϕi,Nϕj,N))+ai,Nsinϕi,N),
(A1)

with the Nth equation being

0=ωNKNj=1NaN,jsinϕN,j=ωN+KNj=1NaN,jsinϕj,N.
(A2)

Since solving f1==fN1=fN=0 is equivalent to solving f1==fN1=i=1Nfi=0, we can replace (A2) by the corresponding sum which simplifies to

0=i=1Nωi
(A3)

using the assumption that the graph is undirected, i.e., ai,j = aj,i, and the anti-symmetry property of the sine function, i.e., sin(x)=sin(x). Since (A3) is satisfied by assumption (or else there are no equilibria), this analysis yields a system of N − 1 equations presented in (A1) in N − 1 variables ϕi,N for i = 1,…,N − 1.

Now, we simply define θ̃i=ϕi,N for i = 1,…,N − 1 and θ̃N=0. Since sin(ϕi,N)=sin(θ̃iθ̃N), (A1) reduces to (3) upon relabeling.

In particular, we have shown that, in our case, setting one angle to 0 is equivalent to working with the difference of phases as methods for handling the O(2) symmetry.

1.
Y.
Kuramoto
,
International Symposium on Mathematical Problems in Theoretical Physics
(
Springer
,
1975
), pp.
420
422
.
2.
J. A.
Acebrón
,
L. L.
Bonilla
,
C. J.
Pérez Vicente
,
F.
Ritort
, and
R.
Spigler
,
Rev. Mod. Phys.
77
(
1
),
137
(
2005
).
3.
S. H.
Strogatz
,
Physica D
143
(
1
),
1
20
(
2000
).
4.
F.
Dörfler
and
F.
Bullo
,
Automatica
50
(
6
),
1539
1564
(
2014
).
5.
L.
Casetti
,
M.
Pettini
, and
E. G. D.
Cohen
,
J. Stat. Phys.
111
,
1091
1123
(
2003
).
6.
D.
Mehta
, “
Lattice vs. continuum: Landau gauge fixing and 't Hooft-Polyakov monopoles
,” Ph.D. thesis,
The University of Adelaide, Australasian Digital Theses Program
,
2009
.
7.
D.
Mehta
and
M.
Kastner
,
Ann. Phys.
326
,
1425
1440
(
2011
).
8.
L.
von Smekal
,
D.
Mehta
,
A.
Sternbeck
, and
A. G.
Williams
, in
PoS LAT
2007
, p.
382
.
9.
L.
von Smekal
,
A.
Jorkowski
,
D.
Mehta
, and
A.
Sternbeck
. in
PoS CONFINEMENT
,
2008
, Vol.
8
, p.
048
.
10.
R.
Nerattini
,
M.
Kastner
,
D.
Mehta
, and
L.
Casetti
,
Phys. Rev. E
87
(
3
),
032140
(
2013
).
11.
M.
Kastner
,
Phys. Rev. E
83
(
3
),
031114
(
2011
).
12.
C.
Hughes
,
D.
Mehta
, and
J.-I.
Skullerud
,
Ann. Phys.
331
,
188
215
(
2013
).
13.
D.
Mehta
and
M.
Schröck
,
Phys. Rev. D
89
,
094512
(
2014
).
14.
D.
Mehta
,
A.
Sternbeck
,
L.
von Smekal
, and
A. G.
Williams
, in
PoS, QCD-TNT09
,
2009
, p.
025
.
15.
J.
Baillieul
and
C. I.
Byrnes
,
IEEE Trans. Circuits Syst.
29
(
11
),
724
737
(
1982
).
16.
A.
Araposthatis
,
S.
Sastry
, and
P.
Varaiya
,
Int. J. Electr. Power Energy Syst.
3
(
3
),
115
126
(
1981
).
17.
S. H.
Strogatz
and
R. E.
Mirollo
,
J. Stat. Phys.
63
(
3–4
),
613
635
(
1991
).
18.
M.
Verwoerd
and
O.
Mason
,
SIAM J. Appl. Dyn. Syst.
8
(
1
),
417
453
(
2009
).
19.
M.
Verwoerd
and
O.
Mason
,
SIAM J. Appl. Dyn. Syst.
7
(
1
),
134
160
(
2008
).
20.
F.
Dörfler
and
F.
Bullo
,
SIAM J. Appl. Dyn. Syst.
10
(
3
),
1070
1099
(
2011
).
21.
D.
Aeyels
and
J. A.
Rogge
,
Prog. Theor. Phys.
112
(
6
),
921
942
(
2004
).
22.
R. E.
Mirollo
and
S. H.
Strogatz
,
Physica D
205
(
1–4
),
249
266
(
2005
).
23.
F.
Dörfler
,
M.
Chertkov
, and
F.
Bullo
,
Proc. Natl. Acad. Sci. U.S.A.
110
(
6
),
2005
2010
(
2013
).
24.
R.
Taylor
,
J. Phys. A: Math. Theor.
45
(
5
),
055102
(
2012
).
25.
J.
Ochab
and
P. F.
Góra
, “
Synchronization of coupled oscillators in a local one-dimensional Kuramoto model
,” preprint arXiv:0909.0043 (
2009
).
26.
P. F. C.
Tilles
,
F. F.
Ferreira
, and
H. A.
Cerdeira
,
Phys. Rev. E
83
,
066206
(
2011
).
27.
A. P.
Morgan
and
A. J.
Sommese
,
Appl. Math. Comput.
29
(
2
),
123
160
(
1989
).
28.
A. J.
Sommese
and
C. W.
Wampler
,
The Numerical Solution of Systems of Polynomials Arising in Engineering and Science
(
World Scientific Publishing
,
Hackensack, NJ
,
2005
).
29.
D. J.
Bates
,
J. D.
Hauenstein
,
A. J.
Sommese
, and
C. W.
Wampler
,
Numerically Solving Polynomial Systems with Bertini
(
SIAM
,
2013
), Vol.
25
.
30.
D.
Mehta
,
H.
Nguyen
, and
K.
Turitsyn
, “
Numerical polynomial homotopy continuation method to locate all the power flow solutions
,” preprint arXiv:1408.2732 (
2014
).
31.
S.
Chandra
,
D.
Mehta
, and
A.
Chakrabortty
, “
Exploring the impact of wind penetration on power system equilibrium using a numerical continuation approach
,” preprint arXiv:1409.7844 (
2014
).
32.
D.
Mehta
,
Adv. High Energy Phys.
2011
,
263937
.
33.
M.
Maniatis
and
D.
Mehta
,
Eur. Phys. J. Plus
127
,
91
(
2012
).
34.
M.
Kastner
and
D.
Mehta
,
Phys. Rev. Lett.
107
,
160602
(
2011
).
35.
D.
Mehta
,
Y.-H.
He
, and
J. D.
Hauenstein
,
JHEP
2012
(
18
),
1207
(
2012
).
36.
D.
Mehta
,
J. D.
Hauenstein
, and
M.
Kastner
,
Phys. Rev. E
85
,
061103
(
2012
).
37.
B.
Greene
,
D.
Kagan
,
A.
Masoumi
,
D.
Mehta
,
E. J.
Weinberg
, and
X.
Xiao
,
Phys. Rev. D
88
(
2
),
026005
(
2013
).
38.
D.
Mehta
,
D. A.
Stariolo
, and
M.
Kastner
,
Phys. Rev. E
87
(
5
),
052143
(
2013
).
39.
D.
Martinez-Pedrera
,
D.
Mehta
,
M.
Rummel
, and
A.
Westphal
,
JHEP
2013
(
110
),
1306
(
2013
).
40.
Y.-H.
He
,
D.
Mehta
,
M.
Niemerg
,
M.
Rummel
, and
A.
Valeanu
,
JHEP
2013
(
50
),
1307
(
2013
).
41.
D. J.
Bates
,
J. D.
Hauenstein
,
A. J.
Sommese
, and
C. W.
Wampler
, “
Bertini: Software for numerical algebraic geometry
,” available at http://bertini.nd.edu.
42.
J. D.
Hauenstein
and
F.
Sottile
,
ACM Trans. Math. Software
38
(
4
),
28
(
2012
).
43.
D. A.
Wiley
,
S. H.
Strogatz
, and
M.
Girvan
,
Chaos
16
,
015103
(
2006
).
44.
S. H.
Strogatz
and
R. E.
Mirollo
,
Physica D
31
,
143
(
1988
).
45.
G. B.
Ermentrout
,
J. Math. Biol.
23
,
55
(
1985
).