\documentstyle[graphicx,multicol,epsf,aps]{revtex}
%\documentstyle[epsf]{ioplppt}
\begin{document}
\title{Mean Field Theory of Polynuclear Surface Growth}
\author{E.~Ben-Naim$^1$, A.~R.~Bishop$^1$, I.~Daruka$^{1,2}$,
  and P.~L.~Krapivsky$^3$}
\address{$^1$Theoretical Division and Center for Nonlinear Studies, 
  Los Alamos National Laboratory, Los Alamos, NM 87545}
\address{$^2$Department of Physics, University of Notre Dame, Notre
  Dame, IN 46556}
\address{$^3$Center for Polymer Studies and Department of Physics,
  Boston University, Boston, MA 02215}
\maketitle 
\begin{abstract} 
  
  We study statistical properties of a continuum model of polynuclear
  surface growth on an infinite substrate.  We develop a
  self-consistent mean-field theory which is solved to deduce the
  growth velocity and the extremal behavior of the coverage.
  Numerical simulations show that this theory gives an improved
  approximation for the coverage compare to the previous linear
  recursion relations approach.  Furthermore, these two approximations
  provide useful upper and lower bounds for a number of
  characteristics including the coverage, growth velocity, and the
  width exponent.

\noindent{PACS numbers:  02.50.Ey, 05.40.+j,  82.20.Mj, 82.60.Nh}
\end{abstract}

\begin{multicols}{2}

\section{Introduction}

Kinetics of surface growth is a fascinating field that has been the
subject of intense current research\cite{god,meakin,zhang,barabasi}.
It is well established that as the surface grows its morphology
remains scale invariant, and for example, fluctuations in the
interface height exhibit an asymptotic scaling behavior. While the
understanding of growth on one-dimensional (1D) substrates is rather
comprehensive, in the physical case of two-dimensional substrates current
theoretical understanding remains incomplete\cite{zhang,barabasi}.

In this study we focus on an appealingly simple yet non-trivial
surface growth problem, the so-called polynuclear growth (PNG)
model\cite{frank,kas,gilmer}.  The PNG model describes the
evolution of islands that nucleate at random on top of previously
nucleated islands and grow in the radial direction. This model is
appropriate for describing situations where there is a competition
between growth along the step edges and growth due to nucleation, as
is the case in polymer crystal growth\cite{keller,bassett}.

The submonolayer version of the PNG model is identical to the
classical Kolmogorov-Avrami-Johnson-Mehl (KAJM) nucleation-and-growth
process\cite{Kolm,Avrami,Johnson} where exact results for the coverage
and more detailed statistical properties are possible
\cite{sek,Sek,bradley,ep}.  Using the exact solution for the KAJM
coverage, an approximate Linear Recursion Relations (LRR) approach to
the PNG process in arbitrary dimension was suggested
\cite{kas,Evans,be,new}.  

In one dimension, the PNG model is equivalent to an ensemble of
``kinks'' (down steps) and ``antikinks'' (up steps) moving
ballistically and annihilating upon collision.  Steady states on
finite 1D substrates were obtained analytically using the fact that
kinks (antikinks) are uncorrelated\cite{frank,ben,gold}.  In the
steady state, the width scales with system size $L$ according to
$w(L)\sim \sqrt{L}$, i.e. the roughness exponent is $\chi=1/2$.
However, the non-equilibrium behavior and especially the asymptotic
time dependence remains an open problem, despite a number of studies
\cite{be,new,ben,gold,gil,krug,kr}.

Our goal is to develop a self-consistent approach that can be viewed
as a mean field theory (MFT) of the PNG process, and to compare it
with the previous LRR approximation as well as with numerical
simulations.  We will show that MFT offers a better description for
the time dependent coverage, and that the two approaches, when
combined, provide upper and lower bounds for the coverage and the
growth velocity.

The rest of this paper is organized as follows.  In the next section
we define the PNG model. In Sec. III, we develop a self-consistent
mean field approximation.  We find that possible growth velocities are
bounded from below, $v\geq v_{\rm min}$, and show that the minimal
velocity is actually selected.  We also solve for the coverage profile
in the tail regions.  In Sec.~IV, we review the LRR approach and
derive the coverage profile analytically.  In Sec.~V, we present
simulation results. Conclusions are given in Sec.~VI.
 

\section{THE PNG MODEL}

The PNG is a continuous time, off-lattice model which contains both
deterministic and stochastic events.  Starting from a flat uniform
$d$-dimensional substrate at time $t=0$, seeds of negligible size
nucleate at random times and random positions at a constant rate per
unit area. Once created, islands grow laterally at a constant velocity
in the radial direction.  When two islands in the same layer meet they
coalesce, and the joint perimeter continues growing in the corresponding
radial direction.  Meanwhile, nucleation continuously generates
additional layers on top of previously nucleated layers.  Clearly, there
are no overhangs in this model, a feature that considerably simplifies
the analysis.  Another important simplification in the PNG model is that
the nucleation rate is uniform in time as well as in space, i.e., it is
independent of the local surface structure. Without loss of generality,
we set the nucleation rate and the radial growth velocity to unity. This
can be achieved by an appropriate rescaling of space and time.

