Microscopic Patterns in the 2D Phase-Field-Crystal Model
Gabriel Martine-La Boissoniere, Rustum Choksi, Jean-Philippe Lessard
MMicroscopic Patterns in the 2DPhase-Field-Crystal Model
Gabriel Martine-La Boissoni`ere ∗ Rustum Choksi † Jean-Philippe Lessard ‡ Department of Mathematics and Statistics, McGill University,Montr´eal, QC, Canada
Abstract
Using the recently developed theory of rigorously validated numer-ics, we address the Phase-Field-Crystal (PFC) model at the microscopic(atomistic) level. We show the existence of critical points and local min-imizers associated with “classical” candidates, grain boundaries, and lo-calized patterns. We further address the dynamical relationships betweenthe observed patterns for fixed parameters and across parameter space,then formulate several conjectures on the dynamical connections (or or-bits) between steady states.
The Phase-Field-Crystal (PFC) model introduced in [1] is a gradient systemcapable of modeling a variety of solid-state phenomena. In its simplest form,the PFC energy can be written as E [ ψ ] = − (cid:90) Ω (cid:0) ∇ ψ + ψ (cid:1) + 14 (cid:0) ψ − β (cid:1) defined on phase-fields ψ ∈ H (Ω) satisfying the phase constraint ¯ ψ = − (cid:90) Ω ψ = 1 | Ω | (cid:90) Ω ψ . The parameter β represents inverse temperature such that β = 0 models maxi-mum disorder. Coupled with this energy is its conservative H − gradient flow ∗ [email protected] † [email protected] ‡ [email protected] a r X i v : . [ m a t h . NA ] F e b igure 1: Left: Details of a grain boundary appearing in a PFC simulation(taken from [3]). Right: Grain boundary network from a PFC simulation (takenfrom [4]). Within each grain is a hexagonal lattice of atoms with a particularorientation.which entails the sixth-order PFC equation ψ t = ∇ (cid:16)(cid:0) ∇ + 1 (cid:1) ψ + ψ − βψ (cid:17) . Note that the PFC model shares its energy with the Swift-Hohenberg equa-tion [2], which is simply the L gradient flow of E . From linear stability analysisapplied to single Fourier mode Ansatz, we find three main candidate global mini-mizers that divide parameter space, see the appendices. In the hexagonal latticeregime, 2D-simulations of the PDE starting with random noise quickly produceatoms that arrange into small patches of hexagonal lattices with random ori-entations. These patches grow and interact with each other, forming grains ofhexagonal lattices of atoms with a particular orientation. The morphology andevolution of these grains have features resembling those in polycrystalline ma-terials (cf. Figure 1). In particular, it has recently been shown that statisticsof many of experimentally observed (universal) grain boundary distributionsare accurately captured by data amassed from simulations of this simple PFCequation [5, 4]. While here we will mostly work with this vanilla PFC formula-tion, we note that a family of PFC-like equations can be derived from Density-Functional-Theory [6] to obtain more complicated models capable of simulatingeutectic and dendritic solidification [7] and graphene structures [8, 9].In this article, we address the PFC model and its steady states at the “mi-croscopic” level - the local atomic arrangement. We believe that such an in-vestigation of microscopic pattern-formation capabilities of PFC is not only ofmathematical interest but is also necessary to construct “designer” models forpolycrystalline behaviour. For example, varying the parameters in the energylead to more complicated states than simple lamellar and hexagonal. Theseinclude 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 toprove any theorem on the exact nature of steady states, local and global min-imizers to this diffuse interface problem. What exists in the physics literature2s numerical simulations, standard linear stability analysis, and Ansatz-drivenenergy comparisons. The recently developed theory of rigorously validated nu-merics (cf. [11, 12, 13, 14, 15]) now provides a powerful new tool to bridge whatcan be observed numerically with rigorous statements on pattern morphology.In a nutshell this approach can be summarized as follows: Given an approximatesteady state, we use the Contraction Mapping Theorem to imply the existenceand local uniqueness of an exact steady state within a controlled distance ofthe approximation. This notion of closeness is strong enough to imply furtheruseful results, including closeness in energy and stability results. In this paperwe use this new approach to address the following aspects of the PFC model: • Are the “classical” candidates obtained from linear stability analysis closeto actual local minimizers? • Are the stable yet complicated patterns observed numerically indeed crit-ical points in the PFC energy landscape? For example, are grain bound-aries steady states or simply metastable states? • What are the dynamical relationships between the observed patterns forfixed parameters and across parameter space?Based upon our results we formulate several conjectures on the connections (ororbits) between steady states. Taken as a whole, our work presents the firststep into a rigorous analysis of the rich PFC energy landscape.The outline of this paper is as follows. We first setup the PFC equation inFourier space and discuss the application of the framework of rigorous compu-tations. We then verify the existence of important steady states of the PFCequation, including localized patterns and grain boundaries. With these statesin hand, we address the energy landscape of PFC with a discussion on con-jectures for connections (or connecting orbits) between steady states. Finally,we presents results in one-parameter numerical continuation to outline someinteresting features of the bifurcation diagram of PFC.
We begin by writing the equation ψ t = 0 in Fourier space to obtain a coupledsystem of equations for the Fourier coefficients of steady states. We will beslightly more general and consider functionals of the form E [ ψ ] = − (cid:90) Ω
12 ( Kψ ) + 14 ( ψ − β ) where K is a linear differential operator K acting on elements of a suitablefunction space. In particular, K = (cid:40) ∇ + 1 for the basic “one-mode” PFC model( ∇ + 1)( ∇ + q ) for the “two-mode” PFC model [16]3here q is the secondary wavelength of two-mode PFC. Taking the H − gradientflow of E , we obtain the PFC-like equation ψ t = ∇ (cid:0)(cid:0) K − β (cid:1) ψ + ψ (cid:1) .For simplicity, we let Ω be the rectangular domain [0 , L x ] × [0 , L y ] withperiodic boundary conditions. We let L x = 4 π √ N x , L y = 4 πN y where N x , N y ∈ N are the number of atoms lined up in the x, y -axes. The mainparameters of the problem are then ( ¯ ψ , β ) and the domain size is given by( N x , N y ).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 ofequations of the form ( a α ) t = F α ( a ) thanks to orthogonality. The steady statesmay then be found numerically by solving F ( a ) = 0 up to some truncationorder 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 bebroken. The simplest way to do so in this context is to also enforce Neumannboundary conditions. It is convenient to write a α = a α ,α so that the symmetryand reality conditions become a | α | , | α | ∈ R .This choice allows us to simplify a complex Fourier series into the cosineexpansion ψ ( x, y ) = (cid:88) α ∈ Z a α exp (cid:18) πi α xL x (cid:19) exp (cid:18) πi α yL y (cid:19) = (cid:88) α ∈ N W α a α cos (cid:18) πα L x x (cid:19) cos (cid:18) πα L y y (cid:19) where W is a weight matrix defined by W α = α = (0 , α = 0 , α (cid:54) = 0 or α (cid:54) = 0 , α = 04 otherwise . The Fourier coefficients of ∇ ψ are given by the elementwise product L α a α where L α = − (cid:32)(cid:18) πα L x (cid:19) + (cid:18) πα L y (cid:19) (cid:33) is the Fourier representation of the Laplacian. Inserting these expressions intothe PFC equation and equating Fourier modes, we obtain( a α ) t = F α ( a ) = L α ( γ α a α + ( a ∗ a ∗ a ) α )where ∗ denotes the discrete convolution and the linear terms combining K and β are γ α = (cid:40) ( L α + 1) − β for PFC( L α + 1) (cid:0) L α + q (cid:1) − β for two-mode PFC . ,
0) Fourier component picks out the average phase so it isfixed to ¯ ψ : this is consistent with ( a , ) t = 0 thanks to L , = 0. To keep trackof the phase constraint directly in F , we replace its first trivial component by F , = a , − ¯ ψ , resulting in: F α ( a ) = (cid:40) a , − ¯ ψ if α = (0 , L α ( γ α a α + ( a ∗ a ∗ a ) α ) otherwise . The operator F then represents the PFC dynamics in the sense that itszeros correspond to steady states of the PFC equation. A numerical advantageof the reduced expansion is that we effectively only have to compute a quarterof the full Fourier series. Obviously, this means we are not treating PFC infull generality over H and will have to address this later. As an aside, theequivalent F for Swift-Hohenberg is simply − ( γ α a α + ( a ∗ a ∗ a ) α ) hence its(0 ,
0) entry is nonzero and average phase is not conserved.
We present a brief overview of the recent framework of rigorously validatednumerics 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 suitableinverse to the derivative DF ( a ). On the one hand, if T is a contraction on aclosed ball, the contraction mapping theorem gives the existence and uniquenessof 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 ¯ a for which F (¯ a ) ≈ T is a contractionaround ¯ a , then we immediately have the existence of an exact steady state (cid:101) a close to ¯ a in an appropriate metric. This relationship is made clear by the radiipolynomial theorem , so-called for reasons that will become clear shortly. Toillustrate 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 applicationto 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 operatorsfrom X to Y , B ( X ) = B ( X, X ) and B r ( a ) ⊂ X for the open ball of radius r around a ∈ X . Theorem 1.
Consider Banach spaces
X, Y , a point ¯ a ∈ X and let A † ∈ B ( X, Y ) , A ∈ B ( Y, X ) . Suppose F : X → Y is Fr´echet differentiable on X and A is injective. In addition, suppose || AF (¯ a ) || X ≤ Y || I − AA † || B ( X ) ≤ Z || A ( DF (¯ a ) − A † ) || B ( X ) ≤ Z || A ( DF ( b ) − DF (¯ a )) || B ( X ) ≤ Z ( r ) r ∀ b ∈ B r (¯ a )5 here Y , Z , Z are positive constants and Z is a positive polynomial in r > .Construct the radii polynomial p ( r ) = Z ( r ) r − (1 − Z − Z ) r + Y . (1) If p ( r ) < for some r > , then there exists a unique (cid:101) a ∈ B r (¯ a ) for which F ( (cid:101) a ) = 0 . The proof of this formulation is given in appendix B, where we show acorrespondence between the sign of the radii polynomial and the contractionconstant of T : if r can be found, T is a contraction and the Newton iterationstarting at ¯ a must converge to some (cid:101) a . This proves not only the existence of theexact steady states but also gives control on its location in X with respect toa known point. In practice, one finds an interval [ r ∗ , r ∗ ] of radii for which p ( r )is negative; r ∗ > maximum distance between ¯ a and (cid:101) a while r ∗ > r ∗ gives the minimum distance between ¯ a and another zero of F . The zeros of F must therefore be isolated for consistency.Each bound may be understood intuitively: Y being small indicates that ¯ a is a good approximation of (cid:101) a while Z being small indicates that A † is a goodapproximation for DF (¯ a ), and so on. These bounds may be simplified ana-lytically but must necessarily be computed numerically. Therefore, we ensurethat our numerical computations go in the same direction as the required in-equalities by using interval arithmetic [23], a formalized approach to deal withnumerical errors. We used the interval arithmetic package INTLAB for MAT-LAB, see [24, 25], to ensure that the radii polynomial approach is numericallyrigorous.This approach allows us to prepare numerical tools that can both find can-didate steady states and compute the radii r ∗ , r ∗ if they exist. If so, we imme-diately have a proof that this candidate provides a good handle on an actualsteady state of the PFC equation. Let us now apply these ideas to PFC by first computing DF and the Newtonoperator. Let σ represent the differentiation indices applied to F α . The deriva-tive of F , is 1 if σ = (0 ,
0) and 0 otherwise, so we use the Kronecker deltanotation to write ∂ a σ F , = δ σ δ σ . For other values of α , the linear terms similarly give ∂ a σ ( L α γ α a α ) = L α γ α δ σ − α δ σ − α . The derivative of the nonlinear triple convolution can be computed by dif-ferentiating with respect to all four a α identified by symmetry. This algebraic6omputation is somewhat tedious but the result can be written succinctly as ∂ a σ ( a ∗ a ∗ a ) α = 3 W σ (cid:16) ( a ∗ a ) | α + σ | , | α + σ | + ( a ∗ a ) | α + σ | , | α − σ | + ( a ∗ a ) | α − σ | , | α + σ | + ( a ∗ a ) | α − σ | , | α − σ | (cid:17) so that the full derivative of F is:[ DF ] σ,α ( a ) = ( ∂ a σ F σ )( a )= (cid:40) δ σ δ σ if α = (0 , L α ( γ α δ σ − α δ σ − α + ∂ a σ ( a ∗ a ∗ a ) α ) otherwise .a, F and the convolutions may be viewed as infinite matrices whose “top-left” 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 truncatedto order M such that a σ = 0 whenever either σ or σ is greater than M . Thisresults in the ( M + 1) matrices a ( M ) , F ( M ) while the derivative becomes the( M + 1) DF ( M ) . Note that the k -convolution of a ( M ) has support kM by definition.We now introduce the Banach space framework. Let ν > (cid:96) ν ( Z )as the space of sequences a α with finite norm || a || ,ν = (cid:88) α ∈ Z | a α | ν | α | = (cid:88) α ∈ Z | a α | ν | α | + | α | . The restriction of (cid:96) ν ( Z ) using the symmetry condition is X = (cid:8) a ∈ (cid:96) ν ( Z ) (cid:12)(cid:12) a α = a | α | , | α | (cid:9) over which the norm simplifies to || a || ,ν = (cid:88) α ∈ N W α | a α | ν | α | = (cid:88) α ∈ N | a α | ν α where ν α is a weight matrix that forces the fast exponential decay of the Fouriercoefficients. The space ( X, || · || ,ν ) can easily be shown to be Banach and the2D discrete convolution forms a Banach algebra over it, immediate results fromthe triangle inequality and the fact that ν > a, (cid:101) a ∈ X have the same meaning as before, with ¯ a = 0 outside of U = { , , ..., M } thanks to the truncation. Let G = DF (¯ a ) ( M ) and denote by A ( M ) the numerical inverse of G . We define approximate operators A † , A as A † α,σ = G α,σ if α, σ ∈ UL α γ α if α = σ, α ∈ N \ U A α,σ = A ( M ) α,σ if α, σ ∈ UL − α γ − α if α = σ, α ∈ N \ U G or its inverse paired withthe linear terms L α γ α as the main “diagonal” of the second block. If G is an7nvertible matrix, 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 operatorproducts against A are bounded on X thanks to the fast decay of L − α γ − α .Thus, we say that A “lifts” the regularity of the other operators back to X ,allowing statements such as T : X → X or ADF (¯ a ) ∈ B ( X ).We show in appendix C how to simplify the bounds into expressions thatcan be evaluated numerically. This allows us to write down the radii polynomial p ( r ) = Z ( r ) r − (1 − Z − Z ) r + Y , noting that Z ( r ) = Z (0)2 + Z (1)2 r , hence thepolynomial is cubic with non-negative coefficients except for maybe the linearterm. We have p (0) > p (cid:48) (0) = Z + Z − p ( r ) → ∞ for large r . Asa consequence, if p is strictly negative for some positive r , there must existexactly two strictly positive roots r ∗ < r ∗ defining the interval where the proofis applicable. When this is satisfied, the radii polynomial theorem gives that There exists an exact solution (cid:101) a of F ( a ) = 0 in B r ∗ (¯ a ). This solution is unique in B r ∗ (¯ a ).Thus, when the radii polynomial is computed using interval arithmetic andhas exactly two real non-negative roots, the zero computed numerically with theNewton iteration is close to an actual steady state of the PFC equation. Note theimportant fact that the ball is in X so a priori, only the Fourier coefficients arecontrolled. Thanks to ν > in value between the phase fields corresponding to ¯ a and (cid:101) a is at most r ∗ . Further, we show in appendix E that the stability of (cid:101) a in X is controlledby the eigenvalues of G . It is important to observe that this matrix will alwayshave a positive eigenvalue because of the trivial condition F , = a , − ¯ ψ . Thisis not indicative of instability in the context of the H − gradient flow because a , is fixed. We shall see later that this unstable direction can be used tocompute a branch of solutions in parameter continuation. For now however,we call the number of positive eigenvalues, minus 1, the Morse index of (cid:101) a ,indicating how many unstable directions are available to a given steady statefor fixed parameters.The procedure to numerically investigate the steady states of the PFC equa-tion is as follows: • Starting from a given initial condition, the Newton iteration is run untilit converges up to numerical precision. • Then, the radii polynomial of the numerical guess is computed and itsroots are tested. The numerical method will fail if G is almost singular, so this is the case in practice. If the proof succeeds, we can characterize an exact steady state in value, inenergy and compute its stability in X . The parameters ( M, ν ) can be ad-justed until the proof succeeds with a trade-off between the computationaleffort and closeness in X . We now have a complete framework for finding verified steady states alongwith their energetic and stability properties. This allows us to understand thebehavior of the PFC system for a given choice of ( ¯ ψ, β ), with three importantcaveats: • We cannot guarantee that we have found all steady states and therefore the global minimizer. Indeed, we may only hope to cover a reasonableportion 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 be-tween two steady states. In particular, large domains and large β increasethe contribution of high frequency Fourier modes, hence the truncationorder can become large even for domains containing only 100 atoms. Thislimits our results to small domains so our analysis is “small scale” in na-ture. • The Neumann boundary conditions restrict us to a “quadrant” of H .While the existence of a steady state, the energy bound and instability obviously extend to H , stability does not as there may be unstable di-rections in the other three Fourier series that are missed by the currentmethod.For the last point, we sometimes observe that translational shifts have adifferent 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 . The candidate global minimizers (constant, stripes, atoms and donuts states)introduced in appendix A have trivial Fourier coefficients by construction, givenby Constant: a , = ¯ ψ Stripes: a , = ¯ ψ a , N y = 12 A s Hexagonal: a , = ¯ ψ a N x ,N y = 12 A h a , N y = 12 A h 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 startingat the atoms state b constructed for ( ¯ ψ, β ) = (0 . , . N x , N y ) = (4 , M = 20. The Newton iteration was used to obtain ¯ b for which the radiipolynomial was tested with ν = 1 .
05, resulting in r ∗ = 1 . · − and r ∗ =6 . · − . The (cid:96) ν distance between b and ¯ b is 1 . · − , indeed smaller than r ∗ .The difference b − ¯ b is mainly captured by new Fourier modes: we find thatthe main Fourier coefficients b N x ,N y = b , N y = − . · − differ by 1 . · − while the largest new Fourier modes are b , = b , = − . · − . Moreover,the distance in the (numerical) sup norm between the two phase fields is ap-proximately 4 . · − which is again smaller than the (cid:96) ν distance, consistentwith the L ∞ bound.This approach was repeated for the other candidates and for a few otherchoices of the PFC parameters in the hexagonal regime, with truncation ad-justed to β . The results are presented in Table 1, showing that such simple can-didates capture well the leading behavior. Note that the agreement decreaseswith increasing β : compare the size of || a − ¯ a || ,ν to || ¯ a || ,ν . The Newton iteration can detect new steady states regardless of stability asit is based on criticality instead of minimality. This allows us to find steadystates that are observed only momentarily or even locally during a PFC sim-ulation. Table 2 presents a few of the 28 distinct steady states found for( ¯ ψ, β ) = (0 . , . N x , N y ) = (8 , ν = 1 .
05 and M = 40. Starting atrandom initial coefficient matrices, the Newton iteration converges in 15 to 50steps. The four main ansatz were also explicitly tested, as only the atoms statecould be reached from random initial conditions.Note that the energy of the exact steady states can be compared from Ta-ble 2: for instance, the energy of the exact atoms state is bounded away fromthe others so it is guaranteed to be the best candidate global minimizer out ofthe observed steady states at the current parameter values.The second and third states presented in the table clearly display two grainsof the same orientation but with boundary atoms meeting “head-to-head.” Thisis essentially an intermediate in the grains slipping on one another that is stabi-lized by the restrictions of the boundary conditions. Such states then representa grain boundary that is stable, at least in X . When PFC simulations [26] areinitialized at these states, the flow appears to be stable for thousands of stepsthen suddenly goes to the hexagonal lattice, meaning there are unstable direc-tions in the rest of H . Nevertheless, the fact remains that grain boundariescan be steady states . 10able 1: Data for selected values of ( ¯ ψ, β ) on the exact steady states (cid:101) a nearthe numerical approximation ¯ a , obtained from the original candidate a . M =20 , ,
40 for each parameter set respectively. The Morse index was verified in X . We write < (cid:15) when the number was numerically computed as 0. E denotesthe energy of the constant state.Ansatz ( ¯ ψ, β ) || ¯ a || ,ν || a − ¯ a || ,ν r ∗ r ∗ E [¯ a ] − E | E [¯ a ] − E [ (cid:101) a ] | MorseindexConstant (0 . , . . < (cid:15) . · − . · − . · − . · − . , .
5) 0 . < (cid:15) . · − . · − . · − . · − . , .
0) 0 . < (cid:15) . · − . · − < (cid:15) . · − . , . . . · − . · − . · − − . · − . · − . , .
5) 1 . . · − . · − . · − − . · − . · − . , .
0) 1 . . · − . · − . · − − . · − . · − . , . . . · − . · − . · − − . · − . · − . , .
5) 1 . . · − . · − . · − − . · − . · − . , .
0) 2 . . · − . · − . · − − . · − . · − . , . . . · − . · − . · − − . · − . · − . , .
5) 0 . . · − . · − . · − − . · − . · − . , .
0) 1 . . · − . · − . · − − . · − . · − ψ, β ) = (0 . , . N x , N y ) = (8 , r ∗ r ∗ E [¯ a ] − E | E [¯ a ] − E [ (cid:101) a ] | Morse index Count1 . · − . · − − . · − . · − . · − . · − − . · − . · − . · − . · − − . · − . · − . · − . · − − . · − . · − . · − . · − − . · − . · − Table 3 presents some steady states found for ( ¯ ψ, β ) = (0 . , . N x , N y ) =(7 , ν = 1 .
01 and M = 60. In this regime, localized or coexistence patternsare observed in PFC simulations, some of which we can confirm to be steadystates: note in particular the existence of a “single atom” state. We see herethat the global minimizer cannot be of the four main ansatz. We observe two atoms states with different amplitudes and stability, highlighting the fact thatthe “linear” candidate is no longer appropriate as β increases and nonlineareffects begin to dominate the energy.Similar results have been obtained previously for a version of Swift-Hohenbergwith broken ψ → − ψ symmetry, see [27, 28].12able 3: Data on steady states for ( ¯ ψ, β ) = (0 . , .
6) and ( N x , N y ) = (7 , r ∗ r ∗ E [¯ a ] − E | E [¯ a ] − E [ (cid:101) a ] | Morse index9 . · − . · − − . · − . · − . · − . · − − . · − . · − . · − . · − − . · − . · − . · − . · − − . · − . · − . · − . · − < (cid:15) . · − . · − . · − . · − . · − β regime Table 4 shows a selection of steady states found in the large β regime, ( ¯ ψ, β ) =(2 . , . N x , N y ) = (4 , ν = 1 .
01 and M = 65. In this regime, the micro-scopic organization is lost as constant patches of phase form, with value closeto ±√ β , meaning that the double well term of the PFC functional dominatesthe oscillation term. 13able 4: Data on steady states for ( ¯ ψ, β ) = (2 . , .
0) and ( N x , N y ) = (4 , r ∗ r ∗ E [¯ a ] − E | E [¯ a ] − E [ (cid:101) a ] | Morse index1 . · − . · − − . . · − . · − . · − − . . · − . · − . · − − . . · − . · − . · − − . . · − The framework allows us to construct a “rigorous” phase diagram for PFC.Here the adjective “rigorous” does not mean that we have identified the groundstate; but rather that the respective candidate state has been rigorously verifiedin 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 globalminimizer. For now, we restrict ourselves to proving that one of the steadystates near the known candidate minimizers has lower energy than all other known steady states at given points . Further, our attempt is somewhat limitedby the small domains we can access. Nevertheless, this construction is usefuland does indicate rigorously where the candidates cannot be global minimizers.Our approach is as follows: we discretize the ( ¯ ψ, β ) parameter space tosome desired accuracy and for each point, we test the four ansatz and severalother candidates obtained from random initial coefficients. When one of the fouransatz 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 , ν = 1 . M = 20. At each point, 30 trials of the Newton iteration were triedand verified rigorously. Note that the points below β = ¯ ψ could have beenskipped since the constant state is known to be the global minimizer in thatregime [10]. This diagram matches the one obtained in the appendices withlinear stability analysis.Fig. 2 (b) shows the phase diagram near ( ¯ ψ, β ) = (0 . , .
6) where localizedpatterns have been observed. The domain is the same size but M = 30 toaccommodate the larger β . At each point, 15 trials were tried and verified,leading to points that have lower energy than the atoms or constant states.14 a) (b) Figure 2: Phase diagram for small parameter values (a) and for the localizedpatterns regime. (b) All points are prepared by rigorously verifying that theexact steady state around the ansatz have lower energy than all other observedsteady states, up to translational shifts. Colored regions are filled in to guidethe eye. The curves show the condition for the energy of the basic ansatz to beequal.This indeed shows the existence of a region where localized patterns are moreenergetically favourable. This region gives an estimate of the full coexistenceregion that ultimately cannot be made explicit without more refined techniques.
As a final example, Table 5 shows three verified steady states for two-mode PFCwith q = 1 / √
2, ( ¯ ψ, β ) = (0 . , . N x , N y ) = (12 , ν = 1 .
01 and M = 64.Note that here, L x = 2 √ πN x and L y = 2 √ πL y to fit the symmetry of thesquare lattice. The second state shows two grains slipping on each other; incontrast, especially to the result for hexagonal lattices, the third state is a grainboundary with non-zero misorientation. Here, the rectangular domains withNeumann boundary conditions can support the geometry of the square latticeat 0 ◦ and ◦ rotations, so we can observe their coexistence. Since this resultcan be extended to larger domains by simple tiling operations, we conclude thatstraight grain boundaries can be steady states even in infinite domains whereboundary conditions cannot “help” stabilizing such defects.Moreover, this grain boundary was observed to be (numerically) stable intwo-mode PFC simulations in the sense that small random perturbations of thephase field always converged back to the grain boundary state. This is nota rigorous proof of stability in H , but it gives a good indication that grainboundaries are likely to be stable features in the PFC model.15able 5: Data on steady states for ( ¯ ψ, β ) = (0 . , . N x , N y ) = (12 , q = 1 / √ E is the energy of the constantstate for two-mode PFC. E [¯ a ] is listed for comparison purposes but it is notrigorously bounded.Visualization r ∗ r ∗ E [¯ a ] − E Morse index2 . · − . · − − . · − . · − . · − − . · − . · − . · − − . · − Suppose Ψ , Ψ represent two steady states, we say that there is a connection(or a connecting orbit) from Ψ to Ψ if there exists a solution ψ ( t ) with theproperty that lim t →−∞ ψ ( t ) = Ψ and lim t → + ∞ ψ ( t ) = Ψ . More precisely, the connecting orbit leaves the unstable manifold of Ψ andends up in the stable manifold of Ψ . Since the PFC equation is a gradientflow, there cannot exist non-trivial homoclinic connections so there is a natural hierarchy of steady states expressed through heteroclinic connections. Thisconcept is extremely useful to “visualize” the energy landscape.States with Morse index 0 are stable (for fixed parameters) and are thus atthe bottom of the hierarchy. Those states with Morse index 1 have one unstabledirection, so there are two distinct perturbations that lead away from the state.For states with Morse index 2, two unstable directions span infinitely manysuch perturbations, and so on. To detect connections, we propose to initializea PFC flow near an unstable steady state offset by such perturbations. If theflow becomes close enough to another known steady states, we stop and proposea conjectured connection between the two steady states. This procedure oftenallows us to find unknown steady states: when the flow stagnates, the Newtoniteration can be run and often converges in very few steps to a steady state thatcan 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 PFC16cheme detailed in [26].While we cannot for the moment prove such claims because “parameteriz-ing” the infinite dimensional stable manifold of the unstable steady states ishighly non-trivial, we are aware of some preliminary work in this direction [29].That said, computer-assisted proofs of connecting orbits from saddle points toasymptotically stable steady states in parabolic PDEs are starting to appear[30, 31, 32].We first consider the standard parameters ( ¯ ψ, β ) = (0 . , . N x , N y ) = (2 , X to simplify the visualization. Wefind seven steady states: both possible translations of the atoms, stripes anddonuts state, and the trivial constant state. Following the program describedabove, we can construct the “connection diagram” shown in Fig. 3 (a) withthe 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 decaysinto the appropriately shifted hexagonal lattices, but this is a slow process asthe sine modes must grow out of numerical noise. This clearly shows that ourmethod cannot be used to guarantee stability in H because it cannot dependon 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 givenby the coefficients a , N y and a N x ,N y . We place the constant state at the originand plot radial lines along linear combination of the unstable directions. The linelength corresponds to the number of PFC steps needed to approach the targetsteady states. In addition, we can color the points along the line as a function ofenergy to indicate energetic relationships. A variant would be to show the energyas the z -component of a surface; essentially giving an indirect visualization of theenergy landscape through 2D unstable manifolds. In particular, this diagramclarifies the relationships between the steady states. For instance, the stripesstates are formed by adding the a , N y mode to the constant state while thedonuts are combinations of the atoms and stripes states.We now consider the localized patterns regime to illustrate these ideas withstates of high Morse index. We do not attempt to build a higher dimensionalvisualization, but simply attempt to recover the “pathways” between the highlyunstable hexagonal lattice with Morse index 20 towards stable steady states.This is visualized in the connection diagram of Fig. 4 (a) which includes afew states of Table 3. In (b), we plot the energy along the PFC flow startingfrom the index 2 state; this plot can be thought of as one of the rays in adiagram 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) ofatoms. In this process, we could not verify that these intermediates are steadystates since the Newton iteration always converged to the endpoint; we thensuppose 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 theMorse index 1 states in Fig. 3 (b). Indeed, unless “trivial” combinations of the17 a)(b) Figure 3: Connection diagram (a) where arrows represent likely connections; theconstant state is connected to all others. The vertical axis gives the ordering inenergy while the numbers give the Morse index. (b) Energy visualization withrespect to the unstable directions of the constant state. This diagram illustrateshow the unstable directions combine to transform the constant state into otherlower energy states. The unstable directions serve as the main axes and thelines represent different initial perturbations. The length of the lines indicatethe number of PFC steps before the flows becomes close to the connecting steadystates. Colors represent energy (red for high and blue for low energy).18 a) (b)
Figure 4: Connection diagram (a) where arrows represent a few of the con-nections found. The two hexagonal lattice states differ in their amplitude andstability. The vertical axis roughly indicates the energy while the numbers givethe Morse indices. We could not obtain (nor disprove) a connection to the singleatom state, indicated with the question mark. The connection labeled with astar is broken down in the energy plot to the right (b). These states appearto be metastable intermediates where the energy gradient becomes small andthe evolution slows down considerably. The blue curve shows the energy as afunction of time in arbitrary units, highlighting momentaneous “flats” in theevolution. 19nstable direction happen to go to an unstable state, we are unlikely to findsuch connections numerically. Similarly, our attempts to find a perturbationthat connects the starting lattice to the single atom state were unfruitful.
A verified steady state (cid:101) a for some parameter ( ¯ ψ, β ) is usually part of a familyof steady states representing a “phase” of matter. In fact, the candidate mini-mizers defined in appendix A as functions of ( ¯ ψ, β ) approximate such families,or branches in the bifurcation diagram . In this context, we can construct suchbranches by starting at a known steady state, vary ¯ ψ and find the closest steadystate 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-arclengthcontinuation [34] in ¯ ψ . Note that the unstable direction that is to followed isprecisely given by the one corresponding to the “fixed” a , = ¯ ψ condition andthis is one of the reasons that we chose to enforce this directly in the formulationof F . As a possible extension, 2D manifolds can be constructed in 2-parametercontinuation 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 . , . ψ . The domain is kept smallwith ( 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 continuethe branches through the folds at ¯ ψ = ± (cid:112) / β . This branch intersects thecheckers state at β = 15 ¯ ψ and the constant and stripes states at β = 3 ¯ ψ . Theenergy plot (b) clearly shows that the donuts state is the “proper” hexagonallattice 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. 20 a) (b) Figure 5: Continuation (bifurcation) diagram showing the L norm of the phase(a) and the energy offset by E (b) as functions of ¯ ψ . The dots represent thestarting points at ( ¯ ψ, β ) = (0 . , . β . In par-ticular, Fig. 6 shows the atoms/donuts branch and the single atom branch inthe localized patterns regime near ( ¯ ψ, β ) = (0 . , .
6) with ( N x , N y ) = (7 , L norm and (b) shows the energy of the phase fieldas functions of ¯ ψ . The hexagonal lattice traces out its usual figure-eight pat-tern while the single atom (and other localized states in general) traces out acomplicated looping path. Such branches illustrate the “snaking” phenomenonpreviously observed in modified Swift-Hohenberg equations that support suchlocalized patterns, see [27] for example. We observe that the path loops on it-self in one direction as the single atom evolves into a localized pattern with 9, 7then 4 atoms before looping back with a 90 ◦ rotation. In the other direction, thebranch moves towards the transition between the hexagonal and constant stateswhere it again loops back. This computation is difficult because the truncationmust remain large and the pseudo-arclength step size must remain small; if thestep size is larger than 0 . a) (b) Figure 6: Continuation (bifurcation) diagram showing the L norm of the phase(a) and the energy offset by E (b) as functions of ¯ ψ . The inset in (a) shows thenorm of ψ − ¯ ψ to better illustrate the snaking phenomenon. Both the hexagonaland single atom branches appear to loop on themselves. We surveyed the basic properties of the PFC equation as a dynamical systemin the framework of rigorous numerics. Thanks to an application of the radiipolynomial approach, we were able to verify the existence of steady states closeto numerically computed approximations. This provided us with important verified information as to the behavior of the energy landscape, especially interms of energetic relationships between steady states. We were also able toprovide partial stability results with the caveat that they only applied to thecosine Fourier series. The Morse indices given were lower bounds in H - thusthose steady states with Morse index higher than 0 must be unstable in H .Such ideas were applied in various regimes of the PFC equation to verify thatcertain important patterns are steady states (as opposed to metastable interme-diates) including single atoms, other localized patterns and grain boundaries. Inparticular, we showed that two-mode PFC supports a non-zero misorientationgrain boundary steady state that we expect to be stable. We also showed theconstruction of the phase diagram with our fully nonlinear approach.Finally, we used such results to further investigate the energy landscapethrough connections or orbits and through parameter continuation. Connec-tions reveal the energetic and dynamical relationships between steady states,highlighting the behavior of unstable patterns as they reach states with lowerenergy. Continuation is especially useful to understand how the important statesevolve across parameter space, highlighting the surprising behavior of the hexag-onal lattice patterns and the snaking behavior of localized patterns.Our work suggests several interesting directions for future work. On onehand, our connection results could be made rigorous with a technique to proveorbits from unstable to stable manifolds. This is a complicated problem because22he stable manifold is infinite dimensional and special techniques must be ap-plied to properly parameterize its “dominant” submanifold. On the other, ourcontinuation results could also be made rigorous or extended to 2-parametercontinuation to reveal more interesting behavior. Alternatively, parameter con-tinuation could be applied to the domain size, for example to investigate prob-lems in elasticity. References [1] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, “Modeling elas-ticity in crystal growth,”
Physical Review Letters , vol. 88, p. 245701, 2002.[2] J. Swift and P. C. Hohenberg, “Hydrodynamic fluctuations at the convec-tive instability,”
Physical Review A , vol. 15, pp. 319–328, 1977.[3] G. Martine La Boissoni`ere and R. Choksi, “Atom based grain extractionand measurement of geometric properties,”
Modelling and Simulation inMaterials Science and Engineering , vol. 26, no. 3, p. 035001, 2018.[4] G. Martine La Boissoni`ere, R. Choksi, K. Barmak, and S. Esedo¯glu,“Statistics of grain growth: experiment versus the Phase-Field-Crystal andMullins models,”
Materialia , p. 100280, 2019.[5] R. Backofen, K. Barmak, K. E. Elder, and A. Voigt, “Capturing the com-plex physics behind universal grain size distributions in thin metallic films,”
Acta Materialia , vol. 64, pp. 72–77, 2014.[6] H. Emmerich, H. L¨owen, R. Wittkowski, T. Gruhn, G. I. T´oth, G. Tegze,and L. Gr´an´asy, “Phase-field-crystal models for condensed matter dynam-ics on atomic length and diffusive time scales: an overview,”
Advances inPhysics , vol. 61, no. 6, pp. 665–743, 2012.[7] M. Greenwood, N. Ofori-Opoku, J. Rottler, and N. Provatas, “Modelingstructural transformations in binary alloys with phase field crystals,”
Phys-ical Review B , vol. 84, no. 6, p. 064104, 2011.[8] M. Seymour and N. Provatas, “Structural phase field crystal approach formodeling graphene and other two-dimensional structures,”
Physical ReviewB , vol. 93, p. 035447, 2016.[9] P. Hirvonen, M. M. Ervasti, Z. Fan, M. Jalalvand, M. Seymour, S. M. V.Allaei, N. Provatas, A. Harju, K. R. Elder, and T. Ala-Nissila, “Multiscalemodeling of polycrystalline graphene: a comparison of structure and de-fect energies of realistic samples from phase field crystal models,”
PhysicalReview B , vol. 94, no. 3, p. 035414, 2016.2310] D. Shirokoff, R. Choksi, and J.-C. Nave, “Sufficient conditions for globalminimality of metastable states in a class of non-convex functionals: a sim-ple approach via quadratic lower bounds,”
Journal of Nonlinear Science ,vol. 25, no. 3, pp. 539–582, 2015.[11] H. Koch, A. Schenkel, and P. Wittwer, “Computer-assisted proofs in anal-ysis and programming in logic: a case study,”
SIAM Review , vol. 38, no. 4,pp. 565–604, 1996.[12] M. T. Nakao, “Numerical verification methods for solutions of ordinaryand partial differential equations,”
Numerical Functional Analysis and Op-timization , vol. 22, no. 3-4, pp. 321–356, 2001.[13] W. Tucker,
Validated numerics: a short introduction to rigorous computa-tions . Princeton University Press, 2011.[14] J. B. van den Berg and J. P. Lessard, “Rigorous numerics in dynamics,”
Notices of the American Mathematical Society , vol. 62, no. 9, pp. 1057–1061, 2015.[15] J. G´omez-Serrano, “Computer-assisted proofs in PDE: a survey,”
SeMAJournal , pp. 1–26, 2018.[16] K. A. Wu, A. Adland, and A. Karma, “Phase-field-crystal model for FCCordering,”
Physical Review E , vol. 81, no. 6, p. 061601, 2010.[17] S. Day, J.-P. Lessard, and K. Mischaikow, “Validated continuation for equi-libria of PDEs,”
SIAM Journal on Numerical Analysis , vol. 45, no. 4,pp. 1398–1424, 2007.[18] A. Hungria, J. P. Lessard, and J. D. Mireles James, “Rigorous numerics foranalytic solutions of differential equations: the radii polynomial approach,”
Mathematics of Computation , vol. 85, no. 299, pp. 1427–1459, 2016.[19] I. Bal´azs, J. B. van den Berg, J. Courtois, J. Dud´as, J. P. Lessard, A. V¨or¨os-Kiss, J. F. Williams, and X. Y. Yin, “Computer-assisted proofs for radiallysymmetric solutions of PDEs,”
Journal of Computational Dynamics , vol. 5,no. 1-2, pp. 61–80, 2018.[20] J. B. van den Berg, “Introduction to rigorous numerics in dynamics: gen-eral functional analytic setup and an example that forces chaos,”
RigorousNumerics in Dynamics, Proceedings of Symposia in Applied Mathematics ,vol. 74, pp. 1–25, 2017.[21] J. B. van den Berg and J. F. Williams, “Validation of the bifurcation di-agram in the 2D Ohta-Kawasaki problem,”
Nonlinearity , vol. 30, no. 4,p. 1584, 2017.[22] J. B. van den Berg and J. F. Williams, “Rigorously computing symmet-ric stationary states of the Ohta-Kawasaki problem in three dimensions,”
SIAM J. Math. Anal. , vol. 51, no. 1, pp. 131–158, 2017=9.2423] R. E. Moore,
Interval analysis , vol. 4. Prentice-Hall, 1966.[24] S. M. Rump, “INTLAB - interval laboratory,” in
Developments in ReliableComputing , pp. 77–104, Springer, 1999.[25] G. I. Hargreaves, “Interval analysis in MATLAB,”
Numerical Algorithms ,no. 2009.1, 2002.[26] M. Elsey and B. Wirth, “A simple and efficient scheme for phase field crys-tal simulation,”
ESAIM: Mathematical Modelling and Numerical Analysis ,vol. 47, no. 5, pp. 1413–1432, 2013.[27] D. J. B. Lloyd, B. Sandstede, D. Avitabile, and A. R. Champneys, “Lo-calized hexagon patterns of the planar Swift-Hohenberg equation,”
SIAMJournal on Applied Dynamical Systems , vol. 7, no. 3, pp. 1049–1100, 2008.[28] J. B. van den Berg, A. Deschˆenes, J. P. Lessard, and J. D. Mireles James,“Stationary coexistence of hexagons and rolls via rigorous computations,”
SIAM Journal on Applied Dynamical Systems , vol. 14, no. 2, pp. 942–979,2015.[29] J. van den Berg, J. Jaquette, and J. Mireles James, “Validated numericalapproximation of stable manifolds for parabolic partial differential equa-tions.” Preprint, 2020.[30] J. Cyranka and T. Wanner, “Computer-assisted proof of heteroclinic con-nections in the one-dimensional Ohta-Kawasaki Model,”
SIAM J. Appl.Dyn. Syst. , vol. 17, no. 1, pp. 694–731, 2018.[31] C. Reinhardt and J. D. Mireles James, “Fourier-Taylor parameterizationof unstable manifolds for parabolic partial differential equations: formal-ism, implementation and rigorous validation,”
Indag. Math. (N.S.) , vol. 30,no. 1, pp. 39–80, 2019.[32] J. Jaquette, J.-P. Lessard, and A. Takayasu, “Global dynamics in noncon-servative nonlinear Schr¨odinger equations.” Preprint, 2020.[33] J. B. van den Berg, J. P. Lessard, and K. Mischaikow, “Global smooth solu-tion curves using rigorous branch following,”
Mathematics of Computation ,vol. 79, no. 271, pp. 1565–1584, 2010.[34] H. B. Keller,
Lectures on numerical methods in bifurcation problems , vol. 79of
Tata Institute of Fundamental Research Lectures on Mathematics andPhysics . Published for the Tata Institute of Fundamental Research, Bom-bay, 1987. With notes by A. K. Nandakumaran and Mythily Ramaswamy.[35] M. Gameiro, J. P. Lessard, and A. Pugliese, “Computation of smooth man-ifolds via rigorous multi-parameter continuation in infinite dimensions,”
Foundations of Computational Mathematics , vol. 16, no. 2, pp. 531–575,2016. 25 ppendices
A PFC ansatz in 2D
PFC simulations can be classified according to a small number of regimes oransatz that represent the (expected) global minimizer. The choice of such can-didates is motivated by numerical experiments but can be obtained analyticallyfrom techniques such as linear stability analysis. Consider a periodic “singleFourier mode” phase field of the form ψ ( x, y ) = ¯ ψ + A cos( y ) + A cos (cid:32) √ x − y (cid:33) + A cos (cid:32) √ x + 12 y (cid:33) on the rectangular domain [0 , π/ √ × [0 , π ]. Inserting this ansatz into thePFC energy yields an expression E [ A , A , A ] that can be optimized in thethree amplitudes. This procedure yields three main classes of states that arewell-known in the PFC literature. • The constant state A = A = A = 0. • The stripes state A = A = 0 and A = √ (cid:112) β − ψ . • The hexagonal lattice states A = A = A = − ψ ± √ (cid:113) β − ¯ ψ .In addition, we also find a “checkers state” where A and A = A are given asmore complicated expressions. The two hexagonal lattices differ in energy: thepositive 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 isnever optimal while the atoms state is more optimal than the donuts state for¯ ψ > β = 3 ¯ ψ , the constant, stripes and donuts states are all ψ ( x, y ) = ¯ ψ . Similarly, the donuts and atoms states merge at β = 12 / ψ . Wecan also compute when states have the same energy. Such behavior occurs ontransition curves of the form β = α ¯ ψ .We can construct the phase diagram in Fig. 7 by labeling with the expectedglobal minimizer. We show in the main text that this “linear” description ofPFC 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 andany fixed point (cid:101) a of T is a zero of F because A is injective. The derivative of T is bounded and also Fr´echet differentiable since we have for any x ∈ X || DT ( x ) || B ( X ) = || I − ADF ( x ) || B ( X ) ≤ || I || B ( X ) + || ADF ( x ) || B ( X ) < ∞ a) (b) Figure 7: Phase diagram (a) and detail (b) constructed by comparing theoptimal energy of the four ansatz. The curves β = 3 ¯ ψ and β = 12 / ψ correspond respectively to the curves on which the amplitude of the stripesand lattice states become complex. The stripes-atoms (blue-yellow) transitioncurve is β ≈ .
22 ¯ ψ while the constant-atoms (red-yellow) transition curve is β = 37 /
15 ¯ ψ . Donuts are never optimal for ¯ ψ >
0. Checkers are never optimal.as DF is the (bounded linear) Fr´echet derivative of F . Now suppose p ( r ) < r >
0, then the radii polynomial in the main text gives Z ( r ) r + Z + Z < p ( r ) /r − Y /r < Y is positive.Let a ∈ B r (¯ a ), we can use the previous inequality to bound || DT ( a ) || B ( X ) = || I − AA † + AA † − ADF (¯ a ) + ADF (¯ a ) − ADF ( a ) || B ( X ) ≤ || I − AA † || B ( X ) + || A ( A † − DF (¯ a )) || B ( X ) + || A ( DF (¯ a ) − DF ( a )) || B ( X ) ≤ Z + Z + Z ( r ) r . Pairing this with the mean value inequality for T , || T ( a ) − ¯ a || X = || T ( a ) − T (¯ a ) + T (¯ a ) − ¯ a || X ≤ sup z ∈ B r (¯ a ) || DT ( z ) || B ( X ) || a − ¯ a || X + || AF (¯ a ) || X ≤ ( Z + Z + Z ( r ) r ) r + Y = p ( r ) + r < r hence T maps B r (¯ a ) to its interior thanks to the strict inequality. Similarlyfor x, y ∈ B r (¯ a ), || T ( x ) − T ( y ) || X ≤ sup z ∈ B r (¯ a ) || DT ( z ) || B ( X ) || x − y || X ≤ ( Z + Z + Z ( r ) r ) || x − y || X < || x − y || X T : B r (¯ a ) → B r (¯ a ) is a contraction with constant κ = ( Z + Z + Z ( r ) r ) < C Computation of the radii polynomial bounds
In the following calculations, we will use usual results such as || Qb || ,ν ≤ || Q || B ( (cid:96) ν ) || b || ,ν , || Q || B ( (cid:96) ν ) = sup || b || ,ν =1 || Qb || ,ν , and the following proposition to compute the norm of operators on X : Proposition 1.
Let Q be an operator such that Q α,σ = c α δ σ − α δ σ − α when-ever α, σ / ∈ U = { , , ..., M } , then || Q || B ( (cid:96) ν ) ≤ max α ∈ U (cid:40) ν α (cid:88) σ ∈ U | Q α,σ | ν σ (cid:41) + sup α/ ∈ U | c α | . Proof.
Let b ∈ X , then Qb can be decomposed as the action of the first finiteblock onto b σ for σ ∈ U and infinitely many diagonal terms c σ b σ for σ / ∈ U . Thenorm of Q can then be written as the sum of two disjoint positive sums usingthe triangle inequality: || Qb || ,ν = (cid:88) σ ∈ N | ( Qb ) σ | ν σ = (cid:88) σ ∈ U (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) α ∈ U Q α,σ b α (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ν σ + (cid:88) σ / ∈ U | c σ b σ | ν σ ≤ (cid:88) σ ∈ U (cid:88) α ∈ U | Q α,σ || b α | ν σ + (cid:88) σ / ∈ U | c σ || b σ | ν σ ≤ (cid:88) α ∈ U (cid:32) ν α (cid:88) σ ∈ U | Q α,σ | ν σ (cid:33) | b α | ν α + (cid:88) α/ ∈ U | c σ || b α | ν α The second term is bounded by the trivial bound C = sup α/ ∈ U | c α | , extractingthe norm of b over N \ U . Similarly, the first term is bounded bymax α ∈ U (cid:40) ν α (cid:88) σ ∈ U | Q α,σ | ν σ (cid:41) (cid:88) α ∈ U | b α | ν α = K (cid:88) α ∈ U | b α | ν α ≤ K || b || ,ν . The norm of Q in B ( (cid:96) ν ) 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 acton different subspaces of (cid:96) ν . The estimate also allows us to compute the normof finite tensors by letting the c α vanish. Let us now apply these results tocompute the four necessary bounds. 28 .1 Z bound This bound is the easiest since by construction, A and A † are approximateinverses up to numerical inversion errors. We then have Z = || I − AA † || B ( (cid:96) ν ) = || I ( M ) − A ( M ) G || B ( (cid:96) ν ) which can be evaluated using proposition 1. C.2 Y bound We must compute AF (¯ a ) so that most terms will be given by the finite product A ( M ) F ( M ) (¯ a ). There still remain some non-zero convolution coefficients in thefull F (¯ a ): since ¯ a ( M ) has M + 1 coefficients in each dimension, L α (¯ a ∗ ¯ a ∗ ¯ a ) α will have 3 M + 1 non-zero coefficients in each dimension. These are multipliedby the appropriate L − α γ − α , resulting in || AF (¯ a ) || ,ν ≤ || A ( M ) F ( M ) (¯ a ) || ,ν + (cid:88) α ∈{ , ,..., M } \ U (cid:12)(cid:12)(cid:12)(cid:12) (¯ a ∗ ¯ a ∗ ¯ a ) α γ α (cid:12)(cid:12)(cid:12)(cid:12) ν α . C.3 Z bound To compute the Z bound, let b, h ∈ X with || h || ,ν = 1 and consider first theeffect of DF ( b ) on h ,( DF ( b ) h ) α = dds F α ( b + sh ) (cid:12)(cid:12)(cid:12)(cid:12) s =0 = L α ( γ α h α + 3( b ∗ b ∗ h ) α ) . Fix r > b = ¯ a + R where || R || ,ν ≤ r , we then have that(( DF ( b ) − DF (¯ a )) h ) α = 3 L α ( b ∗ b ∗ h − ¯ a ∗ ¯ a ∗ h ) α = 3 L α ((2¯ a ∗ R + R ∗ R ) ∗ h ) α . Note that the initial factor L α can instead be represented by the diagonaloperator defined by Λ α,σ = L α δ α − σ δ α − σ . Then, using the fact that theconvolution on X is a Banach algebra, || A ( DF ( b ) − DF (¯ a )) || B ( (cid:96) ν ) ≤ || A Λ || B ( (cid:96) ν ) || a ∗ R + R ∗ R || ,ν || h || ,ν ≤ || A Λ || B ( (cid:96) ν ) (2 || ¯ a || ,ν + r ) r where the norm of A Λ is computed using proposition 1 with the bound Γ =max( γ − ,M +1 , γ − M +1 , ) on the diagonal terms. This computation works for any r >
0, so we have Z ( r ) = 6 || A Λ || B ( (cid:96) ν ) || ¯ a || ,ν + 3 || A Λ || B ( (cid:96) ν ) r = Z (0)2 + Z (1)2 r . This can be done since γ α is monotone decreasing as long as either α or α can dominatethe − β term. For PFC, γ α = ( L α + 1) − β so a sufficiently large M ensures that M/L x or M/L y dominates β . This condition can be checked numerically. .4 Z bound For the final bound, we now consider the action of A † on the same vector h :( A † h ) α = h α if α = (0 , (cid:80) σ G α,σ h σ if α ∈ U \ { (0 , } L α γ α h α otherwiseLet η 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:(( DF (¯ a ) − A † ) h ) α = α = (0 , L α (¯ a ∗ ¯ a ∗ η ) α if α ∈ U \ { (0 , } L α (¯ a ∗ ¯ a ∗ h ) α otherwiseConsider now the action of A on the difference above. The first block of thedifference will be multiplied by the inverse of G while the tail will be multipliedby the appropriate L − α γ − α , thus || A ( DF (¯ a ) − A † ) h || ,ν ≤ (cid:88) σ ∈ U (cid:88) α ∈ U | A α,σ L α (¯ a ∗ ¯ a ∗ η ) α | ν σ +3 (cid:88) σ ∈ Z \ U (cid:12)(cid:12)(cid:12)(cid:12) (¯ a ∗ ¯ a ∗ h ) σ γ σ (cid:12)(cid:12)(cid:12)(cid:12) ν σ using L , = 0 and the triangle inequality for the first term.Let φ ∈ X be such that | (¯ a ∗ ¯ a ∗ η ) α | ≤ φ α whenever α ∈ U and 0 otherwise.Using the Banach algebra property and the bound Γ ≥ | γ σ | − to bound theinfinite sum, we have || A ( DF (¯ a ) − A † ) || B ( (cid:96) ν ) ≤ (cid:88) σ ∈ U | ( A Λ φ ) σ | ν σ + 3Γ || ¯ a || ,ν = 3 || A Λ φ || ,ν + 3Γ || ¯ a || ,ν = Z which can be computed numerically once the (finitely many) φ α have beenobtained. To compute them, we now shift for a moment to Z and extend allvectors appropriately. Now let q = ¯ a ∗ ¯ a , then | ( q ∗ η ) α | ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) σ ∈ Z q α − σ η σ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:88) σ ∈ V ( α ) ∩ W (cid:18) | q α − σ | ν | σ | (cid:19) | h σ | ν | σ | ≤ (cid:88) σ ∈ V ( α ) ∩ W (cid:32) sup τ ∈ V ( α ) ∩ W | q α − τ | ν | τ | (cid:33) | h σ | ν | σ | ≤ sup σ ∈ V ( α ) ∩ W | q α − σ | ν | σ | || h || ,ν where V ( α ) , W ⊂ Z are the regions over which q α − σ and η are non-zero re-spectively. Since q τ = 0 whenever either | τ | or | τ | is larger than 2 M and we30nly need | α | ∈ U , we obtain the overestimate that V ( α ) ⊂ {− M, ..., M } .Further, η σ must vanish for | σ | ∈ U so φ α = max σ ∈{− M,..., M } \{− M,...,M } | (¯ a ∗ ¯ a ) | α − σ | , | α − σ | | ν | σ | for α ∈ U (and 0 otherwise), which completes the computation of Z . D Energy computation and real space norms
The notion of closeness between (cid:101) a and ¯ a extends to their energies. For simplicity,we only handle the basic PFC energy E [ ψ ] = − (cid:90) Ω
12 ( ∇ ψ + ψ ) + 14 ( ψ − β ) . We will write without relabeling E [ ψ ( a ) ] = E [ a ] when ψ ( a ) is the phase fieldcorresponding to the Fourier coefficients a . Taking the average of a Fourierseries returns its constant mode such that E [ a ] = (cid:18)
12 ( La + a ) ∗ ( La + a ) + 14 ( a ∗ a − βδ α δ α ) ∗ ( a ∗ a − βδ α δ α ) (cid:19) , . In the context of the radii polynomial approach, let t = (cid:101) a − ¯ a such that E [ (cid:101) a ] − E [¯ a ] = E [¯ a + t ] − E [¯ a ]= (cid:18)
12 (2 t + Lt ) ∗ ( Lt ) + (¯ a + L ¯ a ) ∗ ( Lt ) + ((1 − β )¯ a + ¯ a ∗ ¯ a ∗ ¯ a + L ¯ a ) ∗ t + 1 − β t ∗ t + 32 ¯ a ∗ ¯ a ∗ t ∗ t + ¯ a ∗ t ∗ t ∗ t + 14 t ∗ t ∗ t ∗ t (cid:19) , . We can simplify (¯ a + L ¯ a ) ∗ ( Lt ) = (cid:0) L ¯ a + L ¯ a (cid:1) ∗ t by integrating by parts andusing the periodic boundary conditions; for example,( a ∗ ( Lb )) , = − (cid:90) Ω ψ ( a ) ∇ ψ ( b ) = − (cid:90) Ω ∇ ψ ( a ) ψ ( b ) = (( La ) ∗ b ) , . We now use the fact that || t || ,ν < r ∗ . In addition, we can overestimate | a , | ≤ || a || ,ν and use the Banach algebra property to obtain the followingbound: | E [ (cid:101) a ] − E [¯ a ] | ≤ (cid:12)(cid:12)(cid:12) ((2 t + Lt ) ∗ ( Lt )) , (cid:12)(cid:12)(cid:12) + (cid:0) | − β | || ¯ a || ,ν + || ¯ a || ,ν + 2 || L ¯ a || ,ν + || L ¯ a || ,ν (cid:1) r ∗ + 12 (cid:0) | − β | + 3 || ¯ a || ,ν (cid:1) r ∗ + || ¯ a || ,ν r ∗ + 14 r ∗ t has been left as a convolution because it is necessaryto control the growth of L α with the r ∗ bound directly. To do so, we have | S | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) α ∈ Z (2 t α + L α t α )( L − α t − α ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) α ∈ N W α (2 L α + L α ) t α (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) which is another way to obtain the previous integration by parts result. Now, the (cid:96) ν norm of t is bounded by r ∗ ; each member of the sum satisfies the inequality W α | t α | ν | α | < r ∗ . Overestimating W α ≥
1, we can write | S | < (cid:88) α ∈ N W α (2 | L α | + L α ) (cid:18) r ∗ W α ν | α | (cid:19) < r ∗ (cid:88) α ∈ N (2 | L α | + L α ) ρ | α | where ρ = 1 /ν <
1. To evaluate this sum, we can compute the polynomialgeometric series ∞ (cid:88) j =0 ρ j = 11 − ρ , ∞ (cid:88) j =0 j ρ j = ρ + ρ (1 − ρ ) , ∞ (cid:88) j =0 j ρ j = ρ + 11 ρ + 11 ρ + ρ (1 − ρ ) which all converge for ρ <
1. Note that the sums can be evaluated by differen-tiating (cid:80) ρ jx with respect to x = 1. The terms in L α can then be expandedand written in such a fashion and assuming that | S | is finite, the sums can besplit and separated. We have S = (cid:88) α ∈ N | L α | ρ | α | = (cid:88) j,k ∈ N (cid:32)(cid:18) πL x (cid:19) j + (cid:18) πL y (cid:19) k (cid:33) ρ j ρ k = (cid:32)(cid:18) πL x (cid:19) + (cid:18) πL y (cid:19) (cid:33) (cid:88) j,k ∈ N j ρ j ρ k = | L , | (cid:88) j ∈ N j ρ j (cid:32)(cid:88) k ∈ N ρ k (cid:33) = | L , | ρ + ρ (1 − ρ ) · − ρ = | L , | ρ + ρ (1 − ρ ) and similarly, S = (cid:88) α ∈ N L α ρ α = ( L , + L , ) ρ + 11 ρ + 11 ρ + ρ (1 − ρ ) + 2 | L , L , | ρ + 2 ρ + ρ (1 − ρ ) . (2)Putting everything together, we arrive at the bound | E [ (cid:101) a ] − E [¯ a ] | ≤ (cid:0) | − β | || ¯ a || ,ν + || ¯ a || ,ν + 2 || L ¯ a || ,ν + || L ¯ a || ,ν (cid:1) r ∗ + 12 (cid:0) S + S + | − β | + 3 || ¯ a || ,ν (cid:1) r ∗ + || ¯ a || ,ν r ∗ + 14 r ∗ ν because of its influence on r ∗ and the growth of S and S . In principle, onecould 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 [¯ a ] computations are addedto the energy bound, both computed with interval arithmetic, we obtain aninterval that is guaranteed to contain the energy of (cid:101) a itself. In particular, thisallows us to prove which of two steady states is more optimal strictly fromnumerical computations.The previous computations illustrate some techniques that allow us to esti-mate the norm of ψ ( t ) (the phase field corresponding to (cid:101) a − ¯ a ) in terms of r ∗ .For example, || ψ ( t ) || ∞ = sup x ∈ Ω | ψ ( t ) ( x ) | = sup x ∈ Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) α ∈ N W α t α cos (cid:18) πα L x x (cid:19) cos (cid:18) πα L y y (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:88) α ∈ N W α | t α | ≤ (cid:88) α ∈ N W α | t α | ν | α | = || t || ,ν < r ∗ (3)provides a pointwise estimate on the value of the exact steady state. Further, asimple calculation shows that Parseval’s identity holds on X ; i.e. || ψ ( a ) || L (Ω) = | Ω | (cid:80) α ∈ N W α a α . We can then bound the L norm of derivatives, for instance ||∇ ψ ( t ) || L = | Ω | (cid:88) α ∈ N L α t α ≤ S | Ω | r ∗ using Eq. (2). Similar results can be built for the lower norms, thus providingan estimate of the form || ψ ( t ) || H ≤ C (Ω , ν ) r ∗ for some constant C that could be computed if necessary. Note the implicitdependence 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 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 (cid:101) a , we can characterize itsstability in X . This is powerful because even linear results are mostly limited totrivial states but a major limitation is that this does not transfer to H because X is restricted to the cosine series.Suppose we have a steady state (cid:101) a ∈ B r (¯ a ) for a verified radius r , then sta-bility is controlled by the spectrum of DF ( (cid:101) a ). This spectrum is real because The energy intervals must be disjoint for such statements to hold. Otherwise, the proofparameters must be improved to tighten the intervals.
33e are in the context of a gradient flow; this can be seen directly from the defi-nition of A † and DF which are symmetric on interchanging indices. Assumingthere are no zero eigenvalues, the positive and negative ones define the unstable and stable manifolds respectively. A steady state with only strictly negativeeigenvalues is said to be stable.While only the approximation A † is known in practice, it has the same signature as DF ( (cid:101) a ) itself; i.e. they have exactly as many strictly positive orstrictly negative eigenvalues. We compute the spectrum of A † in two parts.The ( M + 1) eigenvalues of the finite block G = DF ( M ) (¯ a ) can be computednumerically and verified using interval arithmetic routines. Most eigenvaluescan be unequivocally assigned a sign, but some may be identically 0 or closerto 0 than the available precision. Stability cannot be ascertained in such cases,but assuming ¯ a was verified using the radii polynomial approach, G must besufficiently 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 α (cid:54) = (0 ,
0) and γ is strictly positiveas long as M is sufficiently large. We then have σ ( A † ) = σ ( G ) ∪ { L α γ α } α ∈ N \ U which can be split into a finite number of positive eigenvalues and infinitelymany negative eigenvalues, assuming there are no small eigenvalues. The equiv-alence with σ ( DF ( (cid:101) a )) follows from a homotopy argument. Let H s = (1 − s ) A † + sDF ( (cid:101) a ) for s ∈ [0 , || I − AH s || B ( (cid:96) ν ) = || I − AA † − sA ( DF (¯ a ) − A † ) + sA ( DF (¯ a ) − DF ( (cid:101) a )) || B ( (cid:96) ν ) ≤ || I − AA † || B ( (cid:96) ν ) + s || A ( DF (¯ a ) − A † ) || B ( (cid:96) ν ) + s || A ( DF (¯ a ) − DF ( (cid:101) a )) || B ( (cid:96) ν ) ≤ Z + sZ + sZ ( r ) r ≤ Z + Z + Z ( r ) r < I − AH s is a boundedoperator 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 ( (cid:101) a ) have the same signature and this is in fact trueover the ball of radius r ∗ around ¯ a . This verification is numerically costly for large truncation order; when
M >
40, we onlyverify the eigenvalues that are larger than some arbitrary lower bound, say −1