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