 # An efficient method for computing stationary states of phase field crystal models

Computing stationary states is an important topic for phase field crystal (PFC) models. Great efforts have been made for energy dissipation of the numerical schemes when using gradient flows. However, it is always time-consuming due to the requirement of small effective time steps. In this paper, we propose an adaptive accelerated proximal gradient method for finding the stationary states of PFC models. The energy dissipation is guaranteed and the convergence property is established for the discretized energy functional. Moreover, the connections between generalized proximal operator with classical (semi-)implicit and explicit schemes for gradient flow are given. Extensive numerical experiments, including two three dimensional periodic crystals in Landau-Brazovskii (LB) model and a two dimensional quasicrystal in Lifshitz-Petrich (LP) model, demonstrate that our approach has adaptive time steps which lead to significant acceleration over semi-implicit methods for computing complex structures. Furthermore, our result reveals a deep physical mechanism of the simple LB model via which the sigma phase is first discovered.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The phase field crystal (PFC) model is an important approach to describe many physical processes and material properties, such as the formation of ordered structures, nucleation process, crystal growth, elastic and plastic deformations of the lattice, dislocations, etc  chen2002phase ; provatas2010phase . More concretely, let the order parameter function be , the PFC model can be expressed by a free energy functional:

 E[ϕ(r);Θ]=G[ϕ(r);Θ]+F[ϕ(r);Θ], (1)

where are the physical parameters, is the interaction energy with polynomial type or log-type formulation and is the bulk energy that contains higher-order linear operators to form ordered structures brazovskii1975phase ; lifshitz1997theoretical ; swift1977hydrodynamic . A typical interaction potential function for a bounded domain is

 G[ϕ]=1|Ω|∫Ω[m∏j=1(Δ+q2j)ϕ]2dr,  m∈N (2)

which can be used to describe the pattern formation of periodic crystals, quasicrystals and multi-polynary crystals.

In order to understand the theory of PFC models as well as predict and guide experiments, it requires to find stationary states and construct phase diagrams of the energy functional (1). Mathematically, denote to be a feasible space, one should solve the minimization problem

 minϕ∈VE[ϕ(r);Θ], (3)

with different physical parameters , which brings the tremendous computational burden. Therefore, within appropriate spatial discretization, the goal of this paper is to develop an efficient and robust numerical method for solving (3) with guaranteed convergence.

Most existing numerical methods for computing the stationary states of PFC model can be classified into two categories. The first class of numerical methods solves the steady nonlinear Euler-Lagrange equations of (

3) through different spatial discretization approaches. The second class of numerical methods has been designed via the formulation of nonlinear gradient flow equations. In these numerical PDE approaches, the time-dependent nonlinear gradient flows are discretized in space via different numerical methods. In these time discretized approaches, great efforts have been made to keep the energy dissipation which is crucial for convergence. Typical energy stable schemes to the gradient flows include convex splitting and stabilized factor methods, and recently developed invariant energy quadrature, and scalar auxiliary variable approaches for a modified energy shen2018scalar . It is noted that the gradient flow approach is able to describe the quasi-equilibrium behavior of PFC systems. Numerically, the gradient flow is discretized in both space and time domain via different discretization techniques and the stationary state is obtained with a proper choice of initial data.

Under an appropriate spatial discretization scheme, the infinite dimensional problem (3) can be formulated as a minimization problem over a finite dimensional space. Thus, there may exist alternative numerical methods that can converge to the steady states quickly by using modern optimization techniques. Similar ideas have been shown success in computing steady states of the Bose-Einstein condensate wu2017regularized and the calculation of density functional theory ulbrich2015proximal ; liu2015analysis . In the PFC models, the discretized energy is nonlinear and non-convex which consists of two parts: bulk energy and interaction energy. Motivated by the semi-implicit scheme and the accelerated proximal gradient (APG) method beck2009fast ; tseng2008accelerated

