Exploring the Potential Energy Landscape of the Thomson Problem via Newton Homotopies

Problem via Newton Homotopies Dhagash Mehta,1, 2, a) Tianran Chen,3, b) John W. R. Morgan,4, c) and David J. Wales4, d) 1)Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA. 2)Centre for the Subatomic Structure of Matter, Department of Physics, School of Physical Sciences, University of Adelaide, Adelaide, South Australia 5005, Australia. 3)Department of Mathematics, Michigan State University, East Lansing, MI 48823, USA. 4)University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, UK.


I. INTRODUCTION
Potential energy landscape methods have recently attracted a great deal of research activity in the chemical physics community, with applications for diverse systems, such as metallic clusters, biomolecules, structural glass formers, coarse-grained models of soft matter, thermodynamics, etc. 1,2 Here the potential energy landscape is understood as the surface defined by the potential function of the given system, which is a multivariate nonlinear function V (x) with x = (x 1 , . . . , x n ).
Stationary points (SPs) of the landscape, defined as solutions of the equation ∇V (x) = (∂V /∂x 1 , . . . , ∂V /∂x n ) = 0, are of particular interest as they can be used in global a) Electronic mail: dmehta@nd.edu b) Electronic mail: chentia1@msu.edu c) Electronic mail: jwrm2@cam.ac.uk d) Electronic mail: dw34@cam.ac.uk optimization [3][4][5] and thermodynamic sampling to overcome broken ergodicity. [6][7][8][9] Moreover, SPs can also be employed in the kinetic transition network approach to rare event dynamics. [10][11][12][13][14][15][16][17][18] The stability of the SPs can be analyzed by computing the eigenvalues of the Hessian matrix of V (x), H(V (x)) ij = ∂ 2 V /∂x i ∂x j , at the SPs. An SP is termed degenerate if H(V (x)) at this point has any zero eigenvalues beyond those corresponding to overall translation and rotation, and nondegenerate otherwise. In either case, the number of negative eigenvalues of H(V (x)) corresponds to the Morse-index or simply the index.
Traditionally, minima and transition states (SPs with index 1) have played the most important role in chemical reaction rate theory. In many-body physics, the influence of SPs with higher index saddles in governing the phase transitions has been widely investigated. 2 Recently, phase space geometry and dynamics associated with SPs of higher index have also been considered. [19][20][21] Although analytical solutions for statonary points are unlikely to be tractable in general, computational methods such as the Newton-Raphson approach can rescue the situation. However, most numerical methods have shortcomings of their own, and are not guaranteed to find all the solutions. In particular, the Newton-Raphson method breaks down at singular solutions; 22,23 and the gradientsquare minimization method 24,25 may give numerous spurious solutions, which are not physically relevant. [26][27][28] Other more sophisticated approaches based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, 29,30 and eigenvector-following 26,28 coupled with singleand double-ended searches, 31-38 may be better alternatives.
For potentials having polynomial-like nonlinearity, one can employ symbolic algebraic geometry methods based on the Gröbner basis technique, 39,40 the numerical polynomial homotopy continuation (NPHC) method, [41][42][43][44] or interval based methods, 45 all of which guarantee finding all the complex solutions. Here, the NPHC method has a practical advantage due to its parallelizability, and it has been used to analyse various potential energy landscapes in recent work. [46][47][48][49][50][51][52][53][54][55][56][57][58][59] For a general nonlinear system with nonpolynomial nonlinearity, only the interval based method guarantees finding all solutions, but the method scales badly even for moderately sized systems. Recently, a few specialized homotopies have been claimed to find all isolated solutions for such nonlinear systems, 60,61 but rigorous analysis is still outstanding. As long as the potential is bounded from above and below, an inversion-relaxation method 62 may also find many solutions, provided that one is able to find all the maxima of the potential.
Though not yet directly applied in molecular science, derivative-free soft computing methods have also proved to be efficient in finding many solutions of nonlinear systems of equations. [63][64][65] However, they have yet to be tested for larger systems.
In all of the above methods, certifying if the numerical solution indeed corresponds to an ac-tual solution of the given system may also turn out to be an important issue. [66][67][68] We will address this aspect of the problem in future work.
In the present contribution, which is a follow up to our earlier Communication, 23 we describe in detail a previously known but rather underutilized method known as Newton homotopy. We also explain our implementation of the method specific to exploring the potential energy landscapes of atomic clusters. In addition, we have obtained novel results for the potential energy landscapes of the Thomson problem. Here our motivation is provided by Smale's 7th problem in the celebrated list of eighteen unsolved mathematical problems for the 21st century. 69 This connection suggests that progress on understanding the energy landscapes of the Thomson problem could provide the foundation for fundamental insight into the general issue of global optimisation, as briefly discussed in the Conclusions.
In §II, we first review the construction of the Newton homotopy for locating SPs of a potential function. §II A highlights one key advantage of the Newton homotopy method, namely its exceptional ability in locating degenerate SPs, where methods like Newton-Raphson and gradient-square minimization are relatively inefficient. In §III we briefly review the Thomson problem. Finally, in §III A we present results for the application of the Newton homotopy method to the Thomson problem, including some new results. More detailed discussion of the technical aspects of the Newton homotopy method are included in Appendix §A.

