# A Thermodynamically-Consistent Phase Field Crystal Model of Solidification
with Heat Flux

Cheng Wang Department of Mathematics, The University of Massachusetts, North
Dartmouth, MA 02747<EMAIL_ADDRESS>Steven M. Wise Corresponding author:
Department of Mathematics, The University of Tennessee, Knoxville, TN 37996
<EMAIL_ADDRESS>
###### Abstract

In this paper we describe a new model for solidification with heat flux using
the phase field crystal (PFC) framework. The equations are thermodynamically
consistent, in the sense that the time rate of change of the entropy density
is positive in the bulk and at the boundaries of the domain of interest. The
resulting model consists of two equations, a heat-like equation and a mass-
conservation equation that describes how the atom density changes in time and
space. The model is simple, yet it can properly capture the variation in the
free energy landscape as the temperature is changed. We describe to construct
a temperature-atom-density phase diagram using this energy landscape, and we
give a simple demonstration of solidification using the model.

Keywords: Phase field crystal, classical density functional theory, entropy
production, heat transfer, solidification, melting.

AMS subject classification: 80A22 35K35, 35K55, 49J40

##  1 Introduction

The phase field crystal (PFC) model was introduced in [13, 14] as continuum
description of solidification in a unary material. It was formulated as a mass
conservative version of the classical Swift-Hohenberg equation, but, later,
the model was re-derived, via certain reasonable simplifications, from the
dynamical density functional theory (DDFT) [15]. In particular, assuming a
constant, uniform temperature field $T$, one expresses the Helmholtz free
energy density via

$F=\int_{\Omega}\left\\{f(T,\rho)+\frac{T\kappa_{f,\rho,o}}{2}(\rho-\rho_{o})\mathcal{C}(\rho-\rho_{o})\right\\}d\boldsymbol{x},$

where $\Omega$ is some spatial domain of interest; $\rho:\Omega\to[0,\infty)$
is the number density field of the unary material in $\Omega$; the constant
$\rho_{o}>0$ is a reference density; $\kappa_{f,\rho,o}>0$ is a positive
constant; $f$ is the homogeneous Helmholtz free energy density; and
$\mathcal{C}$ is a symmetric, potentially nonlocal, two-point correlation
operator. The free energy density, $f$, is often taken to satisfy an ideal
“gas” model:

$f(T,\rho)=\rho
k_{B}T\ln\left(\frac{\rho}{\rho_{o}}\right)-k_{B}T\left(\rho-\rho_{o}\right),$

where $k_{B}$ is the Bolzmann constant. Often, one makes a (Taylor) polynomial
approximation of the logarithmic term about the reference density to make the
model more tractable. However, it is the singular nature of the logarithmic
term that guarantees the positivity of the solutions, and this is an important
feature in the numerical and PDE analyses. At constant temperature, one can
argue that the dynamics of the model should satisfy a diffusion-dominated mass
conservation equation of the form

$\dot{\rho}=-\nabla\cdot\boldsymbol{J},\quad\boldsymbol{J}=-M\rho\nabla\mu,$

where $\boldsymbol{J}$ is the diffusion flux; $M>0$ is a mobility; and $\mu$
is the chemical potential:

$\mu:=\delta_{\rho}F=k_{B}T\log(\rho)+T\kappa_{f,\rho,o}\mathcal{C}(\rho-\rho_{o}),$

where we have assumed, for simplicity, that the boundary conditions are
periodic. As a consequence of these assumptions, the total free energy is
dissipated as the system evolves toward equilibrium, and the dissipation rate
is

$\dot{F}=-M\int_{\Omega}\rho|\nabla\mu|^{2}\,d\boldsymbol{x}\leq 0.$

Of course, it would be necessary to justify the property that $\rho>0$ (or at
least $\rho\geq 0$) point-wise for the model to make sense. Numerical analyses
of similar gradient flow models, that is, models that have logarithmic energy
potentials, have been performed in [5, 10, 11, 12, 18, 19, 22, 26].

The PFC modeling framework has a couple of basic, distinctive features. First,
the solutions to the PFC-type models exhibit (at least) two distinct phases.
One is a spatially oscillatory phase, which is identified with the solid
phase, and the other is a spatially uniform phase, which is usually identified
as the liquid (or gas) phase. The peaks of the solutions in the oscillatory
phase are interpreted as the “locations” of the atoms, and typically, one can
choose $\mathcal{C}$ so that the peaks are arranged in a desired crystal
structure [20]. Second, PFC models operate at atomic length scales but
diffusive time scales. Thus, the framework can capture long-time phenomena.

In this paper, we will devise a new model for solidification and melting using
the phase field crystal framework. In particular, we will not assume that the
temperature is uniform in space and time. For an adiabatically isolated
system, this requires that the global entropy is increasing as the system
moves towards equilibrium, while, at the same time, the total energy and the
number density of particles should be conserved. From a modeling perspective,
one has some freedom in choosing the equations so that these properties hold.
In contrast with the case of uniform temperature, the Helmholtz free energy
can increase or decrease in time. While important to the model formulation, it
alone does not dictate the drive of the system toward equilibrium.

While our model is, to our knowledge, new and distinct, an earlier effort has
made some important contributions towards adding heat flux to the PFC
framework for melting/solidification. Specifically, in the paper by Kocher and
Provatas [17], the authors add temperature variation in the study of rapid
solidification. In fact, the framework that they derive is quite general, and,
as they claim, seems to apply to several physical settings. However, the
equations that they derive are ultimately different from those found herein.
Our model is simpler in the description of the internal energy and the latent
heat, but it still captures the most important effects of temperature
variation. In particular, the heat equation that we derive is essentially
linear. Ultimately, Kocher and Provatas introduced a couple of physically
reasonable approximations in the derivation of their working model. In
particular, they have introduced a smoothing operation to the internal energy
density. The result of these approximations is that their working model is no
longer provably entropy non-increasing. Our model, on the other hand, retains
this important property.

In this brief paper, we will construct and demonstrate the key features of
this thermal-PFC model. We plan more thorough numerical, mathematical, and
physical investigations in future works. The paper is arranged as follows. In
Section 2, we derive the model using the framework outlined in classical phase
field community in the 1990s. In Section 3, we show how to construct phase
diagrams from the model. Finally, we use the free energy landscape and the
phase diagram information to do some very simple computations in Section 4.

##  2 Derivation

The derivation here follows the ideas in the now classic papers by Charach and
Fife [4], Wang  _et al._ [23], and Wheeler, McFadden, and Boettinger [24].
These papers described thermodynamically consistent phase field models of
solidification with heat flux, that is, with a non-uniform and variable
temperature field. The main difference between these classical solidification
models and and PFC modeling framework is that, in the latter, it is a non-
trivial task to identify the equilibrium phases. Indeed, in the PFC framework,
great care must be taken in identifying the equilibrium liquid and solid
states via the free energy landscape.

