Lattice Landau Gauge and Algebraic Geometry
Abstract
Finding the global minimum of a multivariate function efficiently is a fundamental yet difficult problem in many branches of theoretical physics and chemistry. However, we observe that there are many physical systems for which the extremizing equations have polynomiallike nonlinearity. This allows the use of Algebraic Geometry techniques to solve these equations completely. The global minimum can then straightforwardly be found by the second derivative test. As a warmup example, here we study lattice Landau gauge for compact U(1) and propose two methods to solve the corresponding gaugefixing equations. In a first step, we obtain all Gribov copies on one and two dimensional lattices. For simple systems their number can already be of the order of thousands. We anticipate that the computational and numerical algebraic geometry methods employed have farreaching implications beyond the simple but illustrating examples discussed here.
Lattice Landau Gauge and Algebraic Geometry
\FullConferenceInternational Workshop on QCD Green’s Functions,
Confinement, and Phenomenology  QCDTNT09
September 07  11 2009
ECT Trento, Italy
1 Introduction
Finding the global minimum of a multivariate function is one of the most important tasks in statistical mechanics, condensed matter theory, lattice gauge theories and theoretical chemistry. Most of the usual methods to minimize a multivariate function are based on the NewtonRaphson method, where a start solution is guessed and then is refined by successive iterations. There are several efficient refinements of this minimization procedure, e.g., simulated annealing. However, they are often not successful as they can get trapped in a local minimum instead of the global one. This problem introduces an error of an unknown order.
A typical example is the Landau gauge on the lattice: in the continuum the corresponding Landau gauge fixing condition is solved to fix the gauge, whereas on the lattice, usually the socalled standard lattice Landau gaugefixing (LLG) functional (whose first derivatives with respect to the gauge parameters are the corresponding gaugefixing conditions on the lattice) is numerically minimized [10, 11, 12, 13]. We take the example of the LLG functional for compact U(), which is given by
(1) 
where the variables sit on the th lattice site of a dimensional lattice grid, is the unit vector in direction, and the variables sit on the links between the th and th lattice sites. Here the variables are related to the gauge potential in the corresponding continuum gauge theory and the variables are the gauge transformations with respect to which is minimized. The special case when all are zero is called the trivial orbit. We will also frequently use random choices for the link angles which we then refer to as random orbits. The Hessian matrix of this functional is the FaddeevPopov operator and its determinant is known as the FaddeevPopov determinant [18] on the lattice.
In order to compare the results from lattice Landau gauge to those from functional methods in the continuum such as studies of DysonSchwinger equations [8], one should in principle find all stationary points of under variation with respect to the gauge variables [3, 4, 6, 2], that is all solutions to the corresponding LLG equations on the lattice,
(2) 
The solutions to these equations are the lattice analogues of Gribov
copies [7].
The reason why this procedure fails is that the signs of all Gribov copies exactly cancel one another, always yielding an exact zero for their sum which computes the vanishing Euler characteristic of the lattice gauge group manifold in the standard implementation of the Landau gauge. This zero sum is the origin for the famous Neuberger problem of standard BRST formulations on the lattice [16, 17], which is the statement that the expectation value of any gaugeinvariant observable is of indefinite form in such a formulation.
The problem can be overcome by stereographically projecting the lattice gauge group onto a manifold whose Euler characteristic is nonzero [6]. This provides a welldefined lattice BRST formulation on the projected manifold and it can serve as a nonperturbative definition of BRST symmetry in the continuum limit.
The Gribov copies of compact U() are of course lattice artifacts. However, the gauge group is a direct product of odddimensional and compact manifolds (circles in this case), as it is for every SU() gauge theory on the lattice also. This is sufficient to conclude that it also shares the same problem with SU() theories. Moreover, it turns out that the Neuberger problem in SU() is avoided when that of its maximal Abelian subgroup is [9], because the coset manifold consist of even dimensional spheres whose Euler characteristic is 2. It is therefore obviously important to understand the Gribov copy and Neuberger problems in compact U() gauge models.
At the same time, for the trivial orbit the functional in Eq. (1) also represents the Hamiltonian of an XY model which stands alongside the Ising and Heisenberg models as one of the most intensely studied systems in statistical mechanics and condensed matter theory. For a random orbit, it is that of a random phase XY model (RPXYM). These have been widely used as simple models of classical superconductivity, to study aspects of high superconductors, to describe the XY magnet with random DzyaloshinskiMoriya interactions or a positional disordered Josephson junction array, to name a few. In either case, the calculations (related to Domain Wall Renormalization Group studies) boil down to obtaining the global minimum of the Hamiltonian with respect to the variables for a given set of ’s (see, e.g., [25]). In Ref. [26], it was shown in lower dimensional models that one can first perform a duality transformation then followed by a numerical minimization so that a minimization algorithm has to search only within the space of local minima. Recently, Akino and Kosterlitz [25] have used this method to obtain the global minimum of a RPXYM Hamiltonian. This method, while giving more confidence in the final results with less numerical effort, still relies on numerical algorithms which are known to fail for larger systems.
The knowledge of all stationary points is important in recent studies of the potential energy landscape of the classical statistical mechanical systems such as the RPXYM and its relation to phase transitions [21, 22, 23]. An effort to numerically classify all Gribov copies of compact U(1) in 2 dimensions was done in [27].
Here, we take up the general problem of finding all solutions to the standard LLG equations for compact U() or the stationary points of the corresponding XY model. The corresponding equations are nonlinear which makes them difficult to deal with analytically or numerically. We observe, however, that this nonlinearity is polynomiallike and the equations can be transformed into systems of polynomial equations. Then we can use Algebraic Geometry techniques to solve these systems and obtain all stationary points of the corresponding function and in particular the global minimum accurately. This procedure is applicable to a wide range of interesting physical problems such as obtaining all stationary points of the classical Hamiltonians of Ising and Heisenberg models; solving classical field equations for pure lattice gauge field theories, Abelian Higgs model, etc.; finding a generic form of SU() matrices; obtaining the vacua of Supersymmetric potentials; and many more (see Ref. [2]). The corresponding gaugefixing equations for the modified LLG, proposed in Refs. [3, 4], for the compact U() case can also be shown to have polynomiallike nonlinearity but their structure is more complicated [2]. Here, we restrict the discussion to the standard LLG.
In the onedimensional lattice case, the LLG equations have been exactly solved for both periodic and antiperiodic boundary conditions elsewhere [2, 3, 4] where the corresponding equations are treated as systems of linear equations. However, for higher dimensional lattices the corresponding equations are genuinely nonlinear and the same method can not be used. We will show that the systems of these nonlinear equations can be viewed as systems of multivariate polynomial equations and therefore can be described by the language of Algebraic Geometry. We will then discuss two promising approaches to solve the polynomial systems: the numerical polynomial homotopy continuation (NPHC) method and the computational Algebraic Geometry approach (in the Appendix) and present our results. Though our goal is to solve the gauge fixing conditions for higher dimensional lattices, we will explain the methods using the onedimensional case because: (1) we already know the exact solutions for this case and so the results from the new approaches will have a precise comparison, (2) the onedimensional case, interpreted as an Algebraic Geometry problem, provides all the essence of the higher dimensional generalizations.
2 The Extremizing Equations as Polynomial Equations
Here, we show that the problem of solving the extremizing equations in terms of the variables can be transformed into that of solving a system of multivariate polynomial equations. We first note that, for a onedimensional lattice having lattice sites, the directional index is irrelevant and is just . Thus, the corresponding extremizing equations are
(3) 
where , for all . These are the gaugefixing equations (2) on the lattice and reproduce the continuum Landau gauge () in the naive continuum limit. These equations are also the steady state equations of a variant of the famous Kuramoto model to study synchronization in mathematical biology [24]. It can be shown that by using antiperiodic boundary conditions on both  and variables the global gauge freedom can be completely fixed [3, 4, 2]. While dealing with periodic boundary conditions one can get rid of the residual gauge freedom by taking, say, to be zero and removing the equation from the system [2]. For simplicity, we take the trivial orbit case (i.e., all ). Using the trigonometric identity, and writing and , we get
(4) 
for all . This is merely a change of notation. However, we add additional equations in the system for each site , namely,
(5) 
for all . Now, the combined system of all and is not just a change of notation but all the and are algebraic variables and the equations are multivariate polynomial equations, i.e., the fact that and are originally and is taken care of by the constraint equations (5). In general, for the onedimensional lattice with lattice sites, we have in total polynomial equations and variables.
3 Numerical Polynomial Homotopy Continuation Method
In the Appendix, a method to exactly solve systems of multivariate equations, called the Groebner basis technique, and its current status are briefly explained. Here we discuss a numerical method to solve a system of multivariate equations, called the numerical polynomial homotopy continuation (NPHC) method. To explain the NPHC method let us consider a system of multivariate polynomial equations, say , where and , that is known to have isolated solutions, e.g., the above mentioned LLG equations after eliminating the global gauge freedom. Now, there is a classical result, called the Classical Bezout Theorem, that asserts that for a system of polynomial equations in variables, for generic values of coefficients, the maximum number of solutions in is , where is the degree of the th polynomial. This bound, called the classical Bezout bound (CBB), is exact for generic values (i.e., roughly speaking, nonzero random values) of coefficients, e.g., for the onedimensional LLG equations with number of lattice sites and with periodic boundary conditions, this number is (because there are polynomials each of which is a degree polynomial). The genericity is welldefined and the interested reader is referred to Ref. [28, 29] for details.
Using the socalled BernsteinKhovanskiiKushnirenko theorem, a tighter bound, which takes the sparsity of the system into account, can be computed. This bound is called the BKK root count or the mixed volume (or stable mixed volume (SMV) [29, 33, 34] in general). The reader is referred to the above references to get a precise definition of the mixed volume and its computation.
Based on either of these bounds on the number of complex solutions, a homotopy can be constructed as
(6) 
where is a random complex number. is a system of polynomial equations with the following properties: (1) the solutions of are known or can be easily obtained. is called the start system and the solutions are called the start solutions, (2) the number of solutions of is equal to the CBB or SMV for , (3) the solution set of for consists of a finite number of smooth paths, called homotopy paths, each parameterized by , and (4) every isolated solution of can be reached by some path originating at a solution of . One can then track all the paths corresponding to each solution of from to and reach . Here, a randomly chosen complex ensures the paths are wellbehaved. By implementing an efficient path tracker algorithm, all isolated solutions of a system of multivariate polynomials system can be obtained.
To illustrate how the method works, we take a univariate polynomial from Ref. [28] to make this discussion selfconsistent, say , pretending that we do not know its solutions, i.e., . We first define a family of problems as
(7) 
where is a parameter. For , we have and at we recover our original problem. The problem of getting all solutions of the original problem now reduces to tracking solutions of from where we know the solutions, i.e., , to . The choice of the piece in Eq. (7), called the start system, should be clear now: this system has the same number of solutions as the CBB of the original problem and is easy to solve. Now, one of the ways to track the paths is to solve the differential equation that is satisfied along all solution paths, say for the th solution path:
(8) 
This equation is called the Davidenko differential equation. Inserting (7) into this equation, we have
(9) 
This initial value problem can be solved numerically (again, pretending that an exact solution is not known) with the initial conditions as and . The other approach is to use Euler’s predictor and Newton’s corrector methods. We do not intend to discuss the actual path tracker algorithm used in practice, but it is important to mention that in the path tracker algorithms used in practice, almost all apparent difficulties have been resolved, such as tracking singular solutions, multiple roots, solutions at infinity, etc. For the sake of completeness, we should also mention here that in the actual path tracker algorithms the homotopy is randomly complexified to avoid singularities by , called the gamma trick, i.e., taking
(10) 
There are several sophisticated computational packages such as PHCpack [30], PHoM [32], Bertini [28] and HOM4PS2 [31] which can be used to solve systems of multivariate polynomial equations. They are all available as freewares from the respective research groups. In the Appendix, we compare the Groebner basis technique and the NPHC method with technical remarks.
3.1 Results
We present the results by classifying the obtained solutions in terms of the number of positive and negative eigenvalues of the corresponding Hessian matrix, or the FaddeevPopov operator, because then we can use the Neuberger zero as a necessary condition for having all the solutions.
Antiperiodic Boundary Conditions
We now explore the simplest nontrivial case in higher dimensional lattices which is the standard LLG functional on a lattice with the trivial orbit and antiperiodic boundary conditions
(i.e., the classical XY model) and also a random orbit case (i.e., randomly chosen ).
For both cases, the CBB was and the SMV was . For the trivial orbit case, the total number of real and complex solutions was out of which there were real solutions, i.e., Gribov copies. Similarly, for the random orbit case, the total number of real and complex solutions was out of which there were real solutions