Below, we concentrate on $S_j(t)$, the uncovered fraction in the $j$th
layer at time $t$.  This important characteristic of multi-layer growth
gives the net exposed fraction of the $j$th layer, $S_j(t)-S_{j-1}(t)$,
and therefore can be used to calculate relevant statistical
properties. The variable $j$ plays the role of the local height $h$,
often used in surface growth studies. In general $\langle f(j)\rangle =
\sum_{j=0}^{\infty} f(j)[S_{j+1}(t)-S_j(t)]$, and in particular the
average height is given by

\begin{equation}
\label{h}
h(t)=\langle j\rangle= \sum_{j=1}^\infty j\left[S_{j+1}(t)-S_j(t)\right].
\end{equation} 
Fluctuations in the height are quantified by the mean square width
$w^2(t)$ defined via
\begin{equation}
\label{w}
w^2(t)=\langle j^2\rangle -\langle j\rangle^2
=\sum_{j=1}^\infty j^2 \left[S_{j+1}(t)-S_j(t)\right]-h^2(t).
\end{equation} 
We expect a linear growth in time for the average height, $h(t)\simeq
vt$, and an algebraic growth of the width, $w(t)\sim t^{\beta}$, with a
priori unknown width exponent $\beta$.  In other words, the uncovered
fraction obeys the following wave-like form
\begin{equation}
\label{form}
S_j(t)=F\left(j-vt\over t^{\beta}\right). 
\end{equation}
The argument reflects the overall shift in the position of the
wave-front, and the multiplicative scale accounts for the algebraic
widening of the front. This expression also serves as a definition of
the width exponent $\beta$, which is equivalent to Eq.~(\ref{w}).
 
Far from the front region, we anticipate the following extremal behavior
of the scaling function $F(z)$:
\begin{equation}
\label{extremal}
F(z)\sim\cases{1-\exp(-z^{\sigma_+})&$z\to\infty$;\cr
               \exp(-|z|^{\sigma_-})&$z\to-\infty$.}
\end{equation}
The exponents $\sigma_{\pm}$ thus characterize relaxation away from
the front region \cite{zhang1}. The exponent $\sigma_+$ which
describes large positive fluctuations in the height can be simply
related to the roughness exponent.  Consider a large positive height
fluctuation, $j=Avt$, with $A\gg 1$.  Such large ``towers'' can be
created only by an anomalously large number, $Avt$, of nucleation
events localized in the same region.  Given the Poisson nature of the
nucleation events, such fluctuations are suppressed exponentially.
Thus, the quantity $1-S_{Avt}(t)$ is estimated by $\exp(-t)$, but
since Eq.~(\ref{extremal}) gives
$\exp\left[-t^{\sigma_+(1-\beta)}\right]$, we conclude that
\begin{equation}
\label{sigp}
\sigma_+={1\over 1-\beta}.
\end{equation}

In this study, we focus on systems with infinite lateral size.  If the
size is finite, it eventually becomes important and the width follows
the scaling form $w(t,L)\sim L^\chi {\cal W}(t/L^z)$.  The infinite
size limit corresponds to $t\ll L^z$ where the size dependence
disappears and therefore $w(t)\sim t^\beta$ with $\beta=\chi/z$
follows\cite{zhang}.

\section{Mean Field Theory}

\subsection{One Dimension}

We start with the PNG model in one dimension where a more comprehensive
analysis is possible.  In this situation, steps nucleate in pairs and
move away from each other with a constant velocity.  The constant
nucleation rate and the growth velocity are set to unity, without loss
of generality.  Hence, the length of an island at time $t$ after birth
equals $2t$.  Consider $f_j(x,t)$, the density of gaps of length $x$ at
time $t$ in the $j$th layer.  This distribution evolves according to
\begin{equation}
\label{fjt}
{\partial f_j(x,t)\over \partial t}=2{\partial f_j(x,t)\over \partial x}
+\gamma_j(t)\left[-xf_j(x,t)+2\int_x^\infty dy f_j(y,t)\right].
\end{equation}
The spatial derivative term describes shrinkage of gaps.  The last two
terms account for changes due to nucleation and thus, are proportional
to the overall nucleation rate at the $j$th layer, $\gamma_j(t)$. The
loss term is proportional to the gap length and the gain term describes
creation of gaps from larger gaps.