which has been successfully applied in image processing and machine learning, we propose an efficient numerical method for calculating the steady states of (

3). As the traditional APG method is proposed for convex problem and its oscillation phenomenon slows down the convergence o2015adaptive ; su2014differential , the restart scheme has been used for accelerating the convergence. Moreover, the numerical speed can be further accelerated by using the line search starting with Barzilai-Borwein steps bao2018coherence . The connection of classical explicit/implicit schemes in gradient flows and proximal gradient methods is also built by defining a generalized proximal operator. Extensive numerical experiments have demonstrated that our approach can quickly reach the vicinity of an optimal solution with moderately accuracy, even for very challenge cases. As a byproduct, our numerical result reaveals a deep physical intension of a simple PFC model, the Landau-Brazovskii (LB) model, by obtaining the sigma phase.

The rest of this paper is organized as follows. Different discretizations of the energy functional via the Fourier pseudospectral approach and the projection method are introduced in section 3. In section 4, we present the gradient type method and the adaptive APG method for solving the discretized minimization problem. The connection between our proposed approach and some existing time discretized schemes in these numerical methods for solving gradient-flow equations has been built in section 5. Numerical results are reported in section 6 to illustrate the efficiency and accuracy of our algorithms. Finally, some concluding remarks are given in section 7.

## 2 Physical models

Two classes of PFC models are considered in the paper. The first one is the Landau-Brazovskii (LB) model which describes periodic structures brazovskii1975phase

. The LB model was introduced to investigate the character of phases and phase transition of periodic crystals. It has been discovered in many different scientific fields, e.g., polymeric materials

shi1996theory . In particular, the energy functional of LB model is

 ELB[ϕ(r)]=1|Ω|∫Ω{ξ22[(Δ+1)ϕ]2+τ2!ϕ2−γ3!ϕ3+14!ϕ4}dr, (4)

where is a real-valued function which measures the order of system in terms of order parameter. is the system volume, is the bare correlation length, is the dimensionless reduced temperature, is phenomenological coefficient. Compared with double-well bulk energy, the cubic term in the LB functional helps us study the first-order phase transition.

The second one is the Lifshitz-Petrich (LP) model that can simulate quasiperiodic structures, such as the bi-frequency excited Faraday wave lifshitz1997theoretical , and the explanation of the stability of soft-matter quasicrystals lifshitz2007soft ; jiang2015stability . Before we present the LP model, an introduction of the average spacial integral, so-called almost periodic integral, is necessary. For a space-filling structure, such as the quasicrystal, the average spacial integral can be defined as

 −∫=limR→∞1|BR|∫BR, (5)

where is the ball centred at origin with radii . Using the above notation, the energy functional of LP model is given by

 ELP[ϕ(r)]=−∫{c2[(Δ+q21)(Δ+q22)ϕ]2+ε2ϕ2−κ3ϕ3+14ϕ4}dr, (6)

is the energy penalty, and are phenomenological coefficients.

Formally, the difference between the LB and LP energy functional is the number of length-scale governed by the differential term. The LB model has a one-length-scale which can be used to study the phase behavior of periodic structures brazovskii1975phase ; zhang2008efficient , while the LP model possesses a two-length-scale that can be used to studied the formation and stability of quasicrystals lifshitz1997theoretical ; lifshitz2007soft ; jiang2015stability ; dotera2014mosaic .

## 3 Discretization of the energy functional

In this section, we introduce different discretization schemes of the energy functionals (4) and (6), and reduce them to finite dimensional minimization problems. Two classes of stationary states are considered. The first class of stationary states is periodic in LB model which can be described in a bounded domain. Thus we can truncate the energy functional from the whole space to a bounded domain with periodic boundary condition. Then we employ Fourier pseudospectral method to discretize LB energy functional. The second class of stationary phases can be quasicrystals in LP model. For these structures, the discretization of the energy functional in a bounded domain results in a significant Diophantine approximation error. In this paper, we apply with the projection method jiang2014numerical , a high dimensional interpretation approach, to discretize the LP energy function (6), which can avoid the Diophantine approximation error.

### 3.1 Fourier pseudospectral discretization

Each of the -dimensional periodic system can be described by a Bravis lattice

 R=d∑j=1ℓjaj,   ℓj∈Z, (7)

where the vector