II. NEWTON HOMOTOPY
This section briefly describes the construction of the Newton homotopy in the restricted context of locating SPs of a given potential function.
Our goal is to locate SPs x = (x 1 , . . . , x n ) ∈ R n of the target potential function V : U → R defined on an open domain U ⊆ R n . That is, we would like to find points x ∈ R n such that ∇V (x) = 0. Among many different approaches for solving this problem, in the present article, we focus on the effective but underutilized method of Newton homotopy. [70][71][72][73][74] The general idea of the homotopy continuation method is to deform the target problem into a closely related problem whose solution is already known. Then, by keeping track of the movement of the solution under the deformation, one can potentially reach the target solution.
In this approach we consider the problem of finding SPs of V (x) within a closely related family of potential functions: for a given point a ∈ R n , we define a family of potential functions smoothly parametrized by a new variable t, given bŷ This family contains the target potential V (x) as a member, since at t = 1,V 1 (x) ≡ V (x). As t varies,V t (x) represents a smooth deformation of the target potential function V (x). This family is constructed so that for one of its members, the problem of finding a SP is trivial: at t = 0, and its gradient vector, with respect to x only, is Clearly, by setting x = a, one obtains the equation ∇V 0 (a) = 0. Therefore x = a is an SP of the functionV 0 : R n → R. We now consider the effect of the deformation (1) on the SP x = a, that is, how the SP x = a ofV 0 relates to SPs of nearby members ofV t as t moves away from t = 0. As a point of departure, we shall assume that for a choice of a ∈ R n , H(V ) at x = a is nonsingular. Since H(V 0 (x)) = H(V (x)), this is equivalent to the assumption that x = a is a nondegenerate SP ofV 0 . This assumption will be made explicit in Assumption 1. One important result from the local theory of smooth real-valued functions is that under a sufficiently small perturbation, nondegenerate SPs remain nondegenerate. Applying this observation to the familyV t [as defined in (1)], as t varies sufficiently close to 0, the SP x = a ofV 0 simply moves accordingly to SPs ofV t while remaining nondegenerate, tracing out a small piece of smooth curve. Stated formally, there is a smooth curve γ 0 ⊂ R n ×R, however short, such that for every point (x, t) ∈ γ 0 , x is a nondegenerate SP of the real-valued func-tionV t for the corresponding t value.
The key step of the homotopy continuation approach hinges on the "continuation" of the small piece of curve γ 0 into a curve that connects the SP x = a ofV 0 to a solution we aim to locate, i.e. an SP of the target potential functionV 1 ≡ V . The possibility for such continuation depends on both the structure of V (x) itself and the choice of the point a ∈ R n . For later reference, we state the condition under which such continuation is possible as the following assumption: (1). We assume that the connected component Under this assumption, let γ be the smooth curve, in C, that contains the point (a, 0). That is, it is the unique curve defined by in the x-t space passing through (a, 0). Starting from this point, one can then numerically trace along this curve, and when a point x = x 1 and t = 1 is reached, the point x 1 is necessarily an SP of the target potential V (x) sinceV 1 ≡ V .
There is a rich selection of sophisticated numerical methods for tracing the curve γ defined by (3). However, within the community of homotopy continuation methods, a special class of methods based on a predictor-corrector scheme has emerged as the method of choice for its superior stability and efficiency. For completeness, we include a technical discussion of the scheme in Appendix §A. It is worth noting that under the smoothness assumption, the smooth curve γ can have "turning points" with respect to the t-direction, as illustrated in Figure 1. Consequently, the variable t cannot be used as a global parameter for the curve. The parametrization by arc-length (Appendix §A) in which the parameter represents the arclength from the starting point x = a is therefore preferred in numerical schemes for tracing the curve γ.

