{'fraction_non_alphanumeric': 0.040229056378355, 'fraction_numerical': 0.012853523389704823, 'mean_word_length': 4.200642742367434, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Microscopic Patterns in the 2D Phase-Field-Crystal Model

Using the recently developed theory of rigorously validated numerics, we address the Phase-Field-Crystal (PFC) model at the microscopic (atomistic) level. We show the existence of critical points and local minimizers associated with"classical"candidates, grain boundaries, and localized patterns. We further address the dynamical relationships between the observed patterns for fixed parameters and across parameter space, then formulate several conjectures on the dynamical connections (or orbits) between steady states.


Introduction
The Phase-Field-Crystal (PFC) model introduced in [1] is a gradient system capable of modeling a variety of solid-state phenomena. In its simplest form, the PFC energy can be written as defined on phase-fields ψ ∈ H 2 (Ω) satisfying the phase constraint The parameter β represents inverse temperature such that β = 0 models maximum disorder. Coupled with this energy is its conservative H −1 gradient flow Figure 1: Left: Details of a grain boundary appearing in a PFC simulation (taken from [3]). Right: Grain boundary network from a PFC simulation (taken from [4]). Within each grain is a hexagonal lattice of atoms with a particular orientation.
which entails the sixth-order PFC equation Note that the PFC model shares its energy with the Swift-Hohenberg equation [2], which is simply the L 2 gradient flow of E. From linear stability analysis applied to single Fourier mode Ansatz, we find three main candidate global minimizers that divide parameter space, see the appendices. In the hexagonal lattice regime, 2D-simulations of the PDE starting with random noise quickly produce atoms that arrange into small patches of hexagonal lattices with random orientations. These patches grow and interact with each other, forming grains of hexagonal lattices of atoms with a particular orientation. The morphology and evolution of these grains have features resembling those in polycrystalline materials (cf. Figure 1). In particular, it has recently been shown that statistics of many of experimentally observed (universal) grain boundary distributions are accurately captured by data amassed from simulations of this simple PFC equation [5,4]. While here we will mostly work with this vanilla PFC formulation, we note that a family of PFC-like equations can be derived from Density-Functional-Theory [6] to obtain more complicated models capable of simulating eutectic and dendritic solidification [7] and graphene structures [8,9]. In this article, we address the PFC model and its steady states at the "microscopic" level -the local atomic arrangement. We believe that such an investigation of microscopic pattern-formation capabilities of PFC is not only of mathematical interest but is also necessary to construct "designer" models for polycrystalline behaviour. For example, varying the parameters in the energy lead to more complicated states than simple lamellar and hexagonal. These include localized patterns in the "glassy regime" -the transition at the liquid (constant) and solid (hexagonal) transitions -and "globules" at large β.
With the exception of the constant (liquid) state (cf. [10]), it is difficult to prove any theorem on the exact nature of steady states, local and global minimizers to this diffuse interface problem. What exists in the physics literature is numerical simulations, standard linear stability analysis, and Ansatz-driven energy comparisons. The recently developed theory of rigorously validated numerics (cf. [11,12,13,14,15]) now provides a powerful new tool to bridge what can be observed numerically with rigorous statements on pattern morphology. In a nutshell this approach can be summarized as follows: Given an approximate steady state, we use the Contraction Mapping Theorem to imply the existence and local uniqueness of an exact steady state within a controlled distance of the approximation. This notion of closeness is strong enough to imply further useful results, including closeness in energy and stability results. In this paper we use this new approach to address the following aspects of the PFC model: • Are the "classical" candidates obtained from linear stability analysis close to actual local minimizers?
• Are the stable yet complicated patterns observed numerically indeed critical points in the PFC energy landscape? For example, are grain boundaries steady states or simply metastable states?
• What are the dynamical relationships between the observed patterns for fixed parameters and across parameter space?
Based upon our results we formulate several conjectures on the connections (or orbits) between steady states. Taken as a whole, our work presents the first step into a rigorous analysis of the rich PFC energy landscape. The outline of this paper is as follows. We first setup the PFC equation in Fourier space and discuss the application of the framework of rigorous computations. We then verify the existence of important steady states of the PFC equation, including localized patterns and grain boundaries. With these states in hand, we address the energy landscape of PFC with a discussion on conjectures for connections (or connecting orbits) between steady states. Finally, we presents results in one-parameter numerical continuation to outline some interesting features of the bifurcation diagram of PFC.

PFC steady states in Fourier space
We begin by writing the equation ψ t = 0 in Fourier space to obtain a coupled system of equations for the Fourier coefficients of steady states. We will be slightly more general and consider functionals of the form where K is a linear differential operator K acting on elements of a suitable function space. In particular, K = ∇ 2 + 1 for the basic "one-mode" PFC model (∇ 2 + 1)(∇ 2 + q 2 ) for the "two-mode" PFC model [16] where q is the secondary wavelength of two-mode PFC. Taking the H −1 gradient flow of E, we obtain the PFC-like equation For simplicity, we let Ω be the rectangular domain [0, L x ] × [0, L y ] with periodic boundary conditions. We let where N x , N y ∈ N are the number of atoms lined up in the x, y-axes. The main parameters of the problem are then (ψ, β) and the domain size is given by Let a α be the Fourier coefficients of ψ and let (a α ) t be the time derivative. Inserting this expansion into the PFC equation results in an infinite system of equations of the form (a α ) t = F α (a) thanks to orthogonality. The steady states may then be found numerically by solving F (a) = 0 up to some truncation order M . We will see later that it is imperative to isolate the zeros of F ; the continuous translational and rotational symmetries of PFC must then be broken. The simplest way to do so in this context is to also enforce Neumann boundary conditions. It is convenient to write a α = a α1,α2 so that the symmetry and reality conditions become a |α1|,|α2| ∈ R.
This choice allows us to simplify a complex Fourier series into the cosine expansion where W is a weight matrix defined by The Fourier coefficients of ∇ 2 ψ are given by the elementwise product L α a α where is the Fourier representation of the Laplacian. Inserting these expressions into the PFC equation and equating Fourier modes, we obtain where * denotes the discrete convolution and the linear terms combining K and β are Note that the (0, 0) Fourier component picks out the average phase so it is fixed toψ: this is consistent with (a 0,0 ) t = 0 thanks to L 0,0 = 0. To keep track of the phase constraint directly in F , we replace its first trivial component by F 0,0 = a 0,0 −ψ, resulting in: The operator F then represents the PFC dynamics in the sense that its zeros correspond to steady states of the PFC equation. A numerical advantage of the reduced expansion is that we effectively only have to compute a quarter of the full Fourier series. Obviously, this means we are not treating PFC in full generality over H 2 and will have to address this later. As an aside, the equivalent F for Swift-Hohenberg is simply −(γ α a α + (a * a * a) α ) hence its (0, 0) entry is nonzero and average phase is not conserved.

Overview of rigorously validated numerics
We present a brief overview of the recent framework of rigorously validated numerics for dynamical systems, see sources including [11,12,13,14] and [15] for a survey of techniques for PDEs.
Consider the Newton-like operator T (a) = a − AF (a), where A is a suitable inverse to the derivative DF (a). On the one hand, if T is a contraction on a closed ball, the contraction mapping theorem gives the existence and uniqueness of a zero of F within this ball. On the other hand, the repeated application of T (allowing A to vary with a) should converge to this fixed point. We can then numerically compute an approximate steady stateā for which F (ā) ≈ 0 up to numerical precision. If in addition we are able to show that T is a contraction aroundā, then we immediately have the existence of an exact steady state a close toā in an appropriate metric. This relationship is made clear by the radii polynomial theorem, so-called for reasons that will become clear shortly. To illustrate the method, we specialize the theorem to the case applicable to PFC, but see [17,18,19,20] for different approaches and [21,22] for an application to Ohta-Kawasaki functional in 2D and 3D, respectively. Given Banach spaces X, Y , we use the notation B(X, Y ) for the space of bounded linear operators from X to Y , B(X) = B(X, X) and B r (a) ⊂ X for the open ball of radius r around a ∈ X.
where Y 0 , Z 0 , Z 1 are positive constants and Z 2 is a positive polynomial in r > 0. Construct the radii polynomial If p(r 0 ) < 0 for some r 0 > 0, then there exists a unique a ∈ B r0 (ā) for which F ( a) = 0.
The proof of this formulation is given in appendix B, where we show a correspondence between the sign of the radii polynomial and the contraction constant of T : if r 0 can be found, T is a contraction and the Newton iteration starting atā must converge to some a. This proves not only the existence of the exact steady states but also gives control on its location in X with respect to a known point. In practice, one finds an interval [r * , r * ] of radii for which p(r) is negative; r * > 0 gives the maximum distance betweenā and a while r * > r * gives the minimum distance betweenā and another zero of F . The zeros of F must therefore be isolated for consistency.
Each bound may be understood intuitively: Y 0 being small indicates thatā is a good approximation of a while Z 1 being small indicates that A † is a good approximation for DF (ā), and so on. These bounds may be simplified analytically but must necessarily be computed numerically. Therefore, we ensure that our numerical computations go in the same direction as the required inequalities by using interval arithmetic [23], a formalized approach to deal with numerical errors. We used the interval arithmetic package INTLAB for MAT-LAB, see [24,25], to ensure that the radii polynomial approach is numerically rigorous.
This approach allows us to prepare numerical tools that can both find candidate steady states and compute the radii r * , r * if they exist. If so, we immediately have a proof that this candidate provides a good handle on an actual steady state of the PFC equation.

Radii polynomial approach for PFC
Let us now apply these ideas to PFC by first computing DF and the Newton operator. Let σ represent the differentiation indices applied to F α . The derivative of F 0,0 is 1 if σ = (0, 0) and 0 otherwise, so we use the Kronecker delta notation to write For other values of α, the linear terms similarly give The derivative of the nonlinear triple convolution can be computed by differentiating with respect to all four a α identified by symmetry. This algebraic computation is somewhat tedious but the result can be written succinctly as so that the full derivative of F is: a, F and the convolutions may be viewed as infinite matrices whose "topleft" entry is the (0, 0) coefficient while the derivative DF is an infinite 4-tensor.
To implement the Newton method numerically, such objects must be truncated to order M such that a σ = 0 whenever either σ 1 or σ 2 is greater than M . This results in the (M + 1) 2 matrices a (M ) , F (M ) while the derivative becomes the (M + 1) 4 4-tensor DF (M ) . Note that the k-convolution of a (M ) has support kM by definition.
We now introduce the Banach space framework. Let ν > 1 and define 1 ν (Z 2 ) as the space of sequences a α with finite norm The restriction of 1 ν (Z 2 ) using the symmetry condition is over which the norm simplifies to where ν α is a weight matrix that forces the fast exponential decay of the Fourier coefficients. The space (X, || · || 1,ν ) can easily be shown to be Banach and the 2D discrete convolution forms a Banach algebra over it, immediate results from the triangle inequality and the fact that ν > 1. Let nowā, a ∈ X have the same meaning as before, withā = 0 outside of U = {0, 1, ..., M } 2 thanks to the truncation. Let G = DF (ā) (M ) and denote by A (M ) the numerical inverse of G. We define approximate operators A † , A as which can be thought of as block tensors containing G or its inverse paired with the linear terms L α γ α as the main "diagonal" of the second block. If G is an invertible matrix, 1 so is A and it is thus injective. The inverse of A is not A † however because A (M ) G ≈ I (M ) only up to numerical inversion errors. Note that F, DF and A † map to a space Y with less regularity than X because of the unbounded L α γ α terms arising from real space derivatives; Y is a space where sequences L α γ α a α have finite norm. However, the operator products against A are bounded on X thanks to the fast decay of L −1 α γ −1 α . Thus, we say that A "lifts" the regularity of the other operators back to X, allowing statements such as T : X → X or ADF (ā) ∈ B(X).
We show in appendix C how to simplify the bounds into expressions that can be evaluated numerically. This allows us to write down the radii polynomial 2 r, hence the polynomial is cubic with non-negative coefficients except for maybe the linear term. We have p(0) > 0, p (0) = Z 0 + Z 1 − 1 and p(r) → ∞ for large r. As a consequence, if p is strictly negative for some positive r, there must exist exactly two strictly positive roots r * < r * defining the interval where the proof is applicable. When this is satisfied, the radii polynomial theorem gives that 1. There exists an exact solution a of F (a) = 0 in B r * (ā).

This solution is unique in
Thus, when the radii polynomial is computed using interval arithmetic and has exactly two real non-negative roots, the zero computed numerically with the Newton iteration is close to an actual steady state of the PFC equation. Note the important fact that the ball is in X so a priori, only the Fourier coefficients are controlled. Thanks to ν > 1 however, we show in appendix D that this control translates into closeness in energy and in real space norms. In particular, the distance in value between the phase fields corresponding toā and a is at most r * .
Further, we show in appendix E that the stability of a in X is controlled by the eigenvalues of G. It is important to observe that this matrix will always have a positive eigenvalue because of the trivial condition F 0,0 = a 0,0 −ψ. This is not indicative of instability in the context of the H −1 gradient flow because a 0,0 is fixed. We shall see later that this unstable direction can be used to compute a branch of solutions in parameter continuation. For now however, we call the number of positive eigenvalues, minus 1, the Morse index of a, indicating how many unstable directions are available to a given steady state for fixed parameters.
The procedure to numerically investigate the steady states of the PFC equation is as follows: • Starting from a given initial condition, the Newton iteration is run until it converges up to numerical precision.
• Then, the radii polynomial of the numerical guess is computed and its roots are tested.
• If the proof succeeds, we can characterize an exact steady state in value, in energy and compute its stability in X. The parameters (M, ν) can be adjusted until the proof succeeds with a trade-off between the computational effort and closeness in X.

Rigorous results on small domains
We now have a complete framework for finding verified steady states along with their energetic and stability properties. This allows us to understand the behavior of the PFC system for a given choice of (ψ, β), with three important caveats: • We cannot guarantee that we have found all steady states and therefore the global minimizer. Indeed, we may only hope to cover a reasonable portion of the underlying space by sampling initial conditions randomly.
• The size of M must be balanced with ν to keep r * as small as possible, keeping in mind that r * is ultimately bounded above by the distance between two steady states. In particular, large domains and large β increase the contribution of high frequency Fourier modes, hence the truncation order can become large even for domains containing only 100 atoms. This limits our results to small domains so our analysis is "small scale" in nature.
• The Neumann boundary conditions restrict us to a "quadrant" of H 2 . While the existence of a steady state, the energy bound and instability obviously extend to H 2 , stability does not as there may be unstable directions in the other three Fourier series that are missed by the current method.
For the last point, we sometimes observe that translational shifts have a different Morse index in X. This is observed for example with the stripes states, see Fig. 3 (a). In this sense, we only provide a lower bound for Morse indices in H 2 .

Verification of the candidate minimizers
The candidate global minimizers (constant, stripes, atoms and donuts states) introduced in appendix A have trivial Fourier coefficients by construction, given by where A s , A h represent amplitudes that optimize the PFC energy calculation. Note that A h differs between the atoms and donuts states.
To illustrate the approach, we first applied the verification program starting at the atoms state b constructed for (ψ, β) = (0.07, 0.025), (N x , N y ) = (4,2) and M = 20. The Newton iteration was used to obtainb for which the radii polynomial was tested with ν = 1.05, resulting in r * = 1.0 · 10 −11 and r * = 6.8 · 10 −3 . The 1 ν distance between b andb is 1.1 · 10 −3 , indeed smaller than r * . The difference b −b is mainly captured by new Fourier modes: we find that the main Fourier coefficients b Nx,Ny = b 0,2Ny = −4.4 · 10 −2 differ by 1.5 · 10 −5 while the largest new Fourier modes are b 8,0 = b 4,5 = −7.4 · 10 −5 . Moreover, the distance in the (numerical) sup norm between the two phase fields is approximately 4.4 · 10 −4 which is again smaller than the 1 ν distance, consistent with the L ∞ bound.
This approach was repeated for the other candidates and for a few other choices of the PFC parameters in the hexagonal regime, with truncation adjusted to β. The results are presented in Table 1, showing that such simple candidates capture well the leading behavior. Note that the agreement decreases with increasing β: compare the size of ||a −ā|| 1,ν to ||ā|| 1,ν .

Steady states in the hexagonal lattice regime
The Newton iteration can detect new steady states regardless of stability as it is based on criticality instead of minimality. This allows us to find steady states that are observed only momentarily or even locally during a PFC simulation. Table 2 presents a few of the 28 distinct steady states found for (ψ, β) = (0.07, 0.025), (N x , N y ) = (8, 5), ν = 1.05 and M = 40. Starting at random initial coefficient matrices, the Newton iteration converges in 15 to 50 steps. The four main ansatz were also explicitly tested, as only the atoms state could be reached from random initial conditions.
Note that the energy of the exact steady states can be compared from Table 2: for instance, the energy of the exact atoms state is bounded away from the others so it is guaranteed to be the best candidate global minimizer out of the observed steady states at the current parameter values.
The second and third states presented in the table clearly display two grains of the same orientation but with boundary atoms meeting "head-to-head." This is essentially an intermediate in the grains slipping on one another that is stabilized by the restrictions of the boundary conditions. Such states then represent a grain boundary that is stable, at least in X. When PFC simulations [26] are initialized at these states, the flow appears to be stable for thousands of steps then suddenly goes to the hexagonal lattice, meaning there are unstable directions in the rest of H 2 . Nevertheless, the fact remains that grain boundaries can be steady states. Table 1: Data for selected values of (ψ, β) on the exact steady states a near the numerical approximationā, obtained from the original candidate a. M = 20, 30, 40 for each parameter set respectively. The Morse index was verified in X. We write < when the number was numerically computed as 0. E 0 denotes the energy of the constant state.  Visualization Table 3 presents some steady states found for (ψ, β) = (0.5, 0.6), (N x , N y ) = (7, 4), ν = 1.01 and M = 60. In this regime, localized or coexistence patterns are observed in PFC simulations, some of which we can confirm to be steady states: note in particular the existence of a "single atom" state. We see here that the global minimizer cannot be of the four main ansatz. We observe two atoms states with different amplitudes and stability, highlighting the fact that the "linear" candidate is no longer appropriate as β increases and nonlinear effects begin to dominate the energy. Similar results have been obtained previously for a version of Swift-Hohenberg with broken ψ → −ψ symmetry, see [27,28]. Table 3: Data on steady states for (ψ, β) = (0.5, 0.6) and (N x , N y ) = (7,4). No count is provided because only a few trials were attempted.  Table 4 shows a selection of steady states found in the large β regime, (ψ, β) = (2.5, 20.0), (N x , N y ) = (4, 2), ν = 1.01 and M = 65. In this regime, the microscopic organization is lost as constant patches of phase form, with value close to ± √ β, meaning that the double well term of the PFC functional dominates the oscillation term. Table 4: Data on steady states for (ψ, β) = (2.5, 20.0) and (N x , N y ) = (4, 2).

Phase diagram with verified steady states
The framework allows us to construct a "rigorous" phase diagram for PFC. Here the adjective "rigorous" does not mean that we have identified the ground state; but rather that the respective candidate state has been rigorously verified in its parameter regime. To this end, one must construct a "patchwork" of (ψ, β) split in regions in which we have a proof that a given state is a global minimizer. For now, we restrict ourselves to proving that one of the steady states near the known candidate minimizers has lower energy than all other known steady states at given points. Further, our attempt is somewhat limited by the small domains we can access. Nevertheless, this construction is useful and does indicate rigorously where the candidates cannot be global minimizers. Our approach is as follows: we discretize the (ψ, β) parameter space to some desired accuracy and for each point, we test the four ansatz and several other candidates obtained from random initial coefficients. When one of the four ansatz has verified lower energy than the others, up to translational symmetries, we label that point accordingly and otherwise leave the point blank. Fig. 2 (a) shows the resulting diagram for small parameter values with (N x , N y ) = (4, 2), ν = 1.01, M = 20. At each point, 30 trials of the Newton iteration were tried and verified rigorously. Note that the points below β =ψ 2 could have been skipped since the constant state is known to be the global minimizer in that regime [10]. This diagram matches the one obtained in the appendices with linear stability analysis. Fig. 2 (b) shows the phase diagram near (ψ, β) = (0.5, 0.6) where localized patterns have been observed. The domain is the same size but M = 30 to accommodate the larger β. At each point, 15 trials were tried and verified, leading to points that have lower energy than the atoms or constant states.  This indeed shows the existence of a region where localized patterns are more energetically favourable. This region gives an estimate of the full coexistence region that ultimately cannot be made explicit without more refined techniques.

Rigorous results for two-mode PFC
As a final example, Table 5 shows three verified steady states for two-mode PFC with q = 1/ √ 2, (ψ, β) = (0.09, 0.025), (N x , N y ) = (12, 4), ν = 1.01 and M = 64. Note that here, L x = 2 √ 2πN x and L y = 2 √ 2πL y to fit the symmetry of the square lattice. The second state shows two grains slipping on each other; in contrast, especially to the result for hexagonal lattices, the third state is a grain boundary with non-zero misorientation. Here, the rectangular domains with Neumann boundary conditions can support the geometry of the square lattice at 0 • and 45 • rotations, so we can observe their coexistence. Since this result can be extended to larger domains by simple tiling operations, we conclude that straight grain boundaries can be steady states even in infinite domains where boundary conditions cannot "help" stabilizing such defects.
Moreover, this grain boundary was observed to be (numerically) stable in two-mode PFC simulations in the sense that small random perturbations of the phase field always converged back to the grain boundary state. This is not a rigorous proof of stability in H 2 , but it gives a good indication that grain boundaries are likely to be stable features in the PFC model.  More precisely, the connecting orbit leaves the unstable manifold of Ψ 1 and ends up in the stable manifold of Ψ 2 . Since the PFC equation is a gradient flow, there cannot exist non-trivial homoclinic connections so there is a natural hierarchy of steady states expressed through heteroclinic connections. This concept is extremely useful to "visualize" the energy landscape.
States with Morse index 0 are stable (for fixed parameters) and are thus at the bottom of the hierarchy. Those states with Morse index 1 have one unstable direction, so there are two distinct perturbations that lead away from the state. For states with Morse index 2, two unstable directions span infinitely many such perturbations, and so on. To detect connections, we propose to initialize a PFC flow near an unstable steady state offset by such perturbations. If the flow becomes close enough to another known steady states, we stop and propose a conjectured connection between the two steady states. This procedure often allows us to find unknown steady states: when the flow stagnates, the Newton iteration can be run and often converges in very few steps to a steady state that can be verified. Alternatively, we could check for inclusion in the target r * ball, but this is a very restrictive criterion that limits our numerical investigation, especially when obtaining connections to unstable states. We use the PFC scheme detailed in [26].
While we cannot for the moment prove such claims because "parameterizing" the infinite dimensional stable manifold of the unstable steady states is highly non-trivial, we are aware of some preliminary work in this direction [29]. That said, computer-assisted proofs of connecting orbits from saddle points to asymptotically stable steady states in parabolic PDEs are starting to appear [30,31,32].
We first consider the standard parameters (ψ, β) = (0.07, 0.025) and use the very small domain (N x , N y ) = (2, 1). This choice is made to ensure that the constant state has Morse index 2 in X to simplify the visualization. We find seven steady states: both possible translations of the atoms, stripes and donuts state, and the trivial constant state. Following the program described above, we can construct the "connection diagram" shown in Fig. 3 (a) with the arrows indicating that a connection was found from a state to the other. Note in particular that the stable stripes state to the right numerically decays into the appropriately shifted hexagonal lattices, but this is a slow process as the sine modes must grow out of numerical noise. This clearly shows that our method cannot be used to guarantee stability in H 2 because it cannot depend on translational shifts.
We also propose a visualization method for such diagrams shown in Fig. 3  (b). Take for example the constant state with its two unstable directions given by the coefficients a 0,2Ny and a Nx,Ny . We place the constant state at the origin and plot radial lines along linear combination of the unstable directions. The line length corresponds to the number of PFC steps needed to approach the target steady states. In addition, we can color the points along the line as a function of energy to indicate energetic relationships. A variant would be to show the energy as the z-component of a surface; essentially giving an indirect visualization of the energy landscape through 2D unstable manifolds. In particular, this diagram clarifies the relationships between the steady states. For instance, the stripes states are formed by adding the a 0,2Ny mode to the constant state while the donuts are combinations of the atoms and stripes states.
We now consider the localized patterns regime to illustrate these ideas with states of high Morse index. We do not attempt to build a higher dimensional visualization, but simply attempt to recover the "pathways" between the highly unstable hexagonal lattice with Morse index 20 towards stable steady states. This is visualized in the connection diagram of Fig. 4 (a) which includes a few states of Table 3. In (b), we plot the energy along the PFC flow starting from the index 2 state; this plot can be thought of as one of the rays in a diagram like Fig. 3 (b). Note that along the flow, the energy decreases in "steps" corresponding to changes in topology, i.e. the formation (or removal) of atoms. In this process, we could not verify that these intermediates are steady states since the Newton iteration always converged to the endpoint; we then suppose they are short-lived "metastable" states.
It is difficult to obtain perturbations that can flow to desired steady states, especially when they are unstable; see how only a few directions reach the Morse index 1 states in Fig. 3 (b). Indeed, unless "trivial" combinations of the   Figure 4: Connection diagram (a) where arrows represent a few of the connections found. The two hexagonal lattice states differ in their amplitude and stability. The vertical axis roughly indicates the energy while the numbers give the Morse indices. We could not obtain (nor disprove) a connection to the single atom state, indicated with the question mark. The connection labeled with a star is broken down in the energy plot to the right (b). These states appear to be metastable intermediates where the energy gradient becomes small and the evolution slows down considerably. The blue curve shows the energy as a function of time in arbitrary units, highlighting momentaneous "flats" in the evolution.
unstable direction happen to go to an unstable state, we are unlikely to find such connections numerically. Similarly, our attempts to find a perturbation that connects the starting lattice to the single atom state were unfruitful.

Parameter continuation for steady states
A verified steady state a for some parameter (ψ, β) is usually part of a family of steady states representing a "phase" of matter. In fact, the candidate minimizers defined in appendix A as functions of (ψ, β) approximate such families, or branches in the bifurcation diagram. In this context, we can construct such branches by starting at a known steady state, varyψ and find the closest steady state at this new parameter value.
Several verified techniques exist for following branches, see [33] and [21] for an application to Ohta-Kawasaki. We use non-verified pseudo-arclength continuation [34] inψ. Note that the unstable direction that is to followed is precisely given by the one corresponding to the "fixed" a 0,0 =ψ condition and this is one of the reasons that we chose to enforce this directly in the formulation of F . As a possible extension, 2D manifolds can be constructed in 2-parameter continuation when both parameters are allowed to vary, see [35]. Fig. 5 shows the norm (a) and offset energy (b) of the main ansatz at (ψ, β) = (0.07, 0.025) are plotted as functions ofψ. The domain is kept small with (N x , N y ) = (2, 1) to keep the bifurcation diagram as simple as possible. The atoms and donuts branches are actually the same since we can continue the branches through the folds atψ = ± 5/12β. This branch intersects the checkers state at β = 15ψ 2 and the constant and stripes states at β = 3ψ 2 . The energy plot (b) clearly shows that the donuts state is the "proper" hexagonal lattice forψ < 0. We note that varying β simply causes the branches to dilate. For example, we expect the 2D hexagonal steady states manifold to be a "conic" figure-eight.  Other "new" branches will appear for larger domains or higher β. In particular, Fig. 6 shows the atoms/donuts branch and the single atom branch in the localized patterns regime near (ψ, β) = (0.5, 0.6) with (N x , N y ) = (7,4). Again, (a) shows the L 2 norm and (b) shows the energy of the phase field as functions ofψ. The hexagonal lattice traces out its usual figure-eight pattern while the single atom (and other localized states in general) traces out a complicated looping path. Such branches illustrate the "snaking" phenomenon previously observed in modified Swift-Hohenberg equations that support such localized patterns, see [27] for example. We observe that the path loops on itself in one direction as the single atom evolves into a localized pattern with 9, 7 then 4 atoms before looping back with a 90 • rotation. In the other direction, the branch moves towards the transition between the hexagonal and constant states where it again loops back. This computation is difficult because the truncation must remain large and the pseudo-arclength step size must remain small; if the step size is larger than 0.0005, the branch breaks away towards the hexagonal lattice solution.

Conclusion
We surveyed the basic properties of the PFC equation as a dynamical system in the framework of rigorous numerics. Thanks to an application of the radii polynomial approach, we were able to verify the existence of steady states close to numerically computed approximations. This provided us with important verified information as to the behavior of the energy landscape, especially in terms of energetic relationships between steady states. We were also able to provide partial stability results with the caveat that they only applied to the cosine Fourier series. The Morse indices given were lower bounds in H 2 -thus those steady states with Morse index higher than 0 must be unstable in H 2 .
Such ideas were applied in various regimes of the PFC equation to verify that certain important patterns are steady states (as opposed to metastable intermediates) including single atoms, other localized patterns and grain boundaries. In particular, we showed that two-mode PFC supports a non-zero misorientation grain boundary steady state that we expect to be stable. We also showed the construction of the phase diagram with our fully nonlinear approach.
Finally, we used such results to further investigate the energy landscape through connections or orbits and through parameter continuation. Connections reveal the energetic and dynamical relationships between steady states, highlighting the behavior of unstable patterns as they reach states with lower energy. Continuation is especially useful to understand how the important states evolve across parameter space, highlighting the surprising behavior of the hexagonal lattice patterns and the snaking behavior of localized patterns.
Our work suggests several interesting directions for future work. On one hand, our connection results could be made rigorous with a technique to prove orbits from unstable to stable manifolds. This is a complicated problem because the stable manifold is infinite dimensional and special techniques must be applied to properly parameterize its "dominant" submanifold. On the other, our continuation results could also be made rigorous or extended to 2-parameter continuation to reveal more interesting behavior. Alternatively, parameter continuation could be applied to the domain size, for example to investigate problems in elasticity.

Appendices A PFC ansatz in 2D
PFC simulations can be classified according to a small number of regimes or ansatz that represent the (expected) global minimizer. The choice of such candidates is motivated by numerical experiments but can be obtained analytically from techniques such as linear stability analysis. Consider a periodic "single Fourier mode" phase field of the form ψ(x, y) =ψ + A 1 cos(y) + A 2 cos Inserting this ansatz into the PFC energy yields an expression E[A 1 , A 2 , A 3 ] that can be optimized in the three amplitudes. This procedure yields three main classes of states that are well-known in the PFC literature.
In addition, we also find a "checkers state" where A 1 and A 2 = A 3 are given as more complicated expressions. The two hexagonal lattices differ in energy: the positive choice is called the "donuts" state while the negative one is the "atoms" state. We can compare the energies directly to show that the checkers state is never optimal while the atoms state is more optimal than the donuts state for ψ > 0. When the coefficients of two candidates are equal, they represent the same regime; for example, at β = 3ψ 2 , the constant, stripes and donuts states are all ψ(x, y) =ψ. Similarly, the donuts and atoms states merge at β = 12/5ψ 2 . We can also compute when states have the same energy. Such behavior occurs on transition curves of the form β = αψ 2 .
We can construct the phase diagram in Fig. 7 by labeling with the expected global minimizer. We show in the main text that this "linear" description of PFC is a good approximation, at least for small β.

B Proof of the radii polynomial theorem
Proof. Consider the Newton operator T (a) = a − AF (a), then T : X → X and any fixed point a of T is a zero of F because A is injective. The derivative of T is bounded and also Fréchet differentiable since we have for any x ∈ X Figure 7: Phase diagram (a) and detail (b) constructed by comparing the optimal energy of the four ansatz. The curves β = 3ψ 2 and β = 12/5ψ 2 correspond respectively to the curves on which the amplitude of the stripes and lattice states become complex. The stripes-atoms (blue-yellow) transition curve is β ≈ 20.22ψ 2 while the constant-atoms (red-yellow) transition curve is β = 37/15ψ 2 . Donuts are never optimal forψ > 0. Checkers are never optimal.
as DF is the (bounded linear) Fréchet derivative of F . Now suppose p(r 0 ) < 0 for some r 0 > 0, then the radii polynomial in the main text gives Let a ∈ B r0 (ā), we can use the previous inequality to bound Pairing this with the mean value inequality for T , hence T maps B r0 (ā) to its interior thanks to the strict inequality. Similarly for x, y ∈ B r0 (ā), is a contraction with constant κ = (Z 0 +Z 1 +Z 2 (r 0 )r 0 ) < 1 and the Banach fixed-point theorem gives the result.

C Computation of the radii polynomial bounds
In the following calculations, we will use usual results such as ||Qb|| 1,ν , and the following proposition to compute the norm of operators on X: Proof. Let b ∈ X, then Qb can be decomposed as the action of the first finite block onto b σ for σ ∈ U and infinitely many diagonal terms c σ b σ for σ / ∈ U . The norm of Q can then be written as the sum of two disjoint positive sums using the triangle inequality: The second term is bounded by the trivial bound C = sup α / ∈U |c α |, extracting the norm of b over N \ U . Similarly, the first term is bounded by The norm of Q in B( 1 ν ) is the supremum of the previous norm over all b with unit norm, therefore the triangle inequality gives the result.
The sharper result max{K, C} can be obtained by noting that the sums act on different subspaces of 1 ν . The estimate also allows us to compute the norm of finite tensors by letting the c α vanish. Let us now apply these results to compute the four necessary bounds.

C.1 Z 0 bound
This bound is the easiest since by construction, A and A † are approximate inverses up to numerical inversion errors. We then have which can be evaluated using proposition 1.

C.2 Y 0 bound
We must compute AF (ā) so that most terms will be given by the finite product A (M ) F (M ) (ā). There still remain some non-zero convolution coefficients in the full F (ā): sinceā (M ) has M + 1 coefficients in each dimension, L α (ā * ā * ā) α will have 3M + 1 non-zero coefficients in each dimension. These are multiplied by the appropriate L −1 α γ −1 α , resulting in

C.3 Z 2 bound
To compute the Z 2 bound, let b, h ∈ X with ||h|| 1,ν = 1 and consider first the effect of DF (b) on h, Fix r > 0 and let b =ā + R where ||R|| 1,ν ≤ r, we then have that Note that the initial factor L α can instead be represented by the diagonal operator defined by Λ α,σ = L α δ α1−σ1 δ α2−σ2 . Then, using the fact that the convolution on X is a Banach algebra, r where the norm of AΛ is computed using proposition 1 with the bound Γ = max(γ −1 0,M +1 , γ −1 M +1,0 ) on the diagonal terms. 2 This computation works for any r > 0, so we have 2 r .

C.4 Z 1 bound
For the final bound, we now consider the action of A † on the same vector h: Let η be the tail of h, i.e. the vector with the same entries as h outside of U and 0 on U . We then have: Consider now the action of A on the difference above. The first block of the difference will be multiplied by the inverse of G while the tail will be multiplied by the appropriate L −1 α γ −1 α , thus using L 0,0 = 0 and the triangle inequality for the first term. Let φ ∈ X be such that |(ā * ā * η) α | ≤ φ α whenever α ∈ U and 0 otherwise. Using the Banach algebra property and the bound Γ ≥ |γ σ | −1 to bound the infinite sum, we have which can be computed numerically once the (finitely many) φ α have been obtained. To compute them, we now shift for a moment to Z 2 and extend all vectors appropriately. Now let q =ā * ā, then where V (α), W ⊂ Z 2 are the regions over which q α−σ and η are non-zero respectively. Since q τ = 0 whenever either |τ 1 | or |τ 2 | is larger than 2M and we only need |α| ∈ U , we obtain the overestimate that V (α) ⊂ {−3M, ..., 3M } 2 . Further, η σ must vanish for |σ| ∈ U so for α ∈ U (and 0 otherwise), which completes the computation of Z 1 .

D Energy computation and real space norms
The notion of closeness between a andā extends to their energies. For simplicity, we only handle the basic PFC energy We will write without relabeling E[ψ (a) ] = E[a] when ψ (a) is the phase field corresponding to the Fourier coefficients a. Taking the average of a Fourier series returns its constant mode such that In the context of the radii polynomial approach, let t = a −ā such that We can simplify (ā + Lā) * (Lt) = Lā + L 2ā * t by integrating by parts and using the periodic boundary conditions; for example, We now use the fact that ||t|| 1,ν < r * . In addition, we can overestimate |a 0,0 | ≤ ||a|| 1,ν and use the Banach algebra property to obtain the following bound: The term strictly in t has been left as a convolution because it is necessary to control the growth of L α with the r * bound directly. To do so, we have which is another way to obtain the previous integration by parts result. Now, the 1 ν norm of t is bounded by r * ; each member of the sum satisfies the inequality W α |t α |ν |α| < r * . Overestimating W α ≥ 1, we can write where ρ = 1/ν 2 < 1. To evaluate this sum, we can compute the polynomial geometric series which all converge for ρ < 1. Note that the sums can be evaluated by differentiating ρ jx with respect to x = 1. The terms in L α can then be expanded and written in such a fashion and assuming that |S| is finite, the sums can be split and separated. We have and similarly, Putting everything together, we arrive at the bound |E[ a] − E[ā]| ≤ |1 − β| ||ā|| 1,ν + ||ā|| 3 1,ν + 2||Lā|| 1,ν + ||L 2ā || 1,ν r * + 1 2 2S 1 + S 2 + |1 − β| + 3||ā|| 2 1,ν r 2 * + ||ā|| 1,ν r 3 * + 1 4 r 4 * which can now be computed numerically. This bound depends strongly on ν because of its influence on r * and the growth of S 1 and S 2 . In principle, one could find an optimal ν that ensures the bound is as small as possible for a given a and a fixed β.
When the numerical errors associated to the E[ā] computations are added to the energy bound, both computed with interval arithmetic, we obtain an interval that is guaranteed to contain the energy of a itself. In particular, this allows us to prove which of two steady states is more optimal 3 strictly from numerical computations.
The previous computations illustrate some techniques that allow us to estimate the norm of ψ (t) (the phase field corresponding to a −ā) in terms of r * . For example, provides a pointwise estimate on the value of the exact steady state. Further, a simple calculation shows that Parseval's identity holds on X; i.e. ||ψ (a) || 2 L 2 (Ω) = |Ω| α∈N W α a 2 α . We can then bound the L 2 norm of derivatives, for instance ||∇ 2 ψ (t) || 2 L 2 = |Ω| α∈N L 2 α t 2 α ≤ S 2 |Ω|r 2 * using Eq. (2). Similar results can be built for the lower norms, thus providing an estimate of the form ||ψ (t) || H 2 ≤ C(Ω, ν)r * for some constant C that could be computed if necessary. Note the implicit dependence on the state itself and (ψ, β) through r * . Combined with the L ∞ bound, this shows that as long as ν > 1, the exact steady state will be in H 2 and can differ from the numerical candidate by at most r * at any point in Ω. The constant C may be large, but it does not affect the pointwise agreement; this is sufficient control for our numerical investigation.

E Stability in X
To complete the analysis of a given steady state a, we can characterize its stability in X. This is powerful because even linear results are mostly limited to trivial states but a major limitation is that this does not transfer to H 2 because X is restricted to the cosine series. Suppose we have a steady state a ∈ B r (ā) for a verified radius r, then stability is controlled by the spectrum of DF ( a). This spectrum is real because we are in the context of a gradient flow; this can be seen directly from the definition of A † and DF which are symmetric on interchanging indices. Assuming there are no zero eigenvalues, the positive and negative ones define the unstable and stable manifolds respectively. A steady state with only strictly negative eigenvalues is said to be stable.
While only the approximation A † is known in practice, it has the same signature as DF ( a) itself; i.e. they have exactly as many strictly positive or strictly negative eigenvalues. We compute the spectrum of A † in two parts. The (M + 1) 2 eigenvalues of the finite block G = DF (M ) (ā) can be computed numerically and verified using interval arithmetic routines. 4 Most eigenvalues can be unequivocally assigned a sign, but some may be identically 0 or closer to 0 than the available precision. Stability cannot be ascertained in such cases, but assumingā was verified using the radii polynomial approach, G must be sufficiently well-conditioned so its eigenvalues cannot be so small.
In the tail, the eigenvalues are simply equal to the diagonal terms L α γ α with α / ∈ U . Thankfully, L is strictly negative for α = (0, 0) and γ is strictly positive as long as M is sufficiently large. We then have ≤ Z 0 + sZ 1 + sZ 2 (r)r ≤ Z 0 + Z 1 + Z 2 (r)r < 1 as in the proof of the radii polynomial approach. Since I − AH s is a bounded operator with norm less than 1, AH s is itself invertible. AH s and thus H s cannot have a zero eigenvalue so that its signature must stay constant for all s. This shows that A † and DF ( a) have the same signature and this is in fact true over the ball of radius r * aroundā.