An efficient and energy-stable scheme for the modified phase-field crystal equation with a strong nonlinear vacancy potential

1.   Introduction

Phase-field models have emerged as a powerful computational approach to modeling and predicting mesoscale morphological and microstructure evolution in materials ^[1], such as crystals ^[2,3], data assimilation ^[4], and nucleation ^[5]. As one example, the modified phase-field crystal (MPFC) equation was adopted to describe diffusive dynamics and elastic interactions ^{[6,7,8,9,10,11]} and extended to the MPFC equation with a strong nonlinear vacancy potential (VMPFC) to include vacancies that exist in real materials when the local density is low ^[12]. We consider the following energy functional ^[12]:

where $\Omega$ is a domain in $\mathbb R^d\; (d = 2, 3)$, $\phi$ is the local atomic density field, $F_{dou} (\phi) = \frac{1}{4} \phi^4 - \frac{\epsilon}{2} \phi^2$ is the double-well potential, $\epsilon > 0$ is the undercooling parameter, $F_{vac} (\phi) = \frac{h_{vac}}{3} (|\phi|^3 - \phi^3)$ is the vacancy potential, and $h_{vac} \gg 1$ is the penalization parameter. For simplicity, we take the periodic boundary condition for the system. Then, the VMPFC equation can be regarded as a pseudo-gradient flow of (1.1) in $H^{-1} (\Omega)$:

where $\alpha, \beta > 0$ are constants, $M > 0$ is a mobility, $f_{dou} (\phi) = F_{dou}' (\phi)$, and $f_{vac} (\phi) = F_{vac}' (\phi)$. When $\alpha = h_{vac} = 0$, the VMPFC equation degenerates back to the classical PFC equation. The initial condition of the VMPFC equation is

The VMPFC equation becomes a mass-conservative equation if $\left(\psi^0 ({\bf x}), {\bf 1} \right)_{L^2} = 0$ is used and the solution of the VMPFC equation decays the following energy functional:

where $(\cdot, \cdot)_{L^2}$ and $(\cdot, \cdot)_{H^{-1}}$ denote the inner products in $L^2 (\Omega)$ and $H^{-1} (\Omega)$, respectively. The $H^{-1}$ inner product is defined as follows: for given $f, g \in H_0$ ($H_0$ is a zero average subspace of a Hilbert space), $(f, g)_{H^{-1}} : = (\nabla v_f, \nabla v_g)_{L^2}$, where $v_f, v_g \in H_0$ are the solutions of the periodic boundary value problems $-\Delta v_f = f, \: -\Delta v_g = g$ in $\Omega$, respectively. And $\left\| \cdot \right\|_{H^{-1}} : = \sqrt{(\cdot, \cdot)_{H^{-1}}}$ denotes the $H^{-1}$ norm.

Various numerical schemes have been proposed to solve the VMPFC equation ^{[13,14,15,16,17,18]}. The scalar auxiliary variable, $u(t) = \sqrt{ \int_\Omega \left(F_{dou}(\phi) + F_{vac}(\phi) \right) d{\bf x} + C_u}$, the multiple scalar auxiliary variables, $u_1(t) = \sqrt{ \int_\Omega F_{dou}(\phi) \: d{\bf x} + C_{u_1}}$ and $u_2(t) = \sqrt{ \int_\Omega F_{vac}(\phi) \: d{\bf x} + C_{u_2}}$, the auxiliary variable, $v({\bf x}, t) = \sqrt{ F_{dou}(\phi) + F_{vac}(\phi) + C_v}$, and the exponential scalar auxiliary variable, $w(t) = \exp \left(\frac{\mathcal{F} (\phi)}{C_w} \right)$, were introduced in ^[13], ^[14,15], ^[16], and ^[17], respectively, to redefine the energy functional and reformulate the VMPFC equation, where $C_u$, $C_{u_1}$, $C_{u_2}$, and $C_v$ are constants that make each radicand positive and $C_w$ is a constant that suppresses the exponential growth. And the second-order backward differentiation formula was used to discretize the reformulated system, and the second-order stabilization term, $S(\phi^{n+1} - 2\phi^n + \phi^{n-1})$, was added to improve energy stability. In ^[18], the linear convex splitting was introduced by truncating $\frac{1}{4} \phi^4 + F_{vac}(\phi)$ by the specific function and the second-order Douglas–Dupont regularization, $- S \Delta t \Delta (\phi^{n+1}-\phi^n)$, was added to improve energy stability. For the balance between accuracy and energy stability, $S$ was chosen to be $h_{vac} O (|\min_{{\bf x} \in \Omega} (\phi({\bf x}, t), 0)|)$ in ^[13] and proportional to $h_{vac}$ in ^{[14,15,16,17,18]}. Thus, it was observed that the convergence constant and energy decay trend are affected by $h_{vac}$.