###  2.1 Basic Assumptions

Set $\Omega\subset\mathbb{R}^{d}$ and let $e,s,f:\Omega\to\mathbb{R}$ denote
the internal energy, entropy, and Helmholtz free energy densities (per unit
volume). The functions $T,\rho:\Omega\to\mathbb{R}$ are the temperature and
the number density (of particles) fields of a unary material occupying the
volume $\Omega$. The densities $e$, $s$, and $f$ are functions only of the
local values of the thermodynamic variables and do not depend upon gradients,
and, for this reason, they are usually called homogenous energy and entropy
densities.

We assume that the total free energy, entropy, and internal energies have the
following forms, respectively:

$F=\int_{\Omega}\left\\{f(T,\rho)+\frac{\kappa_{f,\rho}}{2}(\rho-\rho_{o})\mathcal{C}(\rho-\rho_{o})\right\\}d\boldsymbol{x},$
$S=\int_{\Omega}\left\\{s(e,\rho)-\frac{\kappa_{s,\rho}}{2}(\rho-\rho_{o})\mathcal{C}(\rho-\rho_{o})\right\\}d\boldsymbol{x},$
$E=\int_{\Omega}\left\\{e(s,\rho)+\frac{\kappa_{e,\rho}}{2}(\rho-\rho_{o})\mathcal{C}(\rho-\rho_{o})\right\\}d\boldsymbol{x},$

where $\mathcal{C}$ is a long-range interaction, or correlation operator, and
$\rho_{o}>0$ is a reference density. In this paper, we will assume that
$\mathcal{C}$ is a differential operator of the form

$\mathcal{C}\varrho=L^{-2}\gamma\varrho+2\Delta\varrho+L^{2}\Delta^{2}\varrho,$

where $\gamma$ is a dimensionless parameter, and $L$ is a characteristic
length. In fact, we could instead assume that $\mathcal{C}$ is a nonlocal
operator, as has been done in the classical density functional theory [15,
17]. We require that $F=E-TS$ (globally) and $f=e-Ts$ (locally) , which
implies that

$\kappa_{f,\rho}=\kappa_{e,\rho}+T\kappa_{s,\rho}.$

Our model could be greatly simplified, by choosing $\kappa_{e,\rho}=0$, as was
done in [23, 24], which yields

$\kappa_{f,\rho}=T\kappa_{s,\rho}.$

It is assumed that $\kappa_{f,\rho}$ is linear in temperature:
$\kappa_{f,\rho}=\kappa_{f,\rho,o}T$, where $\kappa_{f,\rho,o}>0$ is a
constant. Thus

$\displaystyle F$
$\displaystyle=\int_{\Omega}\left\\{f(T,\rho)+\frac{\kappa_{f,\rho,o}T}{2}(\rho-\rho_{o})\mathcal{C}(\rho-\rho_{o})\right\\}d\boldsymbol{x},$
$\displaystyle S$
$\displaystyle=\int_{\Omega}\left\\{s(e,\rho)-\frac{\kappa_{f,\rho,o}}{2}(\rho-\rho_{o})\mathcal{C}(\rho-\rho_{o})\right\\}d\boldsymbol{x},$
$\displaystyle E$
$\displaystyle=\int_{\Omega}\left\\{e(s,\rho)\right\\}d\boldsymbol{x}.$

###  2.2 Entropy Production

For the evolution equations, we appeal to conservation laws and entropy
production requirements, assuming those processes are diffusion dominated.
Because energy is conserved locally and, typically, globally, we have the
equation

$\dot{e}=-\nabla\cdot\boldsymbol{J}_{e}.$

It is expected that this will ultimately provide an equation for the
temperature, $T$. Since the number of particles should be conserved locally
and, typically, globally, we employ for $\rho$ a mass conservation equation of
the form

$\dot{\rho}=-\nabla\cdot\boldsymbol{J}_{\rho}.$

Now, we want the entropy to increase locally and globally [7]. To this end, we
calculate the time derivative of the total entropy

$\dot{S}=\int_{\Omega}\left\\{\left(\frac{\partial s}{\partial
e}\right)_{\rho}\dot{e}+\left(\frac{\partial
s}{\partial\rho}\right)_{e}\dot{\rho}\right\\}d\boldsymbol{x}+\int_{\Omega}\kappa_{f,\rho,o}\mathcal{C}(\rho-\rho_{o})\dot{\rho}\,d\boldsymbol{x},$

upon assuming local thermodynamic equilibrium (LTE) boundary conditions for
$\rho$, so that integration-by-parts could be carried out without the
introduction of any boundary integrals. In essence, LTE boundary conditions
ensure that $\mathcal{C}$ is a symmetric operator:
$\int_{\Omega}\varphi\mathcal{C}\varrho\,d\boldsymbol{x}=\int_{\Omega}\varrho\mathcal{C}\varphi\,d\boldsymbol{x}$.

We assume that the internal energy density is, in its most natural form, a
twice continuously differentiable function of $s$ and $\rho$. We further
assume that the system evolves in such a way that it never deviates greatly
from equilibrium, and, consequently, we can use the equilibrium thermodynamic
theory to develop our equations. The first and second laws are encoded in the
Gibbs relation

$de=\left(\frac{\partial e}{\partial s}\right)_{\rho}ds+\left(\frac{\partial
e}{\partial\rho}\right)_{s}d\rho.$

By definition, the temperature satisfies $\left(\frac{\partial e}{\partial
s}\right)_{\rho}=:T>0$. Thus,

$de=Tds+\left(\frac{\partial e}{\partial\rho}\right)_{s}d\rho.$

As we have seen, the Helmholtz free energy density results from a Legendre
transformation of the other densities, namely, $f=e-Ts$. Thus, it follows that

$df=de-Tds-sdT=-sdT+\left(\frac{\partial e}{\partial\rho}\right)_{s}d\rho.$

The natural variables of the free energy density, $f$, are $T$ and $\rho$, as
we have indicated above. This implies the Maxwell relation

$\left(\frac{\partial f}{\partial\rho}\right)_{T}=\left(\frac{\partial
e}{\partial\rho}\right)_{s}.$ (2.1)

We want to choose our diffusion fluxes so that $0\leq\dot{S}$. Observing that

$ds=\frac{1}{T}de-\frac{1}{T}\left(\frac{\partial
e}{\partial\rho}\right)_{s}d\rho,$

it follows that