A. Locating degenerate stationary points
One important advantage of the Newton homotopy method (and homotopy methods in general) is the ability to locate degenerate SPs, where methods like the Newton-Raphson approach may be inefficient. Consider, for example, the real-valued function Its SPs satisfy 22 This system of equations has only one solution in R 2 , namely (0, 0), which is singular. It is shown by Griewank and Osborne (hence, the subscript "GO") in 22 that starting from almost all points in R 2 \ {(0, 0)}, the Newton-Raphson method will diverge. In other words, the Newton-Raphson method is very unlikely to find the SP of V GO . In contrast, in our earlier paper, 23 we obtained the singular solution of this system using the Newton homotopy method. Here, we describe the details of our computation for this system. Applying Newton homotopy (3) to V GO , with starting point (a, b) ∈ R 2 , the problem is transformed to tracing a curve defined by It can be verified that for almost all (a, b) in R 2 , the Newton homotopy H GO , defined above, satisfies the smoothness assumption. Furthermore, for all (a, b) ∈ R 2 such that b > 119a 2 /128, the equation H GO (x, y, t) = (0, 0) defines a smooth curve containing both the starting points (a, b, 1) and (0, 0, 0). In other words, starting from any point (a, b) with b > 119a 2 /128, the Newton homotopy method will reach the degenerate SP (0, 0). We include the proofs in Appendix §B.
For example, choosing the starting point (a, b) = (1, 15/16), the Newton homotopy method obtains the target degenerate SP at (0, 0) without encountering any numerical instability, since the curve defined by the Newton homotopy (6) is actually smooth. Hence the smoothness assumption (Assumption 1) is satisfied. Therefore, numerical curve tracing algorithms can be employed without any instability problems. Figure 2 shows the projection of the smooth curve defined by (6) using the starting point (a, b) = (1, 15/16) onto the x-t space. The target degenerate SP is simply a smooth point on the curve. The degeneracy is reflected by the fact that the curve meets t = 1 at a tangency. However, as long as appropriate curve tracing techniques, such as those described in Appendix §A, are used, such tangency has no impact on the numerical stability of the algorithm. In other words, even though the degenerate SP (0, 0) is a singular solution of the system of equations ∇V GO = 0, and hence infinitely ill-conditioned, it is actually a smooth point if considered as a member of the family of SPs of the family of potentials (1). That is, by treating the problem as a specific instance among a deformation of problems, the degeneracy disappears.
More generally, this result appears to represent typical behavior for Newton homotopy near degenerate SPs, especially when the Hessian matrices at the SPs have corank 1. Indeed, this phenomenon is observed in our numerical experiments with the Thomson problem, described in §III: a large number of degenerate SPs are found by the Newton homotopy without encountering any numerical instability.