In this paper, we present an efficient and energy-stable scheme based on the Crank–Nicolson (CN) formula for the VMPFC equation. In the scheme, $f_{dou} (\phi)$ and $f_{vac} (\phi)$ are treated explicitly, which makes the scheme efficient, and the energy stability is guaranteed by assuming that $f_{dou}' (\phi)$ and $f_{vac}' (\phi)$ are each bounded and by adding two second-order stabilization terms. In particular, by bounding $f_{dou}' (\phi)$ and $f_{vac}' (\phi)$ respectively, we can choose the stabilization parameters independently of $h_{vac}$. As a result, the convergence constant and energy decay trend are not affected by $h_{vac}$. And the scheme can be easily implemented within a few lines of MATLAB code.

The remainder of this paper is organized as follows. We design the numerical approximation and show its mass conservation and energy stability analytically in Section 2 and numerically in Section 3. Conclusions are given in Section 4. The MATLAB code for the numerical approximation is given in the Appendix.

2.   Efficient and energy-stable scheme

We design the numerical approximation for the VMPFC equation based on the CN formula:

where $\delta_t \phi^{n+1} = \phi^{n+1} - \phi^n$, $\phi^{n+ \frac{1}{2}} = \frac{\phi^{n+1} + \phi^n}{2}$, $\phi^{*, n+ \frac{1}{2}} = \frac{3\phi^n-\phi^{n-1}}{2}$, $\delta_t \psi^{n+1} = \psi^{n+1} - \psi^n$, $\psi^{n+ \frac{1}{2}} = \frac{\psi^{n+1} + \psi^n}{2}$, $A, B > 0$ are stabilization parameters, and $\phi^{-1} \equiv \phi^0$.

Theorem 1. The scheme (2.1)-(2.2) with $\left(\psi^0, {\bf 1} \right)_{L^2} = 0$ is mass conserving.

Proof. Suppose that the scheme (2.1)-(2.2) has a solution. From Eq (2.1), we have

by the periodic boundary condition. This gives the relation

With $\left(\psi^0, {\bf 1} \right)_{L^2} = 0$, the relation ensures that

for all $n \geq 0$. □

Before proving the unique solvability of the scheme (2.1)-(2.2), we simplify the scheme as follows:

where $R = \Delta \left(f_{dou}(\phi^{*, n+ \frac{1}{2}}) + f_{vac}(\phi^{*, n+ \frac{1}{2}}) + \frac{1}{2} (1+\Delta)^2 \phi^n + A \Delta t \Delta \phi^n + B (-2\phi^n + \phi^{n-1}) \right)$.

Theorem 2. The scheme (2.1)-(2.2) with $\left(\psi^0, {\bf 1} \right)_{L^2} = 0$ is uniquely solvable for any $\Delta t > 0$.

Proof. We consider the following functional for $\phi$ defined in the constraint space $\left(\phi, {\bf 1} \right)_{L^2} = \left(\phi^n, {\bf 1} \right)_{L^2}$ given by Theorem 1:

It may be shown that $\phi^{n+1}$ is the unique minimizer of $G(\phi)$ if and only if it solves, for any $\varphi$ with $\left(\varphi, {\bf 1} \right)_{L^2} = 0$,

because $G(\phi)$ is strictly convex by

And, Eq (2.4) is true for any $\varphi$ if and only if the given equation holds:

Hence, minimizing the strictly convex functional $G(\phi)$ is equivalent to solving Eq (2.3). □

To prove the energy stability of the scheme (2.1)-(2.2), we assume that there exist constants $L_1, L_2 > 0$ such that

Theorem 3. Assume that (2.5) is satisfied. The scheme (2.1)-(2.2) with $A \geq \frac{M(L_1+L_2)^2}{16\beta}$ and $B \geq \frac{L_1+L_2}{2}$ fulfills the following energy inequality:

where

Proof. By simple calculations,

and

And, for $i = dou$ and $vac$, to handle $f_i (\phi^{*, n+ \frac{1}{2}})$, we expand $F_i(\phi^{n+1})$ and $F_i(\phi^n)$ at $\phi^{*, n+ \frac{1}{2}}$ as

where $\xi^{n+1}$ is a function that is pointwise bounded between $\phi^{n+1}$ and $\phi^{*, n+ \frac{1}{2}}$ and $\xi^n$ is a function that is pointwise bounded between $\phi^n$ and $\phi^{*, n+ \frac{1}{2}}$. Then, we obtain

Using the identities (2.7)–(2.9),

where we have used the assumption (2.5). And

By the definition of $\tilde {\mathcal F}$,

where we have used $\| \delta_t \phi^{n+1} \|_{L^2}^2 = \| \delta_t \phi^{n+1} \|_{H^{-1}} \| \nabla \delta_t \phi^{n+1} \|_{L^2} \leq \frac{1}{2c \Delta t} \| \delta_t \phi^{n+1} \|_{H^{-1}}^2 + \frac{c \Delta t}{2} \| \nabla \delta_t \phi^{n+1} \|_{L^2}^2$ for any $c > 0$ ^[19]. With $c = \frac{M(L_1+L_2)}{4\beta}$, $A \geq \frac{M(L_1+L_2)^2}{16\beta}$, and $B \geq \frac{L_1+L_2}{2}$, we have $\tilde {\mathcal F}^{n+1} - \tilde {\mathcal F}^n \leq 0$. □

3.   Numerical experiments

3.1.   Numerical implementation

The scheme (2.1)-(2.2) can be expressed as follows:

For the first equation with the periodic boundary condition, we can use the Fourier spectral method ^[20,21,22] for achieving efficient computations.

For the energy stability, we replace $F_{dou}(\phi)$ and $F_{vac}(\phi)$, respectively, by

and

where $p, q > 0$ are constants. It is then obvious that there exist $L_1$ and $L_2$ such that (2.5) is satisfied with $f_i$ replaced by $\tilde f_i = \tilde F_i'$ for $i = dou$ and $vac$. Unless otherwise stated, we set $p = 1$, $r = 5$, $q = \frac{r}{h_{vac}}$, $L_1 = 3p^2-\epsilon$, $L_2 = 4h_{vac} q = 4r$, $A = \frac{M(L_1+L_2)^2}{16\beta}$, and $B = \frac{L_1+L_2}{2}$. We note that the stabilization parameters $A$ and $B$ are independent of $h_{vac}$.

3.2.   Numerical verification of the convergence of the scheme and Theorems 1 and 3

First, we verify the convergence of the scheme using

and take $\Omega = [0, 32]^2$, $\Delta x = \Delta y = \frac{1}{3}$, $\alpha = 1$, $\beta = 1$, $M = 1$, and $\epsilon = 0.025$. Figure 1 shows the relative $l_2$-errors of $\phi(x, y, 40)$ with $h_{vac} = 10$, 100, and 1000 for $\Delta t = 2^{-9}, 2^{-8}, \ldots, 1$, where the error is calculated compared to the reference solution using $\Delta t = 2^{-11}$. It is observed that the convergence order is 2 regardless of $h_{vac}$ and the convergence constant is not affected by $h_{vac}$.

Next, we verify Theorems 1 and 3 using