forms the primitive Bravis lattice . The smallest possible periodicity, or named the unit cell, of the system is

 Ω=d∑j=1ζjaj,   ζj∈[0,1). (8)

The associated reciprocal lattice is

 R∗=d∑j=1hjbj,   hj∈Z. (9)

The primitive reciprocal lattice vector satisfies the dual relationship

 aibj=2πδij. (10)

Then the periodic function on the Bravis lattice, i.e., , can be expanded as

 ϕ(r)=∑h∈Zd^ϕ(h)ei(Bh)Tr,   r∈Ω, (11)

where , is invertible. The coefficient, , satisfies

 X:=⎧⎨⎩{^ϕ(h)}h∈Zd:^ϕ(h)∈C, ∑h∈Zd|^ϕ(h)|<∞⎫⎬⎭. (12)

In numerical computations, we need to minimize the LB energy functional (4) in a finite dimensional subspace. More precisely, let , and

 XN:={^ϕ(h)∈X:^ϕ(h)=0, for all |hj|>Nj/2, j=1,2,…,d}. (13)

The number of elements in the set is . The order parameter can be projected into the finite dimensional space , i.e.,

 ϕ(r)≈∑^ϕ(h)∈XN^ϕ(h)ei(Bh)Tr,   r∈Ω. (14)

Due to the orthonormal condition

 1|Ω|∫Ωei(Bh1)Tre−i(Bh2)Trdr=δh1h2, (15)

the LB energy functional can be approximated as

 Eh [^Φ]=ξ22∑h1+h2=0[1−(Bh1)T(Bh2)]2^ϕ(h1)^ϕ(h2)+τ2!∑h1+h2=0^ϕ(h1)^ϕ(h2) −γ3!∑h1+h2+h3=0^ϕ(h1)^ϕ(h2)^ϕ(h3)+14!∑h1+h2+h3+h4=0^ϕ(h1)^ϕ(h2)^ϕ(h3)^ϕ(h4)

where , Fourier coefficient , and

. The convolutions in the above expression can be calculated by Fourier pseudospectral method through the fast Fourier transform (FFT). Therefore, it reduces to a finite dimensional minimization problem:

 min^Φ∈CNEh[^Φ]=Gh[^Φ]+Fh[^Φ] (16)

where and are the discretized interaction and bulk energy, respectively. The gradient of is

 ∇Eh[^Φ]=ξ2Λ^Φ+τ^Φ−γ2F−1N((FN^Φ)2)+16F−1N((FN^Φ)3) (17)

where is a diagonal matrix with entries and is the discretized Fourier transform matrix.

### 3.2 Projection method discretization

For the -dimensional quasicrystals which are the space-filling structures, the spatial integral in the energy functional (1) shall be instead by the almost periodic spacial integral , as defined by Eq. (5). We immediately have the following orthonormal property:

 −∫eik⋅re−ik′⋅rdr=δkk′,   ∀k,k′∈Rd. (18)

For an almost periodic function, the average transformation is

 (19)

and it is well defined katznelson2004anintroduction . In this paper, we carry out the above computation in a higher dimension using the projection method which is based on the fact that a -dimensional quasicrystal can be embedded into an -dimensional periodic structure () hiller1985crystallographic . Using the projection method, the order parameter is

 ϕ(r)=∑h∈Zn^ϕ(h)ei[(P⋅Bh)T⋅r],  r∈Rd, (20)

where is invertible, related to the -dimensional primitive reciprocal lattice and the projection matrix depends on the property of quasicrystals, such as rotational symmetryhiller1985crystallographic . The Fourier coefficient satisfies

 X:={(^ϕ(h))h∈Zn:^ϕ(h)∈C, ∑h∈Zn|^ϕ(h)|<∞}. (21)

Again, in practice, we need to minimize the LP energy functional (4) in a finite dimensional subspace. More precisely, let , and

 XN:={^ϕ(h)∈X:^ϕ(h)=0, for all |hj|>Nj/2, j=1,2,…,n}. (22)

