Exploring Critical Points of Energy Landscapes: From Low-Dimensional Examples to Phase Field Crystal PDEs
Priya Subramanian, Ioannis G. Kevrekidis, Panayotis G. Kevrekidis
EExploring Critical Points of Energy Landscapes: From Low-Dimensional Examples toPhase Field Crystal PDEs
P. Subramanian
Mathematical Institute, University of Oxford, Oxford OX2 6GG, UK
I.G. Kevrekidis
Department of Chemical and Biomolecular Engineering,Johns Hopkins University Whiting School of Engineering andDepartment of Applied Mathematics and Statistics,Johns Hopkins University Whiting School of Engineering
P. G. Kevrekidis
Department of Mathematics and Statistics, University of Massachusetts, Amherst MA 01003-4515, USA andMathematical Institute, University of Oxford, Oxford, OX2 6GG, UK
In the present work we explore the application of a few root-finding methods to a series of pro-totypical examples. The methods we consider include: (a) the so-called continuous-time Nesterov(CTN) flow method; (b) a variant thereof referred to as the squared-operator method (SOM); and(c) the the joint action of each of the above two methods with the so-called deflation method. More“traditional” methods such as Newton’s method (and its variant with deflation) are also broughtto bear. Our toy examples start with a naive one degree-of-freedom (dof) system to provide thelay of the land. Subsequently, we turn to a 2-dof system that is motivated by the reduction ofan infinite-dimensional, phase field crystal (PFC) model of soft matter crystallisation. Once thelandscape of the 2-dof system has been elucidated, we turn to the full PDE model and illustratehow the insights of the low-dimensional examples lead to novel solutions at the PDE level that areof relevance and interest to the full framework of soft matter crystallization.
I. INTRODUCTION AND MOTIVATION
The topic of root finding at the level of stationary states of ordinary and partial differential equations is oneof fairly universal appeal . Especially in the realm of nonlinear waves and pattern formation, a wide variety ofmethods has been proposed including in the recent past , many of which have been summarized in . These typesof methods enable the identification of steady patterns, which are subsequently subject to the evaluation of theirrespective stability properties. When the associated patterns are found to be stable, they suggest natural attractorsof the dynamics (in the case of dissipative dynamical systems) or centers around which the dynamics may oscillate(in conservative dynamical systems). On the other hand, even unstable dynamical solutions (potential maxima orsaddles of the landscape of the system) are interesting in their own right towards characterizing transient dynamicsand transitions between different states.In the past few years, there has been a number of proposals for variations to the well-known gradient flow methodsin order to take advantage of the acceleration that may be provided by the so-called Nesterov method . While theoriginal method was proposed for discrete systems, it has been recently adapted to continuum ones both at the ODElevel and at the PDE level . It has also been shown to be suitable to frame in the realm of Lagrangian and actionextremization principles and has proved to be an efficient method for identifying steady states of even complexmulti-component PDE systems (in connection to spinor atomic Bose-Einstein condensates) in .Our aim in the present work is to expand upon the relevant ideas, using a few simple examples which will provide auseful understanding about how these methods work, especially so in a progressively more complex landscape startingwith 1-dof, progressively moving to 2-dof’s and ultimately connecting the latter with a full PDE. That is to say, ouraim is to explore the relevance of the methods in well-controlled ODE and, subsequently, in (connected to the ODEthrough suitable reductions) PDE examples as a means of identifying some of the advantages and potential caveatsof different methods and as a way to guide the selection of which tool is most well-suited to which example.Finding equilibria of a system of coupled ODEs is simplified if the corresponding conditions for attaining equilibriumreduce to a set of coupled polynomials/Laurent polynomials. In that case, multiple iterative local methods such asNewton’s method or non-iterative global methods can be employed. Local methods need a good initial guess and amethod to either determine the Jacobian or its action on a vector. Alternatively, global methods such as homotopyuse ideas from numerical algebraic geometry to obtain all non-singular solutions (both real and complex) to the setof coupled polynomial equations without any initial guess . Finding equilibria of the PDE system using iterativemethods is subject to the same constraints as described above in the discussion of ODE systems, in that there is a r X i v : . [ n li n . PS ] A ug a need for a “good” initial guess and the need to approximate the action of the Jacobian. One way to generatesuitable initial guesses is to (i) use asymptotic states obtained from the time evolution of the full PDE model, or(ii) generate a reduced order ODE model that reproduces the dynamics of the PDE model in the neighbourhood ofa critical point. Such a reduction is at the heart of a wide variety of techniques at the level of perturbation theory,effective Lagrangians/Hamiltonians, variational methods and many others; see, e.g., the reviews of for just a fewexamples. This second approach utilises appropriate rescaling of the variable(s), key parameters and time in the PDEmodel, to derive effective equations chracterizing the transition at the critical point. Once we obtain a reduced set ofODEs, the different ODE root-finding methods described above can be deployed to obtain the associated equilibriafor the ODE system and then to transcribe those towards initial guesses of the relevant PDE.In the exposition below we start with a very simple 1-dof example where the solution landscape is fully-known, asthe roots are transparent from the very form of the ODE’s right hand side (RHS). This allows us to easily transcribethe CTN method in this context and observe its potential convergence to different types (single, and double) roots.This also simplifies the application of modified methods, such as the method that we will refer to as the square-operator method (SOM) which does effectively use the Jacobian, retaining the form of a flow method but in a newlandscape where the original roots are now landscape minima (and hence each one of them is locally attractive –orfor conservative dynamics, a center–). An additional widely used methodology in the form of the so-called deflationmethod is also brought to bear. This method focuses on eliminating the possibility of the flow to converge to aparticular root. As a byproduct, the use of root finding techniques (such as Newton’s method), in conjunction withthe repeatedly-modified landscape induced by sequentially “removing” from the pool of available attractors differentroots, can offer convergence to numerous distinct roots. We explore how deflation works both with flow methods andwith Newton type methods in our simple examples.Indeed, in all of these methods, we go beyond the informative but rather limited 1-dof example and explore anequally tractable, but more complex 2-dof example. This setting has emerged from a 2-dof reduction (near a suitablecritical point as described above) of an infinite dimensional model of the PFC variety exploring the crystallization ofsoft matter . For the case considered, we find a significant wealth of possible stationary states (9 in total), coveringdifferent types of states: energy minima, energy maxima and saddle points. While naturally root-finding techniquescan identify all 9 of these roots, here we explore the behavior of different types of flow methods, such as CTN, SOM,deflation-based CTN, etc. We examine what converges where, what happens to the landscape under deflation, whatare the basins of attraction of each root, etc. All of these questions reveal the strengths and weaknesses of the differentmethods. Arguably, even more importantly, these set the stage for reconstructing, on the basis of the 2-dof ODEsteady state, viable PDE solutions of the original system. We then explain the effort to capture the correspondingPDE steady states, the comparison of the stability properties of the PDE vs. those of the ODEs and the potentialof the instabilities of the states obtained (based on the reduction) to yield new asymmetric states (via the dynamicalevolution), of not only quasi-crystalline but also other types.Our presentation will be structured as follows: in section II, we will kick off by explaining the different methods inthe simplest 1-dof example. Then, in section III, we will expand considerations to the example at the center of ouranalysis, namely the PFC model of crystallisation in soft matter. The first part of the relevant section will focus onthe more definitive 2-dof considerations: upon performing the relevant reduction of the PDE to the 2-dof model, wewill explore the latter for specific parameter sets at some length. Then, we will return to the PDE to examine whatsolutions we can obtain from the evolution of the reconstructed (based on the 2-dof) PDE solutions. A broad rangeof not only quasi-crystalline, but also rhombic, hexagonal and other asymmetric patterns will be seen to emerge fromthis effort that are of intrinsic interest in their own right for the original PFC system. Finally, in section IV we willsummarize our findings and offer some possibilities for extensions of the present work. II. ONE-DIMENSIONAL EXAMPLE
We consider the dynamics of a system defined as˙ x = f ( x ) = − ( x − x − ( x −
6) = −∇ V ( x ) . (1)Here the overdot indicates a derivative with respect to time and the corresponding energy landscape can then bedefined as V ( x ) = − (cid:82) f ( x ) dx . In order to determine the solution landscape in this simple example, the first methodthat we use is the continuous time Nesterov (CTN) method . In this method, we consider the following ODE (i.e.,the flow) in order to solve the problem f ( x ) = 0, ¨ x + 3 t ˙ x = f ( x ) . (2) FIG. 1: Left: Energy landscape for the one-dimensional example. Right: Continuous time Nesterov (CTN) iteration. Blacklines in both panels is the original system with roots (1,3,3,6). Blue lines indicate the energy landscape and iterations afterremoving one root, i.e., the minimum at x ∗ = 6. This leads to a new minimum being x ∗ = 1. Green line is after removing x ∗ = 1 , x ∗ = 3. One can detect two principal features that distinguish CTN dynamics from the standard gradient descent flow:˙ x = f ( x ). The dynamical system associated with CTN dynamics involves a second order and also non-autonomousODE. As a result, the acceleration vector, rather than the velocity vector, is proportional to f (at least for large t ).Considering a mechanical analogy, the relevant (fictitious) particle x has been given a mass and a time-dependentdissipation, which has the unusual form: (3 /t ) ˙ x . This effective damping on the fictitious particle is tuned to be largeat the initial time, i.e., far from the equilibrium, while it becomes weaker as one (hopefully) approaches the relevantfixed point. This term is therefore responsible for the convergence towards one of the roots of f .We now turn to the specific example of interest. The function f ( x ) has three roots (one of which is double): asingle root x = 1, a double root at x = 3 and another single root at x = 6. Here the equilibrium at x = 6 is the globalminimum of the associated potential energy landscape V ( x ), x = 1 is the global maximum while the double root at x = 3 corresponds to a semi-stable point as shown in Figure 1 left panel. Here the energy landscape V ( x ) is plottedas a function of x in black. The equilibria are identified with red crosses. The right panel of the figure shows thetime evolution of the associated CTN flow. Starting from an initial condition, say, x = 4, we see that the evolution(in black line in the right panel, pink overlay on the black line in the left panel) approaches the energy landscapeminimum at x ∗ = 6. In order to approach any other equilibrium, we need to modify the system. If we explicitlyfactor the known solution out of the system, i.e., define a new system f new = f / ( x − dramatic change, more dramatic in a sense than the deflationconsiderations that will be given below, among other things because this changes the asymptotic behavior of the fieldat large x ; it also alters the stability of x = 1 which is now no longer a maximum, but rather a minimum of thedynamics. This implies that an evolution starting from the same initial condition of x = 4 will now approach thevalue of x ∗ = 1 which is the new minimum. It is interesting to highlight the difference between CTN and a gradientdescent here too. The latter, starting at x = 4 would get attracted to the semi-stable point x ∗ = 3 without being ableto reach x ∗ = 1. Yet, the CTN enables an oscillatory motion (through its acceleration term) which, in turn, allowsthe system to bypass the semi-stable point and reach the, now stable, x ∗ = 1.We can carry this further and remove both roots x ∗ = 1 and x ∗ = 6 in the above rather crude way. Then wefind that the CTN evolution is unable to asymptote to any solution (not shown here). Once again the oscillatorycomponent of the method reveals the effective metastability of the root, rendering the method unable to converge toit. When the double root nature of this point is taken into account (i.e., one factor of x − x ∗ = 3, which of course has by now become a stable fixed point of thehighly modified dynamics.From the above, it becomes clear that CTN can only converge to the stable equilibria of the flow in this one-dimensional example. In order to converge to unstable equilibria, a drastic modification of the relevant landscape(like the division by the monomial factor involving the already identified root, here the x −
6) is needed. Thusrendering x = 1 a minimum, it is also possible converge to the latter root. Semi-stable fixed points appear not to bewell-suited to CTN given its oscillatory nature and, thus, its exploration of the unstable direction of the semi-stableequilibrium. A natural query that arises is whether one can modify the method to be able to converge (for suitableinitial guesses) to all equilibria. To address this issue, we now discuss the squared operator method (SOM).A modification that allows us to convert all the equilibria of a given system f ( x ) into minima is the SOM. Here we FIG. 2: Energy landscape for the one-dimensional example and the related SOM with deflation. The purple curve is the energyassociated with the squared-operator system with deflation, i.e., G ( x ), where the potential energy obtained by the (negative)anti-derivative thereof. Red squares indicate the equilibria for the squared-operator dynamics. The black line is the energyassociated with f ( x ) (and its associated V eff ) shown here multiplied by a factor of 100 (to bring it to the same scale as theother curve) where red crosses are the equilibria. define a new system g ( x ) such that g ( x ) = − f (cid:48) ( x ) f ( x ) . (3)All equilibria of the f ( x ) system are now minima in this new system g ( x ) (as g (cid:48) ( x ∗ ) < f ( x ∗ ) = 0). Defining thecorresponding g ( x ) for the system defined in Eq. (1), we calculate the associated energy landscape. In this squaredoperator system we find that there are 5 equilibria at x = 1 , . , , . ,
6, with the minima at x = 1 , , x = 1 . , .
15. Note that the maxima are spurious in connection to the original root finding problem. They are“topologically” generated separatrices for the dynamics of the squared operator, created by the fact that we are nowminimizing an effective potential energy V eff = f ( x ) . In this setting, the roots that we previously obtained (andonly the roots of the original problem) are the global minima of the novel potential energy landscape V eff . Hence,between these minima, it is topologically necessary that (at nonvanishing values of V eff and thus of f ), maxima of V eff will arise. These maxima now naturally also partition the space between the corresponding minima. In particular,initial conditions in the range −∞ < x < . x ∗ = 1, 1 . < x < .
15 reach x ∗ = 3 and x > .
15 reach x ∗ = 6.These ranges are indeed the one-dimensional basins of attractions of each of the asymptotic minima in this system.However, in addition to obtaining different roots starting with initial data in different basin boundaries, an arguablymore desirable scenario would be to obtain ideally all roots starting from a single initial guess in the spirit of deflationmethods. Thus at the next step, we attempt to combine CTN and SOM with deflation methods. Below, we will showthe case example of (the better behaved, between the two, as we saw above) SOM to illustrate the “pathologies”present. We have identified similar pathologies in the CTN case.In order to explore the relevant technique (SOM with deflation), having obtained an equilibrium of the system, sayat x ∗ = 6, we deflate this equilibrium to define a new system G ( x ) as below, G ( x ) = g ( x ) × (cid:18) || x − x ∗ || p + α (cid:19) . (4)Generally we use a factor with squared norm, i.e., p = 2 and a (positive) shift operator α = 1 in the above factor.The corresponding energy landscape is shown in the Figure 2 in purple, with the equilibria of this deflated squaredoperator system shown as red squares. We also show an exaggerated version of the un-deflated operator in black(multiplied by a factor of 100) for comparison. Starting the time-marching from the initial value of x = 5 .
1, thesquared operator dynamics reaches the equilibrium at x ∗ = 3. If however, we start the time-marching from the samepoint, e.g., x = 5 .
2, the evolution still attempts to reach x ∗ = 6, since the deflated squared operator in purple hascontinuously decreasing value for 5 . < x <
6. Normally the stopping criterion in a CTN method would be to requireboth the velocity ˙ x and acceleration ¨ x to vanish. In this case, we observe that as the evolution approaches x ∗ = 6,the slope of the potential is non-zero. To recognise that we have reached an equilibrium, a modified stopping criterionin this case would be to check if g ( x current ) < e −
10. This makes the evolution stop when we reach an equilibrium,even if the velocity ˙ x , in the CTN formulation is non-zero. FIG. 3: Newton with deflation for the system governed by f ( x ) plotted as a function of the index of the Newton iteration N and the current Newton iterate x N . Starting from the same initial guess of x = 1 .
6, we first converge to x = 1 (black line withmarkers). We this solution deflated, we approach x = 3 (blue line with markers). When both these solutions are deflated, werecover x = 6 solution (red line with markers). Such a modification of the stopping criteria is necessitated when combining deflation with a flow method, evenwell-behaved ones as in the SOM case utilized here. The reason is simple, a posteriori. Deflation eradicates one ofthe roots of the original landscape. In its stead, it leaves behind a pole (typically). This pole may form an infinitelydeep potential well which, in turn, will render it impossible to get unstuck from the flow towards what used to be anequilibrium but is now merely a pole. It may be interesting to consider infinity-crossing methods in this regard (see fora recent example ). Nevertheless, the generic feature remains: the depth of the landscape around a former minimummay change, e.g., with its becoming deeper in the form of a pole. Yet, in the dynamics of the flow considered, thepole continues to attract our original initial condition and thus does not allow it to to converge to another root as isthe goal of the deflation technique.A final technique to employ is the well established combination of Newton iteration with deflation, for which weshow the results in Figure 3. Here, the value of the Newton iterate x N is shown as a function of the number of iteration N . Starting from an initial choice of x = 1 . f ( x ), we obtain the equilibrium at x ∗ = 1 (black curve with circle markers). When this equilibrium of the f ( x ) dynamics is deflated in the same fashionas described in Eqn. (4), the Newton iteration is able to converge to another equilibrium at x ∗ = 3 (blue curve withcircle markers). Deflating both of these equilibria, allows the Newton iteration to approach the third equilibrium at x ∗ = 6 (red curve with circle markers). This confirms that the combination of a Newton iteration with deflation yieldsfar more promising results. For suitable initial guesses, the algorithm can escape from the infinite negative (effective)potential energy pole formed at x ∗ = 1 under deflation and “jump” outside of the relevant effective potential wellthrough its discrete steps. Such a possibility enables it to retrieve additional equilibria and ultimately in this exampleobtain all the roots of the problem. Indeed, this does raise the question of whether deflation should be able to obtainall roots of a nonlinear problem in this way starting from a single initial guess. This is a relevant question to whichwe will return in our 2-dof example that follows below. III. PHASE FIELD CRYSTAL (PFC) MODEL FOR CRYSTALLISATION OF SOFT MATTER
The PFC model describes the density distribution of soft matter forming a liquid or solid on the microscopic scaleand captures the associated time evolution over diffusive time scales. The scalar field U ( x , t ) describes the deviation ofdensity at point x and time t about the average value. Resulting from the conservation of mass across the liquid/solidphase transition, the evolution is governed by a partial differential equation that evolves only to steady asymptoticstates. The amount of useful work that can be extracted from this system being held at a constant volume andtemperature, is the associated Helmholtz free energy ( V ) in the system. In forming the free energy V , we choose alinear operator ( L ) and simple polynomial nonlinearities to include for nonlinear saturation (the − U term) and tobreak the U → − U symmetry (the U term). Then, the free energy V can be written in the form, V = − (cid:90) (cid:20) U L U + Q U − U (cid:21) d x . (5)Here, the choice of L is done with the intention to promote two different length scales with wavenumbers, k = q and k = 1 independently. The choice of the linear operator L allows for marginal instability at two wave numbers k = 1 and k = q , with the growth rates of the two length scales determined by two independent parameters µ and ν ,respectively. The resulting growth rate σ ( k ) of a mode with wave number k is given by a eighth-order polynomial: σ ( k ) = k [ µA ( k ) + νB ( k )] q (1 − q ) + σ q (1 − k ) ( q − k ) , (6)where A ( k ) = [ k ( q − − q + 4]( q − k ) q and B ( k ) = [ k (3 q −
1) + 2 q − q ](1 − k ) , which is also shown inFigure 4. As shown in the inset, the peak height at the two chosen wave-lengths is related to the growth rates, whichin this case are µ = ν = 0 .
1. In the rest of this paper, we analyse this case for concreteness. Similar analyses can beperformed for other choices of the parameters.
FIG. 4: Dispersion relation σ ( k ) for the linear operator as a function of the wavelength k . Two peaks are visible at k = q =1 / π/
12) and at k = 1. The height of the peaks in σ ( k ) is related to the two growth rate parameters µ and ν as shown inthe inset. The relationship between the free energy V given above and the associated dynamics can be written in the form, ∂U∂t = −∇ (cid:18) δ V δU (cid:19) = −∇ (cid:0) L U + QU − U (cid:1) . (7)Here U is a measure of density fluctuations around the mean density, the linear operator L chooses two differentwavelengths for pattern formation as mentioned above and the parameter Q quantifies the effect of destabilisingsecond order nonlinearities.This system is a good example to use in this work as previous results indicate that the solution landscape of thismodel is rich with many different types of patterns. The dynamics conserves the mass (i.e., the integral of U over allof space) in the system. This implies that the dynamics will only possess equilibria that are steady in time. Such asolution landscape which does not have (propagating) waves or other time varying equilibria simplifies the analysis ofapplying different techniques to determine equilibria. Another advantage is that the same model also can be used toanalyze crystallisation of soft matter in 3D. Hence, it bears a number of advantages in connection to our goals hereinand provides a ripe testbed also for possible extensions. A. Reduction to coupled ODEs
As explained in the previous section, the growth rates µ and ν dictate the height of the dispersion peak at the twochosen length scales k = q and k = 1. At the onset of pattern formation, where both length scales are promoted, wehave the situation with both µ and ν being marginal at a co-dimension 2 bifurcation. Close to this onset, we canassume that the resulting pattern in the density field U is formed of wave-vectors at these two marginal lengthscales.This then allows us to perform a weakly nonlinear analysis by expanding the density field U , the linear operator modelparameters µ , ν , the quadratic nonlinearity prefactor Q and time t in terms of a small parameter (cid:15) via a multiplescales expansion. For the choice of scalings with U = (cid:15)U , µ = (cid:15) µ , ν = (cid:15) ν , Q = (cid:15)Q and t = (cid:15) − t , we are ableto balance the effects of the parameters, nonlinearities and time at the same order of expansion. Additional detailsregarding the selection of the scales of the expansion can be found in . We further write U in terms of a finitenumber of marginal modes as U = N (cid:88) j =1 (cid:16) z j e i k . x + w j e i q . x (cid:17) + c.c. (8)For the choice of q = (1 /
2) cos( π/ ≈ . N = 6. The left panel of Figure 5 shows a representativeexample of the wave vectors involved in creating a 12-fold quasicrystal. -0.05 0 0.05-0.0500.05 FIG. 5: Left: Wavevectors at two lengthscales, | k | = q = 0 . | k | = 1 that form a dodecagonal quasicrystal.Wavevectors k j have modulus | k | = 1 and q j have | k | = | q | for j = 1 , . . . ,
6. Complex conjugate wave vectors are identfied withbar, i.e., ¯ k = − k . Right: Solution curves for the two coupled polynomials in (11) with the zero contour of the right hand sideof the first equation in (11) in blue and that of the second equation in red. Solutions of both polynomials indicated with blackmarkers: hollow circles for maxima, crosses for saddles and filled dots for minima. Substituting the above expression into Eq. (5), we can write the volume specific free energy v = V / ( (cid:15) Volume) as v = − µz ¯ z − Q ( z ¯ z + w ¯ w + w ¯ z + ¯ w z ) ¯ z − µ j =6 (cid:88) j =2 µz j ¯ z j − (cid:88) j =1 νw j ¯ w j − Q (28 other cubic terms) − (148 quartic terms) . (9)where (similar to ) we have written the contributions involving ¯ z explicitly up to third order. Assuming fullsymmetry among all wave vectors with the same wave number, i.e., z j = z and w j = w for all j , this expressionbecomes v = (cid:18) − µz − νw − Q w − Qw z − Qwz − Q z + 996 w + 24 w z + 60 w z + 24 wz + 996 z (cid:19) . (10)Taking functional derivatives with respect to z and w then gives us the coupled amplitude equations for the evolutionof z and w as a gradient system associated with this free energy:˙ z = − dvdz = 2 (cid:0) µz + 2 Qz + 4 Qwz + 2 Qw − w − w z − wz − z (cid:1) , ˙ w = − q dvdw = 2 q (cid:0) νw + 2 Qw + 4 Qwz + 2 Qz − z − wz − w z − w (cid:1) . (11) -0.05 0 0.05-0.0500.05 0 0.05-0.050 FIG. 6: Left: Evolutions using CTN (in black) and CTN with deflation (green). Evolutions using SOM combined with CTNbut without (cyan) or with (brown) deflation. Right: Zoomed in view of left panel. We see that both CTN and SOM behavesimilarly with and without deflation. The evolution after deflation in both cases still gets trapped to the original equilibriumwhich is now a pole. The system parameters are µ = ν = 0 . Q = 0 . B. ODE Equilibria
In order to obtain the equilibria of the above coupled equations, we convert them to two coupled polynomials, f = µz + 2 Qz + 4 Qwz + 2 Qw − w − w z − wz − z ,f = q (cid:0) νw + 2 Qw + 4 Qwz + 2 Qz − z − wz − w z − w (cid:1) . (12)Here, our emphasis is on the proof of principle of the techniques rather than on an exhaustive parametric study.Therefore, we solve this set of equations for a given value of µ = ν = 0 . Q = 0 . . Given the features of this numerical algebraic geometry software , we are guaranteed to obtain all realequilibria of the system. The right panel in Figure 5 shows the zero contours of the two functions f (in red) and f (in blue). It is worth mentioning in passing that these zero contours are central to suitable methodologies foridentifying roots (that we do not focus on here) such as the reduced gradient following method . Their crossingpoints then identify the nine possible equilibria (marked in terms of r to r in the right panel of Figure 5 and later).By determining the eigenvalues of the linear operator at each solution, we can determine the dynamical stability ofeach of the determined equilibria in the context of the reduced ODE model of Eqs. (11). We indicate the stabilityinformation of each equilibrium in the figures with black filled dot indicating a stable minimum, black circle indicatingan unstable maximum and black cross indicating a saddle point. The interpretation of these eigenvalues within therealm of the original PDE problem will be briefly discussed in the next section.If the system of interest is f ( z, w ) = [ f , f ] T , in addition to exploring CTN methods for f itself, we can also definea related squared operator as in the 1D example as g ( z, w ), which will have as its minima, all of the equilibria of f ( z, w ). The first step to try is to employ the combination of time stepping in terms of CTN on f ( z, w ) or g ( z, w ).Results of this are shown in Figure 6. In the left panel, we see the solution landscape of the system f just as inthe right panel of 5. The right panel of Figure 6 shows a zoomed in view of the left panel. The overlay of bluefilled circles at every equilibrium indicates that in the associated squared operator system of g ( z, w ), all the originalequilibria are now minima. Moreover, this panel shows that there are no additional minima that have been createdin the transformation. While it is guaranteed in this case, by construction of the effective landscape (associated with f ) that only the roots of the original system will constitute global minima of the associated squared-operator system,in principle additional nonzero minima, as well as separatrices between them, may arise in general. Starting from theinitial conditions of ( z, w ) = (0 . , − .
01) (indicated with a red cross), CTN on f shown as black curve approachesthe minimum denoted as r . Starting from the same initial condition, CTN on g (shown as cyan curve) approachesthe much closer equilibrium of r which, recall, is a minimum for the SOM, but an energy maximum for the originalsystem. This is naturally in line with energetic considerations in these flow-based systems.The next step that we considered (in analogy with our 1-dof exposition) was to try to include deflation in conjunctionwith the CTN method. When we try to deflate r for the CTN on f , the new dynamics (shown as green curve) seemsto follow the un-deflated black curve (i.e., a similar trajectory as before) before stopping. This occurs due to the fact -0.1 -0.05 0 0.05 0.1-0.1-0.0500.050.1 -0.1 -0.05 0 0.05 0.1-0.1-0.0500.050.1 FIG. 7: Left: Asymptotic results of the dynamics stemming from evolving the CTN system with initial conditions on the ( z, w )plane and vanishing initial speeds. The space of possible initial conditions in the ( z, w ) plane has been split into four regions,each corresponding to one of the 4 stable equilibria of the original CTN system. The stable equilibria are shown as blackfilled circles for reference. Right: the same view but now from combining SOM with CTN. Here we observe the non-trivialpartitioning of initial conditions in the plane of ( z, w ) evolving to one of the 9 equilibria of the original CTN system. Stableequilibria of the SOM-CTN system shown as black filled circles (with the same naming convention as in Fig. 5( b ) and later).System parameters are µ = ν = 0 . Q = 0 . that in the deflated dynamics, the earlier energy minimum has become a pole where the particle can slide towardsindefinitely negative effective potential energy. While this landscape interpretation may no longer be formally possible(as the system may not be, by construction, a gradient one), nevertheless, it is used here by analogy with the 1-dof case(where such an interpretation can always be brought to bear). The numerical time stepping algorithm accordinglyadjusts the step-size and takes smaller and smaller steps till the rate of increase exceeds the limits set by the smallestallowed stepsize and the algorithm stops. Nevertheless, even if it did continue, the above analysis suggests that itwould not reach the desired result (i.e., an alternative root). A similar problem persists for the CTN on g (i.e.,combined with SOM). After deflation of r as well, the system follows the brown curve, which is a nearby curve tothe one without deflation (cyan). Hence we see that the pole problem that we encountered in the 1D example, whencombining flow methods with deflation, persists, as may be natural to expect.Having encountered problems upon starting from a single initial condition and using deflation combined with flowmethods to recover all equilibria in this system, we try to use flow methods starting at different initial conditions toidentify the basins of attraction of each minimum in the system. By the term ‘basins of attraction’, in the case of thepresent non-autonomous systems, we will mean the regions leading to the same asymptotic outcome of the dynamicsstarting from arbitrary initial conditions within the plane of ( z, w ), i.e., within the chosen symmetry subspace, alongwith vanishing initial speeds. Results of this analysis are in Figure 7. The left panel shows the basins of attraction ofthe minima for the 2-dof vector field f where a point is colored depending on the final equilibria that a flow methodstarting from that point would converge to. Minima of the system are marked with black crosses. It is interesting tocompare this figure with the right panel of Figure 5. The latter encompasses the saddle points of the energy landscape,while in the left panel of Figure 7 it is only possible to converge to the energy minima.The right panel in Figure 7 shows the basins of attraction for the squared operator system g ( z, w ) along with theminima for g shown as black crosses. This result shows us that using flow methods for the squared operator dynamicsallows to obtain all the equilibria of f , upon suitable scanning of the two-dimensional space of associated initialconditions of the SOM flow method. Naturally, there is a region containing the immediate neighborhood of each(global) minimum that results in convergence to it. However, both figures appear to be somewhat more complex inthat there exist regions of one color within a domain of the two-dimensional space that is predominantly of a differentcolor. This presumably has to do with the presence of some of the unstable nodes and saddle points emerging in themodified energy landscape V eff = (1 / f + f ), as well as with the effectively 4-dimensional nature of the CTNsystems (due to its involving accelerations in each of the two dynamical variables ( z, w )).Next we obtain the equilibria of the system f ( z, w ) using Newton’s method coupled to deflation in the 2-dof system.The results of this analysis are summarised in Fig. 8. Choosing ( z, w ) = ( − . , − . r at ( − . , − . r at ( − . , . -0.05 0 0.05-0.0500.05 FIG. 8: Newton with deflation for the 2D system plotted as a function of the two unknown amplitudes z and w . Startingfrom the same initial guess at ( z, w ) = ( − . , − .
03) (shown as a black cross) and sucessively using deflation, we are able toconverge to 7 of the 9 equilibria in this system (shown as black squares). both of these solutions, we see that successive Newton iterations follow the dark red line in Fig. 8 and converge to thesaddle at r at (0 . , − . r , r , r and r . We see that in this case, the additional equilibria r and r were not approachable usingNewton iteration with deflations of previously determined equilibria. However, choosing other initial conditions, e.g.,( − . , − . in an infinite-dimensional (and hence much richer) problem, there are no a priori guarantees that deflation from agiven initial guess will converge to all solutions. Indeed, even in our simple 2-dof setup, the method is not able toconverge to all solutions. Hence, it is relevant to seek a broader representation of the space of initial guesses in orderto capture the widest possible solution span. C. Extension to the Original PFC Problem
The next question to ask is how much of the information from the ODE system remains relevant to the dynamicsof the (considered numerical discretization of the) full PDE system. In addition to seeking to identify the relevantconfigurations at the level of the PDE, we will also seek to connect the “reduced stability” of the 2-dof systemwith the full PDE stability at the level of the original PFC system bearing infinite degrees of freedom. The twomain assumptions that allow for the reduction from a PDE to an ODE system as derived in section III A are, ( i )the dynamics is assumed to be in the neighbourhood of a critical point, here a co-dimension 2 bifurcation and ( ii )a specific symmetry is assumed in the form of marginal modes, here wave vectors corresponding to a dodecagonalquasicrystal which are chosen in Eqn. (8). This implies that the predictions from the ODE system should matchthe PDE dynamics close to the co-dimension 2 bifurcation when the PDE dynamics is restricted to the dodecagonalsymmetry subspace.With respect to relaxing the first assumption, previous work such as has pointed out that the occurrence of otherglobal bifurcations close to the critical point of interest can restrict the region in parameter space where the ODEsystem reproduces the dynamics of the PDE system. However, in a generic case with no other global bifurcations inthe neighbourhood of the critical point in phase space, the predictions of the reduced ODE system can match wellthe PDE over a wide range of parameters around the critical point. We discuss how the current system falls in thisgeneric category by determining related PDE equilibria to each determined ODE equilbrium.Firstly, we deploy a Newton method to check the existence of equilibria for the full PDE system corresponding tothe ODE solutions. In this process, it is relevant to bear in mind that we can recover associated equilibria for thePDE even if they are dynamically unstable. The results of our Newton iterations are shown in Figure 9 with colored1 -0.05 0 0.05-0.0500.05 FIG. 9: Comparison of Newton iteration for the PDE system with equilibria of the ODE system as initial guess. Each Newtoniteration is shown as a colored circle. The black square is the equilibrium of the PDE system related to r of the ODE system,which is discussed in more detail in Figure 10. System parameters are µ = ν = 0 . Q = 0 . circles indicating converged results at each Newton iteration starting from different ODE equilibria are indicated withdifferently colored circles. The observation is that for the case with µ = 0 . ν = 0 .
1, all of the ODE solutionscan be recovered as solutions to the full PDE system using iterative methods for the full PDE. The reconstructionof the corresponding ODE solutions at the level of PDE waveforms via Eq. (8) allows to construct suitable, rapidlyconverging (at the level of the Newton iteration) initial guesses.Secondly, we compare the stability information for the ODE equilibria with an estimate for their stability in thePDE system from evolutions of the PDE system. Let us consider the example of the ODE equilibria at r which is aminimum with eigenvalues: λ = − .
061 and λ = − . λ is along the vector ( z, w ) = ( − . , − and adding it to the state r . During evolution ofthe PDE system from this initial condition, we observe that the perturbation decays and we recover the equilibrium r as the asymptotic state. We track the decay of the perturbation as a function of time in log scale as shown in theleft panel of Fig. 10. We see that in this case, the slope of the decay in log scale is − . λ . Similarly we also confirm that positive eigenvalues predicted from the ODE systemmatch the linear growth rate of perturbations along the corresponding eigenvector. For example, we consider thesaddle at r with positive eigenvalue λ = 0 . . λ P DE = 5 .
69. The eigenvalues calculated for the PDE equilibrium related to r areshown in the right panel of Figure 10. We note that there are corresponding eigenvalues in the PDE solution at theeigenvalues predicted by the ODE solution, however we see that the PDE solution is highly unstable with multiplepositive real eigenvalues (shown as black crosses). Therefore, we expect flow methods for the PDE to not recover theODE solution in the asymptotic limit for generic initial perturbations that do not reinforce any symmetry.In the next step, we examine the asymptotic results of the PDE dynamics when we consider initial conditions thatfall within the dodecagonal symmetry subspace with different values for ( z, w ) similar to the investigation in the leftpanel of Fig. 7 for the CTN dynamics of the ODE system. For the PDE system, we track the dynamical evolutionby using pseudospectral methods to discretise in Fourier space and by employing an exponential time differencingmethod to step in time; see for details. Starting from different initial guesses for ( z, w ), taken from a grid in the2 -2 0 2 4 6-0.5-0.2500.250.50 20 40 60 80-24-23-22-21-20-19-18 FIG. 10: Left panel: Evolution of a perturbation along the eigenvector corresponding to the largest eigenvalue λ = − . r . Slope of this evolution in log scale − . λ , Right panel: Comparison of eigenvalues from the ODE (in red circles) and corresponding PDE solution(in black crosses) for r . System parameters are µ = ν = 0 . Q = 0 . range ( − . , .
1) for each z and w , we identify the asymptotic states that the PDE system reaches. Figure 11 showsthe types of asymptotic states that are reached when starting from different initial conditions for ( z, w ). The zero/flatstate is when all perturbations decay and the system reaches a quiescent state. When starting with an initial guess of( z, w ) = (0 , w ≈ .
05, irrespective of the value of z . The relevant asymptoticstate has the value ( z, w ) = ( − . , . r which is marked bya black square in Figure 9. All initial conditions in the ( z, w ) plane that evolve to this state are indicated by themustard colour region in Figure 11. Here we note that states with small amplitudes, i.e., less than 10 − , at peaks thatarise from higher order combinations of the original wave vectors are still encompassed within this mustard coloredregion. It is interesting that one of our ODE equilibria is also a relevant potential attractor of the PDE dynamicswhen starting in a quasicrystal symmetry subspace.As discussed using the results in Figs. 9 and 10, equilibria of the ODE system generically persist as equilibria forthe associated PDE system for parameters close to the original critical point and looking for equilibria with a chosensymmetry. However, the stability characteristics of these equilibria in the PDE system are not straightforward todetermine a priori. It may often be the case that the unstable ODE equilibria might be associated with unstablePDE equilibria. Indeed, this has happened for all the unstable ODE equilibria (saddles and maxima) that we havepreviously considered in the 2-dof ODE system. A more complex situation might be when we have a stable ODEequilibrium, i.e., r which is a minimum of the reduced energy landscape and yet the associated PDE equilibrium ishighly unstable with multiple eigenvalues that have a positive real part. Since PDE evolutions that are initialised fromany initial condition are not restricted to stay within a specific symmetry subspace, many other complex patternsbecome possible even when we start the evolution within the ( z, w ) plane. We explore such patterns in further detailin the rest of this section.For chosen evolutions that start from near one of the ODE equilibria, we plot the initial condition, the final stateand a plot that indicates how much a current state during evolution is like the initial pattern at wavenumber k = 1and at wavenumber k = q , i.e., we measure the distance between the current PDE state and its projection on thereduced space of wavenumber k and on that of wavenumber q . The distance from these projections measures howclose (or far) the solution of the PDE lies in comparison to the plane of our 2-dof reduction. In the language ofthese distances and projections, for the stable PDE state discussed above, we start on the 2-dof plane and we stayon it, resulting in a stable PDE solution. However, it is also possible that states may be unstable (or even stable)at the 2-dof level, allowing small perturbations in the direction transverse to the ( z, w ) plane to take advantage ofthe infinite degrees of freedom of the PDE to diverge from the plane and result in a different solution. In that light,such a plot allows us to visualise how the evolution proceeds to redistribute energy between the length scales (andassociated wavevectors) and will allow us to explain our previous observation that initial states with large values of3 -0.1 -0.05 0 0.05 0.10.10.050-0.05-0.1 zero/flatz-rhombw-rhombz-hexw-hexpattern-4i8opattern-12i6opattern-8i10oSymm QCSquaresymm QC FIG. 11: Dynamical outcomes of different initial conditions reconstructed through the projection in the ( z, w ) plane for thePDE system yielding different types of patterns.
FIG. 12: Left panel: Tracking the component of a PDE evolution at the k = 1 and k = q wavelengths. Time stepping startedfrom amplitudes close to r z = − .
01 and w = 0 .
07 (shown as red circle in figure), Middle panel: Prescribedinitial condition and Right panel: Asymptotic state. w ≈ .
05, irrespective of the value of z all converge to the PDE equilibria associated with r (shown as a black filledcircle), while others as we will see below converge to other states.The left panel of Fig. 12 shows the evolution when we start from a point close to r , with ( z, w ) = ( − . , . U field shown in the middle panel, has only components along k = 1 and k = q . This can be inferred from the values of the coordinates of the points which are (1,1,0) in theleft panel. The value 1 here refers to a normalised value of the projection of the current state with respect to thewavevectors at k = 1 and k = q respectively. As time progresses, the evolution moves away from this initial point(identified with a red circle) but stays in the ( z, w ) plane. The corresponding asymptotic state, shown in the rightpanel, looks very similar to that of the initial condition. The only change is in slight difference in the value of theamplitudes ( z, w ) from ( − . , .
07) to ( − . , . r . The reduction in both the ( z, w )amplitudes from the initial value is reflected in the trajectory of the evolution in the || U · ze ikx || and || U · we iqx || plane.When starting from an initial guess ( z, w ) = (0 . , .
05) in Figure 13, which is close to the equilibrium r (a saddlefor the ODE system), we observe the PDE system to evolve to a different looking quasicrystal pattern. This initialpoint is within the mustard colored basin of attraction of the PDE equilibrium at r (as the classification ignorespeaks with amplitudes less than 10 − in the asymptotic state). As the phase relation between ( z, w ) at r is notthe same as that of the chosen initial condition, we observe that the evolution causes the z amplitude to decrease tonegative values, while the w amplitude value increases slightly from the initial guess value. Furthermore, we look atthe Fourier spectra of the initial and final states in Fig. 14 plotted in log scale to highlight the differences at smallvalues. Here we see that the initial state has all its amplitudes in one of the modes at either of the two lengthscales.The evolution of the phase between ( z, w ) adds non-zero peaks at other lengthscales, resulting in the Fourier spectraof the final state having components other than at k = 1 and k = q . This also explains the result in the first panel ofFig. 13, where the path tracked ends at an asymptotic state that is out of the || U · ze ikx || and || U · we iqx || plane.4 -0.50000.5 0.50.51 11.5 11.5 FIG. 13: Left panel: Tracking the component of a PDE evolution at the k = 1 and k = q wavelengths. Time stepping startedfrom amplitudes close to r z = 0 .
04 and w = 0 .
05 (shown as red circle in figure), Middle panel: Prescribedinitial condition and Right panel: Asymptotic state. -15-10-50 -15-10-50
FIG. 14: Left panel: Fourier spectra of the prescribed initial condition in Fig. 13 in log scale and Right panel: Fourier spectraof the asymptotic state also in log scale.
We then explore in Figure 15 an evolution starting from a point close to the r equilibrium (which is a minimumfor the PDE system) with initial values of ( z, w ) = (0 . , − . z, w ) plane for a while as it changes the relative phase between z and w , after which there is a sharp changeto damp out all the w modes resulting in a z − hexagon pattern as the asymptotic state. This is consistent with theinitial guess being in the basin of attraction of z − hex in Figure 11. It is also in line with our discussion about theplane of ( z, w ) not being an invariant plane constraining the PDE dynamics, but rather the possibility of the PDE topotentially move towards other equilibrium states outside of the relevant plane. FIG. 15: Left panel: Tracking the component of a PDE evolution at the k = 1 and k = q wavelengths. Time stepping startedfrom amplitudes close to r z = 0 .
06 and w = − .
01 (shown as red circle in figure), Middle panel: Prescribedinitial condition and Right panel: Asymptotic state.
As a final figure, we also include examples of the other types of pattern in the PDE system along with their powerspectrum in Figure 16. More details about their symmetry subspaces and the number of unknown amplitudes thatneed to be determined using ODE reductions to determine them, are listed in Table I. The key observation out of5
TABLE I: Types of equilibria observed as asymptotic states in the PDE system
Type of equilibrium ( z , z , z , z , z , z ; w , w , w , w , w , w ) No. of unknownamplitudesFlat (0,0,0,0,0,0; 0,0,0,0,0,0) 0 z -stripes (z,0,0,0,0,0; 0,0,0,0,0,0) 1 w -stripes (0,0,0,0,0,0; w,0,0,0,0,0) 1 z -hexagons (z,z,z,0,0,0; 0,0,0,0,0,0) 1 w -hexagons (0,0,0,0,0,0; w,w,w,0,0,0) 1 z -rhomb (0,0,0,z,z,0; w,0,0,0,0,0) 2 w -rhomb (0,0,0,z,0,0; 0,0,w,w,0,0) 2pattern-12i6o (z,0,z,0,z,0; w,w,w,w,w,w) 2symm-QC (z,z,z,z,z,z;w,w,w,w,w,w) 2pattern-4i8o ( z , z ,0, z ,0, z ; 0,0, w , w ,0,0) 4pattern-8i10o ( z ,0, z , z , z , z ; w , w ,0, w , w ,0) 4Squaresymm QC ( z , z , z , z , z , z ; w , w , w , w , w , w ) 4 this is that the original reduction even via its unstable solutions enables the exploration of interesting dynamicsof the associated full PDE model. Indeed, our 2-dof plane was the basis for unraveling a wide variety of PDE-level solutions with different (hexagonal, rhombic, quasi-crystal and other complex pattern) symmetries. Hence, itprovides a methodology that combined with flow, Newton and deflation methods allows an in-depth exploration ofstates accessible to the PDE-problem of interest. IV. CONCLUSIONS AND FUTURE CHALLENGES
In the present work, our scope was to develop some simple examples, suitable for the testing of different types ofmethods (the continuous-time Nesterov method, the square-operator method, the deflation technique, and of coursethe widely used Newton’s method). As such we started with a very simple and transparent in its roots, yet sufficientlyinstructive 1-degree-of-freedom example. There, we saw that flow methods like CTN can only lead to local minima.We modified these then to obtain all roots as local minima using the square-operator method. We appreciated thatneither of these two methods was suitable for utilization in conjunction with the deflation toolbox. When we turnedto Newton’s method plus deflation, we were able to obtain all associated roots.We then extended consideration to a PDE problem, namely a PFC model for crystallization which is an infinitedegree of freedom system. Yet, a reduction can be used at the level of unknown amplitudes at a finite number ofwavevectors at two sets of wave numbers that ultimately projects the problem to a 2-degree-of-freedom setting. Insuch an ODE system, it is possible to obtain all equilibria. It is possible to flow into the stable ones (CTN), flow evento the unstable ones (SOM) and to deflate with Newton to obtain almost all (but not all) starting from the sameinitial guess. The stability of the different equilibria and the corresponding landscape within the two-dimensionalprojection was explained, along with the basins of attraction of different methods.Then, we turned to the full PDE computation. Here, the utilization of Newton iterations enabled the identificationof all the fixed points previously found at the level of the ODE reduction. Working with the same two-dimensionalplane of possible initial conditions, but now at the level of the PDE, the dynamical evolution of the model and itspotential outcomes were evaluated for arbitrary initial conditions reconstructed from the two-dimensional projection.This allowed us to identify numerous novel PFC patterns and to observe how the dynamics may stay within ordepart from the plane of the 2-dof reduction. The eigenvalues at the ODE and PDE levels were compared andsimilarities/disparities were discussed.There are numerous important directions to consider for future studies along this vein. In many of the flow problems,a nontrivial issue that arises is that of running into infinities either because of the inherent dynamics of the problemor because of the effective landscape created, e.g., by poles induced by deflation. It would be especially relevant toattempt to generate globally convergent methods that could bypass this type of issue. Indeed, at the level of algebraicnumerical tools, such as Bertini, there exist suitable complex plane tricks that allow for such properties. Comparingmethods that are used in that latter context, such as the so-called Davydenko method, with methods here such as thesquare operator method is also of interest in its own right. It is relevant to note here that if such methods bypassinginfinities could be easily built-in, then it may be possible to enable flow methods with deflation to work efficientlytowards identifying all possible roots of a problem. Notice also that deflation here is done “isotropically” with respectto different degrees of freedom. Yet, it is not inconceivable to envision a variant of the method which is anisotropic6in its handling of variables. Numerous among these relevant directions are currently under consideration and will bereported in future publications. Acknowledgements.
This material is based upon work supported by the US National Science Foundation underGrants No. PHY-1602994 and DMS-1809074 (PGK). PS acknowledges helpful discussions with C´edric Beaume andAlastair Rucklidge along with support from a Hooke Research fellowship at the Mathematical Institute of the Uni-versity of Oxford. PGK also acknowledges support from the Leverhulme Trust via a Visiting Fellowship and thanksthe Mathematical Institute of the University of Oxford for its hospitality during part of this work. C.T. Kelley,
Iterative methods for linear and nonlinear equations , SIAM (Philadelphia, 1995). D.J. Bates, J.D. Hauenstein, A.J. Sommese, C.W. Wampler,
Numerically solving polynomial systems with Bertini , SIAM(Philadelphia, 2013). M.J. Ablowitz, Z.H. Musslimani, Opt. Lett. (2005) 2140-2142 J. Yang, T.I. Lakoba, Stud. Appl. Math. (2008) 265-292. J Yang,
Nonlinear Waves in Integrable and Nonintegrable Systems , SIAM (Philadelphia, 2010). Y Nesterov, Sov. Math. Doklady, (1983) 367-372 S Weije, S Boyd, EJ Candes, J. Mach. Learn. Res. , 1 (2016). C.B. Ward, N. Whitaker, I.G. Kevrekidis, P.G. Kevrekidis, arXiv:1710.05047. A Wibisono, A Wilson, M Jordan, Proc. Nat. Acad. Sci. , E7351 (2016). C.-M. Schmied, T. Gasenzer, M.K. Oberthaler, P.G. Kevrekidis, Comm. Nonlinear Sci. , 105050 (2020). A. J. Sommese, C. W. Wampler,
The Numerical Solution of Systems of Polynomials arising in Engineering and Science ,World Scientific (Singapore, 2005). Yu.S. Kivshar and G.P. Agrawal,
Optical solitons: from fibers to photonic crystals , Academic Press (San Diego, 2003). B.A. Malomed, Progr. Opt. , 69 (2002). Yu.S. Kivshar, B.A. Malomed, Rev. Mod. Phys. , 763 (1989). P. E. Farrell, A. Birkisson, and S. W. Funke, SIAM J. Sci. Comp. , 2026 (2016). P. E. Farrell, C. H. L. Beentjes and ´A. Birkisson, arXiv:1603.00809 (2016). E. G. Charalampidis, P. G. Kevrekidis, and P. E. Farrell, Commun. Nonlinear Sci. Numer. Simulat. , 482 (2018). E.G. Charalampidis, N. Boull´e, P.E. Farrell, P.G. Kevrekidis, Commun. Nonlinear Sci. Numer. Simulat. , 105255 (2020). P. Subramanian, A.J. Archer, E. Knobloch, A.M. Rucklidge, Phys. Rev. Lett. , 075501 (2016). P. Subramanian, A.J. Archer, E. Knobloch, A.M. Rucklidge, New J. Phys. , 122002 (2018). P.G. Kevrekidis, C.I. Siettos, Y.G. Kevrekidis, Nature Comms. , 1562 (2016). A. M. Rucklidge, M. Silber, and A. C. Skeldon Phys. Rev. Lett. , 074504 (2012). W. Quapp, J. Theor. Comp. Chem. , 385 (2003). R. W. Wittenberg and P. Holmes, Physica D , 1 (1997). -0.200.20.40.6 -0.200.20.40.6 -0.200.20.40.6 -0.4-0.200.20.40.60.010.020.030.040.050.06