$\left(\frac{\partial s}{\partial
e}\right)_{\rho}=\frac{1}{T},\quad\left(\frac{\partial
s}{\partial\rho}\right)_{e}=-\frac{1}{T}\left(\frac{\partial
e}{\partial\rho}\right)_{s}.$

Thus, we can write

$\dot{S}=\int_{\Omega}\left\\{-\frac{1}{T}\nabla\cdot\boldsymbol{J}_{e}+\frac{1}{T}\left(\frac{\partial
e}{\partial\rho}\right)_{s}\nabla\cdot\boldsymbol{J}_{\rho}\right\\}d\boldsymbol{x}+\int_{\Omega}\kappa_{f,\rho,o}\mathcal{C}(\rho-\rho_{o})\nabla\cdot\boldsymbol{J}_{\rho}d\boldsymbol{x}.$

Using the Maxwell relation (2.1), we may write an even more convenient form
for the change in entropy:

$\dot{S}=\int_{\Omega}\left\\{-\frac{1}{T}\nabla\cdot\boldsymbol{J}_{e}+\frac{1}{T}\left(\frac{\partial
f}{\partial\rho}\right)_{T}\nabla\cdot\boldsymbol{J}_{\rho}\right\\}d\boldsymbol{x}+\int_{\Omega}\kappa_{f,\rho,o}\mathcal{C}(\rho-\rho_{o})\nabla\cdot\boldsymbol{J}_{\rho}d\boldsymbol{x}.$

Next, we define the generalized chemical potential

$\mu:=\frac{1}{T}\delta_{\rho}F=\frac{1}{T}\left(\frac{\partial
f}{\partial\rho}\right)_{T}+\kappa_{f,\rho,o}\mathcal{C}(\rho-\rho_{o}).$

Consequently,

$\dot{S}=\int_{\Omega}\left\\{\nabla\left(\frac{1}{T}\right)\cdot\boldsymbol{J}_{e}-\nabla\mu\cdot\boldsymbol{J}_{\rho}\right\\}d\boldsymbol{x}+\int_{\partial\Omega}\left\\{-\frac{1}{T}\boldsymbol{J}_{e}\cdot\boldsymbol{n}+\mu\boldsymbol{J}_{\rho}\cdot\boldsymbol{n}\right\\}d\boldsymbol{a}.$

To get local and global entropy production, we make the following constitutive
choices: for the fluxes,

$\displaystyle\boldsymbol{J}_{e}=M_{T}\nabla\left(\frac{1}{T}\right),\quad\boldsymbol{J}_{\rho}=-M_{\rho}\nabla\mu.$

For the boundary conditions, we take non-negative entropy production
conditions:

$\boldsymbol{J}_{e}\cdot\boldsymbol{n}=-\frac{\beta_{T,\partial\Omega}}{T},\quad\boldsymbol{J}_{\rho}\cdot\boldsymbol{n}=\beta_{\rho,\partial\Omega}\mu\quad\mbox{on
$\partial\Omega$},$

where $\beta_{T,\partial\Omega},\beta_{\rho,\partial\Omega}\geq 0$. These
boundary conditions allow for heat and mass to flow through the outer
boundary, but only in an entropy non-decreasing manner. For adiabatically
insulated materials (meaning no mass or heat is exchanged between $\Omega$ and
the outside world), we assume that
$\beta_{T,\partial\Omega}=\beta_{\rho,\partial\Omega}=0$. Another common,
outer boundary condition results from a constraint on the temperature:

$T=T_{\star}<T_{M}\quad\mbox{on $\partial\Omega$},$

known as the undercooling condition, where $T_{M}$ is the melting temperature.
For this condition, we do not have control of the global entropy production.

For the (general) entropy non-decreasing boundary conditions, the entropy
production rate becomes

$\dot{S}=\int_{\Omega}\left\\{M_{T}\left|\nabla\left(\frac{1}{T}\right)\right|^{2}+M_{\rho}\left|\nabla\mu\right|^{2}\right\\}d\boldsymbol{x}+\int_{\partial\Omega}\left\\{\frac{\beta_{T,\partial\Omega}}{T^{2}}+\beta_{\rho,\partial\Omega}\mu^{2}\right\\}d\boldsymbol{a}\geq
0,$ (2.2)

and the evolution equations are

$\displaystyle\dot{e}$
$\displaystyle=-\nabla\cdot\left(M_{T}\nabla\left(\frac{1}{T}\right)\right),$
$\displaystyle\dot{\rho}$
$\displaystyle=\nabla\cdot\left(M_{\rho}\nabla\mu\right),$ $\displaystyle\mu$
$\displaystyle=\frac{1}{T}\left(\frac{\partial
f}{\partial\rho}\right)_{T}+\kappa_{f,\rho,o}\mathcal{C}(\rho-\rho_{o}).$

###  2.3 The Internal and Free Energy Densities

The only remaining issue is to specify the internal energy density $e$ so that
we can compute its time derivative. To do so, it helps to express the internal
energy density, $e$, as a function of $T$ and $\rho$. These are more natural
variables for the problem. To find $e=\hat{e}(T,\rho)$, first note that

$\frac{f}{T}=\frac{e}{T}-s,$

and, hence, it follows that

$d\left(\frac{f}{T}\right)=-\frac{e}{T^{2}}dT+\frac{1}{T}de-
ds=-\frac{e}{T^{2}}dT+\frac{1}{T}\left(\frac{\partial
e}{\partial\rho}\right)_{s}d\rho.$

Consequently,

$\left(\frac{\partial\left(\frac{f}{T}\right)}{\partial
T}\right)_{\rho}=-\frac{e}{T^{2}}\implies
e=-T^{2}\left(\frac{\partial\left(\frac{f}{T}\right)}{\partial
T}\right)_{\rho}.$ (2.3)

We will begin with an expression of the form $f=f(T,\rho)$. In fact, for the
PFC model, $f$ is usually modeled by an ideal gas law [15]. To gain some
modeling flexibility, we will use a more general non-ideal gas law:

$f(T,\rho)=\rho
k_{B}T\ln\left(\frac{\rho}{\rho_{o}}\right)-k_{B}T\left(\rho-\rho_{o}\right)+\rho_{o}k_{B}Tg\left(\frac{\rho-\rho_{o}}{\rho_{o}}\right)-\alpha\rho_{o}k_{B}T\ln\left(\frac{T}{T_{o}}\right)-\beta\rho
k_{B}T_{o},$

where $k_{B}$ is Boltzmann’s constant; $\rho_{o}$ is a reference density;
$T_{o}>0$ is the reference temperature; $g$ is a polynomial that measures the
deviation of the free energy density from the ideal gas model, and is
typically non-negative; and, finally, $\alpha$, and $\beta$ are positive,
dimensionless constants. Clearly,

