{'fraction_non_alphanumeric': 0.04097831102907245, 'fraction_numerical': 0.013474850023073373, 'mean_word_length': 4.507242693773825, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Construction of a pathway map on a complicated energy landscape

How do we search for the entire family tree without unwanted random guesses, starting from a high-index and high-energy stationary state on the energy landscape? Here we introduce a general numerical method that constructs the pathway map clearly connecting all stationary points branched from a single parent state. The map guides our understanding of how a physical system moves on the energy landscape. In particular, the method allows us to identify the transition state between energy minima and the energy barrier associated with such a state. As an example, we solve the Landau-de Gennes energy incorporating the Dirichlet boundary conditions to model a liquid crystal confined in square box; we illustrate the basic concepts by examining the multiple stationary solutions and the connected pathway maps of the model.

How do we search for the entire family tree without unwanted random guesses, starting from a high-index and high-energy stationary state on the energy landscape? Here we introduce a general numerical method that constructs the pathway map clearly connecting all stationary points branched from a single parent state. The map guides our understanding of how a physical system moves on the energy landscape. In particular, the method allows us to identify the transition state between energy minima and the energy barrier associated with such a state. As an example, we solve the Landaude Gennes energy incorporating the Dirichlet boundary conditions to model a liquid crystal confined in square box; we illustrate the basic concepts by examining the multiple stationary solutions and the connected pathway maps of the model.
Introduction.-Many problems in natural sciences and engineering share a common theme: the energy (freeenergy, or target-function) landscape as a function of the physical variables in the system displays a multitude of minima [1]. The physical variables here could be the amino-acid locations in the protein folding problem [2][3][4][5][6], the atomic positions of an atomic cluster (such as the Lennard-Jones clusters) [7,8], the continuum field variables in AB-diblock copolymer self assembly [9][10][11][12][13], or even network parameters in a to-be-optimized artificial neural networks [14]. One often makes a comparison of a two-variable (2D) problem with the geometric feature on a landscape: minima are basins, maxima are summits, and the two variables are the coordinates on a plane [see Fig. 1(a)]; of great interest is the locations of minima that lie somewhere on the landscape.
Equally important is the possible kinetic pathways between these energy minima, when a system has enough thermal excitations to overcome the energy barriers between them. The minimal energy barrier that partitions two minima is associated with at least one index-1 saddle point, which is a transition state that the system likely experiences when it moves from one minimum to another [15]. Here the Morse index, m, characterizes the nature of a saddle point: it is an energy maximum in m orthonormal directions, and minimum in the rest [16]. The D-variable energy surface possesses a complicated landscape in high-D; recent two decades have witnessed the progress in developing mathematical procedures to determinate and understand saddle points numerically [17][18][19][20].
For example, in Fig. 1(a), the lowest energy barrier between two minima C 1 and C 2 is located at an index-1 saddle B 1 , not the index-1 saddle B 2 . In replacement of the energy landscape, Fig. 1(a), which is impossible to draw in high-D, here we introduce the more informative concept of a pathway map, Fig. 1(b), which de- picts how lower-index saddle points are connected to the higher ones, and ordered according to the values of the energy surface at the saddle points. With the energy values specified vertically, such a pathway map is more useful than simply connecting minima through index-1 saddles [21]. In a high-D system, the pathway map can be constructed from a parent state (a high-index saddle) all the way down to minima (index-0 solutions). It is hence desirable to design a numerically tractable algorithm that constructs a complete pathway map, which guides us on how the system runs into metastable states and the energetic requirement to overcome energy barriers when the system moves between the states. Assume that the energy surface E(x) is a function of a D-dimensional variable x. All stationary points satisfy the nonlinear equations ∇E(x) = 0, and the Morse indices are determined by examining the Hessian ∇∇E(x) at these points [22]. While these mathematical definitions are elementary, capturing all stationary points is practically difficult in computation. Many algorithms, e.g., the minimax method [23,24] or deflation technique [25,26], are dedicated to solving the nonlinear equations, arXiv:1909.13182v1 [physics.comp-ph] 29 Sep 2019 usually relying on an initial guess which deterministically leads to a stationary point. However, as more solutions are discovered, it becomes increasingly difficult to propose and fine-tune initial guesses, to search for remaining solutions. A large number of optimization methods are also developed to directly find local minima; it is often computationally expensive, as the total number of local minima usually grows exponentially in a high-D problem [1]. Therefore, even from a computational viewpoint, it is essential to establish a numerical procedure to discover the pathway map that systematically connects all stationary points, rooted from a high-index saddle as a parent and descending to minima.
Downward Search.-The question then becomes: how do we establish a pathway map, given a parent saddle of E(x)? Here we introduce an efficient computational method that enables the search for such a map starting from an index-m saddle. The assumption is that the location of the saddle point is at x * and that the m normalized vectors v * i corresponding to the m negative eigenvalues of the Hessian matrix are known at this stage, where i = 1, 2, 3, ..., m.
The essential idea is to choose x(0) = x * ± u as the initial search position for a lower index-k (k < m) saddle with a small that "pushes" the system away from the index-m saddle. The direction of pushing, u, is along a linear combination of (m − k) vectors (whose negative eigenvalues have the smallest magnitudes) selected from the set of {v * }. The other k orthonormal vectors from the set are used as the initial unstable directions for the next index-k saddle. This way, the system can gently roll off the original index-m saddle point in the unstable direction within a controlled procedure. Normally, a pair of index-k saddles can be found this way, corresponding to the ± sign in the initial guess.
The details of this numerical procedure, the so-called k-HiOSD method [27], is documented in the Supplementary Material [28], step-by-step. It computationally determines an index-k saddle point based on x(0) and v i (0) (i = 1, 2, 3, ..., k) as an initial input. By relaxing a pseudo-Langevin dynamics for time-dependent x and v i , the computation eventually arrives at all conditions that need to be satisfied by an index-k saddle, producing x * and v * i for it. Hence, by repeating the above procedure starting from an existing index-m saddle as a parent, we establish a systematic scheme to search for all saddle points branched from this parent and to construct a family tree that eventually terminates at index-0 solutions (minima). This avoids the pitfalls existing in previous numerical schemes to search for an arbitrary solution starting from random initial guesses with no control on both search direction and Morse index of the final numerical outcome. The family tree provides the crucial information on the relationships between saddle points. In particular, the pathway between two lower-index saddle points is clearly mapped out through a higher-index saddle point.
Liquid crystals in 2D confinement-To illustrate the basic concepts, we present here the pathway map of 2D nematic liquid crystals confined in a square box (see illustrations in Fig. 2). The mathematical model we use is the Landau-de Gennes (LdG) free-energy in reduced units [29] E[Q(r)] = dr where a 2 × 2 traceless and symmetric tensor field Q(r), characterizing the orientational ordering of the liquid crystals, is considered. Only two elements of the Q-tensor are independent. The variable r = (x, y) resides in the domains [−1, 1] for x and [−1, 1] for y. The specifics of the LdG model, including the Landau function f b and Dirichlet conditions that require the alignment of the liquid crystal directors at the four boundaries, are given in the Supplementary Material [28]. The system parameter α is directly linked to the area of the confining square. The energy functional in (1) must be minimized with respect to the tensor field Q(r) for a given α. There are two equivalent approaches to finding a stationary solution. The first is to take and solve the Euler-Lagrange equation associated with the unknown field Q(r) [30,31]; the second, closer to the spirit presented here, is to discretize the (x, y) region by a grid system containing G grid points and to treat the two independent elements of Q at these points as individual variables; the stationary solution of the energy in (1) is then obtained simply by equating the gradient of E in the D-dimensional space to zero, where D = 2G is the total number of individual variables [32,33].
When α is sufficiently small, it was shown that the socalled well order-reconstruction solution (WORS) is the only stationary solution and corresponds to the global minimum of (1) [34]. Displayed in Fig. 2(a) is such a state with index-0 (minimum). WORS was suggested as a possible stationary state based on the Onsager model in Ref. [35] and was discovered numerically first based on the LdG model in Ref. [36].
One can show that WORS remains a stationary point for all α [34], but the Morse index of WORS changes [31]. Numerically, by using WORS at a smaller α as an initial guess, we can always obtain WORS at a larger α. Table I      longer a stable state and the system begins to descent to other stable states.
Using our k-HiOSD method, we can exhaustively enumerate all possible states starting from WORS. Figure 2 presents the pathways for a number of typical α, where the defect states specified in the nodes are illustrated on the right panel. Take Fig. 2(c) for example. At α = 15, WORS is an index-2 saddle, which bifurcates to two boundary distortion (BD) states (related by π/2 coordinate rotation) that are index-1 saddles. Both BD states then further relax to two diagonal (D) states (again related by π/2 coordinate rotation) as the global minima. This can be compared to the pathway map in (b) when α = 7, where the index-1 WORS is the direct pattern from which the two D states bifurcate. Assume that we wish to design a physical device to take advantage of the different optical properties of the dual D states by switching between them [37]; the energy barrier that the device has to overcome and the intermediate states between the two D states are different for system size α = 7 (b) and α = 15 (c); the former goes through WORS and the latter through BD.
As α increases further, the pathway map becomes more complicated. At α = 30, WORS is now an index-4 saddle. Two types of unrelated defect states, D and rotated (R) (and their related counterparts by π/2 rotation), emerge as index-0 minima. A different lesson from the above is learnt; the switching between the dual D states (D 1 and D 2 ) must go through, for example, the pathway sequence D 1 →J ne →R e →J se →D 2 , or vice versa, where the subscripts indicate the rotated orientations of the pattern (e=east, w=west, n=north, s=south). The system could be trapped in R as a metastable state because it is an energy minimum. In order to move to D where the global minimum resides, the R state needs to overcome an energy barrier E J − E R , which is available from our calculation. From the Kramers theory [38], we can then estimate the trapping time to find out the relaxation kinetics of the bistable system [39].
The relationships between the stationary states could become quite complicated at α = 50, as Fig. 2(f) demonstrates. WORS is now index-6. Starting from WORS as the parent, we search for other states generation-bygeneration to produce the entire family tree (see an illustration of the dynamic downward pathway sequence: WORS→BD→BD+→I→J−→D in Supplementary Material [28]). In total, 89 distinct solutions can be found and after classification to account for basic symmetry operations, they fit into 17 basic types. Surprisingly, the structure that contains more defect features, J± (the + and − signs of J± represent this state accommodates +1/2 defect point and −1/2 defect point inside, respectively), undertakes the simpler J as the transition state that partitions a D state (global minimum) and an R state (local minimum). The switching between the two ground states D 1 and D 2 must go through pathway sequences D J± R J± D as illustrated in the inset where the orientations of the involved states decide the exact pathway.
One of the technical advantages of our method is to produce the entire family tree under a parent. A direct result is our finding of states N±, M±, S±, and T±, which were not reported in the previous numerical studies of the same LdG model [30] or the Onsager model [40]. Previously, the variety of defect states of these models were obtained by simply solving the necessary condition that the gradient of the energy must be zero, without the systematic inquiry of the sufficient condition of whether the found states are high-index saddle points. The conclusion we draw for the LdG model, that WORS is stable at small α and D is stable at large α with R as the possible metastable state, is consistent with the phase diagram concluded in Ref. [40] based on the Onsager model. All these can be further compared with recent experimental observations where no patterns having high Morse index were actually stabilized [37,[41][42][43][44].
Upward Search.-The downward search strategy guarantees the systematic finding of energy minima as index-0 solutions starting from a given parent, which is already superior than most other optimization numerical methods. On the other hand, if multiple parent states exist, we need to conduct the search on each family tree starting from an individual parent state. For example, starting from A 2 in Fig. 1(a), a family tree different from the one in Fig. 1(b) would be found. The A 1 and A 2 trees may share a common B 1 node but the downward searches establish separated maps.
The key to find complete, interwound pathway maps is to conduct a systematic upward search for the parent states. The k-HiOSD algorithm already embeds a mechanism that allows us to search upward starting from an index-m saddle point to find an index-k saddle where k > m. Again, we assume that all eigenvalues and eigenvectors of the Hessian matrix are known. To do so, the initial condition x(0) = x * ± u perturbs x * along the direction u, which is selected along a linearly combination of (k − m) eigenvectors corresponding to small positive eigenvalues. The initial {v * } set used in the search includes these (k − m) eigenvectors and the original m eigenvectors of the negative eigenvalues. By doing so, the search may branch into other family trees through, e.g. the shared B 1 in Fig. 1(a) to end up at parent A 2 . Then a downward search from A 2 covers the entire pathway map different from the one illustrated in Fig. 1(b).
Take the scenario where D is a known, initial state without the knowledge of existence of other structures in Fig. 2(f). Along the upward pathway indicated by dashed arrows, in a single search we can produce the upward pathway sequence D→J+→H→M+→N+→WORS (see Supplementary Material [28]). Once WORS is arrived at, no other upward searches give rise to a solution with a higher index hence WORS is found as a parent.
Summary.-This Letter introduces the concept of pathway maps starting from parent states, which illustrate the relationships between the stationary states and classifies these states using Morse index. The enabler for such construction is our k-HiOSD numerical algorithm. Using the LdG model for confined liquid crystals as an example, we illustrate the usefulness of such a pathway map in understanding the physical properties of a multisolution problem, which are otherwise unobtainable by other methods.
The pathway map starts with a parent structure (nonzero Morse index) and then relates the entire family completely down to energy minima (zero Morse index). Most importantly, the connection between the energy minima and the transition states are clearly shown on the map, which guides our understanding of a physical system when it moves on the energy landscape. The completeness of the family tree is the direct result of the tight control in using k-HiOSD. Lower-index solutions are produced from a high-index solution along the unstable directions in high dimensions, systematically. No random initial guesses are made in the numerical search. The method also avoids the time-consuming calculation of the eigenproblem of the Hessian matrix hence bypasses the unwanted high computation complicity. The stationary state is found together with an assessment of its Morse index.
The k-HiOSD algorithm can also be used for upward search with a selected direction so that the entire search can navigate up and down on the energy landscape. This gives a potential method for complete determination of all energy minima (hence the global energy minimum) if the number of stationary states is finite. The tight control of the initial conditions avoids unwanted random searches hence saves computational time.
We thank Prof. Yucheng Hu for helpful discussions. This work was supported by the National Natural Science Foundation of China (Grants No. 11622102, No. 11861130351, No. 11421110001, No. 11421101 and No. 21873009). J. Yin acknowledges the support from the Elite Program of Computational and Applied Mathematics for PhD Candidates of Peking University.

I. K-HIOSD METHOD
Here, we explain the essential steps in the high-index optimization-based shrinking dimer method for an indexk saddle point (k-HiOSD method) to find an index-k saddle point based on an initial guess [1]. The energy surface E = E(x) is assumed at least twice-differentiable and x is a D-dimensional vector variable.
An index-k saddle point at x * has the following properties. It is a local maximum in k orthonormal directions and minimum in (D − k) orthonormal directions in the D-dimensional space. These can be understood from the Hessian in a matrix form with the elements where i, j = 1, 2, 3, ..., D, or in the tensor form, It has D eigenvalues: (D − k) positive and k negative.
Let v * i be the normalized eigenvector corresponding to the ith smallest eigenvalue, where i = 1, .., k.
The k-HiOSD method locates such a saddle point by solving the following dynamic problem from an initial guess, where both x and v i are treated as functions of time t so that their time-derivativesẋ andv i are considered. The pseudo-dynamics that enables the numerical search for an index-k saddle point is assumed to satisfy with the orthonormal constraints v i · v j = δ ij . Note that I is a D-rank unit tensor. A line-derivative is numerically taken on the vector field ∇E(x) to represent the projection of the Hessian in direction v, where h is a small number. To produce the example in the text, we used h = 10 −6 . The relaxation parameters, β > 0 and γ > 0, guide the rate of the dynamics and are usually taken to be small parameters, according to the actual problem. For example, in the text, when we deal with the liquid-crystal confinement problem, we let β = 5 × 10 −5 and γ = 10 −5 . To complete the differential equations, the initial conditions (initial guess), x(t = 0), v i (t = 0), are required. The time evolution of the involved quantities can be numerically treated by discretization the t variable of Eq. (3), to be used in the explicit Euler scheme, the Barzilai-Borwein gradient method [2], or any other numerical schemes. One can show that once the k-HiOSD dynamics stably converges to a final solution x * and v * i , i.e., when the left hand sides of the equation set completely vanish to zero, all properties of an index-k saddle point are satisfied [1]. Note that no Hessian matrix and no eigenvalue calculation of it are explicitly needed in this method. The algorithm requires an initial guess on x(t = 0) and v i (t = 0) as the input and yields x * and v * i as the output. A converged set v * i guarantees that their eigenvalues to be negative.
On the other hand, if the above pseudo-dynamics does not converge, then an index-k saddle point does not exist starting from the initial search position x(t = 0). In this case, either a new guess is made or another k is attempted.

II. LANDAU-DE GENNES MODEL FOR CONFINED LIQUID CRYSTALS IN 2D
One of the best known models to describe the liquidcrystal ordering is the Landau-de Gennes (LdG) model, which treats the free energy as a functional of the Qtensor field [3,4]. Minimizing the free energy, one finds the stable liquid-crystal configurations of the system. In the two-dimensional example used in the text, Q is taken to be a 2 × 2 traceless and symmetric tensor [5,6], Q = Q 11 (r) Q 12 (r) Q 12 (r) −Q 11 (r) with scalar components Q 11 (r) and Q 12 (r