III. THE THOMSON PROBLEM
In this section, we describe the Thomson problem and apply the Newton homotopy method to explore the potential energy landscape of the model. Although the Thomson model was proposed as a representation of atomic structure, 79 it has since found widespread application in various other areas of chemical physics. Physically, the problem is to find the minimum energy of a system composed of N electrons constrained to move on the surface of a sphere (for simplicity, of unit radius), considered as classical point charges. Although the global minima for N = 4, 6 and 12 are the expected platonic solids, those for N = 8 and 20 are not. Generalizations of the model have also been used to study surface ordering of liquid metal drops confined in Paul traps, to model clusters of proteins on a shell, colloid particles, fullerene structures for carbon clusters, etc. (see, e.g., [80][81][82][83][84].
From the optimization point of view, the model has become a standard benchmark for global optimization algorithms, and there have been many numerical attempts to find local and global minima. [85][86][87][88][89][90] Deterministic results are only known for N = 2, 3, 4, 5, 6 and 12. 91,92 The potential energy function of the Thomson problem is given by where Our goal is to locate the SPs of V N Th that are configurations of N points on the sphere at which ∇V N Th is perpendicular to the sphere. Stated algebraically, for this constrained optimization problem, we can consider the Lagrangian The goal is then to locate the SPs of L N Th . The problem can be significantly simplified by encoding the spherical constraints implicitly via spherical coordinates. Here we use the 2N variables φ i , θ i for i = 1, . . . , N and the transformation Since V N Th and hence L N Th are invariant under any rotation, every SP is represented by a set with three degrees of rotational freedom. This rotational symmetry can be removed by fixing It is important to note that the above transformation is singular at the "poles" (0, 0, ±1) of the sphere and hence can potentially introduce extraneous or degenerate SPs. However, such extraneous SPs are easily identified, and, as noted in §II A, the Newton homotopy is well suited for handling the degeneracy caused by the transformation.

A. Stationary Points of the Thomson Problem
When applied to the Thomson problem (10) (with spherical coordinates), the Newton homotopy method was able to quickly locate a large number of SPs over a wide range of Morse indices.
Since the transformation to spherical coordinates (9) is singular at the poles, it can potentially introduce extraneous results that are solutions to (10) (in the spherical coordinates), but not SPs of (8). To filter out these points, each solution of (10) found is transformed back into Cartesian coordinates. Recall that at an SP of V N Th , the vector ∇V N Th is perpendicular to the surface of the spheres. Since at this point the normal direction of the spheres is given by y 1 , z 1 , . . . , x N , y N , z N ), the angle between −∇V N Th and η can be used as a criterion for identifying extraneous solutions.
In our application to the Thomson problem, we consider a solution extraneous when the cosine of the angle is more than 10 −5 away from 1 (indicating a nonzero angle), and such points are discarded. Table I shows a selection of SPs of V N Th found by the Newton homotopy method described in §II classified by their indices. 100 random starting points are used for each N value. While not every starting point defines a curve that passes through t = 1 and produces an SP, collectively they are able to locate a large number of SPs over a range of different indices.
Since there is an inherent permutation symmetry in the formulation (7), each geometric configuration of the N charges is represented by multiple SPs of (7) (e.g. (x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ) and (x 2 , y 2 , z 2 , x 1 , y 1 , z 1 ) are distinct points in R 6 but they both represent the same geometric configuration of two charges). Table I indicates the distinct values for V N Th . Table II provides some reference data for sizes N = 4 to 7, which was used as a check.
As noted in §II A, the Newton homotopy method (and homotopy methods in general) have a special advantage in locating degenerate SPs when compared to the Newton-Raphson method and gradient-square minimization. This observation is supported by our numerical experiments with the Thomson problem, for which Newton homotopy was able to locate several degenerate SPs without encountering any numerical instability. Table III lists a selection of SPs of V Th that appear to be degenerate in the sense that at the approximated SP at least one eigenvalue of H(V Th ) is very close to zero. Further analysis would be necessary to distinguish the degenerate SPs from the nondegenerate SPs listed in Table I. We leave this analysis for future work.
Equally significant is the ability of the Newton homotopy to obtain multiple SPs from a single starting point, which is illustrated in Figure 3. This capability is even more pronounced in other cases. For example, Figure 4 shows the t-values along another curve defined by the Newton homotopy when applied to the N = 5 Thomson problem. The curve passes through t = 1 at least 9 times and produces 8 (numeri-cally) nondegenerate SPs and one (numerically) degenerate SP.
Our experiments with the Thomson problem suggest that this result is quite typical, as the majority of the curves that intersect t = 1 would continue and turn around to intersect t = 1 again.
Note that for larger N our implementation finds SPs of indices near the middle of the possible range, that is, close to (2N − 3)/2 in the Thomson problem case. We observed the same type of behaviour in the case of Lennard-Jones clusters with our the Newton homotopy earlier, 23   there are exponentially more SPs in the midrange of indices than with indices closer to 0 or (2N − 3) as already shown with general arguments applicable to most potential energy functions. 28