The number of elements in the set is . Together with (18) and (20), the discretized energy function (6) is

 Eh[^Φ] = (23) +ε2∑h1+h2=0^ϕ(h1)^ϕ(h2)−κ3∑h1+h2+h3=0^ϕ(h1)^ϕ(h2)^ϕ(h3) +14∑h1+h2+h3+h4=0^ϕ(h1)^ϕ(h2)^ϕ(h3)^ϕ(h4),

where , , , , . It is clear that the nonlinear (quadratic, cubic and cross) terms in Eq. (23) are -dimensional convolutions in the reciprocal space. A direct evaluation of these convolution terms is extremely expensive. Instead, these terms are simple multiplication in the -dimensional real space. Again, the efficient pseudospectral approach is applied to calculate these convolutions in Eq. (23) through the -dimensional FFT.

Therefore, it leads to the following finite dimensional minimization problem:

 min^Φ∈CNEh[^Φ]=Gh[^Φ]+Fh[^Φ], (24)

where and are the discretized interaction and bulk energies. The gradient of is

 ∇Eh[^Φ]=ξ2Λ^Φ+τ^Φ−γ2F−1N((FN^Φ)2)+16F−1N((FN^Φ)3) (25)

where is a diagonal matrix with entries . The is the discretized Fourier transform matrix and is the corresponding inverse discretized Fourier transform matrix. In the following, we will neglect the superscript of hat for simplicity.

## 4 The proposed numerical approach

In this section, we first review the classical semi-implicit method and accelerated proximal gradient (APG) method and then propose the adaptive APG method with proved convergence. Finally, the connection of generalized proximal operator with gradient flows approaches is present.

### 4.1 Semi-implicit scheme

The semi-implicit scheme is a simple but useful approach for finding the stationary state based on gradient flows. For example, the Allen-Cahn equation of the discretized energy functional is

 Φt=−∇Gh[Φ]−∇Fh[Φ] (26)

with the periodic condition and is the spurious time variable. Given an initial value and the time step , the semi-implicit scheme is

 1α(Φk+1−Φk)=−∇Gh[Φk+1]−∇Fh[Φk], (27)

where is the approximation of the solution at , i.e., . The semi-implicit scheme satisfies the following energy dissipation property.

###### Theorem 4.1.

Let . Assume that there exists a constant such that the bulk energy satisfies , and the time step length , then the solutions of (27) satisfy

 Eh[Φk+1]≤Eh[Φk],    ∀k≥0. (28)
###### Proof.

From (27), it is easy to know

 Φk+1∈argminΦFh[Φk]+⟨∇Fh[Φk],Φ−Φk⟩+12α∥Φ−Φk∥2+Gh[Φ].

It implies

 Fh[Φk]+Gh[Φk] ≥Fh[Φk]+⟨∇Fh[Φk],Φk+1−Φk⟩+12α∥Φk+1−Φk∥2+Gh[Φk+1] ≥Fh[Φk+1]+Gh[Φk+1]+(12α−L2)∥Φk+1−Φk∥2,

where the last inequality is from the Taylor expansion of and the boundedness constraint on . ∎

Therefore, to satisfy the energy dissipation law, the time step length depends on the Lipschitz constant . In a general PFC model, the universal Lipschitz constant may not exist or be very large in bounded domain which leads to a small time step and slows down the convergence speed. Despite its strict requirements on in theory, the semi-implicit scheme works well in practice which inspires us a further exploration of the semi-implicit scheme. In the following context, we will combine modern optimization approaches and the semi-implicit scheme to obtain a more efficient approach.

### 4.2 Accelerated proximal gradient (APG) method

The classical APG method beck2009fast ; tseng2008accelerated is designed for solving the convex composite problem:

 minx∈H H(x)=g(x)+f(x) (29)

where is the finite dimensional Hilbert space equipped with the inner product , and are both continuously convex and has a Lipschitz constant , i.e.

 ∥∇f(x)−∇f(y)∥≤L∥x−y∥, ∀x,y∈H.

Given initializations and , the APG method consists of the following steps:

 tk=(√4(tk−1)2+1+1)/2, (30a) yk=xk+tk−1−1tk(xk−xk−1), (30b) xk+1=Proxαg(yk−α∇f(yk)), (30c)