Eqs.~(\ref{fjt}) contain yet unknown nucleation rates $\gamma_j(t)$
which will be chosen to satisfy the correct kinetic equations for the
uncovered fractions
\begin{equation}
\label{S}
S_j(t)=\int_0^\infty dx\, x\, f_j(x,t),
\end{equation}
and the gap (or island) densities
\begin{equation}
\label{N}
N_j(t)=\int_0^\infty dx f_j(x,t).
\end{equation}
The uncovered fraction decreases with a rate proportional to the island
density, $\dot S_j(t)=-2N_j(t)$, and by integrating Eqs.~(\ref{fjt}) we
indeed recover this exact equation. The island density changes due to
disappearance of gaps as well as due to nucleation. The total nucleation
rate is proportional to the exposed fraction of the $j$th layer and
thus, $\dot N_j(t)=-2f_j(0,t)+S_j(t)-S_{j-1}(t)$.  On the other hand,
by integrating  Eqs.~(\ref{fjt}) we obtain $\dot N_j(t)=
-2f_j(0,t)+\gamma_j(t)S_j(t)$.  Therefore, the choice
\begin{equation}
\label{gamma}
\gamma_j(t)=1-{S_{j-1}(t)\over S_j(t)}
\end{equation}
guarantees that the gap density evolves according to {\em exact} rate
equation.  The expression of Eq.~(\ref{gamma}) for the nucleation rate
$\gamma_j(t)$ in the $j$th layer is intuitively appealing since the
total unit nucleation rate should be reduced to account for nucleation
events below the $j$th layer.  Taking into account that $S_0(t)\equiv 0$
we find $\gamma_1(t)\equiv 1$ and we notice that Eq.~(\ref{fjt}) for the
first layer agrees with  the exact KAJM equation\cite{sek}.  Thus, the set
of rate equations (\ref{fjt}) with the nucleation rates (\ref{gamma})
provides a self-consistent description of the 1D PNG model. It is exact
for the first layer, and additionally, the first two moments of the gap
density satisfy the correct rate equations.  However, it is a mean-field
description since it assumes a spatially homogeneous nucleation rate
$\gamma_j(t)$.

The gap density is found to be exponential, and the formal solution reads
\begin{equation}
\label{fgjt}
f_j(x,t)=g^2_j(t)\,\exp\left[-g_j(t)x-2\int_0^t d\tau
  g_j(\tau)\right], 
\end{equation}
with $g_j(t)=\int_0^t d\tau\,\gamma_j(\tau)$.  The uncovered fraction
and the island density are evaluated using Eqs.~ (\ref{S}) and (\ref{N}):
\begin{eqnarray}
S_j(t)&=&e^{-2\int_0^t d\tau g_j(\tau)},\nonumber\\
N_j(t)&=&g_j(t)S_j(t).
\end{eqnarray}
Evaluating $d^2\ln S_j(t)/dt^2$ together with Eq.~(\ref{gamma}) and $\dot
g_j(t)=\gamma_j(t)$, leads to an infinite set of recursive 
differential equations for the uncovered fraction
\begin{equation}
\label{sj}
\ddot S_j-\dot S_j^2S_j^{-1}+2(S_j-S_{j-1})=0.
\end{equation}
Eqs.~(\ref{sj}) should be solved subject to the initial conditions
$S_j(0)=1$ and $ \dot S_j(0)=0$ for $j\geq 1$.  The recursive
structure of Eqs.~(\ref{sj}) reflects the fact that kinetics of a
given layer is {\em unrelated} to that of all layers above (thanks to
the absence of overhangs).  Additionally, Eqs.~(\ref{sj}) agree with
the nature of the PNG model implying that kinetics of a given layer
should be {\em directly} coupled only to the previous layer.

Using $S_0(t)\equiv 0$, $S_1(t)$ is determined, then $S_2(t)$, etc.  Of
course, for the first layer, the KAJM nucleation-and-growth results are
reproduced \cite{sek}
\begin{eqnarray}
\label{s1}
S_1(t)&=&e^{-t^2},\\
f_1(x,t)&=&t^2e^{-xt-t^2}\nonumber. 
\end{eqnarray}
It is also possible to solve analytically for the second layer 
\begin{eqnarray}
\label{s2}
S_2(t)&=&\cosh^2 t\,e^{-t^2},\\
f_2(x,t)&=&[t\,\cosh\,t-\sinh\,t]^2\,
e^{-(t-\tanh\,t)x-t^2}. \nonumber
\end{eqnarray}
Using the transformation, $S_j(t)=\exp\left[u_j(t)-t^2\right]$, the
differential equations (\ref{sj}) formally simplify to a directed
version of the Toda equations \cite{toda}, $\ddot
u_j=2\exp[u_{j-1}-u_j]$, with the initial conditions $u_j(0)=\dot
u_j(0)=0$ and the boundary condition $u_1(t)=0$. Despite this
simplification it is not possible to integrate these equations, and we
solve numerically for $S_j(t)$.  Fig.~1 shows how the coverage in a
given layer changes with time.  We see that the coverage quickly relaxes
onto a traveling wave form with a finite width, $S_j(t)\to F(j-vt)$.