$\left(\frac{\partial\left(\frac{f}{T}\right)}{\partial
T}\right)_{\rho}=-\alpha\rho_{o}k_{B}\frac{1}{T}+\beta\rho
k_{B}\frac{T_{o}}{T^{2}}.$

It follows from (2.3), therefore, that the internal energy may be expressed as

$e=\hat{e}(T,\rho)=\alpha\rho_{o}k_{B}T-\beta\rho k_{B}T_{o},$

which is linear in $T$ and $\rho$. Note that we use a hat over $e$, since this
functional form is generally different from the functional form of $e$ in its
natural thermodynamic coordinates $s$ and $\rho$.

For the PFC models, we define the latent heat as

$\mathcal{L}(T):=e(T,\rho_{L})-e(T,\rho_{S})=\beta
k_{B}T_{o}\left(\rho_{S}-\rho_{L}\right),$

where $\rho_{L}$ is the (expected) spatially-uniform equilibrium density in
the liquid phase at temperature $T$, and $\rho_{S}$ is the (expected)
spatially-oscillatory field that characterizes the equilibrium solid phase at
temperature $T$. These two density fields generally depend implicitly on the
temperature. Otherwise, there is no direct dependence on the temperature for
the latent heat. We will derive these fields shortly. In the applications that
we will examine, the following inequality is valid:

$\rho_{L}=\overline{\rho_{L}}<\overline{\rho_{S}},$

where the overline represents the spatial average. In this case, the spatial
average of the latent heat is a positive constant, since $\beta$ is assumed
positive:

$\overline{\mathcal{L}(T)}=\beta
k_{B}T_{o}\left(\overline{\rho_{S}}-\overline{\rho_{L}}\right)>0.$

As we see below, the average densities play an important role in the phase
diagram and, therefore, the dynamics of the model.

###  2.4 The Full Model

Finally, we can find the time derivative of the internal energy density:

$\dot{e}=\alpha\rho_{o}k_{B}\dot{T}-\beta\dot{\rho}k_{B}T_{o}.$

With some other standard choices, namely,

$M_{T}=M_{T,o}T^{2},\quad M_{\rho}=M_{\rho,o}\rho,$

where $M_{T,o}$ and $M_{\rho,o}$ are positive constants, the system becomes

$\displaystyle\alpha\rho_{o}k_{B}\dot{T}-\beta k_{B}T_{o}\dot{\rho}$
$\displaystyle=M_{T,o}\Delta T,$ $\displaystyle\dot{\rho}$
$\displaystyle=M_{\rho,o}\nabla\cdot\left(\rho\nabla\mu\right)$
$\displaystyle\mu$
$\displaystyle=k_{B}\ln\left(\frac{\rho}{\rho_{o}}\right)+k_{B}g^{\prime}\left(\frac{\rho-\rho_{o}}{\rho_{o}}\right)-\beta
k_{B}\frac{T_{o}}{T}+\kappa_{f,\rho,o}\mathcal{C}(\rho-\rho_{o}).$

We can now clearly identify the different terms in the chemical potential:
$k_{B}\ln\left(\frac{\rho}{\rho_{o}}\right)$ corresponds to the ideal gas
term, whereas $k_{B}g^{\prime}\left(\frac{\rho-\rho_{o}}{\rho_{o}}\right)$
takes care of the correction to the ideal gas law. The term $-\beta
k_{B}\frac{T_{o}}{T}$ relates to the latent heat of fusion, and
$\kappa_{f,\rho,o}\mathcal{C}(\rho-\rho_{o})$ is concerned with the long-range
interaction among particles. The boundary conditions are of local
thermodynamic equilibrium (LTE) type, coupled with the entropy non-decreasing
conditions

$M_{T,o}T^{2}\nabla\left(\frac{1}{T}\right)\cdot\boldsymbol{n}=-\frac{\beta_{T,\partial\Omega}}{T},\quad-
M_{\rho,o}\rho\nabla\mu\cdot\boldsymbol{n}=\beta_{\rho,\partial\Omega}\mu\quad\mbox{on
$\partial\Omega$}.$

Recapitulating the entropy and energy densities, all of them as functions of
$\rho$ and $T$, we have

$\displaystyle f=f(T,\rho)$ $\displaystyle=\rho
k_{B}T\ln\left(\frac{\rho}{\rho_{o}}\right)-k_{B}T\left(\rho-\rho_{o}\right)+\rho_{o}k_{B}Tg\left(\frac{\rho-\rho_{o}}{\rho_{o}}\right)-\alpha\rho_{o}k_{B}T\ln\left(\frac{T}{T_{o}}\right)-\beta\rho
k_{B}T_{o},$ $\displaystyle e=\hat{e}(T,\rho)$
$\displaystyle=\alpha\rho_{o}k_{B}T-\beta\rho k_{B}T_{o},$ $\displaystyle
s=\hat{s}(T,\rho)$ $\displaystyle=-\rho
k_{B}\ln\left(\frac{\rho}{\rho_{o}}\right)+k_{B}\left(\rho-\rho_{o}\right)-\rho_{o}k_{B}g\left(\frac{\rho-\rho_{o}}{\rho_{o}}\right)+\alpha\rho_{o}k_{B}\left(1+\ln\left(\frac{T}{T_{o}}\right)\right).$

###  2.5 Non-dimensionalization

By appropriate rescaling, we obtain the following non-dimensional version of
the model:

$\displaystyle\alpha\dot{T}-\beta\dot{\rho}$
$\displaystyle=-M\nabla\cdot\left(T^{2}\nabla\left(\frac{1}{T}\right)\right),$
(2.4) $\displaystyle\dot{\rho}$
$\displaystyle=\nabla\cdot\left(\rho\nabla\mu\right),$ $\displaystyle\mu$
$\displaystyle=\ln\left(\rho\right)+g^{\prime}(\rho-1)-\frac{\beta}{T}+\kappa\gamma(\rho-1)+\kappa\left(2\Delta\rho+\Delta^{2}\rho\right),$

where $T\longrightarrow T/T_{o}$, $\rho\longrightarrow\rho/\rho_{o}$, and $M$,
$\alpha$, $\beta$, $\varepsilon$, and $\kappa$, are positive, non-dimensional
constants. The dimensionless constant $\gamma$ can be positive, negative, or
zero. The dimensionless boundary conditions are of LTE type, plus the
following entropy producing conditions

$-T^{2}\nabla\left(\frac{1}{T}\right)\cdot\boldsymbol{n}=\frac{\beta_{T,\partial\Omega}}{T},\quad-\rho\nabla\mu\cdot\boldsymbol{n}=\beta_{\rho,\partial\Omega}\mu\quad\mbox{on
$\partial\Omega$},$