where and the mapping is defined as

 Proxαg(x)=argminy {g(y)+12α∥y−x∥2}. (31)

It is noted that the proximal map in (31) is well defined as is convex. Moreover, the step size can be set adaptively as long as the following inequality holds:

 H(xk+1)≤Qαk(xk+1,xk)≤H(xk)

where

 Qα(x,y)=f(y)+⟨x−y,∇f(y)⟩+12α∥x−y∥2+g(x).

The APG method has the attractive convergence property as follows.

###### Theorem 4.2 (beck2009fast ).

Let be the sequence generated by the (30a)-(30c) and be the optimal objective value of (29) and be the set of minimizers. For any , we have

 H(xk)−H∗≤2∥x0−x∗∥2α(k+1)2,∀x∗∈X∗. (32)

The discretized energy functional in (16) can be reformulated as form (29) by setting

 f=F[Φ]andg=G[Φ]. (33)

We omit the subscript for simplicity. However, there are two main obstacles when directly applying APG method for solving phase field models as is non-convex and has no universal Lipschitz constant. In this paper, we propose an efficient and convergent numerical algorithm for solving the discretized phase field model (16) by combining APG method with restart techniques.

The restart techniques for the APG method was proposed in o2015adaptive which has shown significant acceleration of the APG method by imposing the decreasing property of the objective value when solving convex problems. Furthermore, another restart strategy called speed restart is developed in su2014differential to ensure the linear convergence of the proposed restart APG method for strongly convex problems. In recent years, the restart techniques have been furtherer applied for solving non-convex composite problems in image processing bao2016image . We introduce the details of the proposed algorithm in the following context.

Let and be the current and previous states respectively and the extrapolation weight . We can obtain a candidate state by

 Ψk+1=ProxαkG(~Φk−αk∇F[~Φk]), (34)

where

 ~Φk=Φk+wk(Φk−Φk−1). (35)

It is noted that the proximal mapping in (34) is well defined as is convex. Different from the APG method, the restart technique is to determine whether we accept the result

as the new estimate

. Inspired by the function value restart condition in o2015adaptive , we choose whenever the following condition holds:

 E[Φk]−E[Ψk+1]≥δ∥Φk−Ψk+1∥2 (36)

for some . If the condition (36) does not hold, we restart the APG by setting . In this case, we have

 Φk+1=ProxαkG(Φk−αk∇F(Φk)). (37)

In fact, the scheme (37) provides an adaptive time step semi-implicit approach when varies. From the continuity of , in (16) and the coercive property of , i.e.

 F(Φ)→+∞,Φ→+∞, (38)

the sub-level set is compact for any . Let be the closed ball that contains . From the smoothness of , is Lipschitz continuous in . Denote to be the Lipschitz constant of in the set , i.e.

 ∥∇F[Φ]−∇F[Ψ]∥≤LM∥Φ−Ψ∥, ∀Φ,Ψ∈M.

Thus, we obtain the the next proposition that shows satisfying the sufficient decrease condition for all .

###### Proposition 4.3.

Given an initial point and the iterates with for . If , then

 E[Φk]−E[Φk+1]≥(1/2αk−LM/2)∥Φk+1−Φk∥2, ∀k=1,2,…. (39)

The proof can be easily obtained from Theorem 4.1. Let is from (37) and , the following sufficient condition

 E[Φk]−E[Φk+1]≥η∥Φk+1−Φk∥2 (40)

holds and thus the (37) is a safe-guard step which ensures energy dissipation. On the other hand, by Proposition 4.3, should be less than which might be very small. Thus, it only allows a small step size which may significantly slow down the convergence, and be always too conservative. By line search technique, we can adaptively estimate the step size which will be introduced as follows.

Estimation of step size . Define: , and . We initialize the search step by the Barzilai-Borwein (BB) method barzilai1988two , i.e.

 β0=⟨sk−1,sk−1⟩⟨sk−1,gk−1⟩orβ0=⟨sk−1,gk−1⟩⟨gk−1,gk−1⟩. (41)

Together with the standard backtracking, we adopt the step size whenever (40) holds. The detailed algorithm of estimation the step size is given in Algorithm 1.

