A Revisit of The Energy Quadratization Method with A Relaxation Technique

Consider a domain $\Omega$, and denote the thermodynamic variable as $\phi$. We introduce the general formulation for a thermodynamically consistent reversible-irreversible PDE models as

where $\mathcal{B}$ is a trace operator, and $\mathcal{G}:=\mathcal{G}_{a}+\mathcal{G}_{s}$ is the mobility operator, $\mathcal{G}_{s}$ is symmetric and positive semi-definite to ensure thermodynamically consistency, $\mathcal{G}_{a}$ is skew-symmetric, and $\frac{\delta\mathcal{E}}{\delta\phi}$ is the variational derivative of $\mathcal{E}$, known as the chemical potential. Then, the triplet $(\phi,\mathcal{G},\mathcal{E})$ uniquely defines a thermodynamically consistent model. One intrinsic property of (2.1) owing to the thermodynamical consistency is the energy dissipation law

where the inner product is defined by $(f,g)=\int_{\Omega}fgd\mathbf{x}$,$\forall f,g\in L^{2}(\Omega)$, and $\dot{\mathcal{E}}_{surf}$ is due to the boundary contribution, and $g_{b}$ is the boundary integrand. When $\mathcal{G}_{a}=0$, (2.1) is a purely dissipative system; while $\mathcal{G}_{s}=0$, it is a purely dispersive system. $\dot{\mathcal{E}}_{surf}$ vanishes only for suitable boundary conditions, which include periodic and certain physical boundary conditions.

Now, we illustrate the energy quadratization (EQ) idea on solving (2.1). Denote the total energy as $\mathcal{E}(\phi)=\int_{\Omega}ed\mathbf{x}$, with $e$ the energy density function. We denote $L_{0}$ as a linear operator that can be separated from $e$. Introduce the auxiliary variable

where $A_{0}>0$ such that $q$ is a well defined real variable. Then we rewrite the energy as

With the EQ approach above, we transform the free energy density into a quadratic one by introducing an auxiliary variable to ”remove” the quadratic gradient term from the energy density. Assuming $q=q(\phi,\nabla\phi)$ and denoting $g(\phi)=\frac{\partial q}{\partial\phi}$ and $\mathbf{G}(\nabla\phi)=\frac{\partial q}{\partial\nabla\phi}$, we reformulate (2.1) into an equivalent form

with the consistent initial condition $q|_{t=0}=\left.\sqrt{2\Big{(}e-\frac{1}{2}|\mathcal{L}_{0}^{\frac{1}{2}}\phi|^{2}+\frac{A_{0}}{|\Omega|}\Big{)}}\right|_{t=0}$. Now, instead of dealing with (2.1) directly, we develop structure-preserving schemes for (2.5).

The advantage of using model (2.5) over model (2.1) is that the energy density is transformed into a quadratic one in (2.5). Denoting $\Psi=\begin{bmatrix}\phi\\ q\end{bmatrix}$, we rewrite (2.5) into a compact from

where $\mathcal{L}$ is a linear operator, and $\mathcal{N}(\Psi)$ is the mobility operator as

Here $\mathcal{N}_{0}=(1,\,\,\,g(\phi)+\mathbf{G}(\nabla\phi)\colon\nabla)$, $\mathcal{L}=\begin{bmatrix}\mathcal{L}_{0}&\\ &1\end{bmatrix},$ and $\mathcal{N}_{0}^{\ast}$ is the adjoint operator of $\mathcal{N}_{0}$. When $\dot{\mathcal{E}}_{surf}=0$, the energy law is

To demonstrate this methodology’s generality, we introduce a few widely-used PDE models that can be cast as special cases.

Allen Cahn equation. First of all, let us consider the Allen-Cahn equation

It can be formulated in an energy variation form of (2.1) with $\mathcal{G}=1$ and $F=\int_{\Omega}\frac{\varepsilon^{2}}{2}|\nabla\phi|^{2}+\frac{1}{4}(\phi^{2}-1)^{2}d\mathbf{x}$. Next, if we introduce $q:=h(\phi)=\frac{1}{\sqrt{2}}(\phi^{2}-1-\gamma_{0})$, we will have $g(\phi):=\frac{\partial q}{\partial\phi}=\sqrt{2}\phi$. Then the model is reformulated into the quadratization form of (2.5) as