Some quantitative characteristics of the traveling wave can be
determined analytically.  For $j-vt\gg 1$, the nonlinear term in
(\ref{sj}) is negligible and Eqs.~(\ref{sj}) become linear.  Thus, we
write
\begin{equation}
\label{ansatz}
1-S_j(t)\sim e^{-\alpha(j-vt)}, \quad j-vt\gg 1, 
\end{equation}
with a yet unknown coefficient $\alpha$.  Substituting  into
Eqs.~(\ref{sj}) gives
\begin{equation}
\label{kol}
v^2=2{e^{\alpha}-1\over \alpha^2}.
\end{equation}
The right hand side has a minimum at $\alpha=1.59362$.  Therefore any
velocity in the interval $[v_{\rm min},\infty)$ with $v_{\rm min}
=1.75735$ is possible.  Our numerical integration shows a velocity
that falls within 0.1\% of $v_{\rm min}$, thereby implying that the
minimal velocity is indeed selected.  Such minimum velocity selection
is ubiquitous and occurs for a wide class of initial 
conditions\cite{br,mur}.


\begin{figure}[t]
\vspace{-.2in}
%\centerline{\epsfxsize=9cm \epsfbox{fig1.eps}}
\centerline{\includegraphics[width=9cm]{fig1}}
\vspace{-.2in}
\noindent{\small {\bf Fig.1}
  The uncovered fraction $S_j(t)$ vs. time for layers $j=$ 20, 40, 60,
  and 80.  Clearly, the coverage follows a traveling wave form.}
\end{figure}

To obtain the asymptotic behavior in the other extreme, $z=j-vt\to
-\infty$, we first note that $S_2(t)\gg S_1(t)$, as follows from
(\ref{s1}) and (\ref{s2}).  We assume this for all layers far behind the
front; we will check this assumption $S_j(t)\gg S_{j-1}(t)$ a
posteriori.  This reduces Eqs.~(\ref{sj}) to $\ddot S_j-\dot
S_j^2S_j^{-1}+2S_j=0$ which is solved to yield
$S_j(t)=\exp(-t^2+A_jt+B_j)$.  The traveling wave form implies that
$S_j(t)$ should be a function of a single variable $z=j-vt$.  This
determines the constants $A_j$ and $B_j$, and we find $S_j(t)=F(z)\sim
\exp\left(-z^2/v^2\right)$.  One can verify that the assumption
$S_j(t)\gg S_{j-1}(t)$ is valid. 


\subsection{Arbitrary Dimension}
 
The above analysis cannot be generalized in a straightforward manner
to $d\ne 1$ since the gap distribution is intrinsically
one-dimensional. However, it is still possible to obtain a mean field
description for the uncovered fraction.  

Consider first a KAJM nucleation-and-growth process where the rate of
nucleation events $\gamma(t)$ is homogeneous in space but time
dependent.  Ignoring overlap between growing disks, the uncovered
fraction $S(t)$ decreases with time according to
\begin{equation}
\label{no}
{dS(t)\over dt}=-d\Omega_d\int_0^t d\tau \gamma(\tau) (t-\tau)^{d-1},
\end{equation}
where $\Omega_d=\pi^{d/2}/\Gamma(1+d/2)$ is the volume of the
$d$-dimensional unit sphere and $d\Omega_d$ is its surface area
\cite{Bender} ($\Gamma$ is the gamma function).  Of course,
Eq.~(\ref{no}) overestimates the decay rate since some of the area is
already covered.  Nevertheless, this may be corrected if $(-dS/dt)$ is
reduced by the uncovered fraction, $S$, 
\begin{equation}
\label{yes}
{dS(t)\over dt}=-d\Omega_d\,S(t)\,\int_0^t\, d\tau\,\gamma(\tau)
(t-\tau)^{d-1}.
\end{equation}
For a constant nucleation rate, $\gamma=1$, Eq.~(\ref{yes}) gives 
\begin{equation}
\label{first}
S(t)=\exp\left[-{\Omega_d t^{d+1}\over d+1}\right].
\end{equation}
Thus, the exact KAJM coverage is recovered.  Furthermore, a
generalization of the KAJM solution to time-dependent nucleation rates
appears to be equivalent to Eq.~(\ref{yes}) \cite{Sek}.

We now return to the PNG model. Note first that the nonlocal in time
integro-differential Eq.~(\ref{yes}) can be converted to a higher
order ordinary differential equation $d^{d+1}\ln
S_j(t)/dt^{d+1}=-d!\,\Omega_d \gamma_j(t)$.  Unlike the 1D case, it is
not possible to derive the nucleation rate self-consistently. However,
assuming a spatially homogeneous nucleation rate \cite{gilmer} implies
the total nucleation rate $\gamma_j(t)S_j(t)=S_j(t)-S_{j-1}(t)$ and
therefore Eq.~(\ref{gamma}).  We thus arrive at the following
generalization of the mean-field Eq.~(\ref{sj}) for the uncovered
fraction
\begin{equation}
\label{dsj}
{d^{d+1} \over dt^{d+1}}\,\ln S_j+d!\,\Omega_d
\left[1-{S_{j-1}(t)\over S_j(t)}\right]=0.
\end{equation}