A solution means that a set of values of ’s and ’s satisfies each of the equations with tolerances . All the solutions come with real and imaginary parts. A solution is a real solution if the imaginary part of each of the variables (i.e., all and ) is less than or equal to the tolerance (below which the number of real solutions does not change, i.e., it is robust for all the cases we consider in this discussion). The original trigonometric equations are satisfied with tolerance after and are transformed back to . All these solutions can be further refined with an arbitrary precision. This is a remarkable success of the method because then these solutions are close to the exact solutions.

The sum of the signs of the FaddeevPopov determinants for all real solutions for the trivial and random orbit cases is , which is numerically zero, yielding the expected Neuberger zero as mentioned in the Introduction.

For the trivial orbit case, there are exactly real solutions which have zero FaddeevPopov determinant with tolerance . These solutions constitute the set of singular loci or the Gribov horizons. They can be further classified in terms of the number of zero eigenvalues of the FaddeevPopov operator at each of the solutions. This amounts to classifying the singular solutions of the polynomial system in terms of their multiplicities using the socalled deflation singularities technique[37]. For both the random orbit cases, there is no Gribov horizon, i.e., all the solutions are nonsingular.

For the trivial orbit case, the remaining real solutions have nonsingular FaddeevPopov determinant, i.e., nonsingular solutions. All the nonsingular solutions of both the trivial and the random orbit cases can be classified by the number of negative eigenvalues (Table 2), which shows the expected twofold symmetry giving rise to the Neuberger zero.
boundary conditions orbit real solutions sum of total number zero det. nonzero det. signs of det. antiperiodic trivial 2968 1152 1816 0 random 2480 0 2480 0 periodic trivial 270 52 218 0 random 224 0 224 0 Table 1: An overview on the number of real solutions for the different cases (orbits and boundary conditions). See Section for the detailed description. Case i 0 1 2 3 4 5 6 7 8 9 antiperiodic b.c., trivial 2 18 216 342 330 330 342 216 18 2 antiperiodic b.c., random 2 58 202 402 576 576 402 202 58 2 periodic b.c., trivial 8 28 38 42 56 38 7 1   periodic b.c., random 2 12 30 61 72 38 8 1   Table 2: Summary of the number of real solutions with negative eigenvalues on the lattice, for the different cases (orbits and boundary conditions).
See Table 1 for the summary.
Periodic Boundary Conditions
After eliminating the constant zero mode for the lattice case with periodic boundary conditions (by taking, say, two of the ’s at the corners of the lattice to be , without losing generality. See Ref. [2] for more details), we have equations and variables. The number of solutions is out of which are nonsingular and are singular. The FaddeevPopov operator is now a dimensional matrix. The sum of signs of the FaddeevPopov determinants is zero. Similarly, for a random orbit the number of real solutions is and again the sum of signs of the FaddeevPopov determinants is . The results are summarized in Table 2.
4 Summary
By taking the standard lattice Landau gaugefixing (LLG) functional (the classical Hamiltonian of the random phase XY model (RPXYM)) as an example, we showed that the nonlinearity of the corresponding gaugefixing equations (extremizing equations) is polynomiallike. With additional constraint equations, the combined system of equations can be treated as a system of polynomial equations with all variables defined over . Though we are interested in only real solutions of these polynomial systems, we used Complex Algebraic Geometry concepts over Real Algebraic Geometry due to the stronger results available in the former. To solve the corresponding polynomial systems exactly, one can use an elegant algorithm, called the Buchberger algorithm (see Appendix A). Though the algorithm works well for the onedimensional lattice, for the lattice, it is beyond the capabilities of a standard desktop machine. However, the method can be useful in finding a parametrization of the SU() matrices [39] and hence can then be used, for example, to extend the method of restricting a path integral of the SU() GeorgiGlashow model to the ’t HooftPolyakov sectors by using the twisted Cperiodic boundary conditions [5] for odd and to precisely define the modified lattice Landau gauge for the SU() case.
For the twodimensional case, we used the Numerical Polynomial Homotopy Continuation (NPHC) method, which gives all solutions of a polynomial system numerically. This method does not suffer from the technical difficulties of the Buchberger algorithm or its variants, and in principle this method can find the solutions for those systems which may be intractable for the former method. By computing the eigenvalues of the Hessian matrix at these solutions, one can obtain the global minimum, up to machine precision, of a multivariate function whose extremizing equations can be translated to polynomial equations.
The antiperiodic boundary conditions system for the trivial and a random orbit can be solved using the NPHC method and that there is an exact cancelation of the signs of the FaddeevPopov determinants in this case, i.e., the Neuberger zero. We classified all the solutions in terms of the number of negative and positive eigenvalues of the corresponding FaddeevPopov operator. Apparently, the number of Gribov copies is orbit dependent, though the Neuberger zero is orbit independent. It would be important to compute the number of solutions for different random orbits for the lattice with the modified LLG using the NPHC method. Efforts to solving the corresponding gaugefixing equations for the modified LLG for compact U() and also the linear covariant gaugefixing equations for compact U(), in the spirit of [41, 42], are in progress. It would also be very interesting to solve the Langevin dynamics equations for the XY model and other models [40] using algebraic geometry methods.
5 Acknowledgement
DM would like to thank the developers of the PHCpack, HOM4PS2, Bertini, PHoM and Singular, and to Hirokazu Anai, Guillaume Moroz, Antonio Montes, Akira Suzuki, Vladimir Gerdt, Daniel Robertz, Anton Ilderton, YangHui He, JM Kosterlitz, JonIvar Skullerud, Teresa Mendes, Attilio Cucchieri and Paul Watts for their help and critical remarks on this work. DM was supported by the Science Foundation of Ireland and by the British Council Researcher Exchange Programme to visit the Imperial College London. Support by the Australian Research Council is also acknowledged.
Appendix A Computational Algebraic Geometry
Since many important and useful results in algebraic geometry are mainly developed for the complex field , we treat our polynomial equations as having complex variables, and in the end we use only real solutions and throw the remaining solutions away (see, e.g., [19]). One can then use an important result from computational algebraic geometry: after specifying an ordering of the monomials of the system (denoted as for the ordering in which is placed prior to and is placed prior to ), one can transform a given system of multivariate polynomial equations, referred to as an ideal, to another one which has the same solutions but is easier to solve. This new system is called a Groebner basis for the given monomial ordering.
There is a well defined procedure to find a Groebner Basis for any given ideal and monomial order, called the Buchberger algorithm. It should be noted that the Buchberger algorithm reduces to Gaussian elimination in the case of linear equations, i.e., it is a generalization of the latter. Similarly it is a generalization of the Euclidean algorithm for the computation of the Greatest Common Divisors of a univariate polynomial. Recently, more efficient variants of the Buchberger algorithm have been developed to obtain a Groebner basis, e.g., F4, F5 and Involution Algorithms [35]. Symbolic computation packages such as Mathematica, Maple, Reduce, etc., have builtin commands to calculate a Groebner basis for a given monomial. Singular, COCOA and McCauley2 are specialized packages for Groebner basis and Computational Algebraic Geometry, available as freeware. MAGMA is also such a specialized package available as a nonfree package. Rather than going into the details of the specifics and technicalities of this algorithm, we dive into the practical applications of the Groebner basis technique relevant to our problem and refer the reader to the above mentioned references for further details.
a.1 The Groebner Basis Technique At Work
Here, we provide a practical example of how the Groebner basis technique can be used for our systems. Firstly, for the polynomials for the onedimensional lattice with three lattice sites, antiperiodic boundary conditions and the trivial orbit case, the corresponding ideal (denoted as the polynomials put between ’’ and ’’) is
(11)  
For the ideal in Eq. (11), a Groebner basis for the socalled lexicographic ordering is
(12)  
As noted earlier, the solutions of this system are the same as the original system. Here, the first equation in is a univariate polynomial in variable and solving it is simple because it can be factorized as giving and . Using this and backsubstitution, all the solutions can be obtained and are as follows:
(13) 
Thus, the solution space of the system , called the variety of and denoted as , is the above mentioned set of isolated points in a dimensional affine space. This is the known result from Refs. [3, 4, 2] in terms of ’s, i.e.,
(14) 
All of these solutions are correctly reproduced in Eq. (13). Similar computation can be performed for the periodic boundary conditions case [2].
For an ideal which is known to have only isolated solutions (called a zerodimensional ideal) with a lexicographic ordering one can always find a Groebner basis in an upper diagonal form, analogous to the Gaussian Elimination method, such that at least one polynomial is univariate with the others having an increasing number of variables.
Interpreting the gauge fixing equations in terms of Algebraic Geometry allowed us to deal with the actual nonlinearity of the equations, rather than treating them as linear equations as done in Refs. [2, 3, 4]. Specifically, in this interpretation the corresponding method does not make any distinction between the equations arising from a onedimensional lattice or those arising from a higher dimensional lattice. In theory, as long as one can obtain a Groebner basis, the equations can be exactly solved. For the onedimensional case, it was easy to go up to a fairly big number of lattice sites on a single machine. However, in general, obtaining a Groebner basis is very difficult due to an algorithmic complexity, known as Exponential Space complexity, which roughly means that the RAM required by the computation blows up exponentially. In particular, on a regular single desktop machine with GB RAM, we could not obtain a Groebner basis for the trivial orbit case (i.e., the classical XY model) and with antiperiodic boundary conditions for a lattice, using Singular3.2. The corresponding system is made of equations, each of degree , in variables. However, a more powerful machine should be able to obtain a Groebner basis for these systems.
Recently, V. Gerdt and Daniel Robertz were able to compute a Groebner Basis for this system over using a Linux machine with AMD Opteron (TM) processor MHz, cores, GB memory and with Magma (V2.14.14) in less than hours with degree reverse lexicographic ordering [36]. This is a very important step towards solving the system. However, since the computation is over not all the solutions over can be obtained.
a.2 Comparison between the Groebner basis technique and the NPHC Method
The NPHC method is strikingly different from the Groebner basis technique in that the algorithm for the former suffers from no known major complexities. Moreover, the path tracking is embarrassingly parallelizable, because all the start solutions can be tracked completely independently of each other. This feature along with the rapid progress towards the improvements of the algorithms makes the NPHC well suited for many physical problems arising in condensed matter theory, lattice QCD, etc.
For most systems of polynomial equations in practice, we do not know the actual number of solutions from the beginning. So even after obtaining all solutions through the homotopy continuation method, we would still prefer to have some kind of verification of the solutions. In the standard LLG fixing case, we used the Neuberger zero as a necessary condition. However, this is certainly not a sufficient condition and further checks may be required for bigger systems. This is called certifying the solutions, which has been recently developed using the socalled Tropical Algebraic Geometry [38]. We anticipate that our results mentioned above will serve as ideal test systems for all such new ideas.
To solve the corresponding equations for the periodic boundary conditions case, another method is to leave the constant zero modes in the system and use the Numerical Algebraic Geometry (NAG) [43], which is the generalization of the NPHC method to positivedimensional ideals.
Recently, a computational Algebraic Geometry approach to find vacuum configurations in string phenomenology (where mathematically the problem reduces to minimization of the socalled superpotentials which exhibit polynomiallike nonlinearity) is proposed in Refs. [19, 20]. A very efficient Mathematica package, called STRINGVACUA, has been developed with an interface with the computer algebra system Singular. This package relies on the Groebner basis technique in solving the extremizing polynomial equations. It is anticipated that incorporating the NPHC method can extend the applicability of the package and the approach beyond its present status.
It should be noted that there are several promising methods such as the Discriminant Variety [44] and finding real solutions out of complex curves [45] from computational real algebraic geometry which are emerging as alternatives to the above mentioned methods in some of the important problems mentioned in the Introduction.
Footnotes
 For recent progress in the continuum, see also [14, 15].
 Note that we do not intend to compare the efficiency of the available packages. However, to give an idea how computationally intensive such a calculation is, we note that the HOM4PS2 package took around and minutes, respectively, to run these systems on a Linux singleprocessor desktop machine.