where $\beta_{T,\partial\Omega}\geq 0$ and $\beta_{\rho,\partial\Omega}\geq
0$. The dimensionless entropy is

$S=\int_{\Omega}\left\\{-\rho\ln\left(\rho\right)+\rho-1-g(\rho-1)+\alpha\left(1+\ln\left(T\right)\right)-\frac{\kappa}{2}\left(\gamma(\rho-1)^{2}-2|\nabla\rho|^{2}+(\Delta\rho)^{2}\right)\right\\}d\boldsymbol{x},$

and the entropy production rate is

$\dot{S}=\int_{\Omega}\left\\{\rho\left|\nabla\mu\right|^{2}+MT^{2}\left|\nabla\left(\frac{1}{T}\right)\right|^{2}\right\\}d\boldsymbol{x}+\int_{\partial\Omega}\left\\{\beta_{\rho,\partial\Omega}\mu^{2}+\frac{\beta_{T,\partial\Omega}}{T^{2}}\right\\}d\boldsymbol{a}\geq
0.$ (2.5)

##  3 Free Energies and Phase Diagrams

Let us examine the Helmholtz free energies of the liquid and solid phases,
which inform the phase diagrams of the material. We say phase diagrams because
each choice of the parameter set will give a different phase diagram.
Understanding the free energy landscape will help us calibrate the melting
temperature, $T_{M}$, as well as the equilibrium values of the fields
$\rho_{S}$ and $\rho_{L}$. In particular, if we take $T_{o}=T_{M}$, at the
melting temperature we should have $T=1.0$. We will show how to adjust
parameters so this corresponds to the correct physical case. We will first
recount the method to approximate the free energy analytically [21], and later
we will explain how to get more accurate calculations numerically.

We will start by considering the non-dimensional free energy at the uniform,
dimensionless temperature $T$:

$\displaystyle
F[T,\rho]=\int_{\Omega}\bigg{\\{}T\rho\ln\left(\rho\right)-T\left(\rho-1\right)+g(\rho-1)-\alpha
T\ln\left(T\right)-\beta\rho+\frac{\kappa
T}{2}\left(\gamma(\rho-1)^{2}-2|\nabla\rho|^{2}+(\Delta\rho)^{2}\right)\bigg{\\}}d\boldsymbol{x}.$

The phase diagram can be constructed by minimizing the free energy, or,
equivalently, maximizing the entropy, as long as the temperature is uniform,
which we assume in this section. Now, we will make a small deviation
approximation. Suppose that

$\rho=1+\psi,\quad|\psi|\ll 1.$

For simplicity, let us assume for now that the ideal gas deviation term, $g$,
is identically zero. (We will take up the more general case in future papers.)
Then, using Taylor’s Theorem,

$T\rho\ln(\rho)-T(\rho-1)-\beta\rho+\frac{\kappa\gamma
T}{2}(\rho-1)^{2}=\frac{\lambda
T}{2}\psi^{2}-\frac{T}{6}\psi^{3}+\frac{T}{12}\psi^{4}-\beta(\psi+1)+\mathcal{O}(\psi^{5}),$
(3.1)

where

$\lambda:=1+\gamma\kappa.$

Thus, the approximate free energy is

$\mathcal{F}[T,\psi]=\int_{\Omega}\left\\{\frac{\lambda
T}{2}\psi^{2}-\frac{T}{6}\psi^{3}+\frac{T}{12}\psi^{4}-\beta(\psi+1)-\alpha
T\log(T)+\frac{\kappa
T}{2}\left(-2|\nabla\psi|^{2}+(\Delta\psi)^{2}\right)\right\\}d\boldsymbol{x}.$
(3.2)

We will construct a 2D phase diagram using the procedure outlined in [21]. The
3D version can be done in an analogous way. See, for example, [21]. In the
crystalline phase, we observe through computations that the density field has
a spatially oscillatory “equilibrium” solution that is, up to rotations,
approximately of the form

$\psi(x,y)\approx\chi(x,y):=A+B\left(\cos\left(\frac{2\pi
x}{p}\right)\cos\left(\frac{2\pi
y}{\sqrt{3}p}\right)+\frac{1}{2}\cos\left(\frac{4\pi
y}{\sqrt{3}p}\right)\right),$ (3.3)

where $A$, $B$, and $p>0$ are parameters to be determined. The approximation
$\chi$ is sometimes called the crystal ansatz. Observe that

$\overline{\chi}=A,$

that is, $A$ is the average density and $B$ is the amplitude of oscillations.
This solution has hexagonal symmetry. Its peaks form a hexagonal Bravais
lattice. (See Figure 6.) Next, fixing $\alpha$, $\beta$, $\kappa$, and
$\lambda$, we define

$\displaystyle G(A,B,p,T)$
$\displaystyle:=\frac{2}{\sqrt{3}p^{2}}\int_{0}^{\frac{p\sqrt{3}}{2}}\int_{0}^{p}\left\\{\frac{\lambda
T}{2}\chi^{2}-\frac{T}{6}\chi^{3}+\frac{T}{12}\chi^{4}-\beta(\chi+1)-\alpha
T\log(T)+\frac{\kappa
T}{2}\left(-2|\nabla\chi|^{2}+\left(\Delta\chi\right)^{2}\right)\right\\}dxdy,$

which represents the free energy at a constant temperature $T$ evaluated at
the approximate solution $\chi$, averaged over the crystal’s unit cell, that
is, the smallest repeat unit of $\chi$.

Taking the derivative with respect to $p$, setting this equal to zero, and
solving for $p$, we have

$\frac{\partial G}{\partial
p}=\frac{4B^{2}\kappa(3p^{2}\pi^{2}-16\pi^{4})T}{3p^{5}}=0\quad\implies\quad
p=p_{\rm eq}:=\frac{4\pi}{\sqrt{3}}.$

This gives

$\displaystyle G(A,B,p_{\rm eq},T)$
$\displaystyle=-\beta(1+A)+\frac{T}{1536}\bigg{(}-256A^{3}-48B^{2}(B+6\kappa-6\lambda)+128A^{4}+45B^{4}$
$\displaystyle\quad+96AB^{2}(-3+B)+96A^{2}(8\lambda+3B^{2})\bigg{)}-\alpha
T\log(T).$

Next, taking the derivative with respect to $B$, we have

$\frac{\partial}{\partial B}G(A,B,p_{\rm
eq},T)=\frac{3BT}{128}\left(-4(4A+B+4\kappa-4\lambda)+(16A^{2}+8AB+5B^{2})\right).$

This cubic equation (with respect to $B$) has the solutions $B=0$, which
represents the liquid state, and the two crystalline solutions

