\def\cp #1 #2 #3 {{\sl Chem.\ Phys.} {\bf #1}, #2 (#3)} \def\jetp #1 #2#3 {{\sl Sov.\ Phys.\ JETP} {\bf #1}, #2 (#3)} \def\jpa #1 #2 #3 {{\sl J. Phys.\ A} {\bf #1}, #2 (#3)} \def\jcp #1 #2 #3 {{\sl J.\ Chem.\ Phys.} {\bf #1}, #2 (#3)} \def\jpc #1 #2 #3 {{\sl J.\ Phys.\ Chem.} {\bf#1}, #2 (#3)} \def\jsp #1 #2 #3 {{\sl J.\ Stat.\ Phys.} {\bf #1}, #2 (#3)} \def\jdep #1 #2 #3 {{\sl J.\ de Physique I} {\bf #1}, #2 (#3)} \def\pra #1 #2 #3 {{\sl Phys.\ Rev.\ A} {\bf #1}, #2 (#3)} \def\prb #1 #2 #3 {{\sl Phys.\ Rev.\ B} {\bf #1}, #2 (#3)} \def\pre #1 #2 #3 {{\sl Phys.\ Rev.\ E} {\bf #1}, #2 (#3)} \def\prl #1 #2 #3 {{\sl Phys.\ Rev.\ Lett.} {\bf #1}, #2 (#3)} \def\rmp #1 #2 #3 {{\sl Rev.\ Mod.\ Phys.} {\bf #1}, #2 (#3)} \def\zpc #1 #2 #3 {{\sl Z. Phys.\ Chem.} {\bf #1}, #2 (#3)} \def\zw #1 #2 #3 {{\sl Z. Wahrsch.\ verw.\ Gebiete} {\bf #1}, #2 (#3)} \def\gtwid{\mathrel{\raise.3ex\hbox{$>$\kern-.75em\lower1ex\hbox{$\sim$}}}} \def\ltwid{\mathrel{\raise.3ex\hbox{$<$\kern-.75em\lower1ex\hbox{$\sim$}}}} \def\eg{{\it e.\ g.}}\def\ie{{\it i.\ e.}}\def\etc{{\it etc.}} \def\vs{{\it vs.}}\def\vv{{\it vice versa}}\def\ea{{\it et al.}} \def\apr{{\it a priori}}\def\apo{{\it a posteriori}} \def\pd#1#2{{\partial #1\over\partial #2}} %partial derivative \def\p2d#1#2{{\partial^2 #1\over\partial #2^2}} %second partial derivative \def\td#1#2{{d #1\over d #2}} %total derivative \def\t2d#1#2{{d^2 #1\over d #2^2}} %second total derivative \def\ith{{$i^{\rm th}$}} \def\jth{{$j^{\rm th}$}} \def\kth{{$k^{\rm th}$}} \def\nth{{$n^{\rm th}$}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \def\pvt{P(D,t)} \def\pvT{P(D,T)} \def\pv0{P(D,0)} \def\pvvt{P(D',t)} \def\pvvT{P(D',T)} \def\pkt{P_k(t)} \def\st{S(t)} \def\ct{c(t)} \def\cT{c(T)} \def\ctt{c(t')} \def\cTT{c(T')} \def\cta{{c\sim t^{-\a}}} \def\dtb{{\dav \sim t^{-\b}}} \def\dav{\langle D\rangle} \def\dm{D_{\rm min}} \def\a{\alpha} \def\b{\beta} \def\m{\mu} \def\d{\delta} \def\e{\epsilon} \def\g{\gamma} \def\G{\Gamma} \def\apb{{\a+\b=1}} \def\CP{\hbox{{$\cal P$}}} \def\erf{\rm erf} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \def\indd{{\int \limits_0^{\infty}dD'}} \def\ind{{\int \limits_0^{\infty}dD}} \def\intt{{\int \limits_0^t dt'}} \def\inTT{{\int \limits_0^T dT'}} \def\lr#1#2#3{\left#1 #2 \right#3} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \documentstyle[multicol,aps]{revtex} \begin{document} \title{Kinetics of Heterogeneous Single-Species Annihilation} \author{P.~L.~Krapivsky, E.~Ben-Naim, and S.~Redner} \address{Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215} \maketitle \begin{abstract} We investigate the kinetics of diffusion-controlled {\it heterogeneous} single-species annihilation, where the diffusivity of each particle may be different. The concentration of the species with the smallest diffusion coefficient has the same time dependence as in {\it homogeneous} single-species annihilation, $A+A \to 0$. However, the concentrations of more mobile species decay as power laws in time, but with non-universal exponents that depend on the ratios of the corresponding diffusivities to that of the least mobile species. We determine these exponents both in a mean-field approximation, which should be valid for spatial dimension $d>2$, and in a phenomenological Smoluchowski theory which is applicable in $d<2$. Our theoretical predictions compare well with both Monte Carlo simulations and with time series expansions. \end{abstract} \medskip P.A.C.S. Numbers: 02.50.-r, 05.40.+j, 82.20.-w, 82.20.Wt \begin{multicols}{2} \section{Introduction} The kinetics of diffusion-controlled single-species annihilation, $A+A \to 0$, when each particle has the same diffusion coefficient is now well understood [1]. For spatial dimension $d>2$, the kinetics may be accounted for by the rate equation which predicts that the density decays as $t^{-1}$ in the long-time limit. For $d\leq2$, various phenomenological approaches predict that the density decays as $t^{-d/2}$, but with logarithmic corrections appearing in $d=2$. Accompanying this relatively slow kinetics is a spatial organization in which the probability of finding particles at small separations is reduced compared to a random distribution. In one dimension, exact solutions, either based on an occupation number formalism [2], or by mapping the reaction onto the kinetic Ising-Glauber model [3], provide definitive results about this spatial organization and the reaction kinetics. Our goal in this paper is to describe the kinetics of {\it heterogeneous} single-species annihilation, which is defined by the reaction scheme $A_i+A_k\buildrel K_{i,j}\over\longrightarrow 0$. Here $A_i$ denotes the \ith\ species and the reaction rate matrix $K_{i,j}$ is a function of the diffusivities of the two reacting species. Although we refer to different species in the context of their diffusivity, the reaction itself is single-species annihilation with distinct rates for different reaction channels. We will consider reactant diffusivities which are drawn from a probability distribution. Such a situation arises naturally when the reactants have different masses. As we shall show, the kinetics of the heterogeneous system is considerably richer than that of the homogeneous analog. Similar behavior was also encountered in heterogeneous single-species annihilation with ballistic particle motion [4], where the kinetics depends in an essential way on the form of the initial distribution of velocities. The present investigation is a natural counterpart of this earlier work for diffusive single-species annihilation. When the number of species is finite (\ie, the density distribution, $P(D,t)$, contains a finite number of discrete peaks), the rate equations predicts that the least mobile species (with diffusivity $\dm$) decays as $t^{-1}$, as in homogeneous annihilation. However, the more mobile species each decay at a faster power-law rate with an associated exponent that depends on the diffusivity ratio between the more mobile and the slowest species. When $P(D,t)$ is continuous but with $\dm>0$, the rate equations again show that the least mobile species predominates in the long time limit, and that the decay of the more mobile species is described by non-universal power-law behavior. However, the detailed form of $P(D,t=0)$ near $\dm$ contributes to a logarithmic prefactor in the decay law. If $\dm=0$, the kinetics strongly depends on the initial conditions. For an initial distribution of diffusivities with a power law tail, $\pv0\sim D^{\m}$, as $D \to 0$, the concentration and the average diffusion coefficient decay as $\cta$ and $\dtb$, respectively, with $\a = (2+2\m)/(3+2\m)$ and $\b = 1-\alpha=1/(3+2\m)$. All of these predictions are expected to apply for $d>2$. For $d\leq 2$, we apply the Smoluchowski theory to account for the reaction kinetics. We first test this phenomenological approach on the ``impurity'' problem, in which there is a small concentration of particles with one diffusion coefficient in a homogeneous background of particles with a different diffusivity. The Smoluchowski theory predicts that the impurity species decays as a non-universal power-law in time. These predictions are in good agreement with numerical results that arise from time-series expansions and Monte Carlo simulations in one dimension. Connections are also made with the exactly-soluble problem of the decay of impurities in the related single-species coalescence process. We then generalize the Smoluchowski approach to both continuous distributions of diffusivities and general spatial dimension $d<2$. The rest of this paper is organized as follows. In Section 2, we derive basic results from the rate equations for both discrete and continuous distributions of diffusivities. In Section 3, we investigate the reaction kinetics for $d\leq2$ by the Smoluchowski approach. This leads to phenomenological rate equations with time-dependent reaction rates. The predictions of this approach are then outlined. In Section 4, we present time series expansion and Monte Carlo simulation results in one dimension to support our theoretical findings. Finally, we give a brief discussion in Section 5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Rate Equations Approach} Consider first the mean-field rate equations for the case of two distinct species, $A_1$ and $A_2$, with respective diffusivities $D_1$ and $D_2$. The two species interact according to the bimolecular processes $A_1+A_1 \to 0$, $A_1+A_2 \to 0$, and $A_2+A_2 \to 0$, with respective rates $K_{11}$, $K_{12}=K_{21}$, and $K_{22}$. The corresponding rate equations are \begin{mathletters} \begin{equation} \td {c_1} t = -K_{11}c_1^2 -K_{12}c_1c_2, \end{equation} \begin{equation} \td {c_2} t = -K_{22}c_2^2 -K_{12}c_1c_2, \end{equation} \end{mathletters} where $c_i$ denotes the concentration of the \ith\ species. Smoluchowski suggested a simple way to relate the reaction rate $K_{ij}$ to the diffusion coefficients $D_i$ and $D_j$, and radii $R_i$ and $R_j$ of the reactants [5]. This derivation is based on considering the low-density limit. In the rest frame of a particle of species $i$, it may be considered as a spherical stationary trap of radius $R$ which is surrounded by cloud of $j$ particles which are captured upon contact with the trap. The reaction rate is identified as the flux of particles of type $j$ to the trap under the boundary conditions of absorption at the surface of the \ith\ particle and a fixed concentration as $r\to\infty$. Upon solving the diffusion equation in three dimensions under these conditions, one straightforwardly finds that this flux equals $K=4\pi DR$ in the long-time limit (see, $\eg$, [6]). If both species perform independent Brownian motions with respective diffusion coefficients $D_i$ and $D_j$, then the appropriate generalization of the reaction rate is $K_{ij}= 4\pi(D_i+D_j)(R_i+R_j)$. Assuming, for simplicity, that $R_1=R_2=R$ and absorbing the numerical factor $8\pi R$ into the overall time scale, we may re-write Eq.~(1) as, \begin{mathletters} \begin{equation} \td {c_1} t = -2D_1c_1^2 -(D_1+D_2)c_1c_2, \end{equation} \begin{equation} \td {c_2} t = -2D_2c_2^2 -(D_1+D_2)c_1c_2. \end{equation} \end{mathletters} To solve these equations, consider first the time dependence of the concentration ratio $\psi=c_2/c_1$. By defining the auxiliary variable $dy=c_1dt$, the rate equation for the concentration ratio is \begin{equation} \td \psi y = -(D_2-D_1) (\psi+\psi^2). \end{equation} Because $\psi$ is a strictly decreasing function of $y$ for $D_2>D_1$, the asymptotic behavior may be estimated by ignoring the quadratic term in Eq.~(3). The corresponding solution is simply \begin{equation} \psi(y) \sim e^{-(D_2-D_1)\int_0^tdt'\,c_1(t')}. \end{equation} Two possibilities for the behavior of $y=\int dt'\,c_1(t')$ as $t\to\infty$ can occur. If $y$ reaches a finite limit as $t\to\infty$, then $c_2(t)\propto c_1(t)$; further, $c_1(t)$ and $c_2(t)$ must decay as $t^{-1}$. It is immediate to show that these two conditions are incompatible with the initial rate equations. On the other hand, if $y\to\infty$ as $t\to\infty$, then $c_2/c_1\to 0$, asymptotically. Thus in Eq.~(2a), $c_1c_2$ may be neglected compared to $c_1^2$ and the resulting approximation gives $c_1(t)\simeq (2D_1t)^{-1}$. Using this in the equation for $\psi$, one ultimately finds \begin{equation} c_2 \simeq C (2D_1t)^{-(1+\d)}, \quad C = {c_2(0) \over c_1(0)} \biggl[{1+c_2(0)/c_1(0) \over c_1(0)}\biggr]^\d \end{equation} for $t \to \infty$, with $\d=(D_2-D_1)/2D_1$. While the less mobile species decays as $t^{-1}$ in the long-time limit, as in homogeneous single-species annihilation, the more mobile species decays non-universally as $t^{-\a}$ with $\a = 1+\d =(D_1+D_2)/2D_1$. It is worth noting that when $\d$ is small, the asymptotic behavior given in Eq.~(5) is reached only at very long times. To demonstrate this, consider, $\eg$, equal initial concentrations of the two species. Then solving Eq.~(3) with $\psi(0)=1$ and $0<\d \ll 1$ gives $\psi(y) = c_2/c_1 \simeq 1+\d\ln[c_1(y)/c_1(0)]$. Substituting this into Eqs.~(2), one finds \begin{equation} c_{1,2}(t) \simeq {1\over 4D_1t}\Bigl[1 \pm {\d \over 2}\ln(D_1t)\Bigr], \end{equation} for $1 \ll t \ll e^{2/\d}$. Consequently, the crossover time between the intermediate asymptotics, Eq.~(6), and the final asymptotics, Eq.~(5), diverges as $e^{2/\d}$ for $\d \to 0$. For a finite number of species, the kinetics is similar to that of the two-species case. Namely the concentration of the least mobile species decays as $t^{-1}$, while the more mobile species exhibit non-universal power-law decays with associated exponents $(D_1+D_k)/2D_1$, where $D_k$ is the diffusion coefficient of the \kth\ species. When the diffusivities of the species are drawn from a continuous distribution, the rate equation for $\pvt$, the concentration of species with diffusivity $D$, becomes the integro-differential equation \begin{eqnarray} \pd \pvt t &=& -\pvt\indd\, (D+D')\,\pvvt \nonumber\\ &=& -\pvt \bigl[D\,\CP_0(t)+\CP_1(t)\bigr]. \end{eqnarray} Here $\CP_k(t)$ denotes the $k^{\rm th}$ moment of the diffusivity distribution $\pvt$, \begin{equation} \CP_k(t) = \ind \,D^k \,\pvt. \end{equation} Note that the zeroth moment of the diffusivity distribution is just the particle density, $\ct = \CP_0(t)$, while the average diffusion coefficient is expressed in terms of the zeroth and first moment by $\dav = \CP_1(t)/\CP_0(t)$. Equations similar to Eq.~(7) describe the kinetics of several irreversible processes, such as diffusion-reaction aggregation [7], ballistic annihilation [4], and ballistic aggregation [8]. For the aggregation problem, a reaction rate of the form $(D+D')$ is known as the ``sum'' kernel, for which the rate equations are exactly soluble [7]. In heterogeneous annihilation, the rate equations are also soluble by elementary analysis. A formal, but implicit solution to Eq.~(7) is \begin{equation} \pvt = \pv0 \exp\Biggl[-D\intt \,\CP_0(t') -\intt \,\CP_1(t')\Biggr]. \end{equation} Further, by integrating Eq.~(7) over $D$, the following differential equation which relates the moments $\CP_0(t)$ and $\CP_1(t)$ is obtained \begin{equation} \td {\CP_0} t = -2\CP_0\CP_1. \end{equation} With the initial condition $\CP_0(t=0)=1$, which fixes the scale of the initial distribution, the solution to Eq.~(10) is \begin{equation} \CP_0(t) = \exp\Biggl[ - 2\intt \,\CP_1(t')\Biggr]. \end{equation} Combining Eqs.~(9) and (11) leads to a simplification in which only the unknown zeroth moment, or particle density $\ct = \CP_0(t)$, appears in the formal solution \begin{equation} \pvt = \pv0 \,\sqrt{\ct}\, \exp\Biggl[-D\intt \,\ctt\Biggr]. \end{equation} A scaling analysis of this solution indicates that two types of behavior can occur: one for initial diffusivity distributions with a finite non-zero lower cutoff, $\dm>0$, and the other for situations where $\dm = 0$. For both cases, we consider only those initial distributions which are homogeneous near $\dm$, namely, $\pv0\sim (D-\dm)^\m$ as $D \to \dm$ with $\m > -1$ for normalizability. This restriction leads to mathematical tractability as well as being a natural illustrative choice. For $\dm = 0$, the average diffusion coefficient clearly decay to zero, and under mild restrictions, the particle concentration also decays to zero. For an initial distribution of the form $\pv0 \sim D^\m$ as $D \to 0$, it is natural to assume power-law decays in the average concentration and average diffusivity: $\cta$ and $\dtb$ for $t \to \infty$. Correspondingly, the time-dependent diffusivity distribution is expected to approach the scaling form, \begin{equation} \pvt \simeq t^{\b-\a} \Phi(Dt^\b). \end{equation} Upon substituting this scaling ansatz into Eq.~(12), one finds two relations which then determine the exponents $\a$ and $\b$. First, for the argument of the exponential in Eq.~(12) to be dimensionless, if follows that $\a+\b=1$. Similarly, for the prefactor to be a function only of the scaling variable $x=Dt^\b$, the relation $\a=2\b(1+\m)$ must hold. These two conditions determine the dependence of the fundamental exponents on $\mu$ , \begin{equation} \a = {2+2\m \over 3+2\m}, \quad \b = {1 \over 3+2\m}, \end{equation} while the scaling function is \begin{equation} \Phi(x) = x^\m e^{-x}. \end{equation} Let us now justify the existence of the scaling ansatz by constructing explicit solutions to Eq.~(12). Consider the initial distribution $\pv0=D^\m e^{-D}/\G(\m+1)$, where the numerical factor normalizes the initial density to unity, for convenience. Substituting this initial distribution into Eq.(12) and then integrating the resulting equation over $D$, one finds the simple differential equation \begin{equation} \td u t = u^{-2\m-2}, \end{equation} for the auxiliary function $u(t)=1+\intt \ctt$. With the initial condition, $u(t=0)=1$, Eq.~(16) is readily solved, from which the concentration is \begin{equation} \ct = \td u t= \Bigl[1+(2\m+3)t\Bigr]^{-{2+2\m \over 3+2\m}}. \end{equation} Substituting this in Eq.~(12), we obtain for the diffusivity distribution \begin{eqnarray} \pvt&=&{D^\m \over \G(\m+1)} \Bigl[1+(2\m+3)t\Bigr]^{-{1+\m \over 3+2\m}} \nonumber\\ &&\exp\Biggl(-D \Bigl[1+(2\m+3)t\Bigr]^{1 \over 3+2\m}\Biggr). \end{eqnarray} As expected, this solution approaches the scaling form of Eq.~(13) asymptotically, with the exponents and scaling function given by Eqs.~(14) and (15), respectively. Consider now initial distributions with $\dm > 0$. For concreteness, we examine the case where \begin{equation} \pv0={(D-\dm)^\m \over \G(\m+1)} e^{-(D-\dm)},\qquad {\rm for\ \ } D>\dm. \end{equation} After substitution of this initial distribution in Eq.~(12), the analog of Eq.~(16) is \begin{equation} \td u t = u^{-2\m-2} e^{-2\dm(u-1)}. \end{equation} Solving this equation in the limit of $t \to \infty$, the asymptotic form for the density is \begin{equation} \ct = {1 \over 2\dm t} \Bigl[1-{2+2\m \over \ln(2\dm t)}+ \ldots \Bigr], \end{equation} while the full distribution of diffusivities has the form \begin{equation} {\pvt \over \pv0} \simeq (2\dm t)^{-{D+\dm \over 2\dm}} \Biggl[{\ln(2\dm t) \over 2\dm}\Biggr]^{{(\m+1)D \over \dm}}. \end{equation} Thus the concentration of the species with diffusivity $D$ decays algebraically with a characteristic non-universal exponent $(D+\dm)/2\dm$, but also with a logarithmic prefactor. This distribution can be expressed in a scaling form by identifying the appropriate scaling variable. For this identification, it is helpful to consider the average diffusion coefficient which has the time dependence $\dav - \dm \sim \dm/\ln(2\dm t)$. This suggests that the scaling variable is \begin{equation} x = {\ln(2\dm t) \over 2\dm}(D-\dm). \end{equation} With this identification and in the scaling limit $t \to \infty$ and $D \to \dm$ but with $x$ finite, the asymptotic distribution Eq.~(22) can be written in the form \begin{equation} \pvt \simeq {1 \over \G(\m+1)}(2\dm t)^{-1}{\ln(2\dm t) \over 2\dm}\Phi(x), \end{equation} with scaling function $\Phi(x)$ again equal to $x^\mu e^{-x}$. Thus for both $\dm = 0$ and $\dm > 0$, the diffusivity distribution approaches the scaling form, given by Eq.~(13) and Eq.~(24), respectively, with characteristic width $t^{-1/(2+2\m)}$ in the former case and $1/\ln(t)$ in the latter. \section{Heterogeneous Annihilation in Low Dimensions} In the diffusion-controlled limit, fluctuations effects govern the long-time kinetics of heterogeneous annihilation for $d\leq 2$, in close analogy with homogeneous single-species annihilation [1]. A method which is ideally suited to account for the kinetics in this regime is the Smoluchowski theory [5]. We adapt this approach to treat the kinetics of heterogeneous annihilation when $d\leq 2$. To describe the method, consider first the simple and illustrative example of a two-species system in one dimension in which a background of identical particles with diffusivity $D$ and density $c$ contains relatively rare impurities of diffusivity $D_I$ at concentration $c_I\ll c$. In the Smoluchowski approach, we first compute particle flux to a ``reference'' absorbing particle due to the rest of the particles which comprise the uniform background. The requisite solution to the background concentration is $c(x,t)= c_\infty{\rm erf}(x/\sqrt{4Dt})$, from which the particle flux at the reference particle is $\phi=c_\infty\sqrt{D/\pi t}$. This is identified as the effective microscopic reaction rate, $\tilde k$. For the two-component system of background and impurity, there are distinct rates associated with reactions between background particles and between the background and impurities. In the limit of low impurity density, we neglect the influence of background-impurity reactions on the background density, as well as reactions among the impurities. Under these restrictions, the rate equations become \begin{eqnarray} \dot c &\cong - 2\tilde k_{BB} c^2 \sim -2 \sqrt{2D\over\pi t} c^2,\nonumber\\ \dot c_I &\cong - 2\tilde k_{BI} c\, c_I \sim -2 \sqrt{D+D_I\over\pi t} c\, c_I, \end{eqnarray} where $\tilde k_{BB}$ and $\tilde k_{BI}$ are the effective rates for background-background and background-impurity reactions. Note that the numerator inside the square root involves the relative diffusivities of the two reacting species. From the first equation, the background concentration vanishes as \begin{mathletters} \begin{equation} c(t)=\sqrt{\pi/32Dt}. \end{equation} This deviates from the exact result [2, 3] for single-species annihilation by the numerical factor $c(t)/c_{\rm exact}(t)=\pi/2$. The crucial feature of the Smoluchowski approach is that by the form of the rate equation for $c_I$, the coefficient of $c(t)$ determines the exponent of the decay of the impurity species. We thereby find \begin{equation} c_I(t)\sim t^{-\sqrt{(1+\e)/8}}, \end{equation} \end{mathletters} with $\e=D_I/D$. As we discuss in the next section, this prediction agrees rather well with numerical results. Intuition for this non-universal behavior can be gained by considering the impurity problem in the case where particles coalesce, rather than annihilate. This situation is considerably simpler than that of annihilation, as the enclosing ``cage'', defined by the nearest neighbors of the impurity, evolves only by diffusion. Thus the many-body problem may be reduced to the three-body problem of the impurity and its two nearest-neighbors. By exploiting simple geometric equivalences, this problem can be transformed to the survival of a single random walker which diffuses within an absorbing two-dimensional wedge whose opening angle depends on $D/D_I$. In this latter problem, the survival probability decays as $t^{-\a}$, with $\a=\pi/(2\cos^{-1}(\e/(1+\e)))$. Here $\theta= \cos^{-1}{[\e/(1+\e)]}$ is the opening angle of the wedge [9]. Such a rigorous mapping does not exist for the impurity problem when particles annihilate. This is an intrinsically a many-body process, since the cage surrounding a given particle can involve distant neighbors. Nevertheless, the Smoluchowski approach is essentially identical for both the annihilation and the coalescence reactions, except for overall factors of 2 in the rate equations. Thus the equivalent of Eq.~(26b) for aggregation is $c_I(t)\sim t^{-\sqrt{(1+\e)/2}}$. While this exponent value deviates considerably from the exact result give above, the $\e$ dependence of the exponent is qualitatively correct. The mechanism underlying the non-universal decay of the impurity is the equivalence to the survival probability of a diffusing particle inside an absorbing interval whose length $L$ grows as $t^\b$. Since the probability density of the particle spreads over a spatial extent of the order of $\sqrt{Dt}$, the survival probability decays exponentially for $\beta<1/2$. However, in the marginal case $\beta=1/2$, \ie, $L\sim A\sqrt{Dt}$, the survival probability $S(t)$ decays algebraically in time, $S(t)\sim t^{-\alpha(A)}$, but with the decay exponent dependent on the {\it dimensionless} parameter $A$. In the reaction process, the size of the cage which surrounds the impurity also grows as $t^{1/2}$ for both annihilation and coalescence, since the overall density is decaying as $t^{-1/2}$. However, the microscopic differences in the two reactions are important because they determine the amplitude $A$ in the growth of the cage. As shown above, the phenomenological Smoluchowski treats annihilation and coalescence on the same footing and thus provides a convenient description of the non-universal behavior. We have also explored the case of small difference between diffusion coefficients, $0<\d \ll 1$. We find that intermediate asymptotic behavior initially occurs before the final asymptotic of Eq.~(26) sets in. In particular, for initially equal concentrations \begin{equation} c_{1,2} \simeq \sqrt{\pi \over 128 D_1t} \Bigl[1 \pm {\d \over 8}\ln(D_1t)\Bigr], \end{equation} for $1 \ll t \ll e^{8/\d}$. The crossover time between intermediate asymptotics Eq.~(27) and final asymptotics Eq.~(26) diverges as $e^{8/\d}$ for $\d \to 0$. Let us now consider heterogeneous annihilation in one dimension with a continuous distribution of particle diffusivities. Since the Smoluchowski approach in one-dimension yields a flux that varies as $\propto \sqrt{D/t}$, one finds a rate equation of the form of Eq.~(7) with the kernel $K(D,D')\propto\sqrt{(D+D')/t}$. The correspondence can be made even closer by introducing the modified time variable $T = 4\sqrt{t/\pi}$ which eliminates the explicit time dependence. The rate equation for $\pvT$ becomes \begin{equation} \pd \pvT T = -\pvT\indd \sqrt{D+D'}\,\pvvT. \end{equation} We are unable to find either a scaling solution to this equation or the exponents. However, presumably exact values of the exponents can be obtained by consideration of a closely related and more tractable model. We replace the kernel $K(D,D') = \sqrt{D+D'}$ by one with the same homogeneity degree, $K(D,D') =\sqrt{D}+\sqrt{D'}$. The two kernels obey the bounds ${\sqrt{D}+\sqrt{D'} \over \sqrt{2}} \le \sqrt{D+D'}\le \sqrt{D}+\sqrt{D'}$. This suggests that the asymptotic behavior from the initial kernel should be identical to that predicted by the simpler sum-root kernel $K(D,D')=\sqrt{D}+\sqrt{D'}$. With this modified kernel, the rate equation becomes \begin{equation} \pd \pvT T = -\pvT\, [\sqrt{D}\,\CP_0(T)+\CP_{1/2}(T)], \end{equation} with $\CP_{1/2}(T) = \ind\, \sqrt{D}\,\pvT$. Repeating the steps employed for the mean-field treatment, the formal exact solution to Eq.~(29) is \begin{equation} \pvT = \pv0 \sqrt{\cT} \exp\Biggl[-D\inTT \,\cTT\Biggr]. \end{equation} Again we substitute the scaling ansatz $P(D,t)\sim t^{\b-\a}\Phi(Dt^\b)$ into the above rate equation with the initial condition $P(D,0)\sim D^{\mu}$ as $D\to 0$. By expressing all factors in a scaling form we find the two exponent relations $2\a+\b=1$ and $\a=2\b(1+\m)$ and hence, \begin{equation} \a = {2+2\m \over 5+4\m}, \quad \b = {1 \over 5+4\m}. \end{equation} Here the exponents $\a$ and $\b$ refer to the time dependence of $\ct$ and $\dav$ in terms of the physical time variable $t$. These results can be straightforwardly generalized to arbitrary $d<2$ within the Smoluchowski approach. A simple calculation shows that the flux to a trap behaves as $D^{d/2}t^{-1+d/2}$ for $d<2$. In the rate equations, this suggests the introduction of the modified time $T \sim (Dt)^{d/2}$. One thereby obtains an equation similar to Eq.~(28), but with the kernel $K(D,D') = (D+D')^{d/2}$. Next we replace this kernel by more tractable form, $K(D,D') = D^{d/2}+D'^{d/2}$. Applying a scaling analysis to the rate equations with this kernel then gives the exponents $\a =2d(1+\m)/[d+4(1+\m)]$ and $\b = d/[d+4(1+\m)]$. In two dimensions, the Smoluchowski approach shows that the reaction kernel is $K(D,D')=4\pi(D+D')/\ln(Dt)$. In this sense, two dimensions is the marginal case which demarcates the regime where the reaction rate is time independent (for $d>2$) from the regime where the reaction rate varies as a power law in time ($d<2$). The corresponding rate equation for $\pvt$ in $d=2$ is \begin{equation} \pd \pvt t = -\pvt\indd {4\pi(D+D')\over \ln[(D+D')t]}\pvvt. \end{equation} We analyze this equation in the same spirit as that employed for $d\ne 2$. First, we approximate the slowly varying logarithmic factor by $\ln[(D+D')t] \simeq \ln[\dav t] \simeq \a \ln(t)$ and then introduce the auxiliary time variable $T \to 4\pi t/\a \ln(t)$, so that Eq.~(32) becomes identical in form to the mean-field rate equation, Eq.~(7). Thus replacing $t$ by $T$, the results of the preceding section apply to the present situation. In particular, at asymptotically large times $\pvt$ approaches the scaling form \begin{equation} \pvt \simeq [t/\ln(t)]^{\b-\a} \Phi(D[t/\ln(t)]^\b), \end{equation} with exponents and the scaling function given by Eqs.~(14) and (15). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Numerical Results} Consider the specific and illustrative example of a two-component system consisting of an impurity particle with diffusivity $D_I$ within a background of identical particles with diffusivity $D$. We first describe the results of a time series analysis for the survival probability of a static impurity. Then we discuss complementary Monte Carlo simulations for both the time dependence of the survival probability and the eventual survival probability in a finite system. Both techniques yield estimates for the decay exponent of the impurity that are in good agreement with the Smoluchowski theory predictions. The time series expansion is a general technique for evaluating numerically the initial terms in the exact power-series expansion in time for various observables that characterize reactive systems [10]. The technique involves the successive application of a suitably-defined evolution operator on the initial state of the system. After each application of the evolution operator, configurations for which the impurity particle has not reacted contribute to the survival probability. To compute the first $n$ coefficients of the time power series, it is sufficient to consider the evolution of an $n$-site ring. The primary limitation of the technique is computer memory, since the number of configurations grows exponentially with the order of the expansion. For the case of the static impurity, we obtained the expansion to order 21, as shown in Table~1. We then apply the Pad\'e analysis method suggested in Ref.~[11] to estimate the asymptotic properties of the survival probability. In Figure 1, the diagonal Pad\'e approximants for the survival probability, $S_{[n,n]}(t)$ are plotted for $n=8,9$, and 10. These approximants are essentially identical for $t\ltwid6$ and we conclude that they accurately describe the true survival probability within this time range. These Pad\'e approximants are then fitted to the power-law form $S(t)\sim t^{-\a}$ using the fitting procedure also suggested in [11]. We thereby estimate the exponent value $\a = 0.38 \pm 0.01$ (Table~2). Monte Carlo simulations were also performed for the impurity problem in one dimension. A typical measurement involved 10 configurations of a periodic chain of $5\times 10^5$ sites with initial concentrations $c_I(0)=0.01$ and $c(0)=0.99$. The data for the time dependence of both $c_I(t)$ and $c(t)$ are quite straight on a double logarithmic plot and exponent estimates may be made based on the local slopes of nearby data points. From the slopes of first-neighbor, second-neighbor, and third-neighbor pairs, we estimate that the background species decays as $t^{-\a_B}$, where $\a_B=0.50\pm 0.01$. Here the error bar indicates the magnitude of the fluctuations in the local slopes in the time regime where the data is most linear (between $20\ltwid t\ltwid 10,000$). Since $\a_B=1/2$, this analysis provides a useful confidence test. For the impurity species, the same analysis method suggests $\alpha= 0.37\pm 0.01$ (Fig.~2). A complementary approach is based on a finite size scaling analysis for the {\it eventual} survival probability of an impurity, $S(L)= \lim_{t\to\infty}S(L,t)$, when starting with a single impurity and the rest of the sites occupied by background particles on finite ring with an odd number of sites, $L$. We postulate that $S(L)\sim L^{-\gamma}$ (clearly $\g=1$ for $D_I=D$), while for short times, the survival probability must be independent of the size of the system, so that $S(L,t)\sim t^{-\a}$. Since the two basic length scales are $\sqrt{Dt}$ and $L$, we assume that the survival probability has the scaling form $S(L,t)\sim t^{-\a}f({\sqrt{Dt}/ L})$. The scaling function $f$ must have the limiting behaviors $f(x)\sim 1$ for $x\to 0$ and $f(x)\sim x^{2\a}$ for $x\to\infty$ to account for the two asymptotic limits. These conditions imply that $\g=2\a$. An advantage of the measurement of the eventual survival probability is that it typically gives more accurate exponent estimates compared to the time-dependent observations. We have performed Monte Carlo simulations of finite single impurity systems until the impurity disappears or a single impurity remains to measure the exponent $\gamma$ directly. In Figure~3, we present the average survival probability for $10^7$ realizations, for systems of size $2^n+1$, with $n=1,\ldots,8$. The data suggests the decay exponent of $\gamma=0.75\pm0.01$, consistent with the series and Monte Carlo estimates of $\alpha=0.38\pm0.01$. Table 2 compares the decay exponents found for the other values of $\epsilon=D_I/D$ to the corresponding exponent given by the Smoluchowski approach. The Smoluchowski approximation yield a remarkably good description of the asymptotic behavior of the impurity decay in the annihilation reaction. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Discussion} We have shown that there is non-universal kinetic behavior in diffusion-controlled {\it heterogeneous} single-species annihilation. Typically, the concentration of the species with the smallest diffusion coefficient has the same time dependence as in {\it homogeneous} single-species annihilation, $A+A \to 0$. However, the concentrations of the more mobile species decay as power laws in time, with non-universal exponents that depend on the ratios of the corresponding diffusivities to that of the least mobile species. These exponents were determined by a mean-field approximation, which should be valid for $d>2$, and by a phenomenological Smoluchowski theory which should apply for $d<2$. Our numerical studies focused on the limiting case where there is an infinitesimal concentration of one species in a homogeneous background of another species. For this impurity problem, the behavior of the cage which surrounds the impurity is the crucial ingredient which determines the kinetics of the impurity. For the coalescence reaction, $A+A\to A$, the cage consists of the two nearest-neighbors of the impurity. These enclosing particles continue to undergo diffusion even if there are reactions with other more distant particles in the system. This allows one to reduce the many-body system to a three-particle system which can then be solved exactly in one dimension. In the annihilation reaction, a nearest-neighbor of the impurity can disappear because of a reaction with the next-nearest background particle. This can cause a relatively large rearrangement in which the cage expands to the next background particle. The Smoluchowski approach turns out to be ideally suited for providing a simple description of the non-universal behavior of the impurity survival probability. Although this approach involves an uncontrolled approximation, the Smoluchowski method provides a simple way to quantify the enclosing cage of the impurity for both the coalescence and annihilation reactions. This is a remarkable success in view of the relatively complex cage dynamics in the annihilation reaction. From the information about the size of the cage, the decay exponent follows directly. In particular, for the annihilation reaction, the exponent predictions compare well with numerical estimates. We thank M.~Bramson for stimulating our interest in this problem. We gratefully acknowledge ARO grant \#DAAH04-93-G-0021, NSF grant \#DMR-9219845, and the Donors of The Petroleum Research Fund, administered by the American Chemical Society, for partial support of this research. \begin{thebibliography}{99} \bibitem{} See \eg, M.~Bramson and D.~Griffeath, \zw 53 183 1980 ; D. C. Torney and H. M. McConnell, \jpc 87 1941 1983 . \bibitem{} A.~A.~Lushnikov, \jetp 64 811 1986 . \bibitem{} Z.~R\'acz, \prl 55 1707 1985 ; J. G. Amar and F. Family, \pra 41 3258 1990 . V. Privman, \jsp 69 629 1992 . \bibitem{} E.~Ben-Naim, S.~Redner, and F.~Leyvraz, \prl 70 1890 1993 . \bibitem{} M.~V.~Smoluchowski, \zpc 92 215 1917 . \bibitem{} S.~Chandrasekhar, \rmp 15 1 1943 . \bibitem{} R.~L.~Drake, in {\it Topics in Current Aerosol Research}, Vol. III, part 2, eds. G.~M.~Hidy and J.~R.~Brock (Oxford, Pergamon, 1972). \bibitem{} E.~Ben-Naim, P.~L.~Krapivsky, and S.~Redner, \pre 50 xxx 1994 . \bibitem{} See \eg, M.~E.~Fisher, \jsp 34 667 1984 ; M.~E.~Fisher and M.~P.~Gelfand, \jsp 53 175 1988 ; D.~ben-Avraham, \jcp 88 941 1988 ; M. Bramson and D. Griffeath in {\it Random walks, Brownian motion, and interacting particle systems: A festschrift in honor of Frank Spitzer}, Eds.\ R. Durrett and H. Kesten (Birkhauser, 1991), and references therein; D. Considine and S. Redner, \jpa 22 1621 1989 . \bibitem{} S.~Song and D.~Poland, \jpa 25 3914 1992 ; I.~Jensen and R.~Dickman \jsp 71 89 1993 . \bibitem{} E.~Ben-Naim and J.~Zhuo, \pre 48 2603 1993 . \end{thebibliography} \section*{Figure Captions} {\bf Figure 1.} Time dependence of the diagonal Pad\'e approximants $S_{[n,n]}(t)$ ($n=8,9,10$) based on the time series for the survival probability in the case of a static impurity. A line of slope $-0.38$ is plotted for reference. {\bf Figure 2.} Monte Carlo simulation results for the survival probability of impurity (squares) and background particles (circles) based on 10 configurations of a one-dimensional chain of 500,000 sites. The initial impurity and background concentrations are 0.01 and 0.99, respectively. {\bf Figure 3.} The eventual survival probability $S(L)$ (circles) of a single static impurity on a one-dimensional ring of $L$ sites in which the rest of the ring is initially filled with diffusing particles. The data shown is based on $10^7$ realizations. A line of slope $-0.75$ is shown for reference. \section*{Table Captions} {\bf Table 1} The first 21 coefficients in the time series for the survival probability of a static impurity, defined by $S(t)=\sum_n c_n t^n$. {\bf Table 2} Comparison of the numerical estimates for the exponent $\gamma$ that characterizes the $L$ dependence of the eventual survival probability with the corresponding predictions from the Smoluchowski theory. Estimates are based on averaging over $10^6$ realizations. \end{multicols} \end{document}