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 phenomena^{1} 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.

## I. INTRODUCTION

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

where *K* is the coupling strength, *N* is the number of oscillators, Ω = (*ω*_{1},…,*ω _{N}*) is the vector of intrinsic natural frequencies, and

*a*

_{i}_{,}

_{j}=

*a*

_{j}_{,}

_{i}∈ {0, 1} is the (

*i*,

*j*) entry of the adjacency matrix of the coupling graph which we assume is undirected. The natural frequencies

*ω*indicate how the system oscillates in the absence of any dissipation or exogenous forces. Without loss of generality, we assume $\u2211i=1N\omega i=0$. In particular, if $\omega \xaf:=\u2211i=1N\omega i/N$, then the transformation to a rotating frame, $\theta i\u21a6\theta i\u2212\omega \xaft$ leads to the rescaled natural frequencies $\omega \u0303i=\omega i\u2212\omega \xaf$ satisfying $\u2211i=1N\omega \u0303i=0$.

_{i}The equilibrium conditions are $d\theta 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,

*θ*= 0, and remove the equation $d\theta Ndt=0$ from the system. The remaining system consists of

_{N}*N*− 1 nonlinear equations in

*N*− 1 angles. Since $\u2211i=1Nd\theta idt=\u2211i=1N\omega 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 *K _{c}*(

*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

*a*

_{i}_{,}

_{j}= 1, one may analytically compute

*K*(

_{c}*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

*K*(

_{c}*N*) are known to be prohibitively difficult for a finite but large oscillator population.

Since we are working with undirected coupling graphs, i.e., *a _{i}*

_{,}

_{j}=

*a*

_{j}_{,}

_{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

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* model^{10} and the *XY* model with long-range interactions^{11} 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 equilibria^{15} as well as some counterintuitive examples to plausible conjectures^{16} 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 *K _{c}*(

*N*) for these systems were also computed, followed by providing an algorithm to compute

*K*(

_{c}*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*>

*K*(

_{c}*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

*K*(

_{c}*N*) (see also Refs. 26, 44, and 45).

## II. ALGEBRAIC GEOMETRY INTERPRETATION AND SETUP

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 $\omega N=\u2212\u2211i=1N\u22121\omega i$, fixing *θ _{N}* = 0, and removing the

*N*th equation, we have

We use the identity $sin(x\u2212y)=sin\u2009x\u2009cos\u2009y\u2212sin\u2009y\u2009cos\u2009x$ and substitute $si=sin\u2009\theta i$ and $ci=cos\u2009\theta 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

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* = [*a _{i}*

_{,}

_{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.

## III. NUMERICAL ALGEBRAIC GEOMETRY METHODS

For a fixed *N*, we can interpret (4) as the system *F _{N}*(

*s*,

*c*;

*K*, Ω,

*A*) = 0 with variables

*s*and

*c*and parameters

*K*, Ω, and

*A*. This leads to using parameter homotopies

^{27}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),\u2026,fm(x))=0$ where $x=(x1,\u2026,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, $\u220fi=1mdi$, the system $G=(G1,\u2026,Gm)=0$, where $Gi(x)=xidi\u22121$ 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)=(1\u2212t)f(x)+\gamma tG(x)$ for a random $\gamma \u2208\u2102$. 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 *F _{N}*, suppose, for simplicity, that we want to study the behavior of the solutions of

*F*= 0 as a function of

_{N}*K*, that is, we also fix Ω and

*A*. In the

*ab initio*phase, we pick a random $K=K0\u2208\u2102$ and compute the solutions

*S*

_{0}of $FN(s,c;K0,\Omega ,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\xb7(1\u2212t)+t\xb7K0,\Omega ,A)=0$ starting at *t* = 1 with solutions *S*_{0}. 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 *s _{i}* and cosine

*c*values, we compute the corresponding angle

_{i}*θ*. The eigenvalues of the Jacobian of (1) are used to analyze stability. The

_{i}*index*is the number of negative eigenvalues so that real solutions with index zero correspond to stable steady states.

## IV. RESULTS FOR COMPLETE GRAPHS

We start with the most prominent and well studied case corresponding to the undirected complete graph, namely, *a _{i}*

_{,}

_{j}= 1. We select Ω = (

*ω*

_{1}…

*ω*) to be

_{N}*N*equidistant numbers, namely, $\omega i=\u22121+(2i\u22121)/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.

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 *K _{c}*(

*N*) as follows. We acquire an upper bound on

*K*(

_{c}*N*) by determining the smallest tested value of

*K*for which a real steady state is found. Similarly, we find a lower bound on

*K*(

_{c}*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

*K*(

_{c}*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

*K*(

_{c}*N*).

In Ref. 16, it is conjectured that if there is a stable equilibrium, then there is a stable equilibrium satisfying $|\theta i\u2212\theta j|<\pi /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 $(\theta 1,\theta 2,\theta 3)\u2248(1.6820,0.8410,0)$. Since *θ*_{1} and *θ*_{3} are coupled with $|\theta 1\u2212\theta 3|>\pi /2$, we reject the conjecture.

For larger *N*, the unique stable equilibria can have angles even more spread out, such as $(\theta 1,\u2026,\theta 15)\u2248(2.0981,1.8867,\u2026,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* = *K _{c}* 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\u2264\u2308 25(N\u22121) \u2309$. Based on these results, we conjecture that this phenomenon occurs in general: for equidistant natural frequencies, if 0 ≤ *j* ≪ *N* 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.

## V. RESULTS FOR CYCLIC GRAPHS

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 landscapes^{16} 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.

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 ≤ *K _{c}*(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}

## VI. RESULTS FOR RANDOM GRAPHS

We now turn our attention to random coupling arrangements. In these computations, we choose each *a _{i}*

_{,}

_{j}as 0 or 1 according to a pre-determined probability

*P*while setting

*a*

_{j}_{,}

_{i}=

*a*

_{i}_{,}

_{j}. In other words, [

*a*

_{i}_{,}

_{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.

Next, we fix *N* = 10 and investigate how *K _{c}*(

*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

*K*(10). Figure 8 shows these results sorted according to the number of cycles in the graph.

_{c}For undirected and connected graphs, a theoretic lower and a conjectured theoretic upper bound are^{4,23}

where $degi=\u2211jai,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.

## VII. CONCLUSION

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 *K _{c}*(

*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

*K*(

_{c}*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.

## ACKNOWLEDGMENTS

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.

### APPENDIX: AN EQUIVALENT FORMULATION

In Sec. I, we remarked that the equilibrium conditions $d\theta 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

*θ*= 0 and removing the

_{N}*N*th equation from the system. We claimed no information was lost in this simplification since $\u2211i\omega 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* = [*a _{i}*

_{,}

_{j}], (1) implies that $\u2211id\theta idt=\u2211i\omega i$. In particular, equilibrium, i.e., $d\theta idt=0$ for

*i*= 1,…,

*N*, can occur only if $\u2211i\omega i=0$.

We proceed to solve $d\theta idt=0$ for *i* = 1,…,*N* by considering $\varphi i,j=\theta i\u2212\theta 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

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

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

with the *N*th equation being

Since solving $f1=\cdots =fN\u22121=fN=0$ is equivalent to solving $f1=\cdots =fN\u22121=\u2211i=1Nfi=0$, we can replace (A2) by the corresponding sum which simplifies to

using the assumption that the graph is undirected, i.e., *a _{i}*

_{,}

_{j}=

*a*

_{j}_{,}

_{i}, and the anti-symmetry property of the sine function, i.e., $sin(\u2212x)=\u2212sin(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 $\varphi i,N$ for

*i*= 1,…,

*N*− 1.

Now, we simply define $\theta \u0303i=\varphi i,N$ for *i* = 1,…,*N* − 1 and $\theta \u0303N=0$. Since $sin(\varphi i,N)=sin(\theta \u0303i\u2212\theta \u0303N)$, (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.