References
 Dhagash Mehta, Ph.D. Thesis 2009, The Uni. of Adelaide, Australasian Digital Theses Program.
 L. von Smekal, D. Mehta, A. Sternbeck and A. G. Williams, PoS LAT2007 (2007) 382 [arXiv:0710.2410 [heplat]].
 L. von Smekal, A. Jorkowski, D. Mehta and A. Sternbeck, PoS CONFINEMENT8 (2008) 048 [arXiv:0812.2992 [hepth]].
 S. Edwards, D. Mehta, A. Rajantie and L. von Smekal, Phys. Rev. D 80 (2009) 065030 [arXiv:0906.5531 [heplat]].
 L. von Smekal, arXiv:0812.0654 [hepth].
 V. N. Gribov, Nucl. Phys. B 139 (1978) 1.
 R. Alkofer and L. von Smekal, Phys. Rept. 353 (2001) 281 [arXiv:hepph/0007355].
 M. Schaden, Phys. Rev. D 59 (1999) 014508 [arXiv:heplat/9805020].
 T. Iritani, H. Suganuma and H. Iida, arXiv:0908.1311 [heplat].
 I. L. Bogolubsky, E. M. Ilgenfritz, M. MullerPreussker and A. Sternbeck, Phys. Lett. B 676 (2009) 69 [arXiv:0901.0736 [heplat]].
 A. Cucchieri and T. Mendes, PoS LAT2007 (2007) 297 [arXiv:0710.0412 [heplat]].
 D. B. Leinweber, J. I. Skullerud, A. G. Williams and C. Parrinello [UKQCD Collaboration], Phys. Rev. D 60 (1999) 094507 [Erratumibid. D 61 (2000) 079901] [arXiv:heplat/9811027].
 A. Ilderton, In the Proceedings of 9th Workshop on NonPerturbative Quantum Chromodynamics, Paris, France, 48 Jun 2007, pp 28 [arXiv:0709.1671 [hepth]].
 A. Ilderton, M. Lavelle and D. McMullan, JHEP 0703 (2007) 044 [arXiv:hepth/0701168].
 H. Neuberger, Phys. Lett. B 175, 69 (1986).
 H. Neuberger, Phys. Lett. B 183, 337 (1987).
 L. D. Faddeev and V. N. Popov, Phys. Lett. B 25 (1967) 29.
 J. Gray et. al., Comput. Phys. Commun. 180 (2009) 107 [arXiv:0801.1508 [hepth]].
 J. Gray, Y. H. He, A. Ilderton and A. Lukas, JHEP 0707 (2007) 023 [arXiv:hepth/0703249].
 Marco Pettini, Geometry and topology in Hamiltonian dynamics and statistical mechanics, Interdisciplinary Applied Mathematics 33. New York, Springer. xvi, 451 p.,(2007).
 M. Kastner, Rev. Mod. Phys. 80 (2008) 167.
 D. J. Wales, Energy Landscapes (Cambridge University Press, Cambridge, England, 2007).
 J. A. Acebron, et. al., Rev. Mod. Phys. 77 (2005) 137.
 N. Akino and J. M. Kosterlitz, Phys. Rev. B 66, 5 (2002) 054536,
 M. NeyNifle and H. J. Hilhorst, Phys. Rev. B 51, 13 (1995) 8357,
 P. de Forcrand and J. E. Hetrick, Nucl. Phys. Proc. Suppl. 42 (1995) 861 [arXiv:heplat/9412044].
 A. J. Sommese and C. W. Wampler, The numerical solution of systems of polynomials arising in Engineering and Science, (2005) World Scientific Publishing Company.
 T. Y. Li, Handbook of numerical analysis, XI, p. 209304 (2003).
 J. Verschelde, ACM Trans. Math. Softw., 25, 2, (1999) 251.
 T. L. Lee, T. Y. Li, and C. H. Tsai, HOM4PS2.0: A software package for solving polynomial systems by the polyhedral homotopy continuation method, Preprint, 2008.
 Gunji Takayuki et. al., Computing 73, 1 (2004) 57.
 T. Y. Li and Xiaoshen Wang, Math. Comp. 65 (1996) 1477.
 B. Huber and B. Sturmfels, Disc. Comp. Geom. 17 (2) (1997) 137.
 V. P. Gerdt and Y. A. Blinkov, Math. Comput. Simul. 45, 56 (1998) 519,
 Private communication with V. P. Gerdt and D. Robertz.
 T. Ojika, S. Watanabe and T. Mitsui, J. of Math. Analysis and App. 2, 96 (1983) 463.
 D. Adrovic and J. Verschelde (2008), eprintarXiv: math/0809.0298
 S. Cacciatori, D. Mehta, J. I. Skullerud, Work in progress.
 G. Aarts, Phys. Rev. Lett. 102 (2009) 131601 [arXiv:0810.2089 [heplat]].
 A. Cucchieri, T. Mendes and E. M. S. Santos, arXiv:0907.4138 [heplat].
 A. Cucchieri, A. Maas, T. Mendes, Comput. Phys. Commun. 180 (2009) 215, arXiv:0806.3124 [heplat].
 A. J. Sommese et. al., Algorithms and Computation in Mathematics, A. Dickenstein and I.Z. Emiris (Eds.), Solving Polynomial Equations: Foundations, Algorithms, and Applications, 14 (2005) 339.
 D. Lazard, F. Rouillier, Technical Report RR5322, INRIA, October 2004
 Y. Lu, D. J.Bates, A. J.Sommese, C. W. Wampler, In Proceedings of the Midwest Algebra, Geometry and Its Interactions Conference, Contemporary Mathematics, 2007.
 W. Hanan, D. Mehta, G. Moroz, S. Pouryahya, Proceeding submitted to Joint Conference of ASCM2009 and MACIS2009, Japan, 1417 Dec., 2009.