The analysis presented in the one-dimensional case applies for
arbitrary dimensions. For example, the transformation
$S_j(t)=\exp\left[u_j(t)-\Omega_d t^{d+1}/(d+1)\right]$ reduces
Eqs.~(\ref{dsj}) to a set of generalized directed Toda equations
$d^{d+1} u_j/ dt^{d+1}=d!\Omega_d\exp[u_{j-1}-u_j]$.  Analysis of
these equations or Eqs.~(\ref{dsj}) reveals that the coverage relaxes 
to a traveling wave with a finite width.  To determine the growth
velocity, we insert the ansatz of Eq.~(\ref{ansatz}) into
Eqs.~(\ref{dsj}) to find
\begin{equation}
\label{dv}
v^{d+1}=d!\Omega_d\,{e^{\alpha}-1\over \alpha^{d+1}}.
\end{equation}
This provides the lower bound for the growth velocity, $v_d\geq
v_d^{\rm min}$.  Again the minimal velocity should be selected, and
for example $v_0=1$, $v_1=1.75735$, $v_2=1.67115$, and $v_d\simeq
\sqrt{2\pi e/d}$ when $d\to\infty$.  Furthermore, in finite dimensions
mean field theory predicts universal exponents $\beta=0$,
$\sigma_+=1$, and $\sigma_-=2$. (Note that the relationship of 
Eq.~(\ref{sigp}) is obeyed).

In the zero-dimension limit, the behavior changes qualitatively.
Indeed, for the $d=0$ case Eqs.~(\ref{dsj}) become {\em linear},
\begin{equation}
\label{lin}
{d S_j\over dt}+S_j=S_{j-1}.  
\end{equation}
Solving (\ref{lin}) recursively yields the uncovered fraction,
$S_j(t)=e^{-t}\sum_{i=0}^{j-1} t^i/i!$.  Alternatively, by treating the
variable $j$ as continuous, this set of linear equations reduces to a
simple convection-diffusion equation. Consequently, the width becomes
diffusive, i.e., $\beta=1/2$.

\section{LINEAR RECURSION RELATIONS}

The Linear Recursion Relations (LRR) approach employs the exact
uncovered fraction $S_1(t)$ in the first layer, provided by the KAJM
solution (\ref{first}) \cite{kas,be,ah,ran}.  Nucleation in the
$(j+1)$st layer proceeds only on the already covered fraction of the
$j$th layer, formed with rate $-dS_j/dt$.  Subsequent covering proceeds
as in the KAJM, so one can anticipate that nucleation events in the time
interval $(\tau,\tau+d\tau)$ should provide contribution
$S_1(t-\tau)\left[-\dot S_j(\tau)\right] d\tau$ to the exposed fraction
$S_{j+1}(t)-S_j(t)$ in the $(j+1)$st layer.  This leads to a recursion
relation between adjacent layers:
\begin{equation}
\label{next}
S_{j+1}(t)=S_j(t)-\int_0^t d\tau\, S_1(t-\tau)\,{dS_j(\tau)\over d\tau}.
\label{tj}
\end{equation}
The first layer coverage (\ref{first}) can be recovered by setting the
substrate coverage appropriately, $dS_0/dt=-\delta(t)$.  

Although MFT and LRR are both recursive as every layer is coupled to
the preceding layer, they differ in that the MFT equations are
nonlinear, while the LRR equations are linear. Nevertheless, when
$d\to 0$, both approximations are identical. Indeed, multiplying
Eq.~(\ref{next}) by $e^t$ and differentiating, one recovers the MFT
equation (\ref{lin}).

