\documentclass[twoside]{article}
\usepackage{amssymb, graphicx, amsmath}
\pagestyle{myheadings}
\markboth{\hfil A coupled Cahn-Hilliard particle system \hfil EJDE--2002/73}
{EJDE--2002/73\hfil Tony Shardlow \hfil}
\begin{document}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
{\sc Electronic Journal of Differential Equations},
Vol. {\bf 2002}(2002), No. 73, pp. 1--21. \newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu (login: ftp)}
\vspace{\bigskipamount} \\
%
A coupled Cahn-Hilliard particle system
%
\thanks{ {\em Mathematics Subject Classifications:} 35K35, 65M60, 60H15.
\hfil\break\indent
{\em Key words:} Boundary Value Problems, Finite Element Methods, \hfil\break\indent
Stochastic Partial Differential Equations.
\hfil\break\indent
\copyright 2002 Southwest Texas State University. \hfil\break\indent
Submitted April 5, 2002. Published August 13, 2002.} }
\date{}
%
\author{Tony Shardlow}
\maketitle
\begin{abstract}
A Cahn-Hilliard equation is coupled to a system of stochastic
differential equations to model a random growth process.
We show the model is well posed and analyze the model asymptotically
(in the limit of the interfacial distance becoming small),
to recover a free boundary problem. A numerical method together
with example solutions is presented.
\end{abstract}
\def\Mean#1{\mathbf{E}\Big[#1\Big]}
\def\mean#1{\mathbf{E}#1}
\def\norm#1{\|#1\|}
\def\order#1{\mathcal{O}(#1)\,} % order big o
\newcommand{\bvec}[1]{\mathbf{#1}}
\newcommand{\hvec}[1]{\hat{\mathbf{#1}}}
\newcommand{\tstep}{{\Delta t}} % delta t
\newtheorem{theorem}{Theorem}[section]
\numberwithin{equation}{section}
\section{Introduction and background}
Stochastic PDEs have been used to model random growth processes since
the introduction of the Kardar-Parisi-Zhang (KPZ)
equation~\cite{kardar:1986}. This equation restricts the topology of
an aggregate in $\mathbb{R}^d$, so that its interface may be represented
as a graph $\mathbb{R} \rightarrow \mathbb{R}^{d-1}$. The KPZ equation describes
the evolution of this graph as a fourth order stochastic PDE in $d-1$
dimensions.
Models based on the evolution of a graph are very restrictive on the
topology. By introducing an extra dimension to the model, arbitrary
topologies may be described by writing a PDE for a phase variable
$u(t,x)\colon \mathbb{R}^+\times \mathbb{R}^d \rightarrow \mathbb{R}$.
Then, the
growth will be that of an aggregate $\{x\in \mathbb{R}^d \colon
u(t,x)\approx u_+\}$, where $u_+$ depends on the model under
consideration, and the boundary of the aggregate will be a level set
$\{u= 0\}$ (say). Such phase fields models have been used to described
pattern formation. One well known example is the Cahn-Hilliard
equation~\cite{elliott:1989} for a region $\Omega$:
\begin{gather*}
u_t=\Delta \theta \\
-\theta=\epsilon^2 \Delta u-f(u),
\end{gather*}
where homogeneous Neumann boundary conditions are placed on $u$ and
$\theta$. The function $f$ is the gradient of a double well potential.
The parameter $\epsilon\ll 1$ and measures the interfacial thickness.
The total phase $\int_\Omega u \;dx$ is constant in this model and no
growth is included.
It is of interest to consider perturbing these equations by
noise~\cite{keblinski:96,keblinski:94}. For example, the equation
\begin{equation}\label{eq:dch}
u_t = -\epsilon^2 \Delta^2 u + \Delta (u-u^3) + I,
\end{equation}
with appropriate boundary conditions. For $I=0$, this is the standard
Cahn-Hilliard equation. The driving term $I$ models the random
deposition onto the aggregate. The function $I$ is expressed in terms
of $\nabla u$, to encourage growth normal to the surface of the
aggregate and to concentrate the effects of the noise at the
interface. Two forms for $I$ are suggested,
\[
I_1:=
c_1| \nabla u| + c_2 \sqrt{ |\nabla u|} \dot{W}(t)
\quad \text{or} \quad
I_2:=
| \nabla u|^2
\Big( c_3 + c_4 \dot{W}(t)
\Big),
\]
where $c_1,\dots,c_4$ are constants and $\dot{W}(t)$ is a
``derivative'' of a space-time Wiener process. Numerical
simulations~\cite{keblinski:94} for $I_1$ suggest behaviour similar
to the Eden model, which is characterised by uniform growth rates. The
$I_2$ model is seen to depend on curvature effects at the interface.
(For examples of perturbing the Cahn-Hilliard by thermal fluctuations,
see ~\cite{daprato:1996b, cook:1970}).
This model includes random deposition of particles, the effect of
surface diffusion, and makes no topological assumptions on the
aggregate. One further important physical feature of these systems is
shadowing, where certain areas are free from deposition (because the
flow of material is blocked). This effect is not included in
equation~\eqref{eq:dch}, but by coupling a second equation, Keblinski et
al.~\cite{keblinski:94} have modelled shadowing.
Defining solutions of the above equation in a mathematically precise
way is very difficult, because of the multiplicative forcing functions
$I$. Weak existence theory by standard Faedo-Galerkin convergence
arguments depends on estimates of the type $\mean{(I, u(t))}$ (where
$\mean{}$ denotes average over realisations of $W(t)$ and
$(\cdot,\cdot)$ is the $L_2$ inner product). Space time white noise
will never satisfy such a property because it does not give a well
behaved process in the mean square sense. A smooth process $W(t)$
would need to be introduced.
Another approach to introducing noise into a PDE is the use of
particle systems. Directly modelling the evolution of the depositions
before they hit the aggregate is a natural technique for introducing
noise to the system. Indeed, Diffusion Limited Aggregation
(see~\cite{dla2}) is a spatially discrete model that uses this technique
and incorporates shadowing and arbitrary topologies in a natural
way.
We propose a model for the interaction of a set of particles evolving
according to It\^o stochastic differential equations (SDEs) with a
Cahn-Hilliard system. The particles represent material that is
deposited onto an aggregate, represented by the field $u$. It is a
mean field theory, in that we neglect thermal fluctuations in the
aggregate. In the context of depositions, it is the fluctuations in
the trajectories of the particles that is responsible for the complex
morphology of the aggregate.
Consider the following coupled Cahn-Hilliard
particle system on a domain $\Omega=[0,L]^d$.
\begin{equation}\label{eq:chp}
\begin{gathered}
u_t= \Delta \theta + \gamma_1
\sum_{i\in {\cal P}(t)} |\nabla u| \delta_{X_i} \\
-\theta= \epsilon^2 \Delta u - f(u),
\end{gathered}
\end{equation}
where $\epsilon\ll 1$ is the interfacial parameter, $\gamma_1>0$ is
the coupling strength, $\delta_X(x)=1$ if $|x-X|\le R$ and $=0$
otherwise. The function $f$ is the gradient of a double well potential
$F$ with minima $u_\pm$, often $f(u)=(u^3-u)$; throughout we suppose
that $f$ is an odd polynomial with positive leading order coefficient.
The particles have radius $R$ with centres $X_i\in \Omega$. Particle
$i$ is said to be alive at time $t$ if $i\in{\cal P}(t)$. Further, we
assume that $i\in {\cal P}(t)$ if and only if $t\in[\tau^s_i,
\tau^e_i]$. Particles are introduced to the system at positions
$X_i^s$ at times $\tau^s_i$ and are annihilated independently at times
$\tau^e_i$ with rate
\[
\frac{\gamma_2}{R^d}\int_\Omega |\nabla u|(x) \delta_{X_i}(x)\;dx.
\]
The position of the particle $X_i$ satisfies the
It\^o SDE
\begin{equation}
\label{eq:sdes}
dX_i = \bvec{\lambda}(X_i) \;dt + \sigma(X_i) \,d\bvec{B}_i(t),
\end{equation}
where $\bvec{\lambda}\in C^\infty(\mathbb{R}^d, \mathbb{R}^d)$ and
$\sigma\in C^\infty(\mathbb{R}^d, \mathbb{R}^{d\times d})$ and $\bvec{B}_i$ are IID
standard Brownian motions in $\mathbb{R}^d$ and with (for simplicity)
periodic boundary conditions on $\Omega$.
The Brownian motions $\bvec{B}_i(t)$ live on a probability space with
measure $\mathbf{P}$ and are adapted to a filtration ${\cal
F}_t$. Expectations with respect to $\mathbf{P}$ are denoted $\mathbf{E}$.
The outline of this paper is as follows. \S\ref{sec:rig} considers the
existence and uniqueness of \eqref{eq:chp}, first deriving an a priori
bound for the solution and then sketching the steps necessary to prove
weak solutions exist. This is sketched as many of the steps are quite
standard. \S\ref{sec:asym} describes some basic properties of
\eqref{eq:chp} in a non-rigorous manner, by using asymptotic analysis
and looking for an approximating free boundary problem.
\S\ref{sec:fem} describes a finite element scheme for this system and
how it has been implemented using the deal.II software
package~\cite{BHK:2001}. The final section \S\ref{sec:res} gives some
numerical simulations.
\section{Existence of solutions}\label{sec:rig}
We now discuss the existence of solutions for the coupled
Cahn-Hilliard system~\eqref{eq:chp}. We show how to derive a priori
bounds for this equation, following a standard argument, see for
example~\cite{temam:1988}. With the a priori bound, standard
Faedo-Galerkin arguments can be applied to prove existence of
solutions to the equations. We build up the argument for the following
deterministic equation, before adding in the random components:
\begin{gather*}
u_t = \Delta \theta + \gamma_1 \sum_{i\in {\cal P}(t)} |\nabla u|
\delta_{X_i} \\
-\theta= \epsilon^2 \Delta u - f(u),
\end{gather*}
subject to homogeneous Neumann boundary conditions on $u$ and $\theta$
on $\partial \Omega$ and where $X_i\colon \mathbb{R}^+ \rightarrow \Omega$
are continuous functions of time and $\delta_{X_i}$ is the indicator
function on a ball of radius $R$ centred on $X_i$. For the present
$X_i$ may be considered to be deterministic; we do not discuss the
random aspect until Theorem~\ref{thm:two}.
Throughout the present section, we use $K$ to denote a generic
constant. We work with the function space
\[
V:= \Big\{ \phi\in H^2(\Omega)\colon \frac{\partial \phi}{\partial
\bvec{n}}=0 \text{ on $\partial \Omega$}\Big\}
\]
and denote by $|\cdot|$ and $(\cdot,\cdot)$ the standard norm and
inner product on $L_2(\Omega)$.
The weak formulation of the above equation is achieved by multiplying
by $v\in V$ and applying the boundary conditions with Green's formula:
\begin{align*}
(v, u_t) = & (v,\Delta \theta )
+ \gamma_1\sum_{i\in {\cal P}(t)} (v, |\nabla u|\delta_{X_i})
= -(\nabla v, \nabla \theta )
+ \gamma_1\sum_{i\in {\cal P}(t)} (v, |\nabla u|\delta_{X_i}).
\end{align*}
Now,
\begin{equation}\label{eq:weakF}
-(\nabla v,\nabla \theta )
= (\Delta v, (-\epsilon^2 \Delta u) )
-(\nabla v, \nabla f(u)).
\end{equation}
Take $v=u$:
\[
\tfrac12 \frac{d}{dt} |u|^2 + \epsilon^2 (\Delta^2 u, u) + ( \nabla u,
f'(u) \nabla u)
= \gamma_1\sum_{i\in {\cal P}(t)} ( |\nabla u(x)| \delta_{X_i}, u).
\]
Write $f(s)=b_p s^{2p-1}+ b_{p-1} s^{2p-3} + \dots+ b_1 s$. The leading
term of $f'(s)$ is $(2p -1) b_{p}s^{2p-2}$ and we can find $c>0$ such that
\[
f'(s) \ge (2p-1) b_{p} s^{2p-2} -c.
\]
Substituting for $f'(s)$ and applying Cauchy-Schwartz, we have
\[
\tfrac12 \frac{d}{dt} |u|^2 + \epsilon^2 |\Delta u|^2
+ (2p-1) b_{p} \int_\Omega u^{2p-2} |\nabla u|^2 \;dx \le c|\nabla u|^2 + \tfrac12 |\nabla u|^2
+ \tfrac12 (\gamma_1 N)^2 |u|^2,
\]
where $N$ denotes the maximum number of particles in ${\cal P}(t)$.
Standard interpolation estimates give $|\nabla u|^2 \le K |u| \;
\norm{u}_{H^2(\Omega)}$. Lemma 4.2~\cite{temam:1988} gives that
the norm $\norm{\cdot}_{H^2(\Omega)}$ is equivalent to $|\Delta u| +
|u|$. Hence, for some constant $K$
\[
(c+\tfrac12) |\nabla u|^2 \le K |u| ( |\Delta u| + |u|) \le \tfrac12
\epsilon^2 |\Delta u|^2 + (K+ K^2/2\epsilon^2) |u|^2.
\]
Finally, we obtain for a different constant $K$
\[
\frac{d}{dt} |u|^2 + \epsilon^2 |\Delta u|^2
+ 2(2p-1) b_{2p} \int_\Omega u^{2p-2} |\nabla u|^2 dx \le K \; |u|^2.% + 1.
\]
By the Gronwall Lemma, we have uniform bounds on $[0,T]$ in
$L^2(\Omega)$ and by integrating we have
\[
\int_0^T \norm{u(s)}^2_{H^2(\Omega)} \;ds \le K.
\]
Arguing further, we define the Lyapunov function
\[
{\cal V}(u) := \tfrac12 \epsilon^2 |\nabla u|^2 + \int_\Omega F(u(x)) \;dx.
\]
Arguing formally for a moment (because the integral above need not be
defined for $u\in V$), we note that
\[
\frac{d}{dt} {\cal V}(u) = (\theta, u_t)
\]
and
\begin{align*}
(\theta, u_t) =& \epsilon^2 (\Delta \theta, \theta) + \gamma_1 (
\theta, \sum_{i\in {\cal P}(t)} |\nabla u| \delta _{X_i})\\
=& -\epsilon^2 |\nabla \theta|^2 + \gamma_1 ( \theta,
\sum_{i\in {\cal P}(t)} |\nabla u| \delta _{X_i}).
\end{align*}
Then,
\begin{align*}
\frac{d}{dt}{\cal V}(u(t)) \le& - \epsilon^2 |\nabla \theta|^2 + \tfrac12
|\theta|^2 +
\tfrac12 \gamma_1^2 N^2 |\nabla u|^2\\
\le & - \epsilon^2 |\nabla \theta|^2 + \tfrac12
(2 |\epsilon^2 \Delta u|^2 + 2 |f(u)|^2) +
\tfrac12 \gamma_1^2 N^2 |\nabla u|^2\\
\le& - \epsilon^2 |\nabla \theta|^2 + K
\norm{u}_{H^2(\Omega)}^2,
\end{align*}
where the last inequality uses the fact that $W^2(\Omega)$ can be
continuously embedded in $C^0(\Omega)$ for $d=2,3$. Thus, ${\cal
V}(u(t))$ is bounded on $[0,T]$ as
$\int_0^T \!\norm{u(s)}^2_{H^2(\Omega)}ds$ is finite. The argument can be made rigorous
by truncating $f$ and showing convergence in the limit of the
truncation.
This time take $v=\Delta^2 u$ in~\eqref{eq:weakF} to gain
\begin{align*}
\tfrac12 \frac{d}{dt} |\Delta u|^2 + \epsilon^2 | \Delta^2 u|^2 =
& ( \Delta f(u), \Delta^2 u) + (\gamma_1 \sum_{i \in {\cal P}(t)} |\nabla u(x)| \delta_{X_i},
\Delta^2 u) \\
=& \frac{1}{\epsilon^2}| \Delta f(u)|^2 + \tfrac 14 \epsilon^2 |\Delta^2
u|^2 + \frac{\gamma_1^2 N^2 }{\epsilon^2} |\nabla u|^2 + \tfrac
14 \epsilon^2 |\Delta^2 u|^2.
\end{align*}
Thus,
\[
\frac{d}{dt} | \Delta u|^2 + \epsilon^2 |\Delta^2 u |^2 \le
\frac{2}{\epsilon^2} | \Delta f(u)|^2 + \frac{\gamma_1^2 N^2
}{\epsilon^2} |\nabla u|^2.
\]
It is proved that (see~\cite{temam:1988}), by restricting
$p=2$ in three dimensions and taking an arbitrary positive integer $p$ in
dimensions $d=1,2$, that
\[
|\Delta f(u)|^2 \le \tfrac12 \epsilon^2 |\Delta^2 u|^2 +K.
\]
Then, we have a differential inequality
\[
\frac{d}{dt} | \Delta u|^2 + \epsilon^2 |\Delta^2 u|^2 \le K,
\]
after using the boundedness of ${\cal V}$. This gives a priori estimates in\\
$L^\infty(0,\infty; H^2 (\Omega))$ and $L^2(0,T; H^4(\Omega))$.
\begin{theorem}\label{thm:one}
Consider dimension $d=1,2,3$. For every $u_0 \in L_2(\Omega)$, the
initial value problem~\eqref{eq:chp} has a unique weak solution $u$ in
\[
L^\infty ( [0,T]; H) \cap L^2([ 0,T] ; V)).
\]
The mapping $u(t)$ is continuous in $t$. Let $p$ be a positive
integer, with $p=2$ (i.e., $f$ has degree three) when $d=3$, arbitrary
for dimension 1 and 2. Choose initial data $u_0 \in V$. Then
\[
u(t) \in C([0,T], V) \cap L^2([ 0,T]; {\cal D}(A)).
\]
\end{theorem}
\paragraph{Proof}
This is a standard Faedo-Galerkin approximation argument. Let
$\phi_i$ denote the eigenfunctions of $A$. The idea is to seek
solutions $u_m$ of the form
\[
u_m(t) = \sum_{i=1}^m g_{im}(t) \phi_i,
\]
satisfying
\[
( \frac{du_m}{dt}, \phi_j)
+ \epsilon^2 (\Delta u_m, \Delta \phi_j)
- (f (u_m), \phi_j) = (\gamma_1\sum_i |\nabla u_m| \delta_{ X_i },
\phi_j),
\; j=1,\dots,m.
\]
This is an ODE and existence of solutions is elementary. Further, the a
priori bounds developed above holds for $u_m$ and allow us to take
limits of $u_m$. Standard functional analytic arguments give
convergence to a solution $u$ having the properties described above.
\hfill$\Box$
\begin{theorem} \label{thm:two} Let $d=2$ or $3$ and let $f$ be cubic
($p=2$). Suppose that $\bvec{\lambda}\in C^\infty(\mathbb{R}^d, \mathbb{R}^d)$
and $\sigma\in C^\infty(\mathbb{R}^d, \mathbb{R}^{d\times d})$ and both
functions are globally Lipschitz. Consider initial data $u_0 \in
V$, initial times $\tau^s_i\in \mathbb{R}^+$, and initial positions
$X_i^s$, for $i=1,\dots,N$ ($N$ the number of particles). Then,
there exists a solution of~\eqref{eq:chp} consisting of the phase
variable $u(t)\in C([0,T], V) \cap L^2(0,T; {\cal D}(A))$, the
conditional densities $p_i(t,y;u)$ on $[\tau_i^s,T)\times
C([0,T],V)$ (for the probability the particle is at $y$ at time $t$ given a
phase trajectory $u$), and the particle trajectories $X_i(t)\colon
[\tau^s_i,\tau^e_i]\rightarrow \mathbb{R}^d$. The solution $(u(t),
p_i(t, y; u), X_i(t))$ is uniquely defined by the following properties:
(i) the phase variable obeys $u(0)=u_0$ and for each $v\in V$,
\begin{align*}
(v,u_t)=&(\nabla v, \nabla \theta) + \gamma_1
\sum_{i\in {\cal P}(t)} (v,|\nabla u|\delta_{X_i(t)}),\\
-(v,\theta) = &\epsilon^2 (\nabla v,\nabla u)- (v,f(u));
\end{align*}
(ii) the particle trajectories $X_i(t)$ are continuous functions
$[\tau_i^s, \tau_i^e]\rightarrow \mathbb{R}^d$ satisfying for $\tau_i^s \le t
\le \tau_i^e$
\[
X_i(t)-X_i(\tau_i^s)= \int_{\tau_i^s}^t \lambda(X_i(s)) \;ds
+\int_{\tau_i^s}^t \sigma(X_i(s)) \,d\bvec{B}_i(s),
\]
for independent Brownian motions $\bvec{B}_i$. The annihilation time
$\tau_i^e=\tau^e_i(u)$, where $\tau_i^e(u)$ are independent random
variables on $[\tau_i^s,T)$ satisfying
\[
\mathbf{P}{ \tau_i^e(u)>t }=\int_\Omega p_i(t,y;u) \; dy,\qquad u\in C([0,T],V).
\]
(iii) $p_i(t,y;u)$ is a delta function at $X_i^s$ at $t=\tau_i^s$ and
for $t>\tau_i^e$ satisfies
\[
\frac{dp_i(t, y; u)}{dt} = {\cal L} p_i(t, y; u) -\frac{\gamma_2}{R^d}
\int_\Omega |\nabla u(x,s)|
\delta_{y}(x)p_i(t,y; u) \; dx,
\]
where ${\cal L}$ is the generator for an SDE with drift $\lambda$
and diffusion $\sigma$ with periodic boundary conditions on $[0,L]^d$.
\end{theorem}
\paragraph{Proof}
The solution is constructed as follows. Let the $\bvec{B}_i(t)$ be
independent standard Brownian motions on a probability space
$(\Omega_1, {\cal F}_1)$. Let $X^{*}_i$ for $i=1,\dots,N$ be the
$(\Omega_1, {\cal F}_1)$ random variables taking values in
$C([\tau^s_i,\infty), \mathbb{R}^d)$ which are solutions of~\eqref{eq:sdes}
subject to $X^{*}_i(\tau^s_i)=X^s_i$. The existence
and uniqueness of such solutions are guaranteed by classical theory
under the smoothness assumptions on $\lambda$ and $\sigma$. Let
$u^{(1)}(t)$ denote the unique weak solution of
equation~\eqref{eq:chp} where the particle variables $X_i$ are
replaced by $X_i^{*}$, given by Theorem~\ref{thm:one}. Define
independent random variables
\[
\tau_i^{}\colon C([0,T],V) \rightarrow [\tau_i^s,\infty)
\]
on a second probability space $(\Omega_2, {\cal F}_2)$ such that
$\mathbf{P}{ \tau_i(u) > t}=\int_\Omega p_i(t,y; u) \;dy$. Let
$j(u)$ minimise $\tau_i(u)$ over $i=1,\dots,N$ (i.e., be the
first particle to be annihilated given a phase trajectory). Now on the
joint probability space $(\Omega_1 \times \Omega_2, {\cal F}_1 \times
{\cal F}_2)$, define the $[0,\infty)$ valued random variable
$\tau^e_{j({u^{(1)}})}:=\tau_{j({u^{(1)}})}$.
% Notice that
% $\tau^e_{j(u)}$ so defined are ${\cal F}_{\tau^e_{j(u)}} \times
% {\cal F}_2$ measurable, because the definition of $P_i(t,u)$ depends
% on $u$ only on the time interval $[0,t)$ (as in~\eqref{}).
Finally, let $u(x,s)=u^{(1)}(x,s)$
%, $P_i(t,u)=P^{(1)}_i(t,u)$,
and
$X_i(s)=X_i^{*}(s)$ for $0\le s\le T^{(1)}:= \tau^e_{j(u^{(1)})}$.
% Then, $u(s,x)$ is an ${\cal F}_{1,s}\times {\cal F}_2$ measureable
% random variable that obeys (i) on $[0,T^{(1)}]$.
To generate solutions over the next time period, let $u^{(2)}=u^{(1)}$
on $[0,T^{(1)}]$. For time $t>T^{(2)}$, let $u^{(2)}$ equal
the weak solution of equations~\eqref{eq:chp} again with particles
$X_i^{*}$ but this time with initial data
\[
u^{(2)}(T^{(1)}, x) =u^{(1)}(T^{(1)}, x).
\]
% Then, define independent random variables $\tau^{(2)}_i$, taking
% values in the functions from $C([0,T],V)$ to $[\tau_i^s,\infty)$, for
% $i\ne j(u)$ such that $\mathbf{P}{ \tau_i^{(2)}(u)> t}= \int_\Omega
% p_i(t,y;u) \; dy$ on
% $[\tau_i^s, T^{(1)}]$ and $=P^{(2)}_i(t,u)$ on $[T^{(1)},\infty)$.
Let $k(u)$ minimise $\tau_i(u)$ over $i\ne j(u)$ and set
$\tau^e_{k({u^{(2)}})}:=\tau_{k({u^{(2)}})}$. Finally, let
$u(x,s)=u^{(2)}(x,s)$, $P_i(t,u)=P^{(2)}_i(t,u)$, and
$X_i(s)=X_i^{*}(s)$ for $T^{(1)}\le s\le T^{(2)}:=
\tau^e_{k(u^{(2)})}$.
This process can be iterated $N$ times to generate solutions up to
the time when all particles have died. A solution is generated on
the time interval $[0,T]$, by solving the Cahn-Hilliard equation
on the interval where all particles are dead. We have illustrated
how the solution for the phase variable $u$, particle positions
$X_i$, and annihilation times $\tau_i^e$ can be constructed, with
the solutions pieced together by restarting the processes at the
annihilation times $\tau_i^e$. The construction specifies $u,
X_i,$ and $P_i$ uniquely on each time interval $[0, T^{(1)}]$,
$(T^{(1)}, T^{(2)}]$, $\dots$, $[T^{(N)}, T]$.
\hfill$\Box$
\section{Asymptotic Analysis}\label{sec:asym}
We present a non-rigorous explanation of the properties of the coupled
Cahn-Hilliard particle system~\eqref{eq:chp}. We show that the total
phase, taken as the sum of that in the PDE $\int_\Omega u(t, x) dx$
and a constant amount for each alive particle, $X_i$ for $i\in{\cal
P}(t)$, is conserved. We perform an asymptotic analysis of the
equations to gain a free boundary problem. The analysis provides
conditions that the probability that phase is transferred from
particle to $u$ only at the boundary $\{u\approx 0\}$ tends to one.
\paragraph{Conservation of Phase}
Let the total phase of the coupled Cahn-Hilliard particle system be
denoted
\begin{equation}
\label{eq:phase}
{\cal W}(t)= \int_\Omega u(t,x) dx +
\sum_{i\in {\cal P}(t)} \frac{\gamma_1}{\gamma_2} R^{d}.
\end{equation}
We show that on average the total phase is constant.
% Let $P_i(t,u)$ be the probability particle $i$ is alive at time $t$
% given a phase trajectory $u^*\in C([0,T],V)$ (see Theorem~\ref{}). Then,
% \[
% \frac{d}{dt} P_i(t,u^*) =
% - \frac{\gamma_2}{R^d} \mean{
% \int_\Omega |\nabla u^*|(x) \delta_{X_i(t)}(x) \D x}.
% \]
% Substitute the solution $u$ and taking expectations
% \[
% \frac{d}{dt} \mean{P_i(t,u)} = \frac{d}{dt} P_i(t)=
% - \frac{\gamma_2}{R^d} \mean{}
% \int_\Omega |\nabla u|(x) \delta_{X_i(t)}(x) \D x,
% \]
% where $P_i(t)$ is the probability that particle is alive at time $t$.
For $t> \tau^s_i$, we have from Theorem~\ref{thm:two}
\[
\frac{dp_i(t,y ;u)}{dt} = {\cal L} p_i(t,y; u)
- \frac{\gamma_2}{R^d } \int_\Omega |\nabla u|(x)
\delta_{y} (x) p_i(t,y ;u) \;dx.
\]
Denote the probability the particle is alive at time $t$ (given a
phase trajectory $u$) by $P_i(t)$ (resp. $P_i(t; u)$), so that
$P_i(t; u)=0$ for $t<\tau_i^s$ and $=\int_\Omega p_i(t,y; u) \;dy$
for $t\ge \tau_i^s$ and $P_i(t)=\mean{P_i(t; u)}$. Then for $t>\tau_i^s$,
\begin{align*}
\frac{dP_i(t ;u)}{dt} =& \int_\Omega {\cal L} p_i(t,y ;u) \;dy -
\frac{\gamma_2}{R^d}
\int_\Omega \int_\Omega |\nabla u|(x) \delta_{y}(x) p_i(t,y ;u)\;dx dy \\
=&
- \frac{\gamma_2}{R^d} \Mean{
\int_\Omega |\nabla u^*|(x) \delta_{X_i(t)}(x) \;dx \Big| u^*=u},
\end{align*}
where we have used the boundary conditions to eliminate the first term.
Now average over the phase variable
\begin{align*}
\frac{d\mean{P_i(t ; u)}}{dt}= \frac{d{P_i(t )}}{dt}
=& - \frac{\gamma_2}{R^d} \mean{}
\int_\Omega |\nabla u|(x) \delta_{X_i(t)}(x) \,d x.
\end{align*}
The equation for $U(t):=\mean{ \int_\Omega u(t,x) \;dx}$ is
\[
dU= 0 dt + \gamma_1 \sum_{i\in{\cal P}(t)} \int_\Omega \mean{|\nabla
u|(x) \delta_{X_i(t)}(x) } \;dx,
\]
by using the boundary conditions in the Cahn-Hilliard equation to
eliminate terms. Adding the terms together, we have for $t \ne \tau_i^s$
\begin{equation}
\label{eq:cons}
\frac{d}{dt}(U(t)+\frac{ \gamma_1}{\gamma_2} R^d \sum_{i} P_i(t))
= 0.
\end{equation}
Thus, we have
that rate of change of the expected value of ${\cal W}(t)$ is zero for
$t \ne \tau^s_i$. When $t = \tau^s_i$ some $i$, further particles are
added to the system and the value of ${\cal W}(t)$ increase by the
number of particles added times $\gamma_1 R^d/\gamma_2$.
\paragraph{Asymptotic expansion} We now perform an asymptotic
expansion in small $\epsilon$ in an effort to gain a related free
boundary problem. The analysis follows~\cite{cag:1989}. Assume that
the domain $\Omega$ can be split up into sets $\Omega_+ \cup \Gamma
\cup \Omega_-$, where $u\approx u_\pm$ on $\Omega_\pm$ and
$\Gamma=\{u=(u_+ + u_-)/2\}$.
First, we perform an outer expansion away from the interface
$\Gamma$. Write
\[
u=u^0+\epsilon u^1+\cdots; \quad \theta=\theta^0+\epsilon
\theta^1+\cdots.
\]
\paragraph{Outer order 1 expansion}
\begin{align*}
u^0_t=&-\Delta \theta^0 + \gamma_1\sum_{i\in {\cal P}(t)}
|\nabla u^0| \delta_{X_i}\\
-\theta^0=&-f(u^0).
\end{align*}
In the outer layer, $u^0\approx u^\pm$ (the minima of the potential
$F$). Thus the order 1 expansion is simplified to
\begin{equation}
\label{eq:asym-pde}
u^0_t= \Delta f(u^0) + \gamma_1\sum_{i\in {\cal P}(t)} |\nabla u^0|
\delta_{X_i};\quad
\theta^0=f(u^0).
\end{equation}
Note that the second equation is invertible for $u^0 \approx u_\pm$
and $u^0$ may be eliminated from these equations. We cannot assume the
driving term appears at lower order. The equations are parabolic
PDEs; this contrasts with the Cahn-Hilliard equation where the
interface motion is velocity order $\epsilon$ and the slow motion in
the first order expansion can be treated as an elliptic problem.
Near to the interface $\Gamma$, the solution $u$ varies from $u_+$ to
$u_-$ over an interface of width $\epsilon$. To analyse this
asymptotically, we introduce new coordinates: let $r$ denote the
signed distance from $\Gamma$ (with $r>0$ denoting points in $\Omega_+$)
and $s$ arc length along the interface $\Gamma$. Change variables to
$(r,s)$ space by setting
\begin{gather*}
|\nabla u|^2= u_r^2 + 2u_r u_s (r_x\; s_x + r_y\;s_y) + u_s^2
|\nabla s|^2 \\
\Delta u = u_{rr} + u_{ss} |\nabla s|^2 + u_r \Delta r + u_s \Delta
s \\
\tfrac{d}{dt}u(x_1,x_2,\dots)= \tfrac{d}{dt}u(r,s) + u_r r_t + u_s s_t.
\end{gather*}
The inner expansion is performed by choosing $z=r/\epsilon$ and expanding
\[
U=U^0+\epsilon U^1+\cdots;
\quad
\Theta=\Theta^0+\epsilon \Theta^1+\cdots.
\]
Then, substituting into the PDE, we have
\begin{gather*}
\begin{aligned}(U_t + &\epsilon^{-1} U_z r_t + U_s s_t) \\
=& (\epsilon^{-2} \Theta_{zz}
+ \Theta_{ss} |\nabla s|^2 + \epsilon^{-1}\Theta _z\ \Delta r + \Theta
_s \Delta s) \\
&+ \gamma_1\sum_{i\in{\cal P}(t)} (\epsilon^{-2} U_z^2 + 2
\epsilon^{-1} U_z U_r(r_x s_x + r_y
s_y) + U_s^2 |\nabla s|^2 )^{1/2} \delta_{X_i}
\end{aligned}\\
-\Theta= \epsilon^2 (\epsilon^{-2} U_{zz} + U_{ss} |\nabla s|^2
+\epsilon^{-1} U_z \Delta r + U_s \Delta s) -f(U).
\end{gather*}
Hence
\begin{gather*}
\begin{aligned}
\Theta_{zz}
+ \epsilon ( - U_z r_t + \Delta r \Theta_z )
-\epsilon^2 (U_t + U_s s_t - (\Theta_{ss}
|\nabla s|^2 + \Theta_s \Delta s ))&\\
+ \gamma_1 \sum_{i\in {\cal P}(t)} (\epsilon^{2} U_z^2 + 2
\epsilon^{3} U_z U_r (r_x s_x + r_y
s_y) + \epsilon^4 U_s^2 |\nabla s|^2 )^{1/2} \delta_{X_i} &=0
\end{aligned}\\
U_{zz}- f(U) + \Theta + \epsilon U_z \Delta r +
\epsilon^2( U_{ss} |\nabla s|^2 +U_s \Delta s )=0.
\end{gather*}
\paragraph{Inner order $1$ expansion:}
\begin{gather*}
\Theta^0_{zz}=0; \\
U_{zz} - f(U)+\Theta^0=0.
\end{gather*}
The only bounded solutions that give an interfacial profile for $U$ are
$\Theta^0=0$. Then the first order term $U^0(s,z)=\psi(z)$, where $\psi$
obeys
\begin{equation}
\label{eq:profile}
\psi_{zz}-f(\psi)=0,
\quad
\psi(\pm\infty)=u^{\pm},\quad \psi(0)=0.
\end{equation}
As $f$ is odd, $\psi(x)=-\psi(-x)$.
\paragraph{Inner order $\epsilon$ expansion:}
\[
\Theta^1_{zz} = U^0_z r_t - \Delta r \Theta_z^0 - \gamma_1
\sum_{i\in {\cal P}(t)}|U^0_z|
\delta_{X_i}.
\]
Using $\Theta^0=0$ and $U^0(z)=\psi(z)$, we can integrate to find for $z>0$
\begin{equation}
\label{eq:monday}
\Big[
\Theta^1_z
\Big]_0^z= r_t \psi(z)
- \gamma_1 \sum_{i\in {\cal P}(t)} \int_0^z \psi'(z')\;
\delta_{X_i}(z' \epsilon,s) \;dz'
+c_1(s,t)
\end{equation}
and again
\begin{equation}
\label{eq:theta1}
\Theta^1(z)
= r_t \Psi(z)
-\gamma_1 \sum_{i\in {\cal P}(t)} P_i(z,s,\epsilon)
+c_1(s,t) z + c_2(s,t),
\end{equation}
where $c_1(s,t)$ and $c_2(s,t)$ are constants of integration,
\[
P_i(s,z,\epsilon):=
\int_0^z \int_0^ {z'} \psi'(z'') \delta_{X_i}(z'' \epsilon, s) \;dz''
\]
and $\Psi(z):=\int_0^z \psi(z)\;dz$. Consider the second equation
\[
U^1_{zz} - f'(U^0) U^1 +\Theta^1 = -\kappa \psi'(z),
\]
where $\kappa=\Delta r$ (the mean curvature of $\Gamma$, evaluated at $z=0$).
Let $\Lambda:= (\partial_z)^2 - f'(\psi(z))$ with domain $L^2(-\infty,
\infty)$; then
\[
\Lambda U^1 = -\Theta^1 - \kappa \psi'(z).
\]
$\Lambda$ is an operator with one eigenfunction $\psi'$ that has zero
eigenvalue (because of~\eqref{eq:profile}). Thus the solvability
condition for $\Lambda \Phi = g$ is orthogonality of $g$ with respect
to $\psi'$. The solvability condition is
\[
\int_{-\infty}^\infty \Theta^1(z) \psi'(z) \; dz + \kappa
\int_{-\infty}^\infty (\psi'(z))^2 \; dz=0.
\]
Now, substitute from~\eqref{eq:theta1}, to get
\begin{multline*}
\int_{-\infty}^\infty \Big(r_t \Psi(z)
-\gamma_1 \sum_{i\in {\cal P}(t)} P_i(z, s, \epsilon) +c_1(s,t) z +
c_2(s,t)\Big ) \psi'(z) \;dz\\
+ \kappa \int_{-\infty}^\infty (\psi'(z))^2 \;dz=0.
\end{multline*}
This implies as the integral of the odd term $z \psi'(z)$ disappears
that
\begin{multline*}
\int_{-\infty}^\infty c_2(s,t) \psi'(z) \;dz \\
= - r_t \int_{-\infty}^\infty \Psi(z)\psi'(z) \;dz
+\gamma_1 \sum_{i\in {\cal P}(t)} \int_{-\infty}^\infty P_i(z,s,\epsilon)
\psi'(z)\;dz-
\kappa A,
\end{multline*}
where $A:=\int_{-\infty}^\infty (\psi'(z))^2 dz$. Thus,
\[
(u_+ - u_-) c_2(s,t)= - r_t \int_{-\infty}^\infty \Psi(z)\psi'(z) \;dz
+ \gamma_1 \sum_{i\in {\cal P}(t)} \int_{-\infty}^\infty P_i(z, s,
\epsilon) \psi'(z) \;dz - \kappa A.
\]
\paragraph{Matching conditions} Consider an interface at position
$Y(t,\epsilon)$ and suppose that $Y(t,\epsilon)=Y^0 + \epsilon
Y^{1,\epsilon} + t Y^{1,t} + \order{t\epsilon+\epsilon^2+t^2}$. By
choosing $Y(0,\epsilon)=Y^0$, we have $Y^{1,\epsilon}=0$. Look for
matching conditions at $z=(x-Y(t))/\epsilon$ by writing
\begin{gather}
\begin{split}\label{eq:rev}
\Theta(z,t,\epsilon)=& \theta(Y (t,\epsilon) + \epsilon z, t,
\epsilon)\\
=& \theta^0(Y^0, t)
+ \epsilon \theta^1(Y^0, t) + z \;\epsilon\; \theta_x^0 ( Y^0, t))
+ \order{ \epsilon^2 + \epsilon^2 z+ \epsilon t}.
\end{split}
\end{gather}
It is easy to show that as $z\rightarrow \infty$ with
$\epsilon z\rightarrow 0$,
\[
\int_0^z \psi'(z') \delta_{X_i}(z'\epsilon,s) \;dz' \rightarrow
\begin{cases}
\psi(z),&\text{ $\norm{X_i - (0,s)}\le R$}\\
0,&\text{otherwise}
\end{cases}.
\]
We introduce $\delta^{d-1}_{X_i}(s)$, which takes value one when the
ball of radius $R$ centred at $X_i$ includes the point $s$ on the
interface $\Gamma$, value zero otherwise.
Taking limits with $\epsilon z\rightarrow 0$ and $z \rightarrow
\pm\infty$, we have from~\eqref{eq:monday} and~\eqref{eq:rev}
\begin{align}
\theta^0_r=&
u_{\pm} \Big(r_t - \gamma_1 \sum_{i\in {\cal P}(t)}\delta_{X_i}^{d-1}\Big)
+ c_1(s,t).
\label{eq:asym-a}
\end{align}
Neglecting the
$\order{\epsilon z}$ terms from~\eqref{eq:theta1}, we have that
\begin{equation}
\label{eq:tuesday}
\Theta^1(z,s,t) = (r_t -\gamma_1 \sum_{i\in {\cal P}(t)}
\delta_{X_i}^{d-1}) \Big(\int_0^z
(\psi(z)-u_+) \;dz' + u_+ z \Big)+c_1(s,t) z + c_2(s,t).
\end{equation}
Taking limits with $\epsilon z\rightarrow 0$, we can pick out the
following from~\eqref{eq:rev},~\eqref{eq:tuesday}, and~\eqref{eq:asym-a}:
\begin{align}
\theta^1= & c_2(s,t) +
\Big( r_t - \gamma_1 \sum_{i\in {\cal P}(t)}\delta_{X_i}^{d-1}\Big)
\int_0^\infty
(\psi(z) - u_+) \;dz. \label{eq:asym-b}
\end{align}
Finally then, substituting for $c_2(s,t)$
\begin{align*}
\theta^1 =&\frac {-1} {u_+ - u_-}
\Big( \Big( r_t - \gamma_1 \sum_{i\in {\cal P}(t)}
\delta_{X_i}^{d-1}\Big)
\int_{-\infty}^{\infty}
\Psi(z)\psi'(z) \;dz
+ \kappa A\Big)\\
& + (r_t - \gamma_1 \sum_{i\in {\cal P}(t)} \delta_{X_i}^{d-1})
\int _0^\infty (\psi(z) - u_+) \;dz.
\end{align*}
Now note that
\[
\int_{-Z}^Z \Psi(z) \psi'(z)\;dz = \Big[ \Psi(z) \psi(z)\Big]_{-Z}^Z -
\int_{-Z}^Z \psi(z)^2 \;dz.
\]
Using $u_+=-u_-$, we have integrating by parts
\begin{align*}
\frac{-1}{u_+ - u_-}&\int_{-\infty}^\infty \Psi(z) \psi'(z) \; dz +
\int_0^\infty \psi(z) - u_+\; dz\\
=&\int_0^\infty \frac{2}{u_+-u_-} \psi^2 - 2 \psi + u_+ \; dz
=\int_0^\infty \Big( \sqrt{u_+}- \frac{\psi}{\sqrt{u_+}}\Big)^2 \; dz
=:B.
\end{align*}
Thus,
\begin{equation}
\label{eq:asym-c}
\theta^1 =- ( r_t - \gamma_1 \sum_{i\in {\cal P}(t)}
\delta_{X_i}^{d-1}) B -\frac{1}{u_+- u_-} \kappa A.
\end{equation}
Define $V$ to be the velocity of the interface $\Gamma$ into
$\Omega_-$, so that
$V=r_t$. Collecting our results, for small $\epsilon$, the dynamics of
$\Gamma$ can be determined by the following free boundary problem:
Inside $\Omega^\pm$ \eqref{eq:asym-pde} gives that $\theta$ obeys the
following PDE
\[ (f^{-1}(\theta))_t = \Delta \theta^0 + \gamma_1
\sum_{i\in{\cal P}(t)}
|\nabla f^{-1} (\theta) |
\delta_{X_i}\qquad
\]
subject to boundary conditions on $\Gamma$ given by~\eqref{eq:asym-c}
\[
\theta=-\frac{1}{u_+-u_-}\epsilon A \kappa - V B \epsilon + B
\epsilon
\gamma_1
\sum_{i\in {\cal P}(t)} \delta_{X_i}^{d-1}.
\]
This is the analogue of the Gibbs-Thompson
relation. The interface
$\Gamma$ moves with velocity given by~\eqref{eq:asym-a}
\begin{equation}
\label{eq:intVel}
V=
\frac1 {u_+ - u_-}\Big[ \frac{d \theta}{d \bvec{n}}\Big]_-^+
+ \frac{\gamma_1}{u_+ - u_-} \sum_{i\in {\cal P}(t)}
\delta_{X_i}^{d-1},
\end{equation}
where $[\frac{\partial \theta}{\partial \bvec{n}}]_-^+$ denotes the jump in $\frac{\partial \theta}{\partial
\mathbf{n}}$ on $\Gamma$.
The particles $X_i$ satisfy the SDE
\[
d X_i = \bvec{\lambda}(X_i) \,d t + \sigma(X_i) \,d\bvec{B}_i(t).
\]
To understand the annihilation of the particles, we perform
asymptotics on the rate of annihilation.
\paragraph{Inner Expansion} Under the condition $\epsilon\ll R$, the time for a particle of radius $R$ with
speed $v$ to cross an interface of width $\epsilon$ is order $R/v$.
Across the interface $u(r)\approx \psi(r/\epsilon)$. Thus, the rate of
annihilation in crossing the interface is
\begin{align*}
\frac{\gamma_2}{R}
\int |\nabla u|(x) \delta_{X_i}(x) \;dx
\approx &\frac{\gamma_2}{R^{d}}
\int \psi'(r/\epsilon) \delta_0((r,s)) \;dr\;ds \\
\rightarrow&
\frac{\gamma_2}{R^{d}} R^d \pi_{d-1}
(u_+ - u_-)
= \frac{\gamma_2}{R} (u_+ - u_-) \pi_{d-1},
\end{align*}
as $\epsilon\rightarrow 0$, where $\pi_d$ denotes the volume of a ball
of radius one in dimension $d$. This holds on a time interval of order
$R / v$. Consequently, the probability of annihilation in crossing
the interface is $1- \exp(-K \gamma_2 (u_+ - u_-) \pi_{d-1}/ v)$, some
constant $K>0$. Thus, if we take limits with $\gamma_2 \rightarrow
\infty$, the probability the particles die when they hit the interface
tends to one.
\paragraph{Outer expansion} We expect $|\nabla u|$ to be of order
$e^{-K/\epsilon}$, some constant $K$, in the interior of
$\Omega_{\pm}$ (for $f(u)=u^3-u$, the solution of~\eqref{eq:profile} is
a $\tanh$ profile). Thus, for a particle $X$ that is an order 1
distance from $\Gamma$,
\[
{\dfrac{\gamma_2}{R^d}}
{\int} {|\nabla u|(x)} \delta_X(x) \; dx \approx \frac{\pi_d R^{d}
\gamma_2}{R^d}
e^{-K/\epsilon}= \frac{\pi_d}{R}\; \gamma_2\; e^{-K/\epsilon},
\]
some constant $K$. Thus we expect the particles to live for an
exponentially long amount of time when moving about the interior of
$\Omega_{\pm}$.
Suppose that the time scale for the annihilation of particles at the
boundary ($R/\gamma_2(u_+-u_-) \pi_{d-1}$) is much less than the time
to cross the interface; that is,
\[
\gamma_2 \gg
\frac{ v}{ (u_+ - u_-) \pi_{d-1}},
%\qquad
%\frac{1}{ \gamma_2 (u_+- u_-) \pi_{d-1}} \ll \frac{\epsilon}{v},
\]
then near $\Gamma$ a particle can be expected to live for a time
$R/\gamma_2 (u_+ - u_-) \pi_{d-1}$. Hence, substitute our asymptotic
expression for the interface velocity $V$,
\[
r\Big(t+\frac{R}{\gamma_2 (u_+- u_-) \pi_{d-1}}\Big)- r(t)=
R\frac{\gamma_1}{\gamma_2}
\frac{1}{ (u_+ - u_-)^2 \pi_{d-1} }
\sum_{i\in {\cal P}(t)} \delta_{X_i}^{d-1}+ \order{\epsilon}.
\]
Thus in the $\gamma_2 =\gamma_1$ large limit, on impact of a particle
at the boundary $\Gamma$, the boundary moves by a distance
$R\gamma_1/\gamma_2 (u_+ - u_-)^2 \pi_{d-1}$. Thus we see the impact
of a particle on the interface causes the interface to move by an amount
linear in the particle radius. Further, from~\eqref{eq:cons}, we have
that each particle carries an amount of phase linear in the particle volume
(i.e., $R^d$).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
\section{A numerical method}\label{sec:fem}
The following numerical method has been implemented for dimension
$d=2$ and results are presented in~\S\ref{sec:res}.
Consider a triangulation $T^h$ of $\Omega$ consisting of rectangular
elements. Associated with $T^h$ is the finite element space $S^h$ of
functions that are bilinear functions on the rectangular elements of
$T^h$. Introduce the notation $(\cdot, \cdot)^h$ to be the
$L^2(\Omega)$ inner product evaluated by numerical quadrature based on
nodal values. The quadrate rule for $(\cdot,\cdot)^h$ is chosen to
give exactly $(\cdot,\cdot)$ for elements of $S^h$. Mass lumping is
denoted by $(\cdot, \cdot)^L$. Choose a basis $\chi_j$ for $S^h$ such
that $\sum \chi_j=1$ and equals one at exactly one nodal value and
equals zero at other nodal values. By mass lumping, we mean
\[
(\chi_k, \chi_j)^L =
\begin{cases}
(1, \chi_j),& k=j;\\
0, & k \ne j.
\end{cases}
\]
If $u$ has coordinates $u_j$ with respect to $\chi_j$ then $ (u,
\chi_j)^L = (1,\chi_j) u_j$. Fix a time step $\tstep>0$, our finite
element method generates approximations $X^n_i, \theta^n, u^n$ to
$X_i(n\tstep)$, $\theta(n\tstep)$, $u(n\tstep)$ as follows: Given
particle positions $X_i^n$ and $(u^n, \theta^n)\in S^h \times
S^h$, find $(u^{n+1}, \theta^{n+1}) \in S^h \times S^h$ such that
\begin{gather*}
\tstep^{-1} ( u^{n+1} - u^n, \chi)^L + (\nabla \theta^{n+1}, \nabla
\chi)^h =
\gamma_1 \sum_{i\in{\cal P}(t)} ( |\nabla u^n| \delta_{X_i^n}, \chi)^h\\
-(\theta^{n+1}, \chi)^L = -\epsilon^2 (\nabla u^{n+1}, \nabla \chi)^h -
(f(u^{n+1}), \chi)^L,
\end{gather*}
for all $\chi \in S^h$, where $u^0$ is take to be an approximation to
$u(0)$ in $S^h$. To generate $X_i^{n+1}$, we use
the explicit Euler method:
\[
X_i^{n+1} - X_i^n = \lambda(X_i^n) \;\tstep + \sigma(X_i^n) B^{n,i}_\tstep,
\]
where $ B^{n,i}_\tstep$ are IID Gaussian random variables with mean
$0$ and variance $I \tstep$. The annihilation time of each particle
is determined after computation of $\bvec{u}^{n+1}$ and $X_i^{n+1}$:
generate IID uniform random variables ${\cal U}^{n,i}$ on $[0,1]$ and
annihilate particle $i$ if
\[
{\cal U}^{n,i} \le \tstep \; \frac{\gamma_2}{R^d} (|\nabla u^{n}|
\delta_{X^{n}_i},1)^h.
\]
% This ensures that the conservation law~\eqref{eq:cons} is upheld by the
% numerical method.
Let
\[
%M_{ii} = \sum_j (\chi_j, \chi_i)^h=(\chi_i,1),
%\quad M_{ij}=0 \text{ for $i\ne j$},\quad
M_{ij}= (\chi_i, \chi_j)^L, \quad
A_{ij}=(\nabla \chi_i, \nabla \chi_j)^h.
\]
Let $\hat{u}^n$ be such that
\[
(\hat{u}^n , \chi)^h=
({u}^n , \chi)^h
+ \gamma_1 \tstep M \sum_{i\in {\cal P}(t)} (|\nabla u^n|
\delta_{X_i^n}, \chi)^h,\qquad \chi\in S^h.
\]
Using $\bvec{u}^n, \hvec{u}^n$ to denote coordinates of $u^n,
\hat{u}^n$ with respect to $\chi_j$, we are left to solve the
following nonlinear system:
\begin{gather*}
\frac{1}{\tstep} M (\bvec{u}^{n+1} - \hvec{u}^n) = -A
\bvec{\theta}^{n+1}\\
-M\bvec{\theta}^{n+1} = -\epsilon^2 A \bvec{u}^{n+1} + c \; M\; \hvec{u}^n -
M \phi( \bvec{u}^{n+1}),
\end{gather*}
where $f(u)= \phi(u) - c\;u$, the difference of two monotone
functions. To solve this system, we apply an iterative method of
Lions and Mercier~\cite{lions:1979} (see also
Copetti~\cite{copetti:1991} and Barrett-Blowey~\cite{barrett:1997} for
applications to phase-field equations).
Combining the two equations, we have
\[
A^{-1} \Big( \frac{M \bvec{u}^{n+1} - M \hvec{u}^n}{\tstep} \Big)
+ \epsilon^2 M^{-1} A \bvec{u}^{n+1} - c \; \hvec{u}^n + \phi(\bvec{u}^{n+1}) -
\lambda \bvec{1}
=0.
\]
Note that $A$ has a one dimensional null space spanned by $1$ as we
work with Neumann conditions and hence the Lagrange multiplier $\lambda$
provides a solution. Let
\[
B(\bvec{y}):= A^{-1} M \frac{y-\hvec{u}^n}\tstep + \epsilon^2 M^{-1} A y
-c\; \;\hvec{u}^n.
\]
Then the equation becomes
\begin{equation}
\label{eq:a}
B(\bvec{u}^{n+1}) + \phi(\bvec{u}^{n+1}) - \lambda \mathbf{1}=0.
\end{equation}
Consider a guess $(\bvec{u}^{n+1}_j, \lambda_j)$ for $\bvec{u}^{n+1}$ and the
Lagrange multiplier $\lambda$. In the following iteration $\mu$ is a
relaxation parameter.
\begin{enumerate}
\item Compute
\[
Z_1:= \bvec{u}^{n+1}_j - \mu B( \bvec{u}_j^{n+1}) + \lambda_j \mu
\mathbf{1}=0
\]
and solve
\begin{equation}
\label{eq:b}
\bvec{u}^{n+1}_{j+\tfrac 12} + \mu \phi(\bvec{u}^{n+1}_{j+\tfrac 12})
= Z_1.
\end{equation}
This gives an approximation $\bvec{u}^{n+1}_{j+1/2}$.
\item Now compute $u_{j+1}^{n+1}$, $\lambda_{j+1}$ such that
\begin{align*}
\bvec{u}^{n+1}_{j+1} + \mu B(\bvec{u}^{n+1}_{j+1}) - \lambda_{j+1}
\mu \mathbf{1}
= & \bvec{u}^{n+\tfrac12}_{j+\tfrac 12} - \mu
\phi(\hvec{u}^n_{j+\tfrac 12})\\
= &2 \bvec{u}_{j+1/2}^{n+1}-Z_1,
\end{align*}
To do this, note
\begin{equation}
\label{eq:c}
\bvec{u}^{n+1}_{j+1} + \mu B(\bvec{u}^{n+1}_{j+1}) - \lambda_{j+1} \mu
\mathbf{1}
= \xi,\quad\xi:=2\bvec{u}^{n+1}_{j+1/2}-Z_1,
\end{equation}
Multiply by $A$,
\[
A \bvec{u}^{n+1}_{j+1} + \mu A B(\bvec{u}^{n+1}_{j+1})
=A \xi
\]
and substituting for $B$
\[
A \bvec{u}_{j+1}^{n+1} +
\mu\Big[M \frac{\bvec{u}^{n+1}_{j+1} -\hvec{u}^n}{\tstep}
+\epsilon^2 A M^{-1} A \bvec{u}^{n+1}_{j+1} -c A \hvec{u}^n\Big]
= A \xi.
\]
Let
\[
A_1= (\tstep^{-1} M +\epsilon^2 A M^{-1}
A), \quad B_0=(\tstep^{-1} M + c\;A)\;\hvec{u}^n.
\]
Then we can find $\bvec{u}_{j+1}^{n+1}$ by solving
\[
A_2 \bvec{u}_{j+1}^{n+1}=Z_2,
\quad
A_2=A+\mu A_1,\quad Z_2 =A \xi + \mu B_0.
\]
and the Lagrange multiplier from~\eqref{eq:c}
\[
\lambda_{j+1}=\frac{1}{\mu (\mathbf{1}, \mathbf{1} )} (\mathbf{1},
\bvec{u}_{j+1}^{n+1}- \xi + \mu B(\bvec{u}_{j+1}^{n+1})).
\]
\end{enumerate}
The iteration gives approximations $\bvec{u}^{n+1}_{j+1}$ and Lyapunov
multipliers $\lambda_j$. It is shown in Copetti~\cite{copetti:1991}
that the sequence $\bvec{u}^{n+1}_{j}$ converges to $\bvec{u}^{n+1}$
and \eqref{eq:b} and \eqref{eq:c}
give that \eqref{eq:a} holds.
The triangulation $S^h$ is initially uniform but is adapted according
to the size of $E_\tau=\norm{\mathbf{1}_\tau}^{-1} (\mathbf{1}_\tau,
|\nabla u|)$, where $\mathbf{1}_\tau$ is the indicator function on the
element $\tau$. An element $\tau$ in the triangulation is
\text{coarsen}ed if $E_\tau
E_{\text{refine}}$. The mesh is adapted every $N_{\text{adapt}}$ time
steps. The maximum and minimum size of a rectangle is restricted.
Further, regularisation of the triangulation is performed by the
software package (Deal.II~\cite{BHK:2001}) itself.
The scheme further implements an elementary adaptive time stepping
algorithm. There are two time scales in the interfacial velocity as
indicated in~\eqref{eq:intVel}. When a particle $X_i$ is an order one
distance from the interface $\Gamma$, the interface moves on a time
scale $\epsilon$ and a large time step $\tstep_{+}$ is used. When the
particles $X_i$ are near the interface, the interfacial dynamics are
faster and to capture this is a smaller time $\tstep_-$ is employed.
The algorithm switches between $\tstep_\pm$ when $E_\tstep$ crosses a
critical value $E_{\tstep}^{\text{critical}}$, where $E_\tstep$ equals
the maximum absolute value of $\gamma_1 \tstep_+ \sum_{i\in {\cal P}(t)}
|\nabla u^n| \delta_{X_i^n}$ over nodal values.
\section{Numerical results}\label{sec:res}
Numerical approximations to~\eqref{eq:dch} are now presented for the method
developed in~\S\ref{sec:fem}. We present approximate solutions for the
following parameter values: the domain is $[-2,2]\times [-2,2]$; the
initial phase $u(x,y)=\tanh( -(y+1.7)/\epsilon)$; that is, a band of
aggregate is place on the bottom of the domain. The interfacial
parameter $\epsilon=0.02$; the coupling $\gamma_1=100$; the annihilation rate
$\gamma_2=50$. Particles are introduced into the system at
the top of the domain, $y=2.0$, at a horizontal position $x$ uniformly
distributed on $[-2,2]$ at times $1/200, 2/200,\dots$. Initially $10$
particles are placed in the domain uniformly over $[-2,2]\times
[-1.2,2]$. The particles fall from the top to the bottom of the domain
according to
\[
dX=
\begin{pmatrix}
\lambda_x \\ -\lambda_y
\end{pmatrix} \,dt
+ \begin{pmatrix}
\sigma_x \\ 0
\end{pmatrix} \,d \beta(t),
\]
where $\lambda_x, \lambda_y, \sigma_x\ge 0$ are varied in the three
examples.
The parameters in the numerical method are chosen as follows: the time
steps $\tstep_+=5\times 10^{-5}$ and $\tstep_- =10^{-5}$. The time step
was changed according to critical value
$E_\tstep^{\text{critical}}=0.2$. The relaxation parameter
$\mu=0.8$. Initially the domain has five levels initially and the
adaptive mesh has between 4 and 8 levels (so boxes of size $0.125^2$
to $0.0078125^2$). The meshes are refined with $(E_\text{refine},
E_\text{coarsen}, N_\text{adapt})= (0.08, 0.03, 20)$.
The figures give the evolution of the $u=0$ contour and the
development of complicated patterns is seen, including the appearance
of overhangs and the development of cavities in the aggregate.
Figure~\ref{bild1} shows the growth of a large overhang into a
cavity. Figure~\ref{bild2} used a large value of the noise parameter
$\sigma_x$. Again a complex pattern develops. Figure~\ref{bild3} uses
the same parameters as in Figure~\ref{bild1}, except for the introduction of a
horizontal drift $\sigma_x\ne 0$. The patterns are less complex than
in~Figures \ref{bild1}-\ref{bild2}.
One feature of the figures is the occurrence of instabilities. In
Figure~\ref{bild2}, an island of $u\approx -1$ appears in the second
time frame in the lower left hand corner. The island is a result of
instability in the numerical solution, rather than the aggregation
dynamics. The island disappears in the next time frame. The particles
can cause the phase to flip between $u\pm 1$ even away from the
interface. This behaviour is not well understood, but is believed to
be the result of numerical approximation.
\begin{figure}[th]
\begin{center}
\includegraphics[width=6cm,height=55mm]{e1_10.eps}
\includegraphics[width=6cm,height=55mm]{e1_20.eps}
\includegraphics[width=6cm,height=55mm]{e1_30.eps}
\end{center}
\caption{ $\lambda_x=0$, $\lambda_y=200$, $\sigma_x=5$ at times
$t=0.0547, 0.11615, 0.17925$.} \label{bild1}
\end{figure}
\begin{figure}[th]
\begin{center}
\includegraphics[width=6cm,height=55mm]{e2_15.eps}
\includegraphics[width=6cm,height=55mm]{e2_30.eps}
\includegraphics[width=6cm,height=55mm]{e2_45.eps}
\includegraphics[width=6cm,height=55mm]{e2_60.eps}
\end{center}
\caption{$\lambda_x=0$, $\lambda_y=200$, $\sigma_x=40$ at times
$t=0.08465$, $0.1693$, $0.25395$, $0.33860$.}
\label{bild2}
\end{figure}
\begin{figure}[th]
\begin{center}
\includegraphics[width=6cm,height=55mm]{e3_10.eps}
\includegraphics[width=6cm,height=55mm]{e3_20.eps}
\includegraphics[width=6cm,height=55mm]{e3_30.eps}
\includegraphics[width=6cm,height=55mm]{e3_40.eps}
\end{center}
\caption{$\lambda_x=50$, $\lambda_y=200$, $\sigma_y=5$ at times
$t=0.0547$, $0.1094$, $0.164$, $0.21880$.} \label{bild3}
\end{figure}
\section{Conclusion}
The paper has described a new model for aggregation, by using a phase
field equation coupled to a system of SDEs. The model very naturally
incorporates features of arbitrary topology and shadowing, as well as
incorporating dynamics within the aggregate itself. We have
demonstrated that the equations can be understood in a rigorous
mathematical way, which contrasts with some of the difficult equations
suggested by other authors. A numerical method is suggested for
solving the equations and a number of examples solutions presented.
The numerical method suffers from instabilities and is also slow (it
takes two days on 1 Ghz Linux box to generate each of the test cases).
Further work should develop the linear algebra and analyse the source
of the instability. The numerical solutions computed exhibit effects
known to happen in aggregation processes, such as fingering, but as
the system is based on a diffuse interfaces, the patterns occur on a
large scale. The model will provide further insight when solutions
are computed for smaller $\epsilon$ and on longer time scales.
\subsection*{Acknowledgements}
I am grateful to the Nuffield Foundation NUF-NAL-00 for their support
of this work.
\begin{thebibliography}{10} \frenchspacing
\bibitem{BHK:2001}
{\sc W.~Bangerth, R.~Hartmann, and G.~Kanschat}, {\em {\tt deal.{I}{I}}
Differential Equations Analysis Library, Technical Reference}, IWR.
\newblock \texttt{http://gaia.iwr.uni}-\texttt{heidelberg.de/\~{}deal/}.
\bibitem{barrett:1997}
{\sc J.~W. Barrett and J.~F. Blowey}, {\em Finite element approximation of a
model for phase separation of a multi-component alloy with non-smooth free
energy}, Numerische Mathematik, (1997), pp.~1--34.
\bibitem{cag:1989}
{\sc G.~Caginalp and P.~C. Fife}, {\em Dynamics of layered interfaces arising
from phase boundaries}, SIAM J. Appl. Math., 48 (1989), pp.~506--518.
\bibitem{cook:1970}
{\sc H.~E. Cook}, {\em Brownian motion in spinodal decomposition}, Acta
Metallurgica, 18 (1970), pp.~297--306.
\bibitem{copetti:1991}
{\sc M.~I.~M. Copetti}, {\em Numerical Analysis of Nonlinear Equations arising
in Phase Transitions and Thermoelasticity}, PhD thesis, University of Sussex,
1991.
\bibitem{daprato:1996b}
{\sc G.~Da~Prato and A.~Debussche}, {\em Stochastic {Cahn}--{Hilliard}
equation}, Nonlinear Analysis, 26 (1996), pp.~241--263.
\bibitem{elliott:1989}
{\sc C.~M. Elliott}, {\em The {Cahn}--{Hilliard} model for the kinetics of
phase separation}, in Mathematical Models for Phase Change Problems,
J.F.Rodrigues, ed., Birkh\"{a}user Verlag, 1989.
\bibitem{kardar:1986}
{\sc M.~Kardar, G.~Parisi, and Y.~Zhang}, {\em Dynamic scaling of growing
interfaces}, Physics Review Letters, 56 (1986), pp.~889--892.
\bibitem{keblinski:94}
{\sc P.~Keblinski, A.~Maritan, F.~Toigo, H.~Koplik, and J.~R. Banavar}, {\em
Dynamics of rough surfaces with an arbitrary topology}, Physical Review E, 49
(1994).
\bibitem{keblinski:96}
{\sc P.~Keblinski, A.~Maritan, F.~Toigo, R.~Messier, and J.~R. Banavar}, {\em
Continuum model for the growth of interfaces}, Physical Review E, 53 (1996),
pp.~759--777.
\bibitem{lions:1979}
{\sc P.~L. Lions and B.~Mercier}, {\em Splitting algorithms for the sum of two
nonlinear operators}, SIAM, J. Numer. Anal., 16 (1979), pp.~964--979.
\bibitem{dla2}
{\sc L.~Sander}, {\em Continuum {DLA}: random fractal growth generated by a
deterministic model}, in Fractals in physics, L.~Pietronero and E.~Tosatti,
eds., North-Holland, 1986, pp.~241--246.
\bibitem{temam:1988}
{\sc R.~Temam}, {\em Infinite Dimensional Dynamical Systems in Mechanics and
Physics}, vol.~68 of Applied Mathematical Sciences, Springer-Verlag, 1988.
\end{thebibliography}
\noindent\textsc{Tony Shardlow}\\
Department of Mathematical Sciences,\\
Science Laboratories, South Road, \\
Durham University, Durham DH1 3LE, UK\\
e-mail: \texttt{Tony.Shardlow@durham.ac.uk}
\end{document}