$B=B_{\pm}:=\frac{2-4A\pm 2\sqrt{1+16A+20(\kappa-\lambda)-16A^{2}}}{5}.$ (3.4)

We will take $B=B_{+}$, since, it turns out, using $B=B_{-}$ results in a
higher free energy for the crystal ansatz.

Figure 1: Plots of the free energies (per unit cell) $\mathcal{F}_{S}$ (red)
and $\mathcal{F}_{L}$ (blue) for the parameters $T=1.0$, $\lambda=1.0$,
$\kappa=0.919285$, $\beta=0.025$, and $\alpha=0.1$. At the melting
temperature, $T=1.0$, the minimum values of $\mathcal{F}_{S}$ and
$\mathcal{F}_{L}$ are equal. The orange line is $-0.025315$.

Figure 2: Plots of the free energies (per unit cell) $\mathcal{F}_{S}$ and
$\mathcal{F}_{L}$ for the parameters $T=0.9$ (left), $T=1.1$ (right), and
$\lambda=1.0$, $\kappa=0.919285$, $\beta=0.025$, and $\alpha=0.1$. When the
temperature is above (below) the melting temperature, $T=1.0$, the liquid
phase has lower (higher) free energy. The orange line is $-0.025315-\alpha
T\log(T)$.

Now, the liquid free energy per unit cell can be obtained simply by setting
$B=0$. Thus we define

$\mathcal{F}_{L}(A,T)=G(A,0,p,T)=\frac{\lambda
T}{2}A^{2}-\frac{T}{6}A^{3}+\frac{T}{12}A^{4}-\beta(A+1)-\alpha T\log(T).$
(3.5)

This expression is exact, assuming that (3.2) is the exact free energy
expression. The crystal phase free energy per unit cell is approximated as

$\mathcal{F}_{S}(A,T)=G(A,B_{+},p_{\rm eq},T).$

Now, the melting temperature, $T=1.0$, should be, by definition, the
temperature at which the minimum values of $\mathcal{F}_{S}(A,T)$ and
$\mathcal{F}_{L}(A,T)$ are equal. Thus, for a given parameter set $\lambda$,
$\kappa$, $\beta$, $\alpha$,

$\mathcal{F}_{S}(A_{S},1)=\mathcal{F}_{L}(A_{L},1),$ (3.6)

where

$A_{L}:=\operatorname{argmin}_{A}\mathcal{F}_{L}(A,1)\quad\mbox{and}\quad
A_{S}:=\operatorname{argmin}_{A}\mathcal{F}_{S}(A,1).$

The parameters can be carefully adjusted to make equality in (3.6) happen. We
will show a couple of examples.

In Figure 1, the free energies (per unit cell) $\mathcal{F}_{S}$ (red) and
$\mathcal{F}_{L}$ (blue) are plotted for the parameters $T=1.0$,
$\lambda=1.0$, $\kappa=0.919285$, $\beta=0.025$, and $\alpha=0.1$. At the
melting temperature, $T=1.0$, the minimum values of $\mathcal{F}_{S}$ and
$\mathcal{F}_{L}$ are equal, as shown. If the temperature is below the melting
temperature, the free energy of the solid phase should be lower. This case is
shown in Figure 2(left). Conversely, if the temperature is above the melting
temperature, the free energy of the liquid phase should be lower than that of
the solid phase. This case is shown in Figure 2(right). Thus the model is
flexible enough so that (i) the melting temperature can be calibrated so that
(3.6) holds at $T=1.0$; (ii) if $T<1.0$, the solid phase has lower free energy
than the liquid phase; and (iii) if $T>1.0$, the solid phase has higher free
energy than the liquid phase. Our latent heat term is simple in this model; it
is linear in $\rho$ ($-\beta\rho$). But it is, seemingly, sophisticated enough
to capture the basic physics near the melting temperature.

To compute phase diagrams, one uses the Maxwell common-tangents construction.
The common tangent approach is appropriate when the number density of
particles is a conserved quantity, as is typical with the PFC modeling
framework. The Maxwell construction method is described in detail in [8] and
is understood widely in the Physics literature, so we will not describe it
here. See also [21].

Figure 3: Plots of the free energies (per unit cell) $\mathcal{F}^{\rm E}_{S}$
(red), $\mathcal{F}_{S}$ (blue), and $\mathcal{F}_{L}$ (yellow) for the
parameters $T=1.0$, $\lambda=1.0$, $\kappa=0.917385$, $\beta=0.025$, and
$\alpha=0.1$. At the (recalibrated) melting temperature, $T=1.0$, the minimum
values of $\mathcal{F}^{\rm E}_{S}$ and $\mathcal{F}_{L}$ are equal. The
purple line is $-0.025315$. Notice that, as expected, the approximate solution
obtained by the ansatz (3.3) yields a higher free energy. We adjusted $\kappa$
lower (relative to the value in Figure 1) in order to recalibrate the
parameters so that $T=1.0$ was the “exact” melting temperature.

Figure 4: Plots of the free energies (per unit cell) $\mathcal{F}^{\rm E}_{S}$
(red), $\mathcal{F}_{S}$ (blue), and $\mathcal{F}_{L}$ (yellow) for the
parameters $T=0.9$ (left), $T=1.1$ (right), and $\lambda=1.0$,
$\kappa=0.917385$, $\beta=0.025$, and $\alpha=0.1$. The uniform temperature is
below (above) the recalibrated melting temperature, and the minimum value of
$\mathcal{F}^{\rm E}_{S}$ is below (above) that of $\mathcal{F}_{L}$. The
purple line is $-0.025315-\alpha T\log(T)$.

We have computed the free energy landscape using the crystal solution ansatz
(3.3), which is only an approximation of the solid phase minimizer. In
particular, the approximated free energy is expected to be too large. How big
of an error is incurred by using the approximation? To find out, let us
describe a more accurate method for computing the free energy landscape. In
particular, we need a more accurate representation of the solid phase free
energy. For illustration purposes, let us suppose that the free energy,
$\mathcal{F}$, as given in (3.2), is the “exact” free energy. The liquid free
energy per unit cell (3.5) is the same as before, because no approximation was
used to obtain it. For the crystalline (solid) state, instead of using the
ansatz (3.3), we will use numerical calculations to find the exact energy (per
unit cell)

$\mathcal{F}^{\rm E}_{S}(A,T),$

where $A:=\overline{\rho-1}$, as a function of $A$, at a uniform fixed
temperature, $T$, and a given set of parameters $\lambda$, $\kappa$, $\beta$,
$\alpha$. The procedure for calculating $\mathcal{F}^{\rm E}_{S}$ is outlined
as follows: fixing the parameters $p$, $A$, $T$, $\lambda$, $\kappa$, $\beta$,
and $\alpha$, we solve the system