and the proposed adaptive APG algorithm is present in Algorithm 2.

#### 4.3.1 Convergence analysis

In this section, we show that our proposed method converges to a steady state of the original energy function. Firstly, we present a useful lemma for our analysis.

###### Lemma 4.4 (Uniformized Kurdyka-Lojasiewicz property bolte2014proximal .).

Let be a compact set and defined in (16) be bounded below. Assume that is constant on . Then, there exist , , and such that for all and all , one has,

 ψ′(E(u)−E(¯u))∥∇E(u)∥≥1, (42)

where and .

###### Proof.

The proof is based on the facts that and satisfy the so called Kurdyka-Lojasiewicz property on bolte2014proximal . ∎

###### Theorem 4.5.

Let defined in (16) be bounded below and be the sequence generated by Algorithm 2. If and , then has the global convergence property, i.e. there exists a point such that and .

###### Proof.

Define

 Pk+1=ProxαkG(Φk−αk∇F[Φk]) (43)

and two sets and . It is noted that for any , we have . Let , then there exists some such that for all as is increasing and is reset to at most every iteration. We show the following properties of the sequence generated by Algorithm 2.
Sufficient decrease property. If , we have

 E[Φk]−E[Φk+1]≥max(1/2αk−LM/2,η)∥Φk−Φk+1∥2 (44)

from Proposition 4.3 and the line search criterion (40). Together with the condition (36), the following sufficient decrease property holds

 E[Φk]−E[Φk+1]≥ρ1∥Φk−Φk+1∥2, ∀k, (45)

where . Since , there exists such that and as . This implies

 ρ1∞∑k=0∥Φk+1−Φk∥2≤E[Φ0]−E∗<+∞ and limk→+∞∥Φk+1−Φk∥=0. (46)

Bounded the gradient. If , by the optimality condition of (34), we have

 −∇F[~Φk]+1αk(~Φk−Φk+1)=∇G[Φk+1].

Thus, and

 ∥∇E[Φk+1]∥ ≤(LM+1/¯α)∥Φk+1−~Φk∥ (47) ≤(LM+1/¯α)(∥Φk+1−Φk∥+wk∥Φk−Φk−1∥),

as where is a bounded set and is the Lipschitz constant of in . If , by the optimality condition of (43), we have

 −∇F[Φk]+1αk(Φk−Φk+1)=∇G[Φk+1].

Then, , then

 ∥∇E[Φk+1]∥≤(LM+1/¯α)∥Φk+1−Φk∥. (48)

Combining (47) with (48), it follows that

 ∥∇E[Φk+1]∥≤ρ2(∥Φk+1−Φk∥+wk∥Φk−Φk−1∥)≤ρ2(∥Φk+1−Φk∥+¯w∥Φk−Φk−1∥), (49)

where .
Subsequence convergence. Since which is compact, there exists a subsequence and such that

 limj→+∞Φkj=Φ∗,limj→+∞E[Φkj]=E[Φ∗]andlimj→+∞∇E[Φkj]=∇E[Φ∗], (50)

where the last two equalities are from the continuity of . Moreover, (46) implies

 limj→+∞∥Φkj−Φkj−1∥=0andlimj→+∞∥Φkj−1−Φkj−2∥=0. (51)

Then, we know from (49).
Finite length property. Let be the set of limiting points of the sequence starting from . By the boundedness of and the fact , it follows that is a non-empty and compact set. Moreover, from (45), we know is constant on , denoted by . If there exists some such that , then we have for all which is from (45). In the following proof, we assume that for all . Therefore, , there exists some such that for all , we have , i.e.

 Φ∈Γη(Φ∗,ϵ)for all Φ∗∈w(Φ0). (52)

Applying lemma 4.4, for all we have

 ψ′(E[Φk]−E∗)∥∇E[Φk]∥≥1.

Form (49), it implies

 ψ′(E[Φk]−E∗)≥1ρ2(∥Φk−Φk−1∥+wk−1∥Φk−1−Φk−2∥). (53)

By the convexity of , we have

 ψ(E[Φk]−E∗)−ψ(E[Φk