Using the Laplace transform, analytical results for the growth
velocity and the interface width have been established
\cite{ran,be,new}.  Below, we give an alternative and simpler
derivation which additionally provides the asymptotic behavior of the
coverage profile.  In the long time limit, it is reasonable to treat
the layer number $j$ as a continuous variable.  Replacing the
difference $S_{j+1}-S_j$ by a partial derivative, the recursive
relations (\ref{tj}) become
\begin{eqnarray}
\label{cov}
{\partial S\over \partial j}&\cong& -\int_0^t d\tau\,S_1(\tau)\,
{\partial S(j,t')\over \partial t'}\bigg|_{t'=t-\tau}\\
&\cong& -{1\over v_d}{\partial S\over \partial t}
+{1\over 2v_d^2}\,{\Gamma\left({d+3\over d+1}\right)\over
\Gamma^2\left({d+2\over d+1}\right)}
{\partial^2 S\over \partial t^2},\nonumber
\end{eqnarray}
with $v_d$ the growth velocity 
\begin{equation}
\label{vdim}
v_d=\left({\Omega_d\over d+1}\right)^{1\over d+1}
\Big /\Gamma\left({d+2\over d+1}\right).
\end{equation}
The second line in (\ref{cov}) has been derived by expanding
$S(j,t-\tau)$ in a Taylor series in $\tau$, keeping only the two
dominant terms of the expansion, and replacing the upper limit in the
integral by $\infty$.  The following growth velocities are found:
$v_0=1$, $v_1=2/\sqrt{\pi}=1.12838$, $v_2=1.13719$, and $v_d\to
\sqrt{2\pi e/d}$ as $d\to\infty$. The velocity is almost constant for
physical dimensions (it varies by less than 4\% in the range $1\le d \le
4$) indicating the weak dimension dependence of this approach.

Changing variables from $(j,t)$ to $(j, \xi=j-v_dt)$ recasts
Eq.~(\ref{cov}) into a diffusion equation
\begin{equation}
{\partial S\over\partial t}=D{\partial^2 S\over\partial \xi^2}. 
\label{diff}
\end{equation}
The constant $D={v_d\over 2}\,\Gamma\left({d+3\over d+1}\right)\Big/
\Gamma^2\left({d+2\over d+1}\right)$ plays the role of a diffusion
coefficient and  controls the width of the interface.  In
obtaining Eq.~(\ref{diff}), $j$ was replaced by $v_dt$.  This is clearly
valid in the scaling limit, $j\to\infty$, $|\xi|\to\infty$, $j\sim \xi^2$.
The initial profile of the uncovered fraction, $S(j,0)$, is a step
function: $S(j,0)=0$ for $j\leq 0$ and $S(j,0)=1$ for $j>0$.  Solving
(\ref{diff}) subject to these initial conditions yields
\begin{equation}
S(j,t)={1\over 2}{\rm Erfc}(-z), \quad 
z={\xi\over\sqrt{4Dt}}={j-v_dt\over \sqrt{4Dt}},
\label{cover}
\end{equation} 
with ${\rm Erfc}(z)={2\over \sqrt{\pi}}\int_z^{\infty} du\, e^{-u^2}$
the error function\cite{Bender}. While the two approximations
generally differ in their asymptotic behavior, they do agree in the
extreme cases of $d=0$ and $d=\infty$.

In summary, the LRR approach predicts a dimension-independent
``diffusive'' width exponent $\beta=1/2$.  The shape of the coverage
profile is symmetric and Gaussian far from the front,
$\sigma_{\pm}=2$.



\section{Comparison with Simulations}

To test the two approximations we simulated the PNG process. The
simulation results presented below are for a one-dimensional chain of
length $L=10^4$ and represent a single realization.  Our work is
different from most previous numerical studies which simulated the PNG
process on a lattice\cite{kr}.  The 1D continuous PNG process is
simple. The initial conditions consist of a flat substrate.  Island
nucleation corresponds to the creation of a kink-antikink pair, while
island merger corresponds to a kink-antikink annihilation. A list of
all distances between kink-antikink pairs moving towards each other is
maintained throughout the simulation. These distances decrease
linearly with time as the motion is ballistic. Thus, the first
collision time is picked, and the annihilating pair is removed from the
system. In parallel, kink-antikink pairs are nucleated with unit rate,
and the distance list is updated accordingly.  This simulation treats
time and space as continuous variables. This allows a direct
calculation of the growth velocity, which should vary from the one
obtained from lattice simulations.  In contrast with previous
simulations\cite{be} of the first few layers, this calculation obtains
the growth velocity to much higher accuracy, as it is orders of
magnitude larger in size as well as temporal range.

\begin{figure}[t]
\vspace{-.24in} 
%\centerline{\epsfxsize=9cm \epsfbox{fig2.eps}}
\centerline{\includegraphics[width=9cm]{fig2}}
\vspace{-.35in} 
\noindent{\small {\bf Fig.2.}  
Uncovered fraction $S_j(t)$ versus $t$ for $j=1,2,3,4$. MFT and
LRR approaches give lower and upper bounds, respectively, for the 
uncovered fraction.}
\end{figure}

\begin{figure}[t]
\vspace{-.24in}
%\centerline{\epsfxsize=9cm \epsfbox{fig3.eps}}
\centerline{\includegraphics[width=9cm]{fig3}}
\vspace{-.35in}
\noindent{\small {\bf Fig.3.}
Short time behavior of the height $h(t)$ and the width $w(t)$. 
MFT is closer to the actual behavior.}
\end{figure}

The time dependence of the uncovered fraction for the first four layers
is shown in Fig.~2.  It is seen that the MFT and the LRR approaches
provide upper and lower bounds, respectively, for the actual PNG
coverage.  Additionally, the MFT provides a better approximation for the
uncovered fraction, $S_j(t)$.  For early times, the height and width
predicted by either approximations are quite close to simulation
results, as shown in Fig.~3.  In fact, both approaches agree to the
first significant order in time, as both Eq.~(\ref{sj}) and
Eq.~(\ref{next}) predict $S_j(t)=1-{2^j\over (2j)!}t^{2j}$. The
disagreement between the two is of the order $t^{2j+2}$. As the two
approximations give upper and lower bounds for the PNG process, we
conclude that this is the leading early time behavior of $S_j(t)$.


However, both approximations become progressively worse at later
times. This is due to the fact that the asymptotic behavior of the width
is predicted incorrectly (see Table 1). As the mean-field theory
predicts an asymptotically flat surface, both the roughness and width
exponents vanish $\chi=\beta=0$.  Fig.~4 shows that at least in one
dimension, the PNG asymptotic behavior belongs to the
Kardar--Parisi--Zhang (KPZ) universality class\cite{kpz}. This is
consistent with the analytical result in 1D that gives $\chi=1/2$.  As
our simulations are continuous in space and time, they enable
measurement of the surface growth velocity.  Although in principle,
there is no reason to expect that the surface growth velocity in an
infinite and a finite system are equal, the numerical velocity $v_{\rm
noneq}=1.41\pm 0.01$, which corresponds to (non-equilibrium) growth on
an infinite substrate, is in very good agreement with the known
analytical equilibrium growth velocity $v_{\rm eq}=\sqrt{2}$\cite{ben}.
Numerical values for $v$, $\beta$ and $\sigma_{\pm}$ are summarized in
Table 1. Two-dimensional simulations give a velocity $v_2\approx 1.4$
\cite{gilmer} compare to the values $v_2=1.671$ (MFT) and $v_2=1.137$
(LRR).  Furthermore, the exponent $\beta$ decreases when the spatial
dimension increases\cite{zhang}.  Thus the MFT improves for higher
dimensions.
\begin{figure}[t]
\vspace{-.35in}
%\centerline{\epsfxsize=9cm \epsfbox{fig4.eps}}
\centerline{\includegraphics[width=9cm]{fig4}}
\vspace{-.25in}
\noindent{\small {\bf Fig.4}
Long time behavior of the width. Early behavior is linear and 
late behavior is in the KPZ universality class $t^{1/3}$.}
\end{figure}

It proves useful to note that the dynamical exponent
for the 1D PNG is $z=3/2$, as it follows from the exponent 
relation $z=\chi/\beta$ together with $\chi=1/2$ and $\beta=1/3$.
The finite size effects are negligible as long as $t\ll L^z$,
so for the system size $L=10^4$ we should guarantee that $t\ll 10^6$.
The time range of our simulations, $t<10^3$, is well within the
above bound.

 
We also examined the extremal behavior of $S_j(t)$ using the
simulations. The scaling prediction (\ref{sigp}) holds as the
simulation data is consistent with the exponent value $\sigma_+=3/2$,
in agreement with the KPZ behavior \cite{zhang1}.  Given that the PNG
uncovered fraction is bounded by the two approximations, the
parameters $v$, $\beta$, and therefore $\sigma_+$ are similarly
bounded (see Table 1.)  Since the scaling argument involves the width,
one can not conclude a priory that the same holds for $\sigma_-$.
Nevertheless, the exponent found in the simulation is quite close to
$2$, or possibly slightly larger. In comparison, numerical studies of
the KPZ equation suggest $\sigma_-\cong2.5$ \cite{kbm}.

\vspace{.1in}
\centerline{
\begin{tabular}{|c|c|c|c|}
\hline &MFT&PNG&LRR\\ 
\hline
$v_1$&1.75735&$1.41\pm0.01$&1.12838\\
$\beta$&0&1/3&1/2\\
$\sigma_+$&1&3/2&2\\
$\sigma_-$&2&$\ge 2$&2\\
\hline
\end{tabular}}
\vspace{.1in}
\noindent{\small {\bf Table 1}: Characteristics of the three 
approaches for the one-dimensional PNG model.}



\section{CONCLUSIONS}
 
We have investigated a continuum model of multi-layer growth, the PNG
model.  We confirmed that the LRR approximation implies
dimension-independent growth exponents.  Moreover, we found that the
full coverage profile approaches a dimension-independent form, thus
implying that this approximation effectively ignores interactions.  It
is therefore not surprising that the LRR gives the width exponent
$\beta=1/2$ which is characteristic for non-interacting systems, like
the growth on a zero-dimensional substrate.  We developed an alternative
self-consistent mean-field approach that provides a better approximation
for the uncovered fraction, and close estimates for the early time
behavior.  Mathematically, the mean-field approach leads to an
interesting system of directed higher order Toda equations.  The
mean-field approximation overestimates layer-layer interactions and
implies a sharp interface, i.e., $\beta=0$.  Additionally, the two
approximate approaches combine to give upper and lower bounds for
statistical properties such as the coverage, the velocity, and the
roughness.

The above mean-field theory should allow computation of space-time
correlation functions and structure functions.  Even more detailed
analysis may be possible for 1D substrates.  A bigger challenge is to
solve the PNG process analytically.  Such a solution will undoubtedly
illuminate the theoretical understanding of non-equilibrium growth.

\vskip 0.1in
\noindent
PLK acknowledges support from NSF and ARO.
\vskip 0.1in
\begin{thebibliography}{99}

\bibitem{god}    
     J.~Krug and H.~Spohn, in 
     {\it Solids Far From Equilibrium}, edited by 
     C.~Godr\`eche (Cambridge University Press, Cambridge, 1992).
\bibitem{meakin}
     P.~Meakin,  {\sl Phys. Rep.} {\bf 235}, 191 (1993).
\bibitem{zhang}
     T.~Halpin-Healy and Y.-C.~Zhang, {\sl Phys. Rep.} {\bf 254}, 
     215 (1995).
\bibitem{barabasi} 
     A.-L.~Barabasi and H.~E.~Stanley, {\it Fractal Concepts in 
     Surface Growth}  (Cambridge University Press, Cambridge, 1995). 
\bibitem{frank}
     F.~C.~Frank, {\sl J. Cryst. Growth} {\bf 22}, 233 (1974).
\bibitem{kas}
     D.~Kashchiev, {\sl J. Cryst. Growth} {\bf 40}, 29 (1977).
\bibitem{gilmer}
     G.~H.~Gilmer, {\sl J. Cryst. Growth} {\bf 49}, 465 (1980).
\bibitem{keller}
     A.~Keller, {\sl Rep. Prog. Phys.} {\bf 31}, 623 (1968). 
\bibitem{bassett}
     D.~C.~Bassett, {\it Principles of Polymer Morphology} 
     (Cambridge University Press, Cambridge, 1995).
\bibitem{Kolm} 
     A.~N.~Kolmogorov, {\sl Bull. Acad. Sci USSR, Phys. Ser.} 
     {\bf 1}, 355 (1937). 
\bibitem{Avrami}
     M.~Avrami, {\sl J. Chem. Phys.} {\bf 7}, 1103 (1939);
     {\bf 8}, 212 (1940); {\bf 9}, 177 (1941). 
\bibitem{Johnson}
     W.~A.~Johnson and P.~A.~Mehl, {\sl Trans. AIMME} {\bf 135}, 
     416 (1939)
\bibitem{sek}
     K.~Sekimoto,  {\sl Physica A} {\bf 135}, 328 (1986).
\bibitem{Sek}
     K.~Sekimoto,  {\sl Int. J. Mod. Phys. B} {\bf 5}, 1843 (1991).
\bibitem{bradley}
     R.~M.~Bradley and P.~N.~Strenski, {\sl Phys. Rev. B} {\bf 40}, 
     8967 (1989). 
\bibitem{ep}
     E.~Ben-Naim and P.~L.~Krapivsky,
     {\sl Phys. Rev. E} {\bf 54}, 3562 (1996). 
\bibitem{Evans} 
     J.~W.~Evans, {\sl Rev. Mod. Phys.} {\bf 65}, 1281 (1993).
\bibitem{be}
     M.~C.~Bartelt and J.~W.~Evans,  {\sl J. Phys. A} {\bf 26}, 
     2743 (1993).
\bibitem{new}
     T.~J.~Newman and A.~Volmer,  {\sl J. Phys. A} {\bf 29}, 
     2285 (1996).
\bibitem{ben}
     C.~H.~Bennett, M.~B\"uttiker, R.~Landauer, and H.~Thomas,  
     {\sl J. Stat. Phys.} {\bf 24}, 419 (1981).
\bibitem{gold}
     N.~Goldenfeld,  {\sl J. Phys. A} {\bf 17}, 
     2807 (1984).
\bibitem{gil}
     W.~van Saarlos and G.~H.~Gilmer,
     {\sl Phys. Rev. B} {\bf 33}, 4927 (1986). 
\bibitem{krug}
     J.~Krug and H.~Spohn,
     {\sl Europhys. Lett.} {\bf 8}, 219 (1989). 
\bibitem{kr}
     J.~Krug, P.~Meakin, and T.~Halpin-Healy, {\sl Phys. Rev. A} 
     {\bf 45}, 638 (1992).
\bibitem{zhang1}
Y.~C.~Zhang, {\sl Phys. Rev. B} {\bf 42}, 4897 (1990).
\bibitem{toda}
     M.~Toda, {\sl J. Phys. Soc. Jpn.} {\bf 22}, 431 (1967). 
\bibitem{br}
     M.~Bramson, {\it Convergence of Solutions of the Kolmogorov
     Equation to Travelling Waves} (American Mathematical Society,
     Providence, R.I., 1983).
\bibitem{mur}
     J.~D.~Murray, {\it Mathematical Biology} 
     (Springer-Verlag, New York, 1989). 
\bibitem{Bender}
     C.~M.~Bender and S.~A.~Orszag, {\it Advanced Mathematical Methods
     for Scientists and Engineers} (McGraw-Hill, Singapore, 1984).
\bibitem{ah}
     R.~D.~Armstrong and J.~A.~Harrison, {\sl
     J. Electrohem. Soc.} {\bf 116}, 328 (1969).
\bibitem{ran}
     S.~K.~Rangarajan, {\sl J. Electroanal. Chem.} {\bf 46}, 125 (1973). 
\bibitem{kpz}
     M.~Kardar, G.~Parisi, and Y.~C.~Zhang, {\sl Phys. Rev. Lett.} 
     {\bf 56}, 889 (1986). 
\bibitem{kbm}J.~M.~Kim, A.~J.~Bray, and M.~A.~Moore, {\sl Phys. Rev. A} 
{\bf 44}, 2345 (1991).

\end{thebibliography}
\end{multicols}
\end{document}