$\dot{\rho}=\Delta\mu,\quad\mu=\delta_{\rho}\mathcal{F}[T,\rho],\quad\mbox{on}\quad\Omega=\Omega_{p}:=[0,p]\times\left[0,\frac{\sqrt{3}}{2}p\right]$
(3.7)

subject to periodic boundary conditions, with the initial conditions given by
the crystal ansatz:

$\rho(x,y,t=0)=1+A+B_{+}\left(\cos\left(\frac{2\pi
x}{p}\right)\cos\left(\frac{2\pi
y}{\sqrt{3}p}\right)+\frac{1}{2}\cos\left(\frac{4\pi
y}{\sqrt{3}p}\right)\right),$

where $B_{+}$ is a given in (3.4). Clearly,

$\overline{\rho(\,\cdot\,,0)-1}=A.$

We solve (3.7) to steady state (equilibrium) and compute the energy of the
solution on the unit cell, obtaining the equilibrium solutions

$\rho_{\infty}\equiv\rho(\,\cdot\,,t=\infty)\quad\mbox{on}\quad[0,p]\times\left[0,\frac{\sqrt{3}}{2}p\right]$

and

$\displaystyle F_{\infty}(A,T,p)$
$\displaystyle:=\int_{0}^{p}\int_{0}^{\sqrt{3}p/2}\bigg{\\{}\frac{\lambda
T}{2}(\rho_{\infty}-1)^{2}-\frac{T}{6}(\rho_{\infty}-1)^{3}+\frac{T}{12}(\rho_{\infty}-1)^{4}-\beta\rho_{\infty}$
$\displaystyle\quad-\alpha T\log(T)+\frac{\kappa
T}{2}\left(-2|\nabla\rho_{\infty}|^{2}+(\Delta\rho_{\infty})^{2}\right)\bigg{\\}}dydx.$

Finally,

$\mathcal{F}^{\rm E}_{S}(A,T):=F_{\infty}(A,T,p^{\rm E}_{\rm eq}),\quad p^{\rm
E}_{\rm eq}:=\operatorname{argmin}_{p}F_{\infty}(A,T,p).$

We use a pseudo-spectral, stabilized implicit-explicit (IMEX) method to
perform the equilibrium field calculations. The spectrally accurate
trapezoidal rule is used to compute the energies. The minimization problem in
$p$ is solved using a derivative-free method [3] over the interval $[p_{\rm
eq}-1,p_{\rm eq}+1]$, where $p_{\rm eq}=\frac{4\pi}{\sqrt{3}}$.

The results of our improved free energy computations are shown in Figures 3,
$T=1.0$; 4(left), $T=1.1$; and 4(right), $T=0.9$. Notice that the approximate
solution yields a larger free energy than the true minimizer, as expected. To
recalibrate the model so that $T=1.0$ occurs when $\mathcal{F}^{\rm
E}_{S}(A_{s}^{E},T)=\mathcal{F}_{L}(A_{L},T)$, we need only to adjust the
value of $\kappa$ lower. The approximate free energy landscape computed using
the ansatz is still quite useful, since the resulting free energies are
reasonably good ballpark estimates that can help tune the parameters. The more
accurate method can then refine and recalibrate the parameters after the
ballpark estimates are obtained.

##  4 Numerical Solution of the Model

To conclude this paper, let us perform a simple computation to give just a
small sample of what this model can do. We use the polynomial free energy
(3.2) as our “exact” free energy; the associated entropy is

$\mathcal{S}=\int_{\Omega}\left\\{-\left(\frac{\lambda}{2}\psi^{2}-\frac{1}{6}\psi^{3}+\frac{1}{12}\psi^{4}\right)+\alpha\left(1+\ln\left(T\right)\right)-\frac{\kappa}{2}\left(-2|\nabla\psi|^{2}+(\Delta\psi)^{2}\right)\right\\}d\boldsymbol{x};$

and the associated internal energy becomes

$\mathcal{E}=\int_{\Omega}\left\\{\alpha
T-\beta(\psi+1)\right\\}d\boldsymbol{x},$

where we have made use of the change of variable

$\psi:=\rho-1.$

The simplified evolution equations are

$\displaystyle\alpha\dot{T}-\beta\dot{\psi}$ $\displaystyle=M\Delta T,$ (4.1)
$\displaystyle\dot{\psi}$ $\displaystyle=\Delta\mu,$ $\displaystyle\mu$
$\displaystyle=\frac{1}{T}\delta_{\psi}\mathcal{F}=\lambda\psi-\frac{1}{2}\psi^{2}+\frac{1}{3}\psi^{3}-\frac{\beta}{T}+\kappa\left(2\Delta\psi+\Delta^{2}\psi\right),$

and we use periodic boundary conditions for simplicity. We leave it to the
reader to show that this system is still entropy producing, locally and
globally. Mass is conserved in our simulation, and therefore, this simplified
setup disallows freezing in the usual sense. This is because the solid and
liquid phases have distinct average densities. If we start out with the
equilibrium values of the solid and liquid states and then drop the
temperature below the melting temperature, the solid may not grow much, if at
all. We would have to add mass to the system in order for the solid state to
grow, since the solid state has a higher average density. We will save such
sophisticated simulations for a future paper, where we will explore the model
further.

For the numerical solution, we use a Fourier pseudo-spectral discretization of
space, coupled with a stabilized, linear, first-order IMEX algorithm for time
discretization. This scheme is not designed to keep the temperature field
positive, and will not guarantee entropy production. More sophisticated
schemes will be developed in the future that guarantee these properties
theoretically. See, for example, the related numerical works of the PFC model
with constant temperature [9, 16, 25], the modified PFC equation [1, 2], the
square PFC equation [6]. We perform a singe test, the results for which are
shown in Figure 6. The parameters for the test are given in the captions of
Figure 5, where we give a cartoon description of the free energy landscape for
the parameters.

We seed the center of the domain with a two-grain crystal that is very near to
its equilibrium state at $T=1.0$, with an average density of $A=0.16$. (See
Figure 5(left).) The liquid phase surrounding the crystal seed is
supersaturated, meaning its average density (at $A=0.14$) is higher than that
of the equilibrium state (just below $A=0.11$, see, again, Figure 5(left)).
The initial temperature is roughly $T=0.97$. This shifts the equilibria
slightly; see Figure 5(right). The new equilibrium average would be computed
using the Maxwell common-tangent construction, of course. The extra mass is
ejected from the liquid and attaches to the solid seed, and the seed grows. As
it does the temperature changes due the the release of latent heat, and we see
a commonly observed Gibbs-Thompson-like temperature jump effect at the
boundary of the crystal.

## Acknowledgements

This work is partially supported by the National Science Foundation (USA)
grants NSF DMS-2012669 (C. Wang) and NSF DMS-1719854, DMS-2012634 (S.M. Wise).
SMW thanks Ken Elder for bringing reference [17] to his attention and for
several discussions about PFC models.

Figure 5: Plots of the free energies (per unit cell) $\mathcal{F}^{\rm E}_{S}$
(red), $\mathcal{F}_{S}$ (blue), and $\mathcal{F}_{L}$ (yellow) for the
parameters $T=1.0$ (left), $T=0.97$ (right), and $\lambda=0.6$,
$\kappa=0.459810$, $\beta=0.06$, and $\alpha=0.1$. The uniform temperature is
equal to (below) the recalibrated melting temperature, and the minimum value
of $\mathcal{F}^{\rm E}_{S}$ is equal to (below) that of $\mathcal{F}_{L}$.
The purple line is $-0.0631798-\alpha T\log(T)$.

$t=5.0$

$t=10$

$t=50$

$t=100$

Figure 6: A unary crystal growing in its supersaturated liquid phase: (left
column) the $\psi$ field and (right column) the temperature field. Parameters
for the test are given in the caption of Figure 5. In addition,
$\Omega=(0,400)^{2}$ and $M=0.1$.

## References

  * [1] A. Baskaran, Z. Hu, J. Lowengrub, C. Wang, S.M. Wise, and P. Zhou.  Energy stable and efficient finite-difference nonlinear multigrid schemes for the modified phase field crystal equation.  J. Comput. Phys., 250:270–292, 2013.
  * [2] A. Baskaran, J. Lowengrub, C. Wang, and S. Wise.  Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation.  SIAM J. Numer. Anal., 51:2851–2873, 2013.
  * [3] R.P. Brent.  Algorithms for Minimization Without Derivatives.  Prentice-Hall, Englewood Cliffs, New Jersey, USA, 1973.
  * [4] C. Charach and P.C. Fife.  On thermodynamically consistent schemes for phase field equations.  Open Systems and Information Dynamics, 5:99–123, 1998.
  * [5] W. Chen, C. Wang, X. Wang, and S.M. Wise.  Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential.  J. Comput. Phys.: X, 3:100031, 2019.
  * [6] K. Cheng, C. Wang, and S.M. Wise.  An energy stable fourier pseudo-spectral numerical scheme for the square phase field crystal equation.  Commun. Comput. Phys., 26:1335–1364, 2019.
  * [7] S.R. de Groot and P. Mazur.  Non-Equilibrium Thermodynamics.  North Holland, Amsterdam, Holland, 1962.
  * [8] R.T. DeHoff.  Thermodynamics in Materials Science.  McGraw-Hill, New York, New Yorkk, USA, 1962.
  * [9] L. Dong, W. Feng, C. Wang, S.M. Wise, and Z. Zhang.  Convergence analysis and numerical implementation of a second order numerical scheme for the three-dimensional phase field crystal equation.  Comput. Math. Appl., 75(6):1912–1928, 2018.
  * [10] L. Dong, C. Wang, S.M. Wise, and Z. Zhang.  A positivity-preserving, energy stable scheme for a ternary Cahn-Hilliard system with the singular interfacial parameters.  J. Comput. Phys., 442:110451, 2021.
  * [11] L. Dong, C. Wang, H. Zhang, and Z. Zhang.  A positivity-preserving, energy stable and convergent numerical scheme for the Cahn-Hilliard equation with a Flory-Huggins-deGennes energy.  Commun. Math. Sci., 17:921–939, 2019.
  * [12] L. Dong, C. Wang, H. Zhang, and Z. Zhang.  A positivity-preserving second-order BDF scheme for the Cahn-Hilliard equation with variable interfacial parameters.  Commun. Comput. Phys., 28:967–998, 2020.
  * [13] K.R. Elder, M. Katakowski, M. Haataja, and M. Grant.  Modeling elasticity in crystal growth.  Phys. Rev. Lett., 88:245701, 2002.
  * [14] K.R. Elder, M. Katakowski, M. Haataja, and M. Grant.  Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals.  Phys. Rev. E, 70:051605, 2004.
  * [15] K.R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant.  Phase-field crystal modeling and classical density functional theory of freezing.  Phys. Rev. B, 77:064107, 2007.
  * [16] Z. Hu, S.M. Wise, C. Wang, and J.S. Lowengrub.  Stable and efficient finite-difference nonlinear-multigrid schemes for the phase-field crystal equation.  J. Comput. Phys., 228:5323–5339, 2009.
  * [17] G. Kocher and N. Provatas.  Thermodensity coupling in phase-field-crystal-type models for the study of rapid crystallization.  Phys. Rev. Materials, 3:053804, 2019.
  * [18] C. Liu, C. Wang, and Y. Wang.  A structure-preserving, operator splitting scheme for reaction-diffusion equations with detailed balance.  J. Comput. Phys., 436:110253, 2021.
  * [19] C. Liu, C. Wang, S.M. Wise, X. Yue, and S. Zhou.  A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst-Planck system.  Math. Comp., 2021.  accepted and published online: https://doi.org/10.1090/mcom/3642.
  * [20] N. Provatas, J.A. Dantzig, B. Athreya, P. Chan, P. Stefanovic, N. Goldenfeld, and K.R. Elder.  Using the phase-field crystal method in the multiscale modeling of microstructure evolution.  JOM, 59:83, 2007.
  * [21] N. Provatas and K. Elder.  Phase-Field Methods in Materials Science and Engineering.  Wiley-VCH Verlag, Weinheim, Germany, 2010.
  * [22] Y. Qian, C. Wang, and S. Zhou.  A positive and energy stable numerical scheme for the Poisson-Nernst-Planck-Cahn-Hilliard equations with steric interactions.  J. Comput. Phys., 426:109908, 2021.
  * [23] S.L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Coriell, R.J. Braun, and G.B. McFadden.  Thermodynamically-consistent phase-field models for solidification.  Physica D, 69:189–200, 1993.
  * [24] A.A. Wheeler, G.B. McFadden, and W.J. Boettinger.  Phase-field model for solidification of a eutectic alloy.  Proc. R. Soc. London, 452:495–525, 1996.
  * [25] S.M. Wise, C. Wang, and J. Lowengrub.  An energy stable and convergent finite-difference scheme for the phase field crystal equation.  SIAM J. Numer. Anal., 47:2269–2288, 2009.
  * [26] M. Yuan, W. Chen, C. Wang, S.M. Wise, and Z. Zhang.  An energy stable finite element scheme for the three-component Cahn-Hilliard-type model for macromolecular microsphere composite hydrogels.  J. Sci. Comput., 87:78, 2021.