IV. PHYSICAL INTERPRETATION OF THE NEWTON HOMOTOPY
This section provides a physical interpretation of the Newton homotopy construction in the context of the Thomson problem, although the same interpretation can be applied to any problem of locating SPs of potential energy landscapes. For simplicity, the natural xyz Cartesian coordinate system is used: For an in-  teger N > 1, let V N Th (x) be the potential defined in (7) where x = (x 1 , y 1 , z 1 , . . . , x N , y N , z N ) ∈ R 3N encodes the positions of the N charges.
As described in §II, the Newton homotopy method starts with a random starting point a ∈ R 3N on the potential energy landscape, that is, a random starting configuration of the N charges on the sphere. Unless x = a is already an SP of V N Th , we have ∇V (a) = 0. For each individual electron, say the kth electron, the combined force exerted by the rest of the charges according to the Coulomb potential on this electron is given by f k := (∂V /∂x k , ∂V /∂y k , ∂V /∂z k ) at x = a.
We then apply an artificial counterbalancing force to each electron so that the net force it ex-  periences is zero. Specifically, by applying − f k to the k-th electron, the net force (of the combined effect of the Coulomb potential and the artificial force) on this electron is exactly zero. Note that the potential energy induced by this force alone is represented by − f k · (x k , y k , z k ) and therefore the combined potential energy of the system of N charges is given bŷ Note that the random starting configuration x = a is indeed an SP of this augmented energy landscape since Next, we keep track of the movement of the SP as the artificial counterbalancing force is gradually weakened. This is done by considering the potential with a parameter t (11) in which the parameter t plays the role of time in the continuous weakening of the artificial coun-terbalancing force. When t reaches 1, the artificial force is completely removed and the augmented potential energyV 1 is identical to the original potential V N Th . Comparing (11) with (1), this process describes exactly the mathematical construction of the Newton homotopy.
x = a fold fold Figure 5. An illustration of a trajectory of the moving SP as the artificial force in (11) is gradually weakened to zero (at t = 1) and then reverse direction (t > 1). Each turning point corresponds to a fold catastrophe that potentially produces a new configuration.
The trajectory of the moving SP defined by (11) can certainly be extended beyond t = 1 and return back to t = 1 again, as shown in Figure 5. §III discussed the great benefit of this feature: a large number of SPs can be found by tracing a single curve. In understanding the behavior of the Newton homotopy beyond t = 1, the point of view of catastrophe theory provides an intuitive interpretation. As the value of t exceeds 1, the artificial force applied in (11) reverses direction. It is hence more convenient to use the parametert := t − 1. Also, if we write f = ∇V N Th (a), then (11) becomeŝ We shall focus on the behavior after t = 1, that is, we start witht = 0. Ast increases, this potential represents the application of a continuously increasing artificial force f onto the system of N charges. Here, we considerVt as a family of real-valued functions smoothly parametrized by the single parametert. One basic consequence of catastrophe theory is that while a generic member of this family is expected to be Morse (having only nondegenerate SPs), this family as a whole is likely to contain some isolated non-Morse members. That is, in the process of tracing the moving SP as the artificial force changes, it is likely that some isolated degenerate SP may be encountered. While in general the local structure near a degenerate SP can be complicated, the celebrated classification theorem of elementary catastrophes states that the "generic type" of degenerate SPs can assume only a few standard forms. In this context, the smoothness assumption (Assumption 1) essentially amounts to the assumption that the only type of degenerate SP that can be encountered is the simplest type: "fold catastrophes" which correspond to the "turning points" where the trajectory of the moving SP changes direction int. Figure 5 illustrates a typical curve defined by the Newton homotopy (11) under this smoothness assumption.
Along the trajectory, after the encounter with a turning point (fold catastrophe), thet-value decreases and ift returns back to 0 (i.e. t = 1) then the artificial force is once again removed and the corresponding x will converge to a different SP of the potential V N Th (x). Stated intuitively, the artificial force pulls the cluster beyond its breaking point and a new geometry emerges.

