\documentstyle[graphicx,multicol,aps,epsf]{revtex} \begin{document} \title{Domain Number Distribution in the Nonequilibrium Ising Model} \author{E.~Ben-Naim$\dag$ and P.~L.~Krapivsky$\ddag$} \address{$\dag$Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM 87545} \address{$\ddag$Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215} \maketitle \begin{abstract} We study domain distributions in the one-dimensional Ising model subject to zero-temperature Glauber and Kawasaki dynamics. The survival probability of a domain, $S(t)\sim t^{-\psi}$, and an unreacted domain, $Q_1(t)\sim t^{-\delta}$, are characterized by two independent nontrivial exponents. We develop an independent interval approximation that provides close estimates for many characteristics of the domain length and number distributions including the scaling exponents. \vskip 0.1cm \noindent {PACS numbers: 02.50.Ey, 05.40.+j, 82.20.Mj} \end{abstract} \begin{multicols}{2} \section{INTRODUCTION} The theory of phase ordering kinetics, or domain coarsening, has undergone a rapid development in recent years\cite{Bray}. It has been established that systems quenched from a homogeneous high-temperature disordered state to a low-temperature multi-phase state do not order instantaneously; instead, domains of equilibrium ordered phases form and grow with time as the system approaches local equilibrium on larger and larger scales. Generally, a scale-invariant morphology is developed at late times, and the network of domains is (statistically) independent of time when lengths are rescaled by a single characteristic length scale $L(t)$, the typical domain size. This length scale exhibits an algebraic growth with time, $L(t)\sim t^{\nu}$. However, it was recently realized that additional scaling laws characterized by nontrivial scaling exponents exist in such systems. Examples for such decay modes are the autocorrelation function, $A(L)\sim L^{-\lambda}$\cite{fisher}, and the fraction of the system still frozen in its initial state, $P_0(t)\sim t^{-\theta}$\cite{dbg2,kbr}. The latter ``persistence'' probability has since been investigated theoretically\cite{dbg2,dbg1,dhp,bfk,ms,mbcs,der,lee} and experimentally\cite{ypms} in spin systems, interacting particles systems\cite{kbr,cardy,krl,bhm,mc}, Lotka-Volterra models \cite{lpe,fk}, breath figures growth\cite{bou}, foams\cite{ld}, and even simple diffusion \cite{msbc,w}. Similar to the domain growth exponent, $\nu$, these additional exponents are sensitive to the nonequilibrium dynamics followed by the system, and thus are fundamentally different from their equilibrium counterparts. Precisely how many independent hidden exponents does a coarsening system possesses remains an open question. In this study, we establish that at least in one-dimension, additional exponents describe the survival probability and other more subtle statistical properties of domains. We examine systems with short-range interactions described by a scalar order parameter, namely the 1D T=0 Ising model\cite{Ising} evolving according to nonconserved Glauber dynamics\cite{glauber} and conserved Kawasaki dynamics\cite{kaw}. This paper is organized as follows. We first define the domain number distribution in Sec.~II. In the following section, we review our results for Glauber spin-flip dynamics where we develop and solve analytically an Independent Interval Approximation (IIA) that assumes no correlations between adjacent domains. The IIA predictions compare well with Monte Carlo simulations by giving a correct description of the domain statistics as well as good estimates for the underlying exponents. In Sec.~IV we show that nontrivial exponents underly the zero temperature limit of the 1D Ising model with Kawasaki spin-exchange dynamics as well. The IIA, when carefully modified to conserved dynamics, turns out to be equally useful in this case. Summary and conclusions are given in Sec.~IV. \section{Domain Number Distribution} Although we focus in this study on the Ising model, the statistical properties of domains we are concerned with are relevant to arbitrary coarsening processes in one spatial dimension. For example, we ask, what is the domain survival probability $S(t)$, i.e., the probability that a domain, initially present at the system at time $t=0$, is still present at time $t$ (see Fig.~1). We will present theoretical and numerical evidence supporting an algebraic long time decay of this survival probability, \begin{equation} \label{st} S(t)\sim t^{-\psi}. \end{equation} Such a behavior is robust, as the exponent $\psi$ is not sensitive to the initial state of the system (provided long ranged correlations are absent). Our results will also strongly suggest that the exponent $\psi$ is nontrivial, i.e., it cannot be extracted from so-far known exponents associated with the Ising model. \begin{figure} \vspace{-.25in} %\centerline{\epsfxsize=7cm \epsfbox{fig1.eps}} \centerline{\includegraphics[width=7cm]{fig1}} \noindent {\small {\bf Fig.~1}. Domain motion in the Ising-Glauber model. Surviving domains are marked by $+$, annihilated domains by $-$. The domain number at a later time is also indicated.} \end{figure} In principle, a surviving domain may undergo coalescence with other similar phase domains. Thus, a natural generalization of the domain survival probability is $Q_m(t)$, the density of domains composed of $m$ original domains (see Fig.~1). This quantity satisfies the initial condition $Q_m(0)=\delta_{m,1}$. The total domain density, $N(t)$, is given by $N(t)=\sum_m Q_m(t)$, while the domain survival probability counts {\em initial} domains that have not shrunk and hence contains the density $Q_m(t)$ with weight $m$ \begin{equation} \label{stdef} S(t)=\sum_m m Q_m(t). \end{equation} The average number of domains contained within a surviving domain $\langle m(t)\rangle=S(t)/N(t)$ grows algebraically according to $\langle m(t)\rangle\sim t^{\nu-\psi}$ with $\nu$ the domain decay exponent, $N(t)\sim t^{-\nu}$. If the behavior of $Q_m(t)$ is truly self-similar, it should follow the scaling form \begin{equation} \label{qmscl} Q_m(t)\simeq t^{\psi-2\nu}{\cal Q}(mt^{\psi-\nu}). \end{equation} The scaling function ${\cal Q}(z)$ exhibits the following extremal behavior \begin{equation} \label{qzscl} {\cal Q}(z)\sim\cases{z^{\sigma}&$z\ll 1$,\cr \exp(-\kappa z^{\mu})&$z\gg 1$.} \end{equation} The small argument tail describes domains that contain a very small number of initial domains. In particular, the quantity \begin{equation} \label{qtscl} Q_1(t)\sim t^{-\delta} \end{equation} is of special interest: It gives the density of domains which avoided merging with their neighboring domains up to time $t$. The inequalities $Q_1(t)\leq \sum_m Q_m(t)\leq \sum_m m Q_m(t)$ lead to the bounds $\psi\leq\nu\leq\delta$. Taking into account that at least one surviving domain surrounds a persistent spin gives $P(t)\leq S(t)$, where $P(t)\sim t^{-\theta}$ is the density of persistent spins. Thus we arrive at another upper bound $\psi\leq\theta$ for the exponent $\psi$. These bounds suggest that the domain decay rate is the slowest in the problem. We shall show below that these bounds are strict for the Ising model and more generally for the $q$-state Potts model. However, for the Potts model with $q\to 1$ or $q\to\infty$, and for a few other models\cite{kb} some of these exponents are equal to each other. A useful relation between the scaling exponents can be obtained by substituting $m=1$ in Eq.~(\ref{qmscl}) \begin{equation} \label{sr} \delta-\nu=(\nu-\psi)(1+\sigma). \end{equation} Thus, among the three exponents $\psi$, $\delta$, and $\sigma$, only two are independent. It is well known that under nonconserved (conserved) dynamics $\nu=1/2$ ($\nu=1/3$)\cite{Bray}. Quite obviously, domains disappear when their size vanishes, and therefore the domain size and number distributions are intimately related. Thus, domain survival properties involve the distribution of domains of size $n$ consisting of $m$ original domains at time $t$, denoted by $P_{n,m}(t)$. The aforementioned number distribution is $Q_m(t)=\sum_n P_{n,m}(t)$, and consequently, the domain survival probability is $S(t)=\sum_{n,m} m P_{n,m}(t)$. As will be seen later, studying the joint size-number distribution requires detailed knowledge of the domain size distribution $P_n(t)=\sum_m P_{n,m}(t)$. This distribution obeys the normalization conditions \begin{equation} \label{norm} 1=\sum_n nP_n(t),\qquad N(t)=\sum_n P_n(t). \end{equation} Length conservation implies the first relation, while the second relation gives the total domain density. Since the average domain length grows as $n\sim t^{\nu}$, the length distribution follows the scaling form \begin{equation} \label{pnscl} P_n(t)\simeq t^{-2\nu}{\cal P}(nt^{-\nu}). \end{equation} All of the above scaling behavior emerges from the approximation detailed below. Furthermore, it is satisfied by the simulation data. In the next section, we develop an approximation scheme that helps elucidate many of the qualitative and quantitative features of the domain size and number distributions. \section{Nonconserved Glauber dynamics} We start with the 1D Ising model subject to $T=0$ Glauber dynamics \cite{glauber}. To examine the role of the number of equilibrium phases we also consider a generalization of the Ising model, the $q$-state Potts model. In higher dimensions, the $q$-state Potts model is relevant to physical situations for $q=2$ (the Ising model) and additionally for $q=3,4,\infty$\cite{wu}. For instance, the $q=\infty$ case describes several cellular structures\cite{st} including polycrystals\cite{poly}, foams\cite{ld}, soap froth\cite{glaz}, and magnetic bubbles\cite{bubbl}. We consider uncorrelated initial conditions where each of the $q$ phases is present with equal density $1/q$. The $T=0$ Glauber-Potts dynamics proceeds by selecting a spin at random and changing its value to that of one of its randomly selected neighbors. Thus, domain walls perform a random walk and upon contact, they annihilate or coalesce, depending on the state of the corresponding domains \cite{bbd,racz,amar}. Identifying a domain wall with a particle, ($A$), and absence of a domain wall with a hole ($0$), one finds the single-species diffusion-reaction process \begin{equation} \label{aa} A0\buildrel {1\over 2} \over \rightleftharpoons 0A,\quad AA\buildrel {1\over q-1} \over \longrightarrow 00,\quad AA\buildrel {q-2\over q-1}\over \longrightarrow A0\ {\rm or }\ 0A. \end{equation} The rates indicate the relative probabilities by which each event occurs. \subsection{Domain Size Distribution} Ignoring correlations between neighboring domains allows us to develop an approximate theory for the time-evolution of the domain distribution. Approximations that are similar in nature proved useful in studies of related reaction-diffusion processes\cite{msbc,bbd,ab}. The joint number distribution requires knowledge of the length distribution and we start by deriving a master equation for $P_n(t)$. Under the assumption that the lengths of neighboring intervals are uncorrelated, we write the following rate equation\cite{kb} \begin{eqnarray} \label{pn} {dP_n\over dt}=&&P_{n-1}+P_{n+1}-2P_n\\ +&&{P_1\over (q-1)N^2}\left[\sum_{i=1}^{n-2} P_i P_{n-1-i} -N(P_n+P_{n-1})\right].\nonumber \end{eqnarray} with $N(t)=\sum_{n}P_n(t)$ the total domain density and the boundary condition $P_0(t)=0$. The first three terms reflect that domain walls perform a random walk with hopping rate set to 1/2 without loss of generality. The last two terms are due to domain annihilation: the convolution term accounts for domain merger and the last term for domain loss. In the $q$-state Potts model collision of domain walls results in annihilation with probability ${1\over q-1}$ or in coalescence with probability ${q-2\over q-1}$. Only annihilation events affect the domain distribution and thus the ${1\over q-1}$ prefactor of the annihilation terms. Using the sum rules of Eq.~(\ref{norm}), one verifies that the total length is conserved and the total domain density decays according to the {\it exact} rate equation \begin{equation} \label{N} {dN\over dt}=-{q\over q-1}P_1. \end{equation} The diffusion term in Eq.~(\ref{pn}) implies $\langle n(t)\rangle\sim t^{1/2}$, and since $\langle n\rangle\sim N^{-1}$ the correct decay exponent $\nu=1/2$ \cite{glauber} is recovered. In the following, we will need to determine the asymptotic prefactor $A$, $N(t) \simeq At^{-1/2}$, $A=\int dx\,{\cal P}(x)$, with the scaling function ${\cal P}(x)$ defined according to Eq.~(\ref{pnscl}). The density rate equation (\ref{N}) implies $P_1\simeq {\cal P}'(0)t^{-3/2}$ with ${\cal P}'(0)={q-1\over 2q}A$. A quantitative analysis of Eq.~(\ref{pn}) may be carried by treating the variable $n$ as continuous. The quantity ${\cal P}(x)$ satisfies \begin{equation} \label{px} {\cal P}''+{1\over 2}(x{\cal P})'+{q-2\over 2q}{\cal P} +{1\over 2qA}\,{\cal P}*{\cal P}=0, \end{equation} where ${\cal P}'\equiv d{\cal P}/dx$ and ${\cal P}*{\cal P}\equiv\int_0^x dy\,{\cal P}(y){\cal P}(x-y)$. The normalized Laplace transform of the scaling function ${\cal P}(x)$, $p(s)=A^{-1}\int_0^{\infty} dx\,e^{-sx}\,{\cal P}(x)$, obeys \begin{equation} \label{riccati} {dp\over ds}={p^2\over qs}+\left(2s+{q-2\over qs}\right)p -{q-1\over qs}, \end{equation} subject to the boundary condition $p(0)=1$. The transformation $p(s)=1-qs^2-qs{d\over ds}\ln y(s)$ reduces the Riccati equation (\ref{riccati}) into the parabolic cylinder equation, \begin{equation} \label{par} {d^2 y\over ds^2}+\left(1+{2\over q}-s^2\right)y=0. \end{equation} The solution to (\ref{par}) reads $y(s)=C_-D_{1/q}(-s\sqrt{2})+C_ +D_{1/q}(s\sqrt{2})$, with $D_{1/q}(x)$ the parabolic cylinder function of order $1/q$ \cite{bo}. The large $s$ behavior of $p(s)$, $p(s)\simeq {q-1\over 2q}s^{-2}$, implies $C_-=0$, and we get \begin{equation} \label{ps} p(s)=1-qs^2-qs{d\over ds}\ln D_{1/q}(s\sqrt{2}). \end{equation} The normalization condition $\sum_n nP_n(t)=1$ can be reduced to $Ap'(0)=-1$. This allows us to determine the constant \begin{equation} \label{A} A={\Gamma\big[1-{1\over 2q}] \over\Gamma\big[{1\over 2}-{1\over 2q}\big]}, \end{equation} where $\Gamma$ denotes the gamma function. In deriving (\ref{A}) we have used the properties \cite{bo} \begin{equation} \label{large} D_c(x)\sim x^c\exp(-x^2/4)[1+{\cal O}(x^{-2})], \end{equation} and \begin{equation} D_c(0)={\pi^{1/2}2^{c/2}\over \Gamma(1/2-c/2)},\quad D_c'(0)=-{\pi^{1/2}2^{(c+1)/2}\over \Gamma(-c/2)}. \end{equation} The value of the constant $A$ predicted by the IIA may be compared to the exact one, $A_{\rm exact} =(1-q^{-1})/\sqrt{\pi}$ \cite{amar} In the extreme cases of $q=1$ and $q=\infty$ the prefactor $A$ is exact. The mismatch is worst (roughly 20\%) for the Ising ($q=2$) case where $A=\Gamma(3/4)/\Gamma(1/4)\cong 0.337989$ while $A_{\rm exact} =(4\pi)^{-1/2}\cong 0.28209$\cite{glauber}. The IIA predicts the correct qualitative behavior of length distribution in the limits of small and large intervals \begin{equation} \label{pxscl} {\cal P}(x)\sim\cases{{A(q-1)\over 2q}x&$x\ll 1$,\cr qA \exp(-\lambda x)&$x\gg 1$.} \end{equation} The linear small size behavior is seen from the large $s$ behavior $p(s)\simeq {(q-1)\over 2q}s^{-2}$. On the other hand, the exponential tail follows from the behavior of the Laplace transform $p(s)\simeq q\lambda/(s+\lambda)$ near its pole at negative $s=-\lambda$, given by the first zero of $D_{1/q}(-\lambda\sqrt{2})=0$. For the Ising case one has $\lambda=0.5409$. This value should be compared with the exact value $\lambda=\zeta(3/2)/4\sqrt{\pi} =0.368468$ obtained by Derrida and Zeitak\cite{dz} and the approximate value $\lambda=0.35783$ obtained by Alemany and ben-Avraham\cite{ab}. \subsection{Domain Size-Number Distribution} We are now in a position to tackle the joint size-number distribution, $P_{n,m}(t)$, which captures both the spatial and ``historical'' characteristics of the coarsening domain mosaic. The corresponding rate equation is a generalization of Eq.~(\ref{pn}) \begin{eqnarray} \label{pnm} &&{dP_{n,m}\over dt}=P_{n-1,m}+P_{n+1,m}-2P_{n,m}\\ &&+{P_1\over (q-1)N^{2}}\left[\sum_{i,j}P_{i,j}P_{n-1-i,m-j} -N(P_{n,m}+P_{n-1,m})\right] \nonumber \end{eqnarray} with the initial condition $P_{n,m}(0)=\delta_{n,1}\delta_{m,1}$ and the boundary condition $P_{0,m}(t)=0$. The variable $m$ is almost mute as it appears in a nontrivial way only in the convolution term. One should verify that this master equation is self-consistent. First, by summing over $m$, we recover Eq.~(\ref{pn}). Second, it implies that the domain survival probability satisfies the exact linear equation ${dS/dt}=-\sum_m m P_{1,m}$. We have not succeeded in solving for the joint distribution. Nevertheless, it is possible to obtain analytically many interesting properties of Eqs.~(\ref{pnm}), including the scaling exponents. Given Eqs.~(\ref{pnm}) are recursive in $m$, one can try to solve for $P_{n,1}(t)$, then for $P_{n,2}(t)$, etc. A solution for the former quantity already allows to determine the scaling exponent $\delta$. Thus let us consider the distribution of domains which have not merged with other domains up to $t$, $R_n(t)\equiv P_{n,1}(t)$. For such domains, the convolution term vanishes and they evolve according to the linear rate equation \begin{equation} \label{rn} {dR_n\over dt}=R_{n-1}+R_{n+1}-2R_n-{P_1\over (q-1)N}(R_n+R_{n-1}) \end{equation} with the initial condition $R_n(0)=\delta_{n,1}$ and the boundary condition $R_0(t)=0$. In the continuum limit we again replace $R_{n-1}+R_{n+1}-2R_n$ by $\partial^2 R/\partial n^2$ and $R_n+R_{n-1}$ by $2R_n$ to find a diffusion-convection equation for $R_n(t)$. The transformation $R_n\to \tilde R_n N^{2/q}$ reduces this equation to the diffusion equation for $\tilde R_n$, which is solved to yield $R_n(t)\simeq N^{2/q}t^{-1}{\cal R}(nt^{-1/2})$, with ${\cal R}(x)=x\exp(-x^2/4)/\sqrt{\pi}$. The total density of unreacted domains is $Q_1(t)=\sum_n R_n\sim t^{-{1\over 2}-{1\over q}}$, which gives the decay exponent \begin{equation} \label{delta} \delta={1\over 2}+{1\over q}. \end{equation} Obtaining the second independent exponent $\psi$ is more involved. The natural approach, i.e., a direct investigation of the domain number distribution $Q_m$, appears to be useless, as it requires knowledge of $P_{1,m}$ and hence the entire $P_{n,m}$. The domain survival probability can be alternatively obtained by considering $U_n(t)=\sum_m mP_{n,m}(t)$. This quantity obeys \begin{eqnarray} \label{un} {dU_n\over dt}&=&U_{n-1}+U_{n+1}-2U_n\\ &+&{P_1\over (q-1)N^2} \left[2\sum_{i=1}^{n-2} U_i P_{n-1-i}-N(U_n+U_{n-1})\right],\nonumber \end{eqnarray} obtained by summing Eqs.~(\ref{pnm}). We write $U_n(t)$ in a scaling form $U_n(t)\simeq t^{-\psi-1/2}{\cal U}(nt^{-1/2})$. Asymptotically, the domain survival probability reads $S(t)\simeq Bt^{-\psi}$ with $B=\int dx\, {\cal U}(x)$. The scaling distribution satisfies \begin{equation} \label{ux} {\cal U}''+{1\over 2}(x{\cal U})'+ \left(\psi-{1\over q}\right){\cal U} +{1\over qA}{\cal U}*{\cal P}=0. \end{equation} The normalized Laplace transform of the scaling function ${\cal U}(x)$, $u(s)=B^{-1}\int_0^\infty dx\,e^{-sx}{\cal U}(x)$, obeys \begin{equation} \label{dus} {du\over ds}=2\left({p(s)+q\psi-1\over qs}+s\right)u -{2\psi\over s}, \end{equation} and $u(0)=1$. In deriving (\ref{dus}) we used the relation ${\cal U}'(0)=B\psi$, found by integration of Eq.~(\ref{ux}), combined with $A=\int dx {\cal P}(x)$. Substituting the explicit expression (\ref{ps}) for $p(s)$ into Eq.~(\ref{dus}), and solving for $u(s)$ yields \begin{equation} \label{us} u(s)=2\psi s^{2\psi}D^{-2}_{1/q}(s\sqrt{2}) \int_s^{\infty}\!\!dr\, r^{-2\psi-1}D_{1/q}^{2}(r\sqrt{2}). \end{equation} This solution is consistent with the anticipated $s\to\infty$ behavior, $u(s)\simeq \psi s^{-2}$. Furthermore, evaluating Eq.~(\ref{us}) near the origin gives $u(s)=1+F(\psi)s^{2\psi}+Cs+\cdots$. Therefore, for $u'(s)$ to be finite near $s=0$, we must have $F(\psi)=0$. Evaluating $F(\psi)$ gives \begin{equation} \label{psi} 0=\int_0^\infty dr\,r^{-2\psi}D_{1/q}(r)D_{1/q}'(r), \end{equation} an eigenvalue problem that can be solved numerically to obtain the exponent $\psi$ (see Table 1). In the most interesting case of integer $q$ the domain decay exponent $\psi$ appears to be irrational, in contrast with $\delta$. It is useful to consider the limiting cases that turn out to be solvable. The $q=\infty$ limit is especially simple\cite{bbd} as only domain walls coalesce but cannot annihilate and thereore Eq.~(\ref{pnm}) is linear and thus exact. Furthermore, the domain size number distribution factorizes, $P_{n,m}=P_n(t)\delta_{m,1}$, since similar phase domains never coalesce and therefore the domain number is trivial, $m=1$. Thus $N(t)=S(t)=Q_1(t)\simeq (\pi t)^{-1/2}$ and $\nu=\psi=\delta=1/2$. Additionally, the scaling function is ${\cal P}(x)=x\exp(-x^2/4)/\sqrt{\pi}$. Before going to the opposite limit $q\to 1$, we first note that the Potts model with arbitrary $q\geq 1$ can be mapped onto the Ising model with magnetization $\mu=2/q-1$. Thus the $q\to 1$ limit corresponds to the vanishing volume fraction of minority domains. Therefore minority domains cannot ``meet'', so majority domains change appreciably only due to coalescence. Thus they never disappear, i.e., $S(t)=1$ and $\psi=0$. A majority domain remains unreacted if both its minority neighbors survive, implying $Q_1(t)=N^2(t)$ and $\delta=2\nu=1$. In the above rate equation description, the diffusion term becomes negligible and the IIA is exact. The scaling functions ${\cal P}(x)$ and ${\cal Q}(z)$ are identical exponential functions. However, unlike to the $q=\infty$ case, the joint distribution is not a product of the single function variables\cite{kb}. The joint size-number distribution still obeys the scaling law, $P_{n,m}(t)\sim t^{-5/4}\Phi(x,y)$, with $\Phi(x,y)=x^{-1/2}\exp(-x-y^2/2x)$. The scaling variables $x$ and $y$ are quite unusual, though\cite{aggr,kb}. Indeed, instead of the naive scaling variables $mt^{-1/2}$ and $nt^{-1/2}$, one has $x=(m+n)(\pi t)^{-1/2}$ and $y=(m-n)(\pi t)^{-1/4}$\cite{kb}. The former scaling variable $x$ is just the sum of the naive scaling variables, while the latter ``diffusive'' scale $y$ is hidden. This suggests that generally for the $q$-state Potts model with $q<\infty$ the joint distribution is not necessarily the product of single variable functions and the scaling, if holds, may be rather different from the naive form with scaling variables $n t^{-\nu}$ and $m t^{\psi-\nu}$. \subsection{Simulation Results} To test the IIA predictions, we performed numerical simulations on a spin chain of size $L=10^7$. Random initial conditions and periodic boundary conditions were used. The simulation data represents an average over 10 different realizations. For the Ising case, we found the exponent values $\psi=0.126(1)$ and $\delta=1.27(2)$ (see Fig. 2). These values should be compared with the IIA predictions of $\psi=0.136612$ and $\delta=1$. The IIA neglects correlations that do build up between neighboring domains and thus is not exact. Furthermore, the effects of the correlations is nontrivial, as one exponent is smaller than predicted while the other is larger. As was the case for the persistence exponent, $\theta$, the domain exponents strongly depend on $q$. Numerical values of the exponents $\psi$ and $\delta$ are summarized in Table 1 for representative values of $q$. As $q$ increases, the approximation improves and eventually becomes exact for the extreme case $q=\infty$. Thus, $\psi$ is overestimated by up to 10\% and $\delta$ is underestimated by up to 25\%. We performed several checks to verify that the asymptotic behaviors of Eqs.~(\ref{st}) and (\ref{qtscl}) are robust. For example, they are independent of the initial domain wall concentration (provided that the correlations in the initial condition are short range). We conclude that $\psi$ and $\delta$ are nontrivial exponent, i.e., they cannot be extracted from the known exponents associated with the Ising-Glauber model. Similar to the persistence exponent, $\theta(q)$, the exponents appear to be irrational except for the limiting cases $q=\infty$ ($\psi=\delta={1\over 2}$, $\sigma=0$) and, maybe, for $q=2$ ($\psi={1\over 8}$, $\delta={5\over 4}$, $\sigma=1$). The numerical simulations also confirm that the distribution function $Q_m(t)$ scales according to Eq.~(\ref {qmscl}). The scaling function ${\cal Q}(z)$, defined in Eq.~(\ref{qzscl}), decays exponentially for large argument ($\mu=1$) and is algebraic for small argument. The scaling relations combined with the simulation values give $\sigma=1.05(5)$. This is consistent with the linear behavior seen in our simulations for $z\ll 1$. Comparing with Eq.~(\ref{pxscl}), we conclude that similar scaling functions underlie the domain number and size distributions in the $q=2$ case. \begin{figure} \vspace{-.40in} \centerline{\includegraphics[width=9cm]{fig2}} %\centerline{\epsfxsize=9cm \epsfbox{fig2.eps}} \vspace{-.40in} \noindent {\small {\bf Fig.~2}. Monte Carlo data for the Ising-Glauber model. The domain survival probability $S(t)$, the domain density $N(t)$, and the density of unreacted domains $Q_1(t)$ are shown (top to bottom). The inset plots the local slope $-d\ln S(t)/d\ln t$. Typically, it is stable over a large temporal range, and thus can be used to find the scaling exponents and to estimate the error, typically of the order $0.001$.} \end{figure} On the other hand, direct numerical integration of Eq.~(\ref{pnm}) reveals a number distribution, $Q_m(t)$, that scales according to Eq.~(\ref{qmscl}), and has an exponential tail in agreement with the simulation results. Moreover, the emerging $S(t)$ falls within $5\%$ of the actual survival probability over a significant temporal range, $t<10^3$. In summary, in addition to predicting the correct scaling behavior, Eq.~(\ref{pnm}) provides a good approximation for many quantitative features of the domain distribution, and in particular, good estimates for the decay exponents. \vspace{.1in} \centerline{ \begin{tabular}{|l|c|c|c|c|c|} \hline &\multicolumn{3}{c|}{MC}&\multicolumn{2}{c|}{Eq.~(\ref{pnm})}\\ \hline $q$&$\psi$&$\delta$&$\sigma$&$\psi$& $\delta$\\ \hline 2 &\,0.126\,&\,1.27\,&\,1.05\,&\,0.136612\,&\,1\,\\ 3 &\,0.213\,&\,0.98\,&\,0.67\,&\,0.231139\,&\,5/6\,\\ 4 &\,0.267\,&\,0.85\,&\,0.50\,&\,0.287602\,&\,3/4\,\\ 8 &\,0.367\,&\,0.665\,&\,0.24\,&\,0.385019\,&\,5/8\,\\ 50 &\,0.476\,&\,0.525\,&\,0.03\,&0.480274 &13/25\\ $\infty$ &\, 1/2\,&\,1/2 \,&\,0\,&\, 1/2\,&1/2\\ \hline \end{tabular}} \vspace{.1in} \noindent{\small {\bf Table 1}: Domain exponents for the $q$-state Potts model in one dimension. Local slopes analysis was applied to the simulation data. The theoretical $\psi$ is from Eq.~(\ref{psi}) and $\delta={1\over 2}+{1\over q}$.} \section{Conserved Kawasaki dynamics} We turn now to applying the above methods to the conserved counterpart, the spin-exchange Kawasaki dynamics\cite{kaw}, which describe spinodal decomposition in binary alloys and phase separation in binary liquids. Although some qualitative features are known\cite{tob,cor,huse}, theoretical understanding of the Ising-Kawasaki model is still far from complete even in one dimension. In the following we limit ourselves to the two-phase Ising case. \subsection{Reduction to Domain Diffusion} We start by formulating the appropriate zero-temperature limit of the Ising-Kawasaki model. Consider a two-phase system, e.g., the Ising model (spins up and down) or a binary alloy (atoms of type $A$ and $B$). At zero temperature energy raising transitions are forbidden and only two moves are allowed: the energy-decreasing ``coarsening'' transitions $ABAB\to AABB$ and the energy-conserving ``diffusion'' transitions $ABAA\to AABA$. This dynamics ultimately drives the system to a frozen configuration consisting of strings of alternating domains each of length $\geq 2$ which could not evolve further\cite{zero}. This ``jamming'' behavior arises from the nonergodic nature of the zero-temperature Kawasaki dynamics and it is a robust one: It is independent of the relative transition rates \cite{cor,anal}, as well as the spatial dimension \cite{zero,inf,bin}. Thus, to recover sensible {\em coarsening}, one must consider the zero-temperature limit. Let us assume that temperature is positive; when it is sufficiently low, $T\ll J$ where $J$ is the exchange coupling, the correlation length $\xi\sim e^{J/T}$ is very large and therefore the system exhibits coarsening as long as the mean domain size is small compared with the correlation length. Below, we focus on this intermediate-time regime where the description is drastically simplified\cite{tob,cor}. We assume that the initial stage has been already completed so that single-spin domains disappeared. Coarsening will occur only when a spin splits off a domain wall and penetrates a neighboring domain (say of size $L$). The splitting process occurs with a very small rate $\exp[-4J/T]$. Then, this spin diffuses inside the domain until it is eventually adsorbed by its boundaries. The corresponding probabilities are well-known from elementary probability theory\cite{feller}. The spin will be absorbed by the boundary from which it was issued with probability $1-1/L$. This spin may also be absorbed by the opposite boundary resulting in a one lattice site hop of the entire domain. Thus, the hopping rate is $L^{-1}\exp[-4J/T]$. Rescaling time, $t\to t\exp[-4J/T]$, the spin diffusion will proceed with a huge rate $\exp[4J/T]$ and it therefore may be treated as instantaneous while domain hops proceed with a finite rate reciprocal to the domain size. Thus, the appropriate zero-temperature limit of the Ising-Kawasaki model is realized by taking the limits of infinite ``physical'' time $t_{\rm phys}\to\infty$, while keeping the modified time $t=t_{\rm phys}\exp[-4J/T]$ finite. Hence, entire domains perform a random walk with rate inversely proportional to their lengths (we ignore an anomaly concerning domains of length 2 as it is irrelevant asymptotically). Heuristically, it may be argued that as diffusion is the primary coarsening mechanism, the following scaling for the average domain size, $L\sim \sqrt{Dt}$, holds. However, since the diffusion coefficient and the domain size are reciprocal, $D\sim L^{-1}$, we obtain $L\sim t^{1/3}$ in agreement with the well-known behavior of systems with conserved scalar order parameter\cite{Bray}. \subsection{Domain Size Distribution} For simplicity, we consider the case where the two equilibrium phases are equivalent, as is the case for random initial conditions. Modifying Eq.~(\ref{pn}) to account for domain diffusion, and assuming neighboring domains are uncorrelated, the domain size distribution evolves according to \begin{eqnarray} \label{pn1} {dP_n\over dt}&=&L^{-1}(P_{n-1}-2P_n+P_{n+1})\\ &+&{P_1\over N^2}\left[\sum_{i+j=n}i^{-1}P_i P_j -N(n^{-1}+L^{-1})P_n\right],\nonumber \end{eqnarray} with the domain density $N$ and the inverse average domain size $L^{-1}$ defined via \begin{equation} \label{nln} N=\sum_{n=1}^\infty P_n, \quad L^{-1}=\langle n^{-1}\rangle ={\sum_{n} n^{-1}P_n\over \sum_{n} P_n}. \end{equation} The diffusion term in Eq.~(\ref{pn1}) accounts for change in $n$ due to hoping of a neighboring interval. The convolution term accounts for gain due to domain merger, and the last two terms represent loss due to domain collision or domain merger. This rate equation conserves the total length, $\sum_n nP_n=1$, and summation over $n$ shows that the total domain density evolves according to $\dot N=-2L^{-1}P_1$. The length distribution scales according to Eq.~(\ref{pnscl}), $P_n(t)\simeq t^{-2\nu}{\cal P}(nt^{-\nu})$, with the correct scaling exponent $\nu=1/3$. By inserting that scaling form into Eq.~(\ref{pn1}) we arrive at an integro-differential equation for ${\cal P}(x)$ which is very cumbersome as it involves yet unknown moments of the distribution. We thus resort to numerical integration of Eq.~(\ref{pn1}). The results compare well with Monte-Carlo simulations of the Ising-Kawasaki model (see Fig.3). For example the estimate for the asymptotic prefactor $A$ (defined via $N(t)\simeq At^{-1/3}$), falls within roughly 5\% of the actual value: $A_{\rm MC}=0.441\pm0.001$ while $A_{\rm IIA}=0.415\pm0.005$. \begin{figure} \vspace{-.30in} %\centerline{\epsfxsize=9cm \epsfbox{fig3.eps}} \centerline{\includegraphics[width=9cm]{fig3}} \vspace{-.30in} \noindent {\small {\bf Fig.~3}. Domain density in the Ising-Kawasaki model. The Monte Carlo simulation data (MC) represents an average over 10 systems of size $10^5$. The IIA was obtained by integrating Eq.~(\ref{pn1}) numerically.} \end{figure} The length distribution emerging from the IIA has the same limiting behavior as in the nonconserved case, i.e., it is linear at small $n$ and exponential at large $n$. While the former agrees with our simulation results, there is a disagreement for the latter. Our data is consistent with a Gaussian tail, i.e., ${\cal P}(x)\sim\exp(-x^2)$, for $x\gg 1$. \subsection{Domain Size-Number Distribution} Given the results in the conserved case, it is natural to study the domain exponents and to examine the usefulness of the IIA approach in the conserved dynamics case. In analogy with Eq.~(\ref{pnm}), the master equation for the domain size-number distribution is written \begin{eqnarray} \label{pnm1} &&{dP_{n,m}\over dt}=L^{-1}(P_{n-1,m}-2P_{n,m}+P_{n+1,m})\\ &+&{P_1\over N^2}\left[\sum_{i+j=n}\sum_{k+l=m}i^{-1}P_{i,k}P_{j,l} -N(n^{-1}+L^{-1})P_{n,m}\right].\nonumber \end{eqnarray} Summing the above equations over $m$, we indeed recover Eq.~(\ref{pn1}) for the length distribution. To determine the exponents, it is again simpler to consider the distributions $R_n(t)$ and $U_n(t)$ instead of the joint distribution $P_{n,m}(t)$. The density of single parent domains, $R_n(t)\equiv P_{n,1}(t)$, evolves according to the linear rate equation similar to Eq.~(\ref{rn}) \begin{eqnarray} \label{rn1} {dR_n\over dt}&=&L^{-1}(R_{n+1}+R_{n-1}-2R_n)\\ &&-{P_1\over N}(n^{-1}+L^{-1})R_n.\nonumber \end{eqnarray} We expect that the distribution $R_n(t)$ scales according to $R_n(t)\sim t^{-\delta-1/3}{\cal R}(nt^{-1/3})$. Integrating this equation we get $Q_1(t)=\sum R_n(t)\sim t^{-\delta}$ with $\delta\cong 0.645$. On the other hand, Monte-Carlo simulations of the domain diffusion process give $\delta\cong 0.705$ (see Fig.~4). The domain survival probability can be found by using the auxiliary function $U_n(t)=\sum_m mP_{n,m}(t)$ which satisfies the analog of Eq.~(\ref{un}) \begin{eqnarray} \label{un1} {dU_n\over dt}&=&L^{-1}(U_{n+1}+U_{n-1}-2U_n)\\ &&+{P_1\over N^2} \left[2\sum_{i+j=n}i^{-1}P_iU_j-N(n^{-1}+L^{-1})U_n\right].\nonumber \end{eqnarray} This distribution $U_n(t)$ should scale according to $U_n(t)\sim t^{-\psi-1/3}{\cal U}(nt^{-1/3})$. The domain survival probability $S(t)=\sum_n U_n(t)$ then decays according to $S(t)\sim t^{-\psi}$. \begin{figure} \vspace{-.30in} %\centerline{\epsfxsize=9cm \epsfbox{fig4.eps}} \centerline{\includegraphics[width=9cm]{fig4}} \vspace{-.30in} \noindent {\small {\bf Fig.~4}. Monte Carlo data for the Ising-Kawasaki model from the same simulation as in Fig.~(3). The domain survival probability $S(t)$, the domain density $N(t)$, and the density of unreacted domains $Q_1(t)$ are shown (top to bottom).} \end{figure} Again, the agreement with the simulation is remarkable. Numerical integration data give an estimate of $\psi\cong 0.147$ while Monte-Carlo simulations (see Fig.~4)) give $\psi\cong 0.130$. We also verified that the scaling relations of Eqs.~(\ref{qmscl}) and (\ref{qzscl}) are satisfied by the IIA as well as the simulation data. We conclude that in the Kawasaki case as well, nontrivial exponents characterize domain statistics. Furthermore, the approximate approach reproduces most qualitative features of the domain size and number distribution, and provides good estimates for the scaling exponents. We now describe how to obtain $\psi$ in the limit where one of the two phases occupies a vanishing volume fraction. Denote by $L_A(t)$ and $L_B(t)$ the average sizes of minority and majority domains, respectively. They both grow as $t^{1/3}$ but remain greatly different throughout the evolution, $L_A\ll L_B$. The domains diffusion rates, $D_A\sim L_A^{-1}$ and $D_B\sim L_B^{-1}$, thus greatly differ as well: $D_A\gg D_B$. In principle, two neighboring minority domains can overtake a separation distance of order $L_B$ and coalesce; this requires the coalescence time $t_c\sim L_B^2/D_A\sim L_AL_B^2$. On the other case, a minority domain can shrink due to diffusion of neighboring majority domains; this requires the shrinking time $t_s\sim L_A^2/D_B\sim L_A^2L_B$. We see that $t_s\ll t_c$, so we should just take into account the shrinking of minority domains. Thus majority domains do not disappear, implying $S(t)=1$ and $\psi=0$. For the minority phase we anticipate $\psi=\nu=1/3$, while for the symmetric case of equal concentrations $\psi\cong 0.130$. This indicates that similar to the nonconserved dynamics case, the domain exponents vary continuously as the volume fraction is varied. \section{Conclusions} In summary, we investigated the one-dimensional Ising model subject to zero temperature Glauber and Kawasaki dynamics. We introduced the domain size-number distribution and showed that it obeys scaling and is characterized by two independent nontrivial decay exponents. Similar to the persistence exponent, these exponents are sensitive to the type of the dynamics and the volume fraction of the (globally or locally) conserved equilibrium phase. We also introduced an approximation which is based on terminating the hierarchy of rate equations describing the domain density. This approximation is very useful in predicting the qualitative nature of the domain distribution as well as estimating important parameters including the scaling exponents. In the proper T=0 limit of the Ising model with Kawasaki dynamics, this approximation is especially useful as very little is known analytically about the domain distribution. It will be interesting to generalize the domain survival concept to higher dimensions. At least for the $q\to\infty$ limit of the Potts model, domains are well defined, and such a generalization is possible. The nonconserved dynamics can indeed be studied in this limit, and the domain exponents $\psi=\theta=\delta=d/2$ (for $d\geq 2$) have been reported\cite{kb}, consistent with simulations\cite{ld} and with experiments on $d=2$ soap froths\cite{tam}. Recently, it was pointed out that coarsening mosaics may be characterized by more than one algebraically growing length scale, and that morphologies consisting of domain and super-domains may exist \cite{lpe}. This, together with the above results suggest that our current understanding of such systems is only partial. Domain statistics indicates that several nontrivial decay laws underlie the evolution of elementary processes such as the nonequilibrium Ising model. These nontrivial exponents do not emerge naturally from studies of traditional quantities such as spatiotemporal correlations. It remains a challenge to find and obtain these underlying ``hidden'' exponents from a more systematic method. It is also intriguing whether an entire hierarchy or a finite number of independent decay modes are present in these systems. \bigskip\noindent This paper is dedicated to Leo Kadanoff on the occasion of his 60th birthday. The research of PLK has been supported by grants from NSF and ARO. \begin{thebibliography}{99} \bibitem{Bray} A.~J.~Bray, {\sl Adv. Phys.} {\bf 43}, 357 (1994). \bibitem{fisher} D.~S.~Fisher and D.~A.~Huse, {\sl Phys.~Rev.~B} {\bf 38}, 373 (1988). \bibitem{dbg2} B.~Derrida, A.~J.~Bray, and C.~Godr\`eche, {\sl J.~Phys.~A} {\bf 27}, L357 (1994). \bibitem{kbr} P.~L.~Krapivsky, E.~Ben-Naim, and S.~Redner, {\sl Phys. Rev. E} {\bf 50}, 2474 (1994). \bibitem{dbg1} A.~J.~Bray, B.~Derrida, and C.~Godr\`eche, {\sl Europhys. Lett.} {\bf 27}, 175 (1994). \bibitem{dhp} B.~Derrida, V.~Hakim, and V.~Pasquier, {\sl Phys. Rev. Lett.} {\bf 75}, 751 (1995); {\sl J. Stat. Phys.} {\bf 85}, 763 (1996). \bibitem{bfk} E.~Ben-Naim, L.~Frachebourg, and P.~L.~Krapivsky, {\sl Phys. Rev. E} {\bf 53}, 3078 (1996). \bibitem{ms} S.~N.~Majumdar and C.~Sire, {\sl Phys. Rev. Lett.} {\bf 77}, 1420 (1996). \bibitem{mbcs} S.~N.~Majumdar, A.~J.~Bray, C.~Cornell, and C.~Sire, {\sl Phys. Rev. Lett.} {\bf 77}, 3704 (1996). \bibitem{der} B.~Derrida, {\sl Phys. Rev. E} {\bf 55}, 3705 (1997). \bibitem{lee} B.~P.~Lee and A.~D.~Rutenberg, {\sl Phys. Rev. Lett.} {\bf 79}, 4842 (1997). \bibitem{ypms} B.~Yurke, A.~N.~Pargellis, S.~N.~Majumdar, and C.~Sire, {\sl Phys. Rev. E} {\bf 56}, 40 (1997). \bibitem{cardy} J.~Cardy, {\sl J.~Phys.~A} {\bf 28}, L19 (1995). \bibitem{krl} P.~L.~Krapivsky, S.~Redner, and F.~Leyvraz, {\sl Phys. Rev. E} {\bf 51}, 3977 (1995). \bibitem{bhm} E.~Ben-Naim, {\sl Phys. Rev. E} {\bf 53}, 1566 (1996); M.~Howard, {\sl J.~Phys.~A} {\bf 29}, 3437 (1996); C.~Monthus, {\sl Phys. Rev. E} {\bf 54}, 4844 (1996). \bibitem{mc} S.~N.~Majumdar and C.~Cornell, {\it cond-mat/9707344}. \bibitem{lpe} L.~Frachebourg, P.~L.~Krapivsky, and E.~Ben-Naim, {\sl Phys. Rev. Lett.} {\bf 77}, 2125 (1996); {\sl Phys. Rev. E} {\bf 54}, 6186 (1996). \bibitem{fk} L.~Frachebourg and P.~L.~Krapivsky, {\sl Phys. Rev. E} {\bf 55}, 252 (1997). \bibitem{bou} S.~N.~Marcos-Martin, D.~Beysens, J.~P.~Bouchaud, C.~Godr\`eche, and I.~Yekutieli, {\sl Physica A} {\bf 214}, 396 (1996). \bibitem{ld} B.~Levitan and E.~Domany, {\sl Europhys. Lett.} {\bf 38}, 485 (1997). \bibitem{msbc} S.~N.~Majumdar, C.~Sire, A.~J.~Bray, and S.~J.~Cornell, {\sl Phys. Rev. Lett.} {\bf 77}, 2867 (1996); B.~Derrida, V.~Hakim, and R.~Zeitak, {\sl Phys. Rev. Lett.} {\bf 77}, 2871 (1996). \bibitem{w} A.~Watson, {\sl Science} {\bf 274}, 919 (1996). \bibitem{Ising} E.~Ising, {\sl Z. Phys.} {\bf 31}, 253 (1925). \bibitem{glauber} R.~J.~Glauber, {\sl J. Math. Phys.} {\bf 4}, 294 (1963). \bibitem{kaw} K.~Kawasaki, {\sl Phys. Rev.} {\bf 145}, 224 (1966). \bibitem{kb} P.~L.~Krapivsky and E.~Ben-Naim, {\sl Phys. Rev. E} {\bf 56}, 3788 (1997). \bibitem{wu} F.~Y.~Wu, {\sl Rev. Mod. Phys.} {\bf 54}, 235 (1982). \bibitem{st} For a review of cellular structures see J.~Stavans, {\sl Rep. Prog. Phys.} {\bf 56}, 733 (1993). \bibitem{poly} S.~K.~Kurtz and F.~M.~A.~Carpay, {\sl J. Appl. Phys.} {\bf 51}, 5725 (1981). \bibitem{glaz} J.~A.~Glazier, M.~P.~Anderson, and G.~S.~Grest, {\sl Philos. Mag. B} {\bf 62}, 615 (1990); H.~Flyvbjerg, {\sl Phys. Rev. E} {\bf 47}, 4037 (1993). \bibitem{bubbl} K.~L.~Babcock and R.~M.~Westervelt, {\sl Phys. Rev. Lett.} {\bf 64}, 2168 (1990). \bibitem{bbd} D.~ben-Avraham, M.~A.~Burscka, and C.~R.~Doering, {\sl J. Stat. Phys.} {\bf 60}, 695 (1990). \bibitem{racz} Z.~Racz, {\sl Phys. Rev. Lett.} {\bf 55}, 1707 (1985). \bibitem{amar} J.~Amar and F.~Family, {\sl Phys. Rev. A} {\bf 41}, 3258 (1990). \bibitem{ab} P.~A.~Alemany and D.~ben-Avraham, {\sl Phys. Lett. A} {\bf 206}, 18 (1995). \bibitem{bo} C.~M.~Bender and S.~A.~Orszag, {\it Advanced Mathematical Methods for Scientists and Engineers} (McGraw-Hill Book Co., Singapore, 1984). \bibitem{dz} B.~Derrida and R.~Zeitak, {\sl Phys. Rev. E} {\bf 54}, 2513 (1996). \bibitem{aggr} E.~Ben-Naim and P.~L.~Krapivsky, {\sl Phys. Rev. E} {\bf 53}, 291 (1996). \bibitem{tob} R.~Cordery, S.~Sarker, and J.~Tobochnik, {\sl Phys. Rev. B} {\bf 24}, 5402 (1981). \bibitem{cor} S.~J.~Cornell, K.~Kaski, and R.~B.~Stinchcombe, {\sl Phys. Rev. B} {\bf 44}, 12263 (1991). \bibitem{huse} S.~N.~Majumdar and D.~A.~Huse, {\sl Phys.~Rev.~E} {\bf 52}, 270 (1995). \bibitem{zero} A.~Levy, S.~Reich, and P.~Meakin, {\sl Phys. Lett. A} {\bf 87}, 248 (1982); P.~Meakin and S.~Reich, {\sl Phys. Lett. A} {\bf 92}, 247 (1982). \bibitem{anal} Y.~Elskens and H.~L.~Frisch, {\sl J. Stat. Phys.} {\bf 48}, 1243 (1987); V.~Privman, {\sl Phys. Rev. Lett.} {\bf 69}, 3686 (1993); P.~L.~Krapivsky, {\sl J. Stat. Phys.} {\bf 74}, 1211 (1994). \bibitem{inf} R.~G.~Palmer and H.~L.~Frisch, {\sl J. Stat. Phys.} {\bf 38}, 867 (1985); S.~N.~Majumdar and C.~Sire, {\sl Phys. Rev. Lett.} {\bf 70}, 4022 (1993). \bibitem{bin} A.~Sadiq and K.~Binder, {\sl J. Stat. Phys.} {\bf 35}, 517 (1984). \bibitem{feller} W.~Feller, {\it An Introduction to Probability Theory, Vol. 1} (Wiley, New York, 1971). \bibitem{tam} W.~Y.~Tam, R.~Zeitak, K.~Y.~Szeto, and J.~Stavans, {\sl Phys. Rev. Lett.} {\bf 78}, 1588 (1997). \end{thebibliography} \end{multicols} \end{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%