\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}