V. CONCLUSIONS
Following our Communication, 23 in this paper we have described a Newton homotopy approach specifically adapted for the potential energy landscapes of atomic clusters. We have shown that compared to the traditional numerical solvers, our method is efficient in finding multiple SPs of any index in the following sense: unlike most other methods, our Newton homotopy implementation finds multiple SPs given a single initial point. Another feature of the homotopy approach that we have exploited in our work is its ability to find degenerate (singular, or non-Morse) SPs. To demonstrate this feature, we first applied our implementation to the celebrated Griewank-Osborne (GO) system, 22 whose only SP is the origin, which is a sin-gular SP of the system. The Newton-Raphson method (and most other numerical approaches) are known to fail to find the SP of this system starting from any point except the SP itself.
In an appendix, we prove that if the system is treated with the Newton homotopy approach, the singular SP of the GO system is straightforwardly obtained.
We have applied our implementation to obtain some new results for the potential energy landscape (PEL) of the Thomson problem. Locating global minima for the Thomson problem has been an active area of research for quite some time. Our goal in the present work is not to obtain exhaustive database for the model, rather it is to give a proof of concept of a method and its implementation for a particular application. Since the number of local minima, and stationary points of higher Hessian index, are expected to increase exponentially with the number of charges, 28,93,94 locating the global minimum becomes increasingly more challenging. However, the rate at which the number of stationary points increases with N is relatively slow for the Thomson problem because the Coulomb potential is longranged, [95][96][97] which effectively smooths the landscape. Nevertheless, locating the true global minimum inevitably becomes more difficult for larger sizes. Our results in the present article suggest that for this model the number of SPs of any index also increases relatively slowly with the number of charges. We also show that in addition to the SPs, the Thomson problem possesses singular SPs, highlighting the richness of the PEL.
From the homotopy continuation point of view, the Thomson problem has gained special attention because it appeared as the 7th problem in Steven Smale's list of eighteen unsolved problems for the 21st century. 69 Leaving the technicalities of the statement of the problem aside, it asks us to algorithmically construct a collection of points whose value of the energy 1 is close (meaning less than O(log(N )) away) to the global minimum value. It has been proved that if a sequence minimizes the Thomson problem, then it asymptotically minimizes the general logarithmic potential, 98 and vice-versa. 99 (see also 100 in conjuction with. 66,101 ) We note that the index of the initial point in the Newton homotopy a as an SP of V 0 seems to determine the index of the first SP the curve obtains, i.e, the index appears to be preserved during the first homotopy run. Since, NH is capable of finding more than one SP along one curve, after the first SP it encounters, the curve may snake back and forth crossing t = 1 several times. In that case, it obtains SPs of alternating parities in terms of the index.
We point out another homotopy based approach called the numerical polynomial homotopy continuation (NPHC) method, which guarantees to find all the complex (including real) SPs for potentials having polynomial-like nonlinearity. 41,42 There, one needs to track as many solution-paths as the estimated upper bound on the number of complex solutions of the system. When the gradient equations of the Thomson problem are transformed to the polynomial form, by clearing out the denominators and using dummy variables to remove the square-root terms, the known upper bounds, such as the Bézout bound, blow up immediately even for small system sizes, making it prohibitively difficult to solve the systems, even though the actual number of real solutions is very small. In such situations, although the Newton homotopy approach is not guaranteed to find all the SPs, it has an advantage in that it quickly finds many SPs, restricting only over the real solutions from the beginning. A more rigorous comparison between the two approaches is underway and will be published elsewhere.
Our implementation of the Newton homotopy method finds SPs of indices close to (2N − 3)/2 when N is moderately large. While this result may stem from the fact that the actual number of SPs at this mid-range of indices is usually the general logarithmic potential of which the Thomson problem is a special case. exponentially greater than the number of SPs of indices near 0 or (2N − 3), in many applications we require SPs of index 0 (minima) and 1 (transition states). An effort to restrict the homotopy search of SPs to a given index is currently underway. In §II, via the construction of theV t , the problem of locating an SP of the target potential V is transformed into the problem of tracing a smooth curve consisting of SPs of the family of potential functionsV t (of which V is a member) in the x-t space. This section briefly outlines the necessary numerical techniques and refers to standard texts 102-104 for comprehensive analysis of this subject. Throughout, we shall assume that the smoothness assumption (Assumption 1) is satisfied. Under this assumption, we have a smooth curve γ passing through the point (a, 0) such that for every point (x, t) ∈ γ, x is an SP of V t . Let H(x, t) := ∇V t (x), then γ is the unique smooth curve consisting of points that satisfy and passes through (a, 0). The goal is to "trace along" γ from the starting point (a, 0) via reliable, efficient, and stable numerical methods. If this smooth curve crosses the plane of R n × R defined by t = 1 at the point (x, 1), then, by continuity, x is an SP of the target potential function V .
There is a wide range of variations on the basic methods. Here we briefly review the method based on Euler's predictor and Generalized Newton's corrector: 105 While in the construction of the deformationV t , it is natural to consider t as the designated parameter, as illustrated by Figure 3 and 4, it is very common (and indeed typical) to have "turning points" along γ where the curve cannot be parametrized by t. Fortunately, the smooth curve γ is always parametrized by arc length. To simplify the notation let y = (x, t) ∈ R n × R and write H(y) = H(x, t). Hence we consider H as a function H : R n+1 → R n where no specific variable is considered as the parameter.
The arc length parameterization of γ is a smooth function y : R + → γ such that y(0) = (a, 0), H(y(s)) = 0, and the tangent vectorẏ always has unit length for all s ∈ R + . These conditions can be characterized by the system of differential equations DH(y(s)) ·ẏ(s) = 0 ẏ(s) = 1 where DH is the n × (n + 1) Jacobian matrix of H : R n+1 → R n . Parameterizations satisfying the above system are clearly not unique: near every point on a smooth curve, there are always two different arc length parameterizations leading in opposite directions. Therefore, to trace along a curve without backtracking, one must be able to determine and maintain a consistent orientation.
As stated in (A2),ẏ(s) is in the null space of the n×(n+1) matrix DH(y(s)), which is of rank n by the smoothness assumption. Therefore, the (n + 1) × (n + 1) square matrix DH(y(s)) y(s) must be nonsingular, and hence its determinant never vanishes for all s ∈ R + . This result, in turn, implies that the sign σ(y) := sgn det DH(y(s)) y(s) never changes along the curve for all s ∈ R + and it can therefore be used to determines the orientation of the parameterization. Once an orientation σ 0 = ±1 is chosen, it must be kept consistent while tracing along the curve to prevent backtracking.
In principle, any ODE solver capable of integrating the above system can be used to trace the curve. However, a "prediction-correction scheme" is generally preferred due to numerical stability concerns. 72,104 In this scheme, one traces along the smooth curve γ via a series of discrete "prediction-correction" steps. In a "prediction" step, starting from a point known to be close to γ, an efficient but potentially less accurate numerical method is used to produce an approximation of a point "one step" further along the curve γ in the sense of the arc-length for some given step size. Such a prediction step is followed by a "correction" step in which Newton-like method is applied to refine the prediction to within a close proximity of the curve γ. One repeats these steps to move forward along γ. While there are many different choices for the "predictors" and the "correctors", the combination of Euler's predictor and the Generalized Newton's corrector provides a simple but effective combination.