and set $\Omega = [0,128]^2$, $\Delta x = \Delta y = 1$, $\alpha = 0.01$, $\beta = 1$, $M = 1$, and $\epsilon = 0.9$, where $rand(x, y)$ is a random number between $-0.01$ and $0.01$ at the grid points. Figure 2(a, b) show the evolution of $\int_\Omega \phi(x, y, t) \: dxdy$ and $\tilde {\mathcal F}(t)$, respectively, with $h_{vac} = 5000$ and 500000 for $\Delta t = 2^{-2}, 2^{-1}, \ldots, 2^6$. As proved by Theorems 1 and 3, the masses are conserved, and the energies decrease monotonically regardless of $h_{vac}$. Moreover, the energy decay trend is not affected by $h_{vac}$. Figure 3(a, b) show the evolution of $\phi(x, y, t)$ with $\Delta t = 2^{-2}$ for $h_{vac} = 5000$ and 500000, respectively.

And, to investigate the effect of $A$ and $B$ on the energy stability, we perform simulations with small $A$ and $B$. Figure 4 shows the evolution of $\tilde {\mathcal F}(t)$ with small $A$ and $B$ for $h_{vac} = 5000$. In Figure 4(a, b), $A = \frac{1}{100} \frac{M(L_1+L_2)^2}{16\beta}$ and $B = \frac{1}{100} \frac{L_1+L_2}{2}$ lead to the increase of the energy because they do not satisfy the stability condition ($A \geq \frac{M(L_1+L_2)^2}{16\beta}$ and $B \geq \frac{L_1+L_2}{2}$, Theorem 3).

3.3.   Influence of the vacancy potential on the evolution

To explore the influence of the vacancy potential on the evolution, we employ

and take $\Omega = [0, 32]^3$, $\Delta x = \Delta y = \Delta z = \frac{1}{2}$, $\Delta t = \frac{1}{4}$, $\alpha = 1$, $\beta = 0.8$, $M = 1$, and $\epsilon = 0.9$, where $rand(x, y, z)$ is a random number between $-0.01$ and $0.01$ at the grid points. Figures 5 and 6 show $\phi(x, y, z, 320)$ without the vacancy potential ($h_{vac} = 0$) and with the vacancy potential ($h_{vac} = 3000$), respectively, for $\bar \phi = 0.06$, 0.11, and 0.16. The evolution of $\tilde {\mathcal F}(t)$ is shown in Figure 7. For $h_{vac} = 0$, no vacancies appear even when $\bar \phi = 0.06$. In contrast, for $h_{vac} = 3000$, there are many vacancies. At the same time, the number of atoms increases with $\bar \phi$. The simulation results are in good agreement with those of ^[12].

4.   Conclusions

We presented the efficient and energy-stable scheme based on the CN formula for the VMPFC equation. In the scheme, $f_{dou} (\phi)$ and $f_{vac} (\phi)$ were treated explicitly, which made the scheme efficient, and the energy stability was guaranteed by assuming that $\max_{\phi \in \mathbb R} |f_{dou}'(\phi)| \leq L_1 = 3p^3-\epsilon$ and $\max_{\phi \in \mathbb R} |f_{vac}'(\phi)| \leq L_2 = 4r$ and by adding $-A \Delta t \Delta \delta_t \phi^{n+1}$ and $B (\delta_t \phi^{n+1} - \delta_t \phi^n)$. We proved that the scheme is energy-stable under $A \geq \frac{M(L_1+L_2)^2}{16\beta}$ and $B \geq \frac{L_1+L_2}{2}$ which are independent of $h_{vac}$. Through numerical experiments, it was observed that (ⅰ) our scheme is second-order accurate in time, mass-conservative, energy-stable, and suitable for simulating vacancies; and (ⅱ) the convergence constant and energy decay trend are not affected by $h_{vac}$.

In future work, we have a plan to carry out a rigorous error analysis for the scheme and to analyze the stability of the scheme.

Use of Generative-AI tools declaration

The author declares they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The corresponding author (H.G. Lee) thanks the reviewers for the constructive and helpful comments on the revision of this article and was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. RS-2022-NR069708).

Conflict of interest

The author declares no conflict of interest in this paper.

Appendix

The MATLAB code for the numerical scheme in 2D (Figure 3(a)) is given as follows.

The MATLAB code for the numerical scheme in 3D (Figure 6) is given as follows.