with the consistent initial condition $q|_{t=0}=h(\phi|_{t=0})$. Its matrix formulation of (2.6) can also be easily obtained, by noticing $\mathcal{L}_{0}=-\varepsilon^{2}\Delta+\gamma_{0}$.

Cahn-Hilliard (CH) equation. In the second example, we consider the Cahn-Hilliard equation given as

where $f(\phi)$ is the bulk potential. In this paper, we choose the double-well potential $f(\phi)=\frac{1}{4}(\phi^{2}-1)^{2}$. Other potentials, such as the Flory-Huggins potential can also be used. It can be reformulated in an energy variational form of (2.1) with the notations $\mathcal{G}=-\Delta$, and $F=\int_{\Omega}\Big{[}\frac{\varepsilon^{2}}{2}|\nabla\phi|^{2}+f(\phi)\Big{]}d\mathbf{x}$. Then, we introduce $q:=h(\phi)=\frac{1}{\sqrt{2}}(\phi^{2}-1-\gamma_{0})$ and $g(\phi):=\frac{\partial q}{\partial\phi}=\sqrt{2}\phi$, leading us to the quadratized form of (2.5) as

with the consistent initial condition $q|_{t=0}=h(\phi|_{t=0})$.

Molecular beam epitaxy (MBE) model. Then, we focus on the molecular beam epitaxy (MBE) model with slope section. The model reads

Similarly, this model can be written in an energy variation form of (2.1) as $\mathcal{G}=M$ and $F=\int_{\Omega}\Big{[}\frac{\varepsilon^{2}}{2}(\Delta\phi)^{2}+\frac{1}{4}(|\nabla\phi|^{2}-1)^{2}\Big{]}d\mathbf{x}$. If we introduce the notations $q:=h(\phi)=\frac{\sqrt{2}}{2}(|\nabla\phi|^{2}-1-\gamma_{0})$, and $\mathbf{G}(\nabla\phi):=\frac{\partial q}{\partial\nabla\phi}=\sqrt{2}\nabla\phi$, the MBE model can be written in the quadratized form as

with a consistent initial condition for $q$.

Phase field crystal (PFC) equation. Next, we investigate the phase field crystal (PFC) equation which reads as

with $a_{0}$ and $b_{0}$ are model parameters. The PFC model can also be written in an energy variational form of (2.1) with $\mathcal{G}=\Delta$ , and $F=\int_{\Omega}\Big{\{}\frac{1}{2}\phi\Big{[}-b_{0}+(a_{0}+\Delta)^{2}\Big{]}\phi+\frac{1}{4}\phi^{4}\Big{\}}d\mathbf{x}$. Then, if we introduce the notations $q:=h(\phi)=\sqrt{2}\sqrt{\frac{1}{4}\phi^{4}-\frac{b_{0}}{2}\phi^{2}-\frac{\gamma_{0}}{2}\phi^{2}+C_{0}}$ and $g(\phi):=\frac{\partial q}{\partial\phi}=\frac{\phi^{3}-(b_{0}+\gamma_{0})\phi}{\sqrt{2}\sqrt{\frac{1}{4}\phi^{4}-\frac{b_{0}}{2}\phi^{2}-\frac{\gamma_{0}}{2}\phi^{2}+C_{0}}}$, we can rewrite the PFC model into the quadratizated form of (2.5) as

We recall the classical EQ numerical algorithms to solve the generalized thermodynamically consistent problems. Mainly, we discretize all the linear terms implicitly and all the nonlinear terms explicitly. We will end up with some energy stable and linear numerical algorithms.

It could be easily verified that the EQ schemes are uniquely solvable and unconditionally energy stable [10]. Meanwhile, the known issue is that energy is modified, since $q$ does not necessarily equal to its original definition of (2.3) after discretization, i.e. $q^{n+1}\neq h(\phi^{n+1})$ in general.

To overcome this inconsistency issue of the modified energy and the original energy in the baseline EQ method above, we come up with the novel relaxed EQ schemes as follows, as an analogy to the two baseline EQ schemes.

Similarly, we can propose the relaxed BDF2 scheme as follows.