Euler's predictor
Starting from a point y 0 ∈ R n+1 on the curve γ, Euler's predictor simply moves one step along the tangent direction of γ at this point for some given step size ∆s. The tangent directionẏ is given by (A2), which can be computed efficiently via the QR-decomposition of DH(y 0 ) since, by construction,ẏ lies in the one-dimensional null space of DH(y 0 ). As a by-product, the QR-decomposition of DH(y 0 ) also allows one to quickly compute the determinant det DH(y 0 ) ∆y and thereby determine the correct sign forẏ.
Withẏ computed, the point y = y 0 + ∆s ·ẏ is produced as the resulting prediction.

Generalized Newton's corrector
The predictionŷ produced by the Euler's predictor is usually inaccurate. If successive prediction steps were to continue using such inaccurate data, numerical errors would quickly accumulate to an unacceptable level. As a result, within the community of numerical homotopy continuation methods, it is a standard practice to have each prediction step followed by a correction step via a series of Newton-like iterations. Note, however, that the defining equation H(y) = 0 of γ has n equations and n + 1 variables (y = (x, t)), so the standard Newton's iteration cannot be used directly. The Generalized Newton's Iteration 105 is therefore required: Define the map N : R n+1 → R n+1 given by N (y) = y − (DH(y)) + H(y) where (DH(y)) + is the Moore-Penrose inverse of the n×(n+1) Jacobian matrix DH(y) of H at the given point y. Starting from the potentially inaccurate predictionŷ produced by the Euler's predictor detailed in the previous section, the map N is iteratively applied toŷ until a given convergence criterion is met. Stated formally, the result of the correction is y = N k (ŷ) = N • N • · · · • N (ŷ) where k ∈ Z + is determined by some convergence criterion.
It is, of course, possible for the Generalized Newton's corrector to converge too slowly or fail to converge at all. This usually means that either the starting pointŷ of the correction is too far from the curve γ, or the defining equation has become ill-conditioned. In both cases, one should repeat the previous prediction step with a smaller step size.
Moreover, for choices (a, b) ∈ R 2 with b > 119a 2 /128, the curve defined by H GO (x, y, t) = 0 starting from (a, b, 0) will always reach the degenerate SP (0, 0) of (4). To verify this result, we can first solve for y in the second equation in the system H GO (x, y, t) = 0 and substitute it into the first equation, resulting in a cubic polynomial equation in x with coefficients involving a,b, and t. Its discriminant ∆ is We are interested in the sign of ∆ as t changes from 1 to 0. Substituting b = r119 128 a 2 , ∆ becomes Note that 2 20 ∆/(1 − t) 2 a 6 is a linear function in t. Straightforward calculation reveals that this function takes negative values at t = 0 and t = 1 for any r > 0. Therefore ∆ < 0 for all t ∈ [0, 1] and r > 0. Recall that b = 119+r 128 a 2 , so for b > 119a 2 /128 and t ∈ [0, 1], the discriminant ∆ is negative indicating that the corresponding polynomial equation has a unique real root for each fixed t. Together with the smoothness assumption, these results imply that the equation H GO (x, y, t) = (0, 0), defines a single smooth curve in R 2 × [0, 1] which necessarily connects the starting point and the target solution (0, 0) of the system (5).