\documentclass[twoside]{article}
\usepackage{amssymb} % font used for R in Real numbers
\pagestyle{myheadings}
\markboth{\hfil Models for phase separation \hfil EJDE--2000/48}
{EJDE--2000/48\hfil Paul C. Fife \hfil}
\begin{document}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
{\sc Electronic Journal of Differential Equations},
Vol.~{\bf 2000}(2000), No.~48, pp.~1--26. \newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu \quad ftp ejde.math.unt.edu (login: ftp)}
\vspace{\bigskipamount} \\
%
Models for phase separation and \\ their mathematics
\thanks{ {\em Mathematics Subject Classifications:} 35K55, 80M35, 74A25.
\hfil\break\indent
{\em Key words:} Cahn-Hilliard equation, phase-field equations,
phase separation, gradient flows.
\hfil\break\indent
\copyright 2000 Southwest Texas State University and University of
North Texas. \hfil\break\indent
Submitted February 10, 2000. Published June 22, 2000. \hfil\break\indent
This paper (except for the last section,
``Postscript'', which was added later) was part of the proceedings of
a workshop held in Katata, Japan in August, 1990 [Nonlinear
Partial Differential Equations with Applications to
Patterns, Waves, and Interfaces, M. Mimura and T. Nishida,
eds.] In its revised form the paper was submitted in June,
1991. Those proceedings were slated to be published, but in
fact never were.
} }
\date{}
%
\author{Paul C. Fife}
\maketitle
\begin{abstract}
The gradient flow approach to the Cahn-Hilliard and phase field
models is developed, and some basic mathematical properties of
the models, especially phase separation phenomena, are reviewed.
\end{abstract}
\newcommand{\pref}[1]{(\ref{#1})}
\newcommand{\bm}[1]{\mbox{\boldmath $#1$}}
\newtheorem{thm}{Theorem}
%\newtheorem{lm}{Lemma}
%\newtheorem{prop}{Proposition}
\section{ Introduction}
Materials exhibiting a fine mixture of phases are a
familiar and well studied phenomenon in many settings. Such
mixtures can arise for a variety of reasons, but the process
upon which I mainly focus here is that of ``spinodal
decomposition'', roughly described as follows. Consider a
homogeneous substance, initially homogeneous, with
controllable thermodynamic properties. Those properties are
adjusted so that the substance enters an extremely unstable
state, such that the unfolding of its instabilities entails
phase separation into a fine-scaled heterogenous mixture of
phases.
The prototypical scenario involves lowering the temperature of a binary
alloy to a level where the alloy can no longer exist in equilibrium in its
homogeneous state. The most basic properties of this important process have
been described with the aid of a model proposed by Cahn and
Hilliard
\cite{CHi,Cah2}, the subject of considerable investigation.
Although originally
proposed in a metallurgical connection, the spontaneous heterogenization
behavior of the Cahn-Hilliard equation has also occasionally captured the
imagination of modelers of pattern formation in other contexts.
Another related model for phase transitions, the phase-field system, has
been proposed and used successfully to gain insight and computational
advantage in studying solid-liquid phase boundaries. One of my objects will
be to sketch a rationale for both models, and to show that the latter, which
is not as well known as the former, possesses a theory of phase separation
which is very much like that of the Cahn-Hilliard model. This analogy is a
mathematical one at the present stage, because the phase field model is no
doubt too simplistic to adequately describe known fine-grained phase
separation events involving the liquid-solid transition.
The two models examined here are also capable of representing a mode of
phase separation quite different from that alluded to above, namely nucleation
and growth. This is described in Sections \ref{s4} and \ref{s9}.
For each model, I present a derivation
(Sections \ref{s2},\ref{s5},\ref{s6}) based on constrained
gradient dynamics for a physically motivated functional. In the first case,
the functional represents a free energy for the system, and in the second
case, the entropy. The constraints are that the total mass (in the first
case) and the total energy (in the second case) be conserved.
I then review (Sections \ref{s3} and \ref{s7}) what is known about the spontaneous phase
separation event for each model. In the case of the phase field equations, I
also describe (Section \ref{8}) some one-dimensional computer simulations of the appearance
and disappearance of fine-mixtures of phases under the action of internal
heating.
Finally, a third model (due to Novick-Cohen and Pego \cite{NP}) is briefly
described in Section \ref{s10}.
\section{ Cahn-Hilliard model} \label{s2}
The Cahn-Hilliard equation
\begin{equation}
\alpha u_t = \Delta\left( -\epsilon^2\Delta u + f(u) \right)~~~(\Delta \equiv \nabla^2) \label{1}
\end{equation}
was proposed in \cite{CHi} and \cite{Cah2} as a simple model
for
the process of phase
separation of a binary alloy at a fixed temperature. The
function
$f(u)$ is of
``bistable type'' with three simple zeros at $u_1$, $u_2$,
and an intermediate value $u_0$. More specifically,
\begin{equation} \begin{array}{l}
f(u) = 0 \mbox{ only at }
u = u_1,~u_0,~u_2.\\
f'(u) > 0,~~u < u_1^* \mbox{ or } u > u_2^*,\\
f'(u) < 0,~~u \in (u_1^*,u_2^*),\end{array} \label{1-a}\end{equation}
where $u_1 < u_1^* < u_0 < u_2^* < u_2$. The three
intervening intervals have names:
$(u_1,u_1^*)$ is metastable interval 1;
$(u_1^*,u_2^*)$ is the spinodal interval;
$(u_2^*,u_2)$ is metastable interval 2.
Here the function $u(x,t)$ represents the concentration
of one of the two metallic components of the alloy. If we
assume the total density is constant for simplicity, the
composition of the mixture may be adequately expressed by
the single function $u$. The parameter $\epsilon$ is an
``interaction length'', small compared to characteristic
dimensions on the laboratory scale, and $\alpha$ is a relaxation
parameter.
If the alloy is contained in a vessel $\Omega$, the equation \pref{1} should be
supplemented with boundary conditions on the boundary
$\partial\Omega$.
These are usually
taken to be
\begin{equation}
\partial_\nu u = \partial_\nu \left( -\epsilon^2\Delta u + f(u) \right) = 0 \mbox{ on } \partial\Omega,
\label{2} \end{equation}
where $\partial_\nu$ denotes differentiation normal to $ \partial\Omega$.
The physical meaning of the
second of these two conditions is that none of the mixture can pass through
the walls of the container. The first condition is the most natural way to
ensure that the total ``free energy''
of the mixture decreases in time (this, as
we shall see below, is a requirement from thermodynamics), when there is no
interaction between the alloy and the containing walls. Since
$\partial_\nu f(u) = f'(u)\partial_\nu u = 0$,
these conditions take the equivalent form
\begin{equation}
\partial_\nu u = \partial_\nu \Delta u = 0 \mbox{ on } \partial\Omega. \label{2'}
\end{equation}
Now suppose a uniform binary mixture at a high temperature and
concentration $u^*$ is suddenly quenched to a given lower temperature. A
commonly occurring situation is that there is a pair of
values of $u$,
say $u_1$ and $u_2$,
and a pair of solid phases, denoted by I and II, say, such that the
two phases can be in equilibrium and in physical contact only if the former
has concentration $u_1$, and the latter has concentration
$u_2$.
It is assumed that
this is the case for the alloy and (quenched) temperature being studied here;
these distinguished concentrations are then represented by the outermost zeros
of the function $f$.
It is further assumed that $u^*$ lies between $u_1$ and
$u_2$,
more specifically in
the spinodal interval shown. It is known that the uniform
mixture $u = u^*$ is
then very unstable, and that the growth of instabilities results in phase
separation, also called spinodal decomposition.
The Cahn-Hilliard model was originally derived using
sound physical concepts to describe this process. However,
an alternate rationale as a gradient flow is possible, as I
will now show. We begin with the observation that the
nonequilibrium mixture $u = u^*$ at the given temperature
will naturally evolve so that some of it becomes a solid in
phase I with $u = u_1$, and the rest becomes phase II with
$u = u_2$. On a simple level of understanding, this
evolving process might suggestively and conveniently be
modeled by proposing that a certain functional ${\cal
F}_0[u]$ decrease in time and approach a minimum. This
functional is
\begin{equation}
{\cal F}_0[u] \equiv \int_\Omega W(u(x))dx, \label{3} \end{equation}
where the function $W(u)$ is smooth and nonnegative
with equal minima at $u = u_1$ and $u_2$.
At this level, any ``double well potential'' $W(u)$
of this type would
do.
Such a framework, based on a time-monotone functional, is especially
natural because thermodynamics itself tells us that for isothermal systems
there should be a ``free energy'' functional which steadily decreases and
approaches a minimum. As desired, the minimum in this case is assumed by
functions $u$ which take on only the two values $u_1$ and
$u_2$.
At the same time
that ${\cal F}_0$ is decreasing,
the total amount of each species in the vessel must
remain equal to the given original amount; this means
that the evolution is
subject to the constraint
\begin{equation}
\int_\Omega u(x,t)dx = \mbox{ constant}. \label{4} \end{equation}
What is desired is a law of evolution for the function
$u(x,t)$ which will ensure that ${\cal F}_0$ decreases while \pref{4}
holds.
The most obvious choice for such a law is to try to make
$u$ evolve in the opposite direction of the gradient of the
functional in \pref{3}, constrained according to \pref{4}, and
in fact this is the approach (but applied to a functional
more general than \pref{3}) which I shall use in motivating
the models discussed in this paper. As explained below,
even using this guideline, there are many possible choices
of evolution law, because of the many Hilbert space settings
available. Unfortunately, it will turn out that if \pref{3}
is used, none of the most standard Hilbert spaces will give
us an evolution law which is both local and well-posed.
However, at least one nonstandard space can be used; that
idea is developed in Section \ref{s10}.
For now, we use a more common device to ameliorate this
unpleasant situation. The device is to add an extra term to
\pref{3} which exacts a penalty for the function $u$'s having
sudden changes with respect to the space variable $x$.
Physically, such a penalty would arise if, in addition to
$u_1$ and $u_2$ being preferred states, there is a tendency
for the material to prefer to be as uniform as possible.
One widely accepted way to impose this penalty is to add a
gradient term to the free energy functional ${\cal F}_0$ to
produce
\begin{equation}
{\cal F}_\epsilon[u] \equiv \int _\Omega \left[ W(u(x)) + \frac12 \epsilon^2 |\nabla u(x) |^2
\right] dx, \label{5} \end{equation}
where $\epsilon$ is a small positive parameter.
Gradient-corrected functionals have been used a great deal in the physics
literature, beginning with van der Waals \cite{W}. See, for
example,
\cite{LL,HHM,HH}, as well as \cite{CHi}.
From here on, for convenience I shall drop the subscript
``$\epsilon$'' from ${\cal F}$.
As mentioned before, it is natural to seek a law of
evolution in the form
\begin{equation}
\frac{\partial u}{\partial t} = -K\mbox{grad}_0{\cal F}[u], \label{6} \end{equation}
for some positive constant or function $K$. The symbol
``grad$_0$'' here denotes a
constrained gradient, namely the gradient on the manifold, in some Hilbert
space, defined by \pref{4}.
Without giving a general discussion of gradients, I say what will be meant
by grad$_0{\cal F}[u]$ in the present context. The domain ${\cal
D}$ of the functional ${\cal F}$ in \pref{5}
will be taken to be smooth enough functions, defined in $\Omega$ and satisfying the
first of the two conditions \pref{2}. Let $M_0$ be the linear
manifold in ${\cal D}$ of
functions $v$ satisfying $\int v(x)dx = 0$. Given any $u
\in {\cal D}$, let $M_u$ be the parallel
manifold $u + M_0$.
Let $H$ be a Hilbert space of functions containing
${\cal D}$. Then for a function
$u_0 \in {\cal D}$, grad$_0{\cal F}_0{u_0}$, if it exists, is a
function $w \in M_0$ such that for any
function $u(t)$ smooth in $t$ (as well as in $x$),
defined for small real $t$ and
taking values in $M_{u_0}$ with $u(0) = u_0$,
\[
\frac{d}{dt}\left.{\cal F}[u]\right|_{t=0} = \left\langle w,
\frac{du}{dt}(0) \right\rangle_H. \]
We can now see the reason for \pref{6}: it guarantees that
the free energy will decrease in time. In fact, for
functions forced to evolve on any such parallel manifold,
\[
\frac{d}{dt}{\cal F}[u] = \left\langle \mbox{grad}_0{\cal F}[u],\frac{\partial u}{\partial
t} \right\rangle = -K \|\mbox{grad}_0{\cal F}[u] \|^2 \le 0. \]
The scalar product and norm here are with respect to the chosen Hilbert space,
and as we shall see, the choice of a physically relevant space is crucial.
Given any Hilbert space $H$ of functions of $x$
containing ${\cal D}$, one can set
out to construct the constrained gradient grad$_0{\cal F}$.
In the case of ${\cal L}_2$, we
obtain a nonlocal operator, which may be rejected on the physical grounds that
the evolution law for the distribution $u$ should be local: there should be no
``action at a distance''. Higher order Hilbert spaces
$H^k,~k > 0$, have the
same difficulty. Choosing $H$ to be the zero-average
subspace of the dual of $H^1$, however, produces a more
reasonable result, namely the Cahn-Hilliard equation
\pref{1}. For convenience we denote this subspace by the
symbol $\bar{H}^{-1}$; it is different from $H^{-1}$.
To understand this claim, recall that the constrained
gradient operator grad$_0{\cal F}$ should
satisfy the following: For any smooth $v(x)$ satisfying
\begin{equation}
\int v(x)dx = 0 \label{7} \end{equation}
and the first of \pref{2}, we should have
\begin{equation}
\frac{d}{dt}\left.{\cal F}[u + tv]\right|_{t=0} =
\left\langle\mbox{grad}_0{\cal F}[u],v \right\rangle. \label{8} \end{equation}
Now \pref{7} holds if and only if the homogeneous Neumann
problem for
\[
\Delta\Phi(x) = v(x)\]
has a (unique) solution $\Phi(x)$ satisfying $\int\Phi(x)\,dx =
0$.
Employing this
expression for $v$ and integrating by parts, we obtain
\[
\frac{d}{dt}\left.{\cal F}[u + t\Delta\Phi] \right|_{t=0} = \int_\Omega\left[
W'(u) - \epsilon^2\Delta u\right] \Delta\Phi\,dx;\]
the boundary term vanishes because of the assumed first
condition in \pref{2}.
A
second integration by parts then yields
\begin{equation}
-\int_\Omega\nabla\left[ W'(u) - \epsilon^2\Delta u\right] \cdot\nabla\Phi\,dx, \label{9}
\end{equation}
where I have used the fact that $\nabla\Phi$ has zero normal component on $\partial\Omega$.
Now let $H = \bar{H}^{-1}$ be as described above, so
that its functions satisfy \pref{7}. The scalar product
in $H$ can be defined by
\begin{equation}
\langle v_1,v_2 \rangle = \langle \nabla\Phi_1,\nabla\Phi_2 \rangle_{{\cal L}_2},
\label{10}\end{equation}
where $\Phi_1$ and $\Phi_2$ are related to $v_1,~v_2 \in
\bar{H}^{-1}$ as indicated above.
The expression in \pref{9} is the ${\cal L}_2$-scalar
product of $
-\nabla \left[ W'(u) - \epsilon^2\Delta u\right] $ with $\nabla\Phi$.
We apply \pref{10} with these two vector
fields playing the roles of $\nabla\Phi_1$ and $\nabla\Phi_2$
respectively. To apply \pref{10}, we require that the first field have zero
normal component on $\partial\Omega$; hence the second condition in
\pref{2}.
We thus obtain,
from \pref{8}, that \pref{9} equals
\begin{eqnarray*}
\left\langle -\nabla\cdot\nabla\left[ W'(u) - \epsilon^2\Delta u \right], \nabla\cdot\nabla\Phi
\right\rangle_H &=& \left\langle -\Delta \left[ W'(u) - \epsilon^2 \Delta u
\right], v \right\rangle_H \\
&=& \langle\mbox{grad}_0{\cal F}[u],v \rangle_H.
\end{eqnarray*}
In the space $H = \bar{H}^{-1}$ we therefore identify
\begin{equation}
\mbox{grad}_0{\cal F}[u] \equiv -\Delta \left[ W'(u) - \epsilon^2\Delta u \right], \label{11} \end{equation}
and specify its domain as those functions in ${\cal D}$ which satisfy \pref{2}.
The ansatz \pref{6} now gives us the law of motion
\[
\frac{\partial u}{\partial t} = K \Delta \left[ W'(u) - \epsilon^2 \Delta u\right], \]
which is the Cahn Hilliard equation \pref{1} when $K$
is chosen to be a (positive)
constant and $f(u) = W'(u)$. In this regard, notice that $f(u)$ has the
appropriate bistable properties if $W(u)$ is smooth, has local minima at
$u = u_1$ and $u_2$, and has only one local maximum between
them.
Throughout the
rest of the paper, I shall assume that $f(u)$ is a cubic, odd with respect to
its middle zero. This assumption is merely to avoid certain technicalities.
It is now clear why the penalty term in \pref{5} is necesssary if one insists
on using this Hilbert space $\bar{H}^{-1}$ to get a
deterministic well-posed evolution
law for $u$. (Similar objections apply when using $H^k$ for
any integer
$k$.) If
we set $\epsilon = 0$, thus reverting to the functional \pref{3},
the resulting evolution
law becomes \pref{1} with $\epsilon = 0$. Since $f$ is not
monotone,
this is a forward-backward
nonlinear diffusion equation with the spinodal interval
representing a regime of ``negative diffusion''. This does not satisfy the
requirement of determinism. Nevertheless, this and similar equations have
been studied.
Further discussions of physical aspects of this and related models can be
found in many papers, such as
\cite{Cah3,Cah4,Cah5,GD,HHM,HH,L1,L2,LS,
Pen,SS}. Some basic mathematical properties of the equation are
found, among other places, in
\cite{ABF,BF1,CGS,CN,EFG,E,EF,EZ1,GM,N,NSe,NS,
NST1,NST2,P,Z}.
\section{Spinodal decomposition} \label{s3}
As mentioned before, the focus here is on the unfolding of the
instabilities of a uniform solution $u = u^* =$ const
(notice that any constant
is a solution of \pref{1}).
First consider problems on the whole space, i.e. for
all $\bm{x} \in \bm{R^n}$. If in fact the domain $\Omega$
is large but finite, the unstable
modes we construct may provide a clue to the qualitative features of the
solution far from the boundary.
We subject the given constant solution to a linear stability analysis by
linearizing \pref{1} about it:
\begin{equation}
\alpha w_t = \Delta\left( -\epsilon^2\Delta w + f'(u^*)w \right), \label{12} \end{equation}
and looking for solutions of the form $w(\bm{x},t) =
\exp{(i\bm{k\cdot x} + \sigma t)}$. Such
solutions are possible if and only if
\begin{equation}
\alpha\sigma = |\bm{k}|^2(\beta(u^*) - \epsilon^2|\bm{k}|^2), \label{13}\end{equation}
where $\beta(u) = -f'(u)$.
This tells us (i) that it is possible to have $\sigma > 0$
(implying
linear
instability) only for $u^*$ in the spinodal interval,
where $\beta = \beta(u^*) > 0$; and
(ii) the fastest growing unstable modes are those with
$|\bm{k}|^2 = \frac12 \beta\epsilon^{-2}$.
These modes correspond to spatial patterns with
characteristic
length scale $\epsilon$;
for example in one dimension the wave length will be
$\lambda = \frac{2\pi}{k} = \sqrt{\frac{8}{\beta}}\,\pi\epsilon$.
Since we are assuming $\epsilon$ is small, it is expected that these unstable modes
lead to a state $u$ with fine-grained variation in $\bm{x}$.
Even if we neglect all modes except the fastest growing ones, namely those
with this given value of $|\bm{k}|^2$,
the pattern which emerges is by no means
unique, nor is it expected to be regular. This is because the spatial
orientation of these modes, as well as their phase, can be quite arbitrary.
Cahn presented a description \cite{Cah4} of what can be expected from this
decomposition. In that paper, he calculated the superposition of many such
modes with randomly chosen orientations and phases. The results indicate
rather irregular, sometimes snake-like patterns.
The situation is simpler if our problem has only one space dimension and
the spatial domain is a finite interval, say $0 \le x \le
1$. Then $k$,
which is a
scalar, is restricted to be an integer multiple of $\pi$.
So there will be only a
finite (but large, of the order $\epsilon^{-1}$)
number of unstable modes, namely those
with
\begin{equation}
k = n\pi < \sqrt{\beta}\,\epsilon^{-1}, \label{14}\end{equation}
and there will typically be a unique fastest growing one.
In this case, we expect that the spinodal decomposition process will lead
to a much more unique patterned state, namely one which is more or less
periodic in $x$, with wave number $k$ close to
$\sqrt{\beta/2}\,\epsilon^{-1}$.
This expectation is
reasonable for small enough random initial perturbations of
the given constant solution. It can be shown that for each
$k$ satisfying \pref{14},
there exists an
exact stationary solution $u_k(x)$
of the nonlinear problem \pref{1} with total mass
$\int u_k(x)dx$ equal to that of the given constant
solution $u = u^*$, and which is
periodic in $x$ with period $2\pi/nk$, $n$ an integer
(I conjecture there is a unique
one with $n = 1$).
If $k$ is chosen to correspond to the maximal $\sigma$ in
\pref{13},
this solution $u_k$
of the nonlinear problem is a candidate for the state which is typically
approached through the spinodal decomposition process. Indeed, what appears
to happen is that given any neighborhood of $u_k$, a solution of \pref{1} with
random initial data close enough to $u^*$ (depending on $\epsilon$) will with high
probability enter that neighborhood at some later time.
Grant \cite{G} has recently proved a similar statement,
and in fact was able to estimate such a probability, given a
distribution of initial data concentrated near $u^*$. His
method is based on establishing the existence of a one-
dimensional pseudo-unstable manifold (the locus of ideal
solutions growing ``purely'' with the maximal growth rate)
emanating from $u^*$ in a suitable functional space, and
finding precise estimates for the shapes of orbits in a
neighborhood of $u^*$ in terms of that manifold. He is then
able to estimate which initial data will enter a given
neighborhood of the omega limit set of points on the
manifold. Finally, it is observed that the omega limit set
consists only of points representing spinodally decomposed
states.
On the other hand, the solution will not usually stay in that
neighborhood, because the exact solution $u_k$ is itself unstable. It will
eventually leave it and further evolve under longer time scales. What we
have, then, is a state characterized by a rapid spatial variation of the
concentration function $u$, which attracts solutions which start very near
constant solutions in the spinodal interval. But this attraction is only to a
neighborhood of itself, and is only temporary; solutions will then slowly
drift far away from that state. Physically, the state with rapid variation of
the function $u$ represents a fine mixture of phases,
because large $u$
(near $u_2$)
represents a concentration for which the alloy is in phase
II, whereas
small $u$
(near $u_1$) represents phase I.
The subsequent dynamics will be characterized by coarsening, in which the
characteristic length scales slowly increase, and the time scales as well. A
discussion of this, together with a study of the time scales associated with
instabilities of stationary solutions of varying coarseness, is found in
\cite{BF1}.
One result proved in this latter paper lends credence to
the hypothesis ``coarser = less unstable'', possibly first
enunciated by Langer \cite{L2}. We examined this hypothesis
with two concepts of ``less unstable'': (i) the growth rates
of the unstable modes are smaller, and (ii) the dimension of
the unstable manifold is smaller. In the case of one space
dimension, it was shown by comparison with the spectrum of
the bistable nonlinear diffusion equation that coarser
stationary exact solutions have unstable manifolds with
smaller dimension. It was also shown that the growth rates
of any such nonconstant solution are individually smaller
than the corresponding growth rates for the unstable
constant solution, making the latter the ``most unstable''
of all. Finally, comparing two solutions, if the coarser
one has wave number which divides that of the finer one and
periodic boundary conditions are imposed, then the same
result about the individual growth rates holds. If boundary
conditions \pref{2} are used instead, then the ratio of the
two wave numbers must be odd.
The method of obtaining these results is based on a
characterization of the spectrum of linearized CH operators
in terms of a max-min principle for a certain Rayleigh
quotient (this quotient was used already by Langer
\cite{L2}). This characterization is independent of
dimension. It also allows comparison of such spectra with
that of the second order bistable diffusion operator.
The later evolution of solutions of \pref{1} leads to
configurations characterized by thin layers with width
$O(\epsilon)$ separating regions in which the function $u$ is
relatively flat. The layers move very slowly; and their
laws of motion form a very interesting chapter in the theory
of \pref{1}, which we shall not explore here. Basic work on
these types of solutions has been done in \cite{ABF} and
\cite{P}.
Solutions with layers (representing interfaces) can be thought of as
carrying the free energy in two different places: in the bulk material away
from the layers, and in the interfaces themselves. For the phase field model
to be discussed below, the same is true of the entropy. Theories of phase
transitions have been constructed which from the start envisage such a
partition of thermodynamic properties into bulk and interface portions; the
phase boundaries are then considered to be infinitely thin. See, for example,
\cite{Gu1,Gu2,Gu3}.
The analogous phenomenon of slow motion of interfaces in
the context of the bistable nonlinear diffusion equation,
which is \pref{1} with the prefactor $\Delta$ replaced by $-1$, has
been the subject of much more research, which again will not
be referenced here. It should be noted at this point that
this latter equation is a special case (when $\lambda = 0$)
of the phase field model to be discussed below.
\section{Nucleation} \label{s4}
It is well known that for binary alloy systems describable by the
Cahn-Hilliard equation, there are two mechanisms by which ``phase
separation'' may occur following the quenching of a mixture to a
temperature at which that mixture is no longer globally stable.
The first mechanism is that of spinodal decomposition,
described above in Section 3. This occurs when, at the new
temperature, the average value $u^*$ of the composition
function $u$ is a number in the spinodal interval (and $\epsilon$ is
small enough).
The second mechanism may occur when $u^*$ lies in one of the two metastable
intervals associated with the function $f$.
In that case, it turns out that the solution $u \equiv u^*$
is a locally stable solution of the CH equation, so that in fact the
free energy (under the
constraint that the average composition be equal to $u^*$)
has a local minimum at that uniform state, but the global minimum under that
constraint is lower. In this case, computer simulations and intuitive
arguments indicate that certain perturbations of the constant solution may
result in the solution being drawn toward the globally stable solution, which
will be one in which the vessel is more or less divided into two connected
regions, corresponding to phase I and phase II.
A typical way in which this latter movement may occur is through
nucleation. The perturbation will then be local but large enough. This local
deviation may initiate a local phase change which grows in extent. During
this process, the two phases will be characterized by $u$
taking values only in
the two metastable intervals. When $u \in$ met.\ int.\ 2, we have
phase II; in the other interval we have phase I.
For definiteness we suppose the average value $u^*$ lies in met. int. 2.
Then the perturbation would cause a small region of phase I to appear,
surrounded by phase II. As the former region grows, the value of $u$ in the
latter region must increase, in order for the average value of $u$ to remain
equal to $u^*$. Eventually, the globally stable state is
reached,
for which the
value of $u$ in the two phases is approximately
(for small $\epsilon$) $u_1$ and $u_2$, and the
amounts of the two phases are adjusted so that the average
value of $u$
is as it
should be. Between the two regions, there is a thin transition layer.
The problem of finding and characterizing this ultimate global minimizer
of ${\cal F}$ has been studied extensively for small $\epsilon$ by
Modica and
others \cite{M1,M2,LM,FT,KS}. In one space dimension it can be characterized
completely for any $\epsilon$ \cite{CGS} (see also \cite{GM}).
Clearly, given a local perturbation of the uniform solution, it may or may
not result in such a global transition as described. In other words, it may
or may not lie in the domain of attraction of the global minimizer. If we had
a one-parameter family of local perturbations, there would typically be a
threshold value of the parameter across which the perturbation would pass into
that domain from the domain of attraction of the original locally stable
uniform (constant) solution.
This suggests that there will be a steady solution of the CH equation
which lies on the boundary between these two domains of attraction.
Consider the case of one space dimension, the domain being the unit
interval $(-1,1)$. And for simplicity let us impose a condition of symmetry on
all our solutions, namely that they be even functions of $x$.
It is known \cite{BF2} that for small $\epsilon$ there is an
unstable exact stationary solution with the prescribed mass,
which is in the form of a localized perturbation of the
constant solution. We shall call it a ``spiked solution''
$w(x)$. It is a solution of
\[
\epsilon^2 w'' + f(w) + c = 0,\]
together with Neumann boundary conditions at $x = 0,1$, with the constant $c$
chosen so the required mass condition $\int w(x)dx = u^*$
is satisfied. The spike is
thin, with width of the order $\epsilon$, so it contributes little to the mass.
Outside the spike, therefore, $w$ must be approximately equal to $u^*$, and
therefore $c \approx -f(u^*)$.
This spiked solution no doubt serves as a threshhold of
some sort, lying on the boundary between the domain of
attraction of the constant solution $u \equiv u^*$ and that of
the global minimizer. If so, it may represent a typical
minimum nucleus required to induce motion toward the
globally stable solution. This motion will apparently
consist of the nucleus growing in extent, while at the same
time the value of $u$ surrounding the nucleus will increase
towards $u_2$.
Analogous assertions are very likely true for radially symmetric regions
in higher space dimensions. Also when the domain is the entire space, the
spiked solution probably differs very little from that for a finite domain,
nor will the subsequent evolution differ greatly. For
further discussion, see \cite{BF2}.
\section{Phase-field models} \label{s5}
Let me still speak about a multiphase material, but this time a pure
material rather than an alloy. The phases will be liquid and solid. The
scenario is no longer such that the temperature is constant, but rather now we
suppose the system to be energy-controlled. In the simplest case, it will be
a closed system with no energy flux entering from the outside. I shall also
speak later of the situation when internal heating provides addition of energy
to the material at a prescribed rate, and fixed temperature boundary
conditions extract energy through the boundary.
As before, we still envisage a parameter, called an
order parameter or rather ``disorder parameter'', which will
vary continuously but somehow describe the phase of the
material. In the Cahn-Hilliard model, that parameter is the
concentration $u$. Order parameters for liquid-solid
mixtures have been defined on the basis of statistical
mechanical models, as in \cite{RY,SN} (also see \cite{CA}).
Such a parameter should be a measure of the local
microscopic order or disorder of the material.
In addition to the order parameter, we allow the temperature to vary
throughout the material; but these two will be the only thermodynamic
variables describing the local state of the system.
The first use of a continuously varying order parameter
concept to model a solidification process with diffuse
interface was apparently that of Cahn \cite{Cah1}, who
proposed a much simplified version of the model described
below. In that paper, Cahn also discussed conditions under
which the motion of phase interfaces may be expected to
proceed in accordance with the resulting diffuse interface
theory, as opposed to the classical mechanism of stepwise
lateral motion.
Halpern, Hohenberg and Ma \cite{HHM,HH} described
several classes of theories involving continuously varying
order parameters. One of their models, ``Model C'', was
simplified by Langer and conveyed to Fix. Fix \cite{F},
Caginalp \cite{C1}, and Langer \cite{L3} were the first to
write about this model. In particular, Caginalp popularized
it among mathematicians in this and later papers. The model
was later independently proposed by Collins and Levine
\cite{CL}. Since these papers, it has been studied
extensively; see for example
\cite{C2,CF1,CF2,CF3,FG1,PF,EFFGM,EZ2}. Most of these
studies have concentrated on phase interfacial dynamics and
statics arising from the model, which I shall not go into.
Many separate studies have been made of the stationary solutions of this
model; I shall not survey them here. The stationary solutions satisfy the
same type of equation as do the stationary solutions of \pref{1}: a second order
nonlinear elliptic equation involving an additive constant parameter which is
to be adjusted so that the solution will satisfy a given integral condition.
\section{Derivation} \label{s6}
The derivation of the phase field model in most of the
above papers (\cite{HHM} excluded) was on the basis of a
gradient-corrected free energy functional ${\cal F}$ as in \pref{5},
but with ${\cal F}$ depending on two functions: the order parameter
$\phi$ and the temperature $T$. The function $\phi$ in those
derivations is analogous to the concentration $u$ in the
Cahn-Hilliard setting. We therefore have ${\cal F}[\phi,T]$. The
dynamics of $\phi$ was postulated to be such that ${\cal F}$ will
decrease on solution paths when $T$ is fixed. Specifically,
$\frac{\partial\phi}{\partial t}$ should be along the direction of the negative
gradient of ${\cal F}$ with respect to its first variable $\phi$. A
second equation coupling $\phi$ and $T$ is an equation
expressing conservation of energy when energy flux through
the material is Fickian, i.e. proportional to the
temperature gradient.
This second equation destroys the desired property that
${\cal F}$ decrease on solution paths, and this should really be
no surprise, since the temperature field in ${\cal F}$ is not
fixed, but rather coupled to $\phi$. In fact the gradient
with respect to $\phi$ alone is not significant. Easy
examples verify that indeed ${\cal F}$ does not always decrease on
solution paths.
Besides this, a second objection to this approach in this setting is that
the free energy is not an appropriate thermodynamic quantity to use in
formulating a dynamical theory based on monotonicity, because temperature is
not constant in time.
In view of this, Penrose and Fife \cite{PF} and others developed,
in a thermodynamically consistent way, a class of models of
phase field type in which the dynamics is based on the
requirement of monotonic increase of entropy.
In outline, the procedure was to develop, in accordance with the
constraints of thermodynamics, a reasonable class of expressions for the bulk
entropy density $s(e,\phi)$, where $e$
is the energy density, in systems for which
the total energy is conserved. Then for the same reasons that a gradient term
was included in the Cahn-Hilliard free energy functional \pref{5}, such a term
involving the gradient of $\phi$ is added to $s$
and the result integrated to form a
total entropy functional
\begin{equation}
S[e,\phi] = \int_\Omega - \left[ \frac12 \epsilon^2 |\nabla\phi|^2 + s(e,\phi) \right]
dx. \label{16} \end{equation}
The dynamics of $e$ and $\phi$ are to be chosen so that on solution paths, $S$
increases while the total energy
\[ E = \int e\,dx \]
remains constant.
Let $\Phi$ denote the pair of functions $\Phi = (e,\phi)$.
As before, the proper
dynamics consists in setting
\begin{equation}
\frac{\partial\Phi}{\partial t} = K\mbox{grad}_0S[\Phi], \label{17} \end{equation}
the gradient here denoting a constrained gradient
in a suitable Hilbert space,
the constraint being
\begin{equation}
E = \mbox{ const.} \label{18} \end{equation}
We have to select the proper Hilbert Space (call it $H$)
for pairs $\Phi$,
and
determine this constrained gradient.
I argue in the same way as in Section \ref{s2}. For
any function $h(x)$ (representing a perturbation to the
energy density) with
\begin{equation}
\int h\,dx = 0, \label{19} \end{equation}
and any function $\psi(x)$, this constrained gradient
satisfies
\begin{equation}
\frac{d}{dt}\left.S[e + th,\phi + t\psi]\right|_{t=0} =
\frac{d}{dt}\left.S[\Phi + t\Psi] \right|_{t=0} =
\langle\mbox{grad}_0S[\Phi],\Psi \rangle_H, \label{20} \end{equation}
where I have used the notation $\Psi = (h,\psi)$. Now \pref{19}
is equivalent to the
condition that $h = \Delta W$ for some unique function W
with zero normal derivative
on the boundary and average zero. Make this replacement and
calculate,
from
\pref{16} and \pref{20},
\[
\langle \mbox{grad}_0S,\Psi \rangle_H = \int_\Omega \left[ s_e\Delta W + (\epsilon^2
\Delta\phi + s_\phi)\psi \right] dx, \]
provided that $\phi$ satisfies a zero Neumann condition on
$M\Omega$.
We henceforth
assume that.
Integrating by parts, we see that the right side becomes
\[
-\langle \nabla(s_e),\nabla W \rangle_{{\cal L}_2} + \langle \epsilon^2 \Delta\phi +
s_\phi,\psi \rangle_{{\cal L}_2}. \]
As before, the first term here can be expressed as
\[
\langle -\Delta s_e,h\rangle_{\bar{H}^{-1}},\]
provided that $s_e$ has zero normal derivative on $\partial\Omega$:
\[
\frac{\partial}{\partial\nu}s_e = 0 \mbox{ on } \partial\Omega. \]
Therefore if we choose as Hilbert space $H =
\bar{H}^{-1} \times {\cal L}_2$ with scalar
product equal to the sum of the scalar products in the two component spaces,
we obtain
\begin{equation}
\langle \mbox{grad}_0S,\Psi\rangle_H = \left\langle \left( -\Delta s_e,\epsilon^2 \Delta
\phi + s_\phi \right),\Psi \right\rangle_H. \label{21} \end{equation}
On the other hand, we could equally as well take the scalar product to be a
weighted sum, say with weights $\frac{1}{\mu}$ and 1 ($\mu >
0$).
In that case, instead of \pref{21} we get
\begin{eqnarray}
\langle \mbox{grad}_0S,\Psi\rangle_H &=& \frac{1}{\mu}\langle
-\mu\Delta s_e,h \rangle_{\bar{H}^{-1}} + \langle \epsilon^2 \Delta\phi
+ s_\phi,\psi \rangle _{{\cal L}_2} \nonumber \\
&=& \langle (-\mu\Delta s_e,\epsilon^2\Delta\phi + s_\phi),\Psi \rangle_H.
\label{21'}
\end{eqnarray}
In this space $H$, we therefore conclude that grad$_0 S$ is the pair of
differential expressions in parentheses on the right of
\pref{21'}.
From \pref{17} we
therefore obtain the pair of equations
\begin{equation}\begin{array}{c}
\frac{\partial e}{\partial t} = -K_1\Delta s_e,\\
\frac{\partial\phi}{\partial t} = K_2 (\epsilon^2\Delta\phi + s_\phi), \end{array} \label{22} \end{equation}
where $K_1 = \mu K$ and $K_2 = K$. The boundary conditions
posited,
in order for this
approach to succeed, are
\begin{equation}
\frac{\partial\phi}{\partial\nu} = \frac{\partial}{\partial\nu}s_e = 0 \mbox{ on }
\partial\Omega. \label{22c} \end{equation}
Given these formulas, the next object will be to indicate reasonable
choices for the function $s(e,\phi)$. Moreover, since we shall want to obtain
equations for the dynamics of $T$ and $\phi$ rather than of
$e$ and $\phi$,
we need to
specify the function $e(T,\phi)$.
As a first step, we postulate a form for the energy density function. We
propose
\begin{equation}
e(\phi,T) = v(T)q(\phi) + w(\phi), \label{23} \end{equation}
for some functions $q,~v$, and $w$ to be chosen. Here the first term on the right
typically represents the kinetic energy density of the material, and the
second is the potential energy. It is reasonable that the potential energy
density should depend on the local degree of order in the material, and the
kinetic energy on the temperature. For this reason also, $v$ should be an
increasing function of $T$ (the absolute temperature).
Given \pref{23}, it now follows from
thermodynamic constraints that the entropy
density must be of the form
\[
s(e,\phi) = q(\phi)y(T(e,\phi)) + s_0(\phi), \]\[
y(T) = \int dv(T)/T, \]
where $s_0$ is arbitrary at this stage.
As an example, consider the choice $v = -\frac{\gamma}{T},~q =
1,~ w = \lambda\phi$. Then the
parameter $\lambda$ will denote a latent heat: the potential energy will be
proportional to the disorder parameter $\phi$
and the proportionality constant is
$\lambda$. (The liquid phase has the greatest disorder,
hence the greatest $\phi$, and
also the greatest potential energy.)
With this choice, we obtain $s(e,\phi) = -\frac{\lambda}{2T^2} +
s_0(\phi) = -\frac{1}{2\lambda}(\lambda\phi - e)^2 + s_0(\phi)$, so that
\[
s_\phi = -\frac{\lambda}{\gamma}(\lambda\phi - e) + s_0'(\phi) = \frac{\lambda v}{\gamma} +
s_0'(\phi);~s_e = - \frac{v}{\gamma},\]
and the equations \pref{22} become
\begin{equation}
\begin{array}{c}
\frac{\partial v}{\partial t} + \lambda\frac{\partial\phi}{\partial t} = \frac{K_1}{\gamma}\Delta v,\\
\frac{\partial\phi}{\partial t} = K_2 \left(\epsilon^2\Delta\phi + \frac{\lambda v}{\gamma} + s_0'(\phi)
\right),\\
\partial_\nu \phi = \partial_\nu v = 0 \mbox{ on } \partial\Omega. \end{array} \label{24} \end{equation}
An approximation to \pref{24}
may be made when we assume $T$ does not deviate too
far from the melting temperature $T_0$. Then we replace
$-\frac{1}{T}$
by $-\frac{1}{T_0} + \frac{T - T_0}{T_0^2}$ and set $c = \gamma
T_0^{-2},~K_3 = c\gamma^{-1}K_1,~s_1(\phi) = s_0(\phi) -
\frac{\lambda\phi}{T_0}$, to obtain
\begin{equation}
\begin{array}{c} c\frac{\partial T}{\partial t} + \lambda\frac{\partial \phi}{\partial t} = K_3\Delta T,\\
\frac{\partial\phi}{\partial t} = K_2 \left( \epsilon^2\Delta\phi + \frac{\lambda c}{\gamma} (T - T_0)
+ s_1'(\phi) \right),\\
\partial_\nu\phi = \partial_\nu T = 0 \mbox{ on } \partial\Omega. \end{array} \label{25} \end{equation}
This is the form of the phase field equations which is
found most often in the mathematical literature, if the
function $-s_1'(\phi)$ is chosen to be like the function
$f(\phi)$ considered before (but now it is a function of $\phi$
rather than $u$). The parameter $c$ is meaningful
as a specific heat.
At this point, one could wish for a more enlightened
approach to selecting the proper functions $s_0$ and
$e(T,\phi)$. In \cite{PF}, a couple of models from
statistical mechanics were considered, in which the disorder
parameter can be defined and the dependence of $e$ and $s$ on
$\phi$ specified in a less ad hoc manner. On the basis of
these examples, it appears that a sound physical argument
can be made for making the following alterations in the
traditional phase field model:
1. The potential energy $w(\phi)$ should not be linear,
but rather a quadratic
function. Thus the latent heat $\lambda(\phi) = w'(\phi)$ should be linear rather than
constant: $\lambda(\phi) = a + b\phi$.
2. The term $s_\phi$ should be equal to $-\frac{\lambda(\phi)}{T} +
s_0'(\phi)$.
Of course, again for $T$
restricted to a neighborhood of $T_0$ this dependence on
$T$
can be linearized,
as in the second of \pref{25};
however the coefficient of $T$ will be a function of $\phi$, which
cannot reasonably be made linear.
What results as a more reasonable model is a pair of
equations like
\pref{25},
except that in both equations $\lambda$ is set equal to a linear
function of $\phi$.
\section{Phase separation for the traditional model} \label{s7}
By traditional model, I mean the equations \pref{25}. For convenience I
rewrite them by setting $u = \frac{\lambda c}{\gamma}(T - T_0)$ and
$f(\phi) = -s_0'(\phi)$ (the bistable function).
A rescaling can always be effected so that $K_3/c = 1$. If
$K_2 = O(\epsilon^{-2})$,
then the usual Stefan problem for phase interface motion can be
obtained as a formal limit as $\epsilon \to 0$ \cite{CF3}.
Make this assumption, and in fact
for simplicity set $K_2 = \epsilon^{-2}$. Finally, set $\Lambda =
\frac{\lambda^2}{\gamma}$.
Doing all this, we
obtain
\begin{equation}\begin{array}{c}
\frac{\partial u}{\partial t} + \Lambda\frac{\partial\phi}{\partial t} = \Delta u, \\
\epsilon^2\frac{\partial\phi}{\partial t} = \epsilon^2 \Delta\phi - f(\phi) + u,\\
\partial_\nu\phi = \partial_\nu u = 0 \mbox{ on }\partial\Omega. \end{array} \label{26}\end{equation}
I shall first approach the question of phase separation from the point of
view of the unfolding of instabilities of a constant solution. To do this one
should have to find what those constant solutions are. From \pref{26} it is
evident that they are pairs $(u_0,\phi_0)$ satisfying
\begin{equation}
u_0 = f(\phi_0). \label{27} \end{equation}
We linearize \pref{26} about this constant solution to obtain a dispersion
relation. This amounts to looking for solutions of \pref{26}
in the form
\[
\phi \approx \phi_0 + \delta\Phi\exp{(\sigma t + i\bm{k}\cdot \bm{x})},~~~u \approx u_0 + \delta
U\exp{(\sigma t + i\bm{k} \cdot \bm{x})}. \]
for small $\delta$, and discarding terms in the resulting
equation of order
$\delta^2$ and
higher. Doing this, we find that the various parameters in this assumed form
for the solution must satisfy (setting $k = |\bm{k}|$) \[
\sigma(U + \Lambda\Phi) = -k^2U,\]\[
\epsilon^2\sigma\Phi = -k^2\epsilon^2\Phi + \beta(\phi_0)\Phi + U,\]
where (as before) $\beta(\phi) \equiv -f'(\phi)$. Let $ \bar{\sigma} =
\epsilon^2\sigma,~\bar{k} = \epsilon k$. Then
\[
U(\bar{\sigma} + \bar{k}^2) + \Lambda\bar{\sigma}\Phi = 0,\]\[
-U + (\bar{\sigma} + \bar{k}^2 - \beta) \Phi = 0. \]
Now drop the bars. The larger of the two values of $\sigma$ allowing a nontrivial
$(U,\Phi)$ is given by
\begin{equation}
\sigma = \frac12 \left[ -2k^2 + \beta - \Lambda + \sqrt{(\beta - \Lambda)^2 + 4\Lambda k^2}
\right]. \label{28}\end{equation}
Let $\eta = \Lambda - \beta,~~p = 2k^2$. Then
\[
\sigma = \frac12 \left[ -\eta - p + \sqrt{\eta^2 + 2\Lambda p} \right].\]
If $\eta > 0~(\Lambda > \beta)$, then $\sigma = 0$ for $k = 0$ and for
$k^2 = \beta$. If $\eta < 0$, then $\sigma = 0$ for $k^2 = \beta$
only. If $\beta > 0,~\sigma > 0$ exactly for $k^2$ between 0 and
$\beta$; and if $\beta \le 0,~ \sigma$ is never positive.
This means, as in the Cahn-Hilliard case, that
$(u_0,\phi_0)$
is linearly
unstable for exactly those values of $\phi_0$ where $\beta(\phi_0)
= -f'(\phi_0) > 0$, i.e.\ for $\phi_0$ in
the ``spinodal'' of the function $f(\phi)$.
For $\beta > 0$, the maximum value of $\sigma$ occurs for
$2\frac{\partial\sigma}{\partial p} = 0$, which is for
\[
p = 2k^2 = \beta \left( 1 - \frac{\beta}{2\Lambda} \right) \equiv 2k_m^2 \equiv p_m.\]
For $\Lambda > \beta$ the dispersion relation \pref{28} clearly has
the same qualitative features as the dispersion relation \pref{13} for the
Cahn-Hilliard equation.
Following the same sort of reasoning as in Section 3, we can
expect the phase separation process in higher dimensions to lead to a
nonunique configuration with order parameter $\phi$ varying in $\bm{x}$ according to a
characteristic length of the order $\epsilon$.
In one space dimension, one might expect most small perturbations of the
constant solution to eventually enter, but not stay, in a neighborhood of an
exact stationary solution with wavelength corresponding to
$k = k_m$. Let us
examine the possible stationary solutions. Clearly, $u =$
const.\ and $\phi$ satisfies
\begin{equation}
\epsilon^2\Delta\phi - f(\phi) = \mbox{ const } = -u. \label{29}\end{equation}
Note that stationary solutions $u(x)$
of the Cahn-Hilliard equation satisfy this
same equation, but generally with a different constant. We can think of the
constants in the two cases as functions of the total mass (for the CH case)
and of the total energy (in the phase field case), if these quantities are
known.
The stability properties of these nonconstant stationary
solutions are relevant for understanding the dynamics of
phase separation and coarsening. They were studied in
\cite{BF1} for both the Cahn-Hilliard and traditional
phase-field models. In particular, comparisons were
made between the unstable eigenvalues of the stationary
solutions in these two contexts, and those in the context of
the bistable model, which is the phase-field model with $\Lambda
= 0$ and $u =$ const.
\section{ Internal heating effects} \label{s8}
In the previous section, I indicated that unstable
constant solutions of the phase field equations develop into
a finely grained decomposed state in a manner quite
analogous to spinodal decomposition for the Cahn-Hilliard
equation. In the present section, I shall argue that the
same type of decomposition occurs when a stable initial
state in solid phase is subjected to internal heating.
Although this is suggestive of the well-known appearance
of mush during the heating of an ordinary solid substance,
in actuality the phase field model is inadequate to describe
this latter event in typical cases, and is probably also
inadequate to describe it in the opposite scenario of a
liquid being cooled internally. The phase field model
simply applies to an imperfection-free medium whose only
field variables are temperature and an order parameter. In
actual cases, structural imperfections and temperature
fluctuations nearly always play a governing role in phase
transitions due to internal heating, and they are not
accounted for in this model. Therefore the following is an
idealization whose practical value will be only realized
when other effects are added.
To put in internal heating, we adjust the previous model
\pref{26} so that the
total energy of the system is no longer assumed constant; rather by some
external means, it is increased at a constant rate. This additional energy is
assumed to be supplied to the system at a rate which is not only uniform in
time, but in space as well. Let $q$ denote this constant volumetric heating
rate. Then in place of \pref{26} we have
\begin{equation}
\begin{array}{c}
\frac{\partial u}{\partial t} + \Lambda\frac{\partial\phi}{\partial t} = \Delta u + q,\\
\epsilon^2\frac{\partial\phi}{\partial t} = \epsilon^2\Delta\phi - f(\phi) + u. \end{array} \label{30}\end{equation}
Numerical experiments \cite{FG2} have been performed for
this model, in which initially the system is entirely solid
at a temperature below $T_0$ (the melting temperature of the
material). Thus at $t = 0,~u \equiv u_0 < 0$
and $\phi \equiv \phi_0$, where $f(\phi_0) = u_0$. Recall that
$f(\phi)$ is a bistable cubic-like function.
In our experiment, the temperature is fixed at the boundary, so that the
boundary conditions, in place of the third equation in
\pref{26},
are taken to be
\begin{equation}
\partial_\nu\phi = 0,~u = u_0 \mbox{ on } \partial\Omega. \label{30'}\end{equation}
Several distinct stages in the subsequent evolution of the system are
apparent. These stages are shown in the graphs reproduced in
\cite{FG2}. Here we shall be content with describing their
qualitative features, as follows.
1. In the first stage the system, at each point in space and time, stays
very near the curve $\Gamma:~u = f(\phi)$. The initial point is
on $\Gamma$,
to the left of
the place where $u = 0$. Near the boundaries of the domain interval, in fact,
the solution remains on $\Gamma$ in all later stages as well.
Significant deviations
occur in the interior.
2. Let $\Gamma_1$ denote the ascending (left hand) portion
of $\Gamma$ for $u > 0$. In analogy with the CH model,
call it the left metastable portion of $\Gamma$
corresponding to a superheated solid. It ends at the local
maximum. At some time in the experiment, there will be a
portion of the material whose state passes beyond this
apex. When this happens, the material fairly quickly
develops an instability, provided the rate of heating $q$ is
small enough. This instability results in two effects: (a)
the function $\phi$ develops large spatial oscillations, and
(b) the temperature function $u$ is attracted to 0.
Actually, the immediate behavior of $u$ when the oscillations
first appear is that it also incurs spatial oscillations;
but these are quickly smoothed out so that $u$ is
approximately constant, and that constant, as was
indicated, quickly approaches zero.
3. These spatial oscillations persist for a time longer
than it took them to develop in the first place. Then they
begin disappearing in the following fashion: one by one, the
``dips'' in the oscillatory profile are filled in; this
results in portions of the profile of $\phi$ which are
relatively flat-topped. There is no particular order in
which they fill in; this part is random, except that the
dips which happen to be the least deep are filled in
first. The resulting configuration is a partially coarsened
mush.
4. Eventually, all the dips are filled in; at this
point, one may say that the mush has disappeared. Then the
$\phi$ profile is relatively flat. Shortly thereafter, the $u$
profile flattens as well. The next stage in the evolution
is characterized by the state in the interior of the
material moving up the right hand ascending branch of $\Gamma$.
After a long enough period of time, the system arrives near
its ultimate steady profile, in which the temperature
distribution is parabolic, and the $\phi$ distribution
exhibits two regions: solid and liquid. In each of these,
the system's state is on an ascending stable portion of $\Gamma$.
In the solid, that portion is to the left of $\Gamma_1$, and
in the liquid, it is to the right of $\Gamma_2$. Between the
two, the function $\phi$ makes sharp transitions between $+1$
and $-1$.
\section{Nucleation in the phase field context} \label{s9}
Since we have drawn considerable analogies between the
CH and PF models, it is natural to look for nucleation
effects in the latter context as well. The physical setting
is indeed a natural one: a supercooled liquid represented by
a constant solution $(u_0,\phi_0)$ of \pref{26}. The point
$(u_0,\phi_0)$ will lie on the right hand ascending metastable
branch $\Gamma_2$ of \pref{26}. In this case, an analysis will
undoubtedly show that this uniform solution is locally
stable, but does not maximize the global entropy. We
conjecture that nucleation phenomena will occur, very much
as described before in Section \ref{s3}.
In other words, for the phase field model there are undoubtedly also two
mechanisms for phase separation of the uniform solution
$(u_0,\phi_0)$:
one is
analogous to spinodal decomposition, for which this point lies on the spinodal
portion $\Gamma_s$, and the other is a nucleation and growth process, for which it
lies on one of the two metastable portions $\Gamma_1$ or $\Gamma_2$.
\section{The Novick-Cohen-Pego alternative} \label{s10}
It is possible to model the relaxation of the functional
${\cal F}_0(0)$
with a
local gradient flow, without adding a van der Waals penalty term as in \pref{5}.
One way to do this was devised by Novick-Cohen and Pego. Their model was
motivated in part by consideration of viscous effects in the flow of the basic
material induced by phase separation, and I shall not go into that aspect.
Their model turns out to be a gradient flow for the
functional ${\cal F}_0$
with
reference to a somewhat nonstandard Hilbert space. To define this space, we
first consider a certain selfadjoint operator $A$,
defined on the space $M \equiv$ the
subspace of ${\cal L}_2(\Omega)$ consisting of
functions with average zero: $\int u(x)dx = 0$.
We define, for some $\nu > 0$,
\[
Au = (\nu\Delta - I)^{-1}\Delta u.\]
This can be considered as the selfadjoint extension, in the space $M$, of this
operator acting on smooth functions satisfying a homogeneous Neumann condition
on $\partial\Omega$. For such a function $u$, $(\nu\Delta - I)^{-1} u$ is
defined to satisfy the same
boundary condition. The operator $A$
is readily seen to preserve the condition
of having average zero, and so the selfadjoint extension does so also.
Finally, $A$ is seen to be positive definite.
In $M$, we define a scalar product
\[
\langle u,v\rangle_M \equiv \left\langle A^{-1}u,v
\right\rangle_{{\cal L}_2}.\]
We may now determine grad${\cal F}_0$ in $M$. For this, note
that
\[
\frac{d}{dt}\left.{\cal F}_0[u + tv]\right|_{t=0} = \langle W'(u),v
\rangle_{{\cal L}_2} = \langle AW'(u),v\rangle_M =
\langle\mbox{grad}{\cal F}_0[u],v\rangle_M,\]
so that grad${\cal F}_0[u] = AW'(u)$, and the gradient dynamics gives rise to the
evolution equation\[
\frac{\partial u}{\partial t} = -KAW'(u) = K(I - \nu\Delta)^{-1}\Delta W'(u),\]
or \begin{equation}
\frac{\partial u}{\partial t} + K\Delta \left( f(u) + \nu\frac{\partial u}{\partial t} \right),
\label{31} \end{equation}
where as before, $f = W'$.
A global existence theory for (31) was obtained in
\cite{NP}, as well as some very surprising stability results
which will not be recounted here. Solutions can approach,
as $t \to \infty$, a much greater variety of steady
states that can those of the Cahn-Hilliard equation.
Constant solutions in the spinodal region are unstable,
as for the CH equation, but their dispersion relation does
not identify a preferred wave number with maximal growth
rate. In fact for large $k$, the growth rates of the
unstable modes are approximately independent of $k$, and a
phase separation can take place with arbitrarily small
characteristic length.
\section*{Acknowledgments} Some of the work reported here was supported by the
National Science Foundation and by the Science and Engineering Research
Council.
\section*{Postscript} \label{ps} Since this paper was written
in 1991, a tremendous amount of work has been done relevant
to its content. I indicate here just a few of the
references; they are at the end of the bibliography below.\\[.1in]
\noindent {\bf Relevant to Section \ref{s2}:} In \cite{CT}, the
notion of gradient flow is exploited to provide insight into
other models. In \cite{Gur} and in other papers by Gurtin
and coworkers, an alternative basis for the derivation of
these models is provided and studied.\\[.1in]
\noindent {\bf Relevant to Section \ref{s3}:} There has been
considerable work on the phenomena of spinodal
decomposition, slow interfacial dynamics, and other aspects
of the dynamics of the CH equation. Some of the relevant
papers are \cite{ABrF}-- \cite{SWan1}. The paper \cite{fife}
is another review article dealing partly with these
issues.\\[.1in]
\noindent {\bf Relevant to Section \ref{s4}:} Higher dimensional
nucleation phenomena were dealt with in \cite{Fusco}.\\[.1in]
\noindent {\bf Relevant to Section \ref{s6}:} See
\cite{ChF}--\cite{Um} for a few of the many papers on
thermodynamic phase-field modeling.
\begin{thebibliography}{}
\bibitem{ABF} N. D. Alikakos, P. W. Bates, and G. Fusco, Slow motion for
the Cahn-Hilliard equation in one space dimension,
J. Differential Equations,
{\it 90}, 81-135 (1991).
\bibitem{BF1} P. W. Bates and P. C. Fife, Spectral comparison
principles related to coarsening for the Cahn-Hilliard and
phase-field equations, Physica {\it D 43}, 335-348 (1990).
\bibitem{BF2} P. W. Bates and P. C. Fife, Dynamics of nucleation for the
Cahn-Hilliard equation, SIAM J.\ Appl.\ Math.\ {\it 53,}
990-1008 (1993).
\bibitem{C1} G. Caginalp, An analysis of a phase field model of
a free boundary, Arch.\ Rat.\ Mech.\ Anal.\ {\it 92}, 205-245
(1986).
\bibitem{C2} G. Caginalp, The role of microscopic anisotropy in the macroscopic
behavior of a phase boundary, Ann.\ Phys.\ {\it 172}, 136-155 (1986).
\bibitem{CF1} G. Caginalp and P. C. Fife, Phase field methods for interfacial
boundaries, Phys.\ Rev.\ {\it B, 33}, 7792-7794 (1986).
\bibitem{CF2} G. Caginalp and P. C. Fife, Higher order phase
field models and detailed anisotropy, Phys.\ Rev.\ {\it B, 34},
4940-4943.
\bibitem{CF3} G. Caginalp and P. C. Fife, Dynamics of layered
interfaces arising from phase boundaries, SIAM
J. Appl.\ Math.\ {\it 48}, 506-518 (1988).
\bibitem{Cah1} J. W. Cahn, Theory of crystal growth and interface motion in
crystalline materials, Acta Metallurgica {\it 48}, 554-562 (1960).
\bibitem{Cah2} J. W. Cahn, On spinodal decomposition, Acta
Metallurgica {\it 9}, 795-801
(1961).
\bibitem{Cah3} J. W. Cahn, On spinodal decomposition in cubic crystals, Acta
Metallurgica {\it 10}, 179-183 (1962).
\bibitem{Cah4} J. W. Cahn, Phase separation by spinodal decomposition in isotropic
systems, J. Chem.\ Physics {\it 42}, 93-99 (1965).
\bibitem{Cah5} J. W. Cahn, Spinodal decomposition, Trans.\
Metallurg.\ Soc.\ of AIME {\it 242} , 166-180 (1968).
\bibitem{CHi} J. W. Cahn and J. E. Hilliard, Free energy of a nonuniform system I.
Interfacial free energy, J. Chem.\ Phys.\ {\it 28}, 258-267 (1958).
\bibitem{CGS} J. Carr, M. E. Gurtin, and M. Slemrod, Structured
phase transitions on a finite interval,
Arch.\ Rat.\ Mech.\ Anal.\ {\it 86}, 317-351 (1984).
\bibitem{CL} J. B. Collins and H. Levine, Diffuse interface
model of diffusion-limited crystal growth, Phys.\ Rev.\ B {\it
31}, 6119-6122 (1985). Erratum loc. cit. {\it 33}, 2020.
\bibitem{CN} E. A. Coutsias and J. C. Neu, The aging of nuclei
in a binary mixture, Physica {\it D 12}, 295-302 (1984).
\bibitem{CA} W. A. Curtin and N. W. Ashcroft,
Weighted-density-functional theory of inhomogeneous liquids
and the freezing transition, Phys.\ Rev.\ {\it A 32},
2909-2919 (1985).
\bibitem{EFFGM} J. C. Eilbeck, P. C. Fife, J. E. Furter,
M. Grinfeld, and M. Mimura, Stationary states associated
with phase separation in a pure material. I--The large
latent heat case, Phys.\ Rev.\ Lett.\ {\it A 139}, 42-46
(1988).
\bibitem{EFG} J. C. Eilbeck, J. E. Furter, and M. Grinfeld, On a
stationary state characterization of transition from
spinodal decomposition to nucleation behaviour in the
Cahn-Hilliard model of phase separation, Phys.\ Lett.\ {\it A
135}, 272 (1989).
\bibitem{E} C. M. Elliott, The Cahn-Hilliard model for the
kinetics of phase separation, in Mathematical Models for
Phase Change Problems, J. F. Rodrigues, ed., pp. 35-73,
Int. Ser. of Numerical Math. Vol. 88, Birkh\"{u}ser Verlag,
Basel (1989).
\bibitem{EF} C. M. Elliott and D. A. French, Numerical studies
of the Cahn-Hilliard equation for phase separation, IMA
Jour. Appl. Math {\it 38}, 97-128 (1987).
\bibitem{EZ1} C. M. Elliott and S. Zheng, On the Cahn-Hilliard
equation, Arch.\ Rat.\ Mech.\ Anal.\ {\it 96}, 339-357 (1986).
\bibitem{EZ2} C. M. Elliott and S. Zheng, Global existence and
stability of solutions to the phase field equations, in:
Free Boundary Problems, K-H. Hoffman and J. Sprekels, eds.,
Int. Ser. of Numerical Math. Vol. 95, Birkh\"{u}ser Verlag,
Basel (1990).
\bibitem{FG1} P. C. Fife and G. S. Gill, The phase-field
description of mushy zones, Physica {\it D 35}, 267-275
(1989).
\bibitem{FG2} P. C. Fife and G. S. Gill, Phase-transition
mechanisms for the phase-field model under internal
heating, Phys.\ Rev.\ A {\it 43}, 843-851 (1991).
\bibitem{F} G. Fix, Phase field methods for free boundary
problems, in: Free Boundary Problems, A. Fasano and
M. Primicerio, eds., 580-589, Pitman, London (1983).
\bibitem{FT} I. Fonseca and L. Tartar, The gradient theory of
phase transitions for systems with two potential wells,
Proc.\ Roy.\ Soc.\ Edinburgh {\it A 111}, 89-102 (1989).
\bibitem{G} C. Grant, Spinodal decomposition for the
Cahn-Hilliard equation, Commun. in Partial Differential
Equations {\it 18}, 453--490 (1993); see also: The dynamics of
phase separation for the Cahn-Hilliard equation, PhD Thesis,
Univ. of Utah, 1991.
\bibitem{GD} J. D. Gunton and M. Droz, Introduction to the
Theory of Metastable and Unstable States, Lecture Notes in
Physics No. 183, Springer-Verlag, Berlin (1983).
\bibitem{Gu1} M. E. Gurtin, On a theory of phase transitions
with interfacial energy, Arch.\ Rat.\ Mech.\ Anal.\ {\it 87 },
187-212 (1985).
\bibitem{Gu2} M. E. Gurtin, On the two-phase Stefan problem with
interfacial energy and entropy, Arch.\ Rat.\ Mech.\ Anal.\ {\it
96}, 199-241 (1986).
\bibitem{Gu3} M. E. Gurtin, On phase transitions with bulk,
interfacial, and boundary energy, Arch.\ Rat.\ Mech.\ Anal.\
{\it 96}, 243-264 (1986).
\bibitem{GM} M. E. Gurtin and H. Matano, On the structure of the
equilibrium phase transition within the gradient theory of
fluids, Quart.\ Appl.\ Math.\ {\it 46} (1988), 301-317.
\bibitem{HHM} B. I. Halperin, P. C. Hohenberg, and S-K. Ma,
Renormalization-group methods for critical dynamics:
I. Recursion relations and effects of energy conservation,
Phys.\ Rev.\ {\it B 10}, 139-153 (1974).
\bibitem{HH} P. C. Hohenberg and B. I. Halperin, Theory of
dynamic critical phenomena, Rev.\ Modern Physics {\it 49},
435-485 (1977).
\bibitem{KS} R. V. Kohn and P. Sternberg, Local minimizers and
singular perturbation, Proc.\ Roy.\ Soc.\ Edin.\ {\it 3,}
69-84 (1989).
\bibitem{LL} L. D. Landau and E. M. Lifshitz, Statistical
Physics, 2nd Ed., Pergamon Press, Oxford (1969)
\bibitem{L1} J. S. Langer, Theory of the condensation point,
Annals of Physics {\it 4 }, 108-157 (1967).
\bibitem{L2} J. S. Langer, Theory of spinodal decomposition in
alloys. Ann.\ Physics {\it 65 }, 53-86 (1971).
\bibitem{L3} J. S. Langer, Models of pattern formation in first-order phase
transitions, pp. 164-186 in Directions in Condensed Matter Physics, World
Science Publishers (1986).
\bibitem{LS} Yu. S. Lipatov and V., V. Shilov, Spinodal
decomposition in polymeric systems, Russ.\ Chem.\ Revs.\ {\it
53 }, 697-710 (1984).
\bibitem{LM} S. Luckhaus and L. Modica, The Gibbs-Thompson
relation within the gradient theory of phase transitions,
Arch.\ Rat.\ Mech.\ Anal.\ {\it 107 }, 71-83 (1989).
\bibitem{M1} L. Modica, The gradient theory of phase transitions
and the minimal interface criterion, Arch.\ Rat.\ Mech.\ Anal.\
{\it 98 }, 123-142 (1987).
\bibitem{M2} L. Modica, Gradient theory of phase transitions
with boundary contact energy, Ann. Inst. H. Poincar\'{e},
Anal. Nonlin\'{e}are {\it 5}, 453-486 (1987).
\bibitem{NS} B. Nicolaenko and B. Scheurer, Low-dimensional
behavior of the pattern- forming Cahn-Hilliard equation,
Trends in the Theory and Practice of Nonlinear Analysis,
V. Lakshmikantham, ed., 323-336, Elsevier (1985).
\bibitem{NST1} B. Nicolaenko, B. Scheurer, and R. Temam,
Inertial manifolds for the Cahn-Hilliard model of phase
transition, in: Ordinary and Partial Differential Equations,
Proc. Dundee Conf. (1986), B. Sleeman and R. J. Jarvis,
eds., Pitman Res. Notes in Math. Vol. 157 , Longman,
N.Y. (1987).
\bibitem{NST2} B. Nicolaenko, B. Scheurer, and R. Temam, Some
global dynamical properties of a class of
pattern-forming equations, Comm.\ Par.\ Diff.\ Eqns.\ {\it 14},
245-297 (1989).
\bibitem{N} A. Novick-Cohen, The nonlinear Cahn-Hilliard equation: transition from
spinodal decomposition to nucleation behavior,
J. Stat.\ Physics {\it 38}, 707-723
(1985).
\bibitem{NP} A. Novick-Cohen and R. L. Pego, Stable patterns in
a viscous diffusion equation, Trans.\ Amer.\ Math.\ Soc.\
{\it 324}, 331-351 (1991).
\bibitem{NSe} A. Novick-Cohen and L. A. Segel, Nonlinear aspects
of the Cahn-Hilliard equation, Physica {\it 10D,} 277-298
(1984).
\bibitem{P} R. L. Pego, Front migration in the nonlinear
Cahn-Hilliard equation, Proc.\ Roy.\ Soc.\ London {\it A 422,}
261-278 (1989).
\bibitem{Pen} O. Penrose, Statistical mechanics and the kinetics
of phase separation, in: Material Instabilities in Continuum
Mechanics and Related Mathematical Problems, J. M. Ball,
ed., 373-394, Clarendon Press, Oxford (1988).
\bibitem{PF} O. Penrose and P. C. Fife, Thermodynamically
consistent models of phase- field type for the kinetics
of phase transitions, Physica {\it D 43}, 44-62 (1990).
\bibitem{RY} T. V. Ramakrishnan and M. Yussouff, First-principle
order-parameter theory of freezing, Phys.\ Rev.\ {\it B 19},
2775-2794 (1979).
\bibitem{SN} S. Sachdev and D. R. Nelson, Order in metallic
glasses and icosahedral crystals, Phys.\ Rev.\ {\it B 32},
4592-4606 (1985).
\bibitem{SS} V. P. Skripov and A. V. Skripov, Spinodal
decomposition (phase transitions via unstable states).
Soviet Physics Uspekhi {\it 22}, 389-410 (1979).
\bibitem{W} J. D. van der Waals, The thermodynamic theory of
capillarity under the hypothesis of a continuous variation
of density (translation of Dutch title):
Konink. Akad. Weten. Amsterdam (sec. 1), Vol. 1, No. 8
(1893). English trans.: J. S. Rowlinson, J. Stat.\ Phys.\
{\it 20}, 197-244 (1979).
\bibitem{Z} S. Zheng, Asymptotic behavior of the solution to the
Cahn-Hilliard equation, Appl.\ Anal.\ {\it 23}, 165-184
(1986).
\bibitem{CT} J. W. Cahn and J. E. Taylor, Surface motion by
surface diffusion, Acta metall.\ mater.\ {\it 42}, 1045-1063
(1994).
\bibitem{Gur} M. E. Gurtin, Generalized Ginzburg-Landau and
Cahn-Hilliard equations based on a microforce balance,
Physica {\it D 92}, 178-192 (1996).
\bibitem{ABrF} N. Alikakos, L. Bronsard, and G. Fusco, Slow
motion in the gradient theory of phase transitions via
energy and spectrum, Calculus
of Variations 6, 39-66 (1998).
\bibitem{AF1} N. Alikakos, G. Fusco, The spectrum of the Cahn-Hilliard
operator for generic interface in higher space dimensions, Indiana
Mathematics Journal 42, 637-674 (1993).
\bibitem{AF2} N. Alikakos, G. Fusco, Slow dynamics for the Cahn-Hilliard
equation in higher space dimensions. Part I: Spectral estimates,
Comm.\ in P.D.E. 19, 1397-1447 (1994).
\bibitem{AF3} N. Alikakos, G. Fusco, Slow dynamics for the
Cahn-Hilliard equation in higher space dimensions: the
motion of bubbles, Arch.\ Rat.\ Mech.\ Anal.
141, 1-61 (1998).
\bibitem{AFK} N. Alikakos, G. Fusco, and M. Kowalczyk,
Finite-dimensional dynamics and interfaces intersecting the
boundary: equilibria and quasi-invariant manifolds, Indiana
Mathematics Journal 45, 1119-1155 (1996).
\bibitem{BX1} P. Bates and J. Xun, Metastable patterns for the
Cahn-Hilliard equation: Part I, J. Differential Equations 111, 421-457 (1995).
\bibitem{BX2} P. Bates and J. Xun, Metastable patterns for the
Cahn-Hilliard equation: Part II, J. Differential Equations 117, 165-216 (1995).
\bibitem{ChenX} Xinfu Chen, Global asymptotic limit of solutions
of the Cahn-Hilliard equation, J. Differential Geometry {\it
44}, 262-311 (1996).
\bibitem{ERD} K. Elder, T. Rogers, and R. Desai, Early stages of
spinodal decomposition for the Cahn-Hilliard-Cook model of
phase separation. Phys.\ Rev.~ B 38, 4725--4739 (1988).
\bibitem{E1} D. Eyre, Coarsening dynamics for solutions of the
Cahn-Hilliard equation in one dimension, preprint.
\bibitem{E2} D. Eyre, Systems of Cahn-Hilliard equations, SIAM
J. Appl.\ Math.\ 53, 1686-1712 (1993).
\bibitem{E3} D. Eyre, Cascades of Spinodal Decomposition in the
Ternary Cahn-Hilliard Equations, to appear
\bibitem{fife} P. C. Fife, Pattern formation in gradient
systems, to appear as a chapter in Handbook for Dynamical
Systems. III: Toward Applications, B. Fiedler and N. Kopell,
eds., Elsevier.
\bibitem{MW1} S.\ Maier-Paape and T.\ Wanner, Spinodal
decomposition for the Cahn-Hilliard equation in higher
dimensions. Part I: Probability and wavelength estimate,
Comm. Math. Physics {\it } 195, 435--464 (1998).
\bibitem{MW2} S.\ Maier-Paape and T.\ Wanner, Spinodal
decomposition for the Cahn-Hilliard equation in higher
dimensions: Nonlinear dynamics, Arch.\ Rat.\ Mech.\ Anal.\
{\it 151}, 187-219 (2000).
\bibitem{SWan} E. Sander and T. Wanner, Monte Carlo simulations
for spinodal decomposition, J. Stat.\ Phys. {\it 95},
925-948 (1999).
\bibitem{SWan1} E. Sander and T. Wanner, Unexpectedly linear
behavior for the Cahn-Hilliard equation, SIAM J. Appl.\
Math.\ {\it 60}, 2182-2202 (2000).
\bibitem{Fusco} P. Bates and G. Fusco, Equilibria with many
nuclei for the Cahn-Hilliard equation, J. Differential
Equations {\it 160}, 283--356 (2000).
\bibitem{ChF} C. Charach and P. Fife, On thermodynamically
consistent schemes for phase field equations, Open Systems
and Information Dynamics {\it 5}, 99-123 (1998).
\bibitem{FrG1} E. Fried and M. E. Gurtin, Continuum theory of
thermally induced phase transitions based on an order
parameter. Physica {\it D 68}, 326-343 (1993).
\bibitem{Um} A. P. Umantsev, Thermodynamic stability of phases
and transition kinetics under adiabatic
conditions. J. Chem.\ PHys.\ {\it 96}, 189-200 (1993).
\end{thebibliography}
\noindent{\sc Paul C. Fife}\\
Mathematics Department \\
University of Utah, 155 S. 1400 E.,\\
Salt Lake City, UT 84112-0090, USA \\
e-mail: fife@math.utah.edu
\end{document}