\documentstyle[aps,multicol]{revtex} \begin{document} \title{Time Series Expansion for Reaction Processes} \author{E.~Ben-Naim and J.~Zhuo} \address{Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215} \maketitle \begin{abstract} The temporal evolution of reactive systems is studied by means of exact enumeration of all states of the system for short times. An algorithm for the evaluation of the power series for basic quantities of interest is presented. In addition, a consistent methd of estimating the asymptotic behavior of the series is suggested. This approach is applied to two reaction schemes, ballistic annihilation and diffusive-driven three species annihilation, where previous Monte Carlo simulations suggest power-law decays for the concentration. In both cases, reasonable estimates for the decay exponents are obtained. \end{abstract} \pacs{PACS numbers:05.20.Dd, 05.40.+j., 02.50.+s, 82.20.Mj.} \begin{multicols}{2} \section{Introduction} There has been considerable recent interest in understanding the kinetics of simple reaction models. Thus far, there are relatively few theoretical methods available to study these reactions and these approaches are usually limited to specific situations. Typically, the transport mechanism of the reactants is specified and the particles interact upon contact. Hence, the use of the representation number formalism is natural and has the advantage that a {\it linear} rate equation emerges [1,2]. From a theoretical viewpoint, exact methods [3] as well as perturbative methodsK [4] can be applied to this rate equation. From a numerical viewpoint, the exact solution to any reaction scheme can be evaluated for the first few time steps. If this can be carried to a sufficiently high order, then one might hope to extract useful information about the long-time kinetic behavior of the reaction. There have been two noteworthy studies which have followed this general approach [5,6]. Dickman and Jensen applied the series expansion formalism to nonequilibrium interacting particle systems, while Song and Poland investigated several generic reaction-diffusion processes. While this latter approach appears to be promising, the series calculation was carried out by hand. Furthermore, the analysis methods used sometimes lacked internal consistency, leading to an uncertain interpretation of the results. In this study, we describe a systematic approach to obtain time power series for general types of reaction processes. This exact enumeration method holds in any spatial dimension; however, it is hard to implement numerically for dimensions higher than one. Therefore our discussion will be restricted to the one-dimensional case only. In section II, the theoretical description of reactions in terms of creation and annihilation operators is presented. We then discuss how to convert the theoretical formalism into a computational algorithm simply by considering the time evolution of small clusters with periodic boundary conditions. In section III, we develop a simple and relatively direct method to determine the asymptotic behavior of basic physical quantities from the results of the time series expansion. We first illustrate the failure of the straightforward ratio method for analyzing these types of series. This leads us to suggest a direct application of Pad\'e approximants as a method for determining the asymptotic behavior. We exploit the exactly soluble case of single species annihilation to demonstrate the relative utility of the various series analysis methods. In section IV, the time series expansion is applied to ballistic annihilation and to three-species diffusive annihilation. In both cases, analysis of the series provides useful estimates for the asymptotic kinetic behavior. \section{Theory and Implementation} In this section we describe how the temporal evolution of a reactive system can be determined in terms of the action of an evolution operator on an initial state. This approach is advantageous since it leads to a linear equation of motion. Thus a formal solution can be written for the state of the system at any given time. The use of the occupation number representation is natural in this formulation and facilitates the numerical evaluation of the coefficients of the exact time power series. In general, a specific site $i$ can be either empty or occupied by one or more particles. For simplicity, we restrict ourself to models with singly occupied sites, although a generalization to multiple occupancy is possible. In addition, we are interested in situations where there are only a few different species. The state of a single site $i$ will take the symbolic values $| \phi_i \rangle=| o_i\rangle, |a_i\rangle,|b_i\rangle$ {\it etc.}, to describe an empty site, an $A$-occupied site or a $B$-occupied site, respectively. A typical configuration of reactants can be constructed as a product of all the site states $ |\phi_i\rangle=\prod_{i} | \phi_i\rangle$. The state of a system may now be written as \begin{equation} | \Psi\rangle=\sum_{\Phi}P(\Phi,t)| \phi_i\rangle, \end{equation} where $P(\Phi,t)=\langle\Phi| \Psi(t)\rangle$ is the probability that the system is in configuration $| \phi_i\rangle$ at time $t$. The state function $| \Psi\rangle$ obeys a linear equation of motion \begin{equation} {\partial\over\partial t}| \Psi\rangle={\cal L}|\Psi\rangle, \end{equation} where the operator ${\cal L}=\sum_{\Phi\Phi'}{\cal L}_{\Phi\Phi'}|\Phi\rangle\langle\Phi'|$ enumerates all possible transitions between two configurations. Thus ${\cal L}_{\Phi\Phi'}$ is the rate of transition from $\Phi'$ to $\Phi$. This evolution operator can be represented in a compact form using creation and annihilation operators. For example, the creation operator for the $A$ species at site $i$, $A_i^{\dagger}$, is non-zero only when applied to the empty state, $A_i^{\dagger}|o_i\rangle=| a_i\rangle$. Similarly, the annihilation operator $A_i$ at site $i$ is non-vanishing when applied to $| a_i\rangle$ only, $A_i| a_i\rangle=| o_i\rangle$. With these elemental operators, it is now possible to express the effect of the diffusion step and the reaction step on an arbitrary configuration of particles. For example, the situation where all the $A$ particles hop to the right can be represented in terms of the operator \begin{equation} {\cal L}_{AO}=D\sum_i (1-A^{\dagger}_iA_{i+1})A_iA^{\dagger}_{i+1}. \end{equation} The positive term accounts for hopping from site $i$ to site $i+1$ with rate $D$, and the negative term corresponding to the loss of probability from site $i$ due to hopping to the right. To be concrete, for the two site configuration $| ao\rangle$, ${\cal L}_{AO}| ao\rangle=D(| oa\rangle-| ao\rangle)$, while ${\cal L}_{AO}$ yields zero when applied to any other two-site configuration. The operator that corresponds to the actual reaction can also be represented in a similar fashion. For example, for the case of mutual annihilation of two reactive species $A$ and $B$ via the reaction $A+B\to 0$, the rate operator for the annihilation process is \begin{equation} {\cal L}_{AB}=2D\sum_i (1-A_i^{\dagger}B^{\dagger}_{i+1})A_i B_{i+1}. \end{equation} The factor of two accounts for the fact that a hop of either the $A$ or the $B$ in the correct direction leads to an annihilation event. This operator gives a non-zero result only when applied to $| ab\rangle$, that is $ {\cal L}_{AB}| ab\rangle=2D(| oo\rangle-| ab\rangle)$. In fact, if only nearest neighbors interact, ${\cal L}$ is just a linear combination of operators of the form ${\cal L}_{AO}$ and ${\cal L}_{AB}$. In other words, a reaction process is well defined by the interaction of a particle with a neighboring empty site and a neighboring particle . Once the initial state of the system is specified, the state of system at any later time is given by the following formal solution of Eq.~(2) \begin{equation} | \Psi(t)\rangle=e^{{\cal L}t}| \Psi(0)\rangle. \end{equation} In most cases, an exact solution is hard to obtain, and one is led to expand the state function as a power series around $t=0$. By defining $| \Psi_j\rangle={\cal L}^j| \Psi(0)\rangle/j!$, the series expansion to the state function takes the form \begin{equation} | \Psi(t)\rangle=\sum_{j=0}^{\infty} t^j| \Psi_j\rangle. \end{equation} Typical quantities of interest, such as the concentration of a given species or particle correlation functions, can be obtained by evaluating expectation values of a suitably-defined operator. For example, the concentration of the $A$-species is given by \begin{equation} c_A(t)={1\over N}\sum_{\Phi}N_A (\Phi)\langle \Phi|\Psi(t)\rangle, \end{equation} where $N_A (\Phi)$ is the number of $A$ particles in configuration $| \phi_i\rangle$. We now discuss how the aforementioned formalism can be developed into a calculational procedure to provide a numerical evaluation of the time series. The elements of the time evolution operator typically represent interactions between nearest neighbors only, and thus the successive application of the operator $j$ times connects sites which are separated by a distance $j$. Hence, to obtain a series to $n^{\rm th}$ order, one can restrict attention to a ring of size $N=n+1$ with periodic boundary conditions. It is useful to note that in the case of nearest neighbors interactions only, the rate operator is well defined by its action on all possible pair configurations. Therefore, one may consider the full evolution operator as a superposition of pair operators and then apply the pair operator to all $N$ pairs of a given configuration. Typically, a pair operator is non-vanishing only for a few pairs and thus a code written for computing series in one process can be converted easily to evaluate the series for another process. In practice, the number of terms in $| \Psi_j\rangle$ increases exponentially in $j$ and it is advantageous to classify the initial state and the rate operator in terms of symmetries to optimize the computation. For spatially homogeneous systems, the initial state as well as the rate operator are invariant under discrete rotations, and it is therefore useful to represent rotationally equivalent configurations by only a single configuration. Moreover, processes may obey a conjugation symmetry. For example, in the case of diffusion-limited $A+B\to 0$, only states with a majority of $A$ particles need be considered. By taking advantage of these various symmetries, we are able to reach higher order in the series enumeration. As an example of the formalism, let us calculate the first three terms in a ballistic annihilation process. Further details of this model will be discussed in Sec.~IV. In this model one species ($A$) can hop to the right only, and the other species ($B$) can hop to the left only. Whenever a particle lands on any occupied site an annihilation event occurs, {\it i.e.}, mutual reaction as well as self-reaction occur. The rate operator takes the form \begin{equation} {\cal L}={\cal L}_{AO}+{\cal L}_{OB}+{\cal L}_{AB}+{{\cal L}_{AA}\over 2}+ {{\cal L}_{BB}\over2}, \end{equation} with the definitions of Eqs.~(3) and (4). Here ${\cal L}_{AA}$ is obtained from ${\cal L}_{AB}$ by replacing $B$ with $A$. Note, however, that self-reaction occurs with half the rate of mutual reaction since a hop of only one of the particles in a same species pair will lead to a reaction. Let us consider an initial state in which each site has the same probability of being occupied by an $A$ or a $B$. Thus using rotational invariance, the initial state can be written as \begin{equation} | \psi_0\rangle=(| aaa\rangle+3| aab\rangle+3| abb\rangle+| bbb\rangle)/8. \end{equation} By applying the rate operator of Eq.~(8) twice, one can directly verify that the next two terms in the time series expansion of the state function are, \begin{eqnarray} | \psi_1\rangle&=3D(| aoo\rangle+| boo\rangle-2| \psi_0\rangle)/2\nonumber\\ | \psi_2\rangle&=9D^2(2| \psi_0\rangle-| aoo\rangle-| boo\rangle)/4. \end{eqnarray} Next, by applying Eq.~(7), we obtain the concentration to second order as, \begin{equation} c_A(t)=c_B(t)=1/2-Dt+3(Dt)^2/2+\cdots\ . \end{equation} Although this exercise is quite simple, the number of terms in the $| \Psi_j\rangle$ grows rapidly in $j$ and the computation is best done systematically by a computer program. \section{Asymptotic analysis of the time expansion} We are now interested in extracting the asymptotic behavior from the short-time power series representation of observables. In previous studies, this extrapolation procedure was accomplished by direct methods, as well as by using various forms of Pad\'e analysis. Unfortunately, the accuracy and consistency of specific methods varies widely from series to series. Thus many previous approaches incorporate some degree of subjectivity in the analysis. This is a feature which one would hope to avoid in the ideal situation. We will propose an extrapolation approach which ameliorates subjective aspects to a large extent. We first consider the analysis of the series for the exactly soluble process $A+A\to 0$ to illustrate a misuse of the ratio method. We also demonstrate how direct application of Pad\'e approximants can serve as a simple yet useful way to determine the asymptotic behavior. The basic question that we shall address is to ascertain the appropriate analysis method for estimating the exponent $\alpha$ in the asymptotic form for the time dependent concentration \begin{equation} c(t)\sim t^{-\alpha} \end{equation} from the first $n+1$ terms in the series expansion \begin{equation} c(t)\approx c_0+c_1t+\cdots+c_nt^n \end{equation} To test the utility of any analysis method, we expand the exact solution for the single species annihilation process to an arbitrarily large order and then analyze the series as if the exact solution was not available. For single species annihilation the rate operator takes the form \begin{equation} {\cal L}={\cal L}_{AO}+{\cal L}_{OA}+{\cal L}_{AA}\, . \end{equation} Using Eqs.~(3) and (4) it is found that the quartic term $A_i^{\dagger} A_i A_{i+1}A^{\dagger}_{i+1}$ vanishes, and an exact solution can be obtained via a transformation to Fermion operators [3,7]. The exact expression for the concentration is \begin{equation} c(t)=I_0(4Dt)\exp(-4Dt), \end{equation} from which the asymptotic behavior is $c(t)\sim t^{-1/2}$. Let us now attempt to determine the exponent $\alpha$ defined in Eq.~(12) by a direct ratio analysis. In this method, one defines $y=1-c(t)/c(0)$, and using Eq.~(12), one obtains $t(y)\sim (1-y)^{-1/\alpha}$. Given the series expansion for $c(t)$, series inversion formulae yield the series expansion for $t(y)$, namely $t=\sum_{i=1}^{n} b_i y^i$. The latter series is now compared to the expansion of $(1-y)^{-1/\alpha}$ and an estimate for the exponent is expressed in terms of the successive ratios $b_i/b_{i-1}$, $\alpha_i=1/\bigl(1-i(1-b_i/b_{i-1})\bigr)$ [6]. Following this recipe we expand Eq.~(15) to the $30^{\rm th}$ order and get the following approximants for $\alpha$ \begin{eqnarray} \alpha_i=&2,\ 1.2,\ .89,\ .73,\ .62,\ .55,\ .50,\ .46,\ .43,\ .40,\ \nonumber\\ &.38,\ .36,\ .35,\ .34,\ .33,\ .33,\ .32,\ .32,\ .32,\ .33,\ \\ &.34,\ .35,\ .37,\ .40,\ .45,\ .54,\ .74,\ 1.40,\ -6.4,\ -.72\nonumber \end{eqnarray} which do not converge to the exact value $\alpha=1/2$. A possible source of the non-convergence may be that the transformation from infinite time to $y=1$ introduced singularities that ``interfere'' with the asymptotic behavior. Similar results occurred for various related processes and in all the 1D processes considered by Song and Poland. We now suggest an alternative and simple method which appears to be better suited for asymptotic analysis than the ratio method. The analysis can be summarized by the following steps: \begin{itemize} \item Representing the original function by different Pad\'e approximants. \item Considering the two approximants which agree over the longest temporal range. \item Fitting the aforementioned approximants to a power-law function and thereby determining the decay exponent. \end{itemize} For the purpose of asymptotic analysis, it is useful to return to a more basic question and consider the best way to approximate the {\it concentration}, given the Taylor series of Eq.~(13). We again consider single-species annihilation and limit ourself to the first 25 series terms, since this is a typical upper limit of what we can attain by the series computation with current resources. A plot of the truncated series for the concentration {\it vs.}\ the exact solution given Eq.~(15) shows that the Taylor series converges to the exact solution only up to $t\cong 1$ (see Fig.~1(a)). If we define a cut-off point $\tau_0$ as the earliest time when the discrepancy between two functions exceeds $\epsilon=0.001$, we obtain $\tau_0=1.1$. Since we are interested in determining the behavior of the concentration in the long-time limit, we see that na\"\i ve use of the Taylor series is of limited practical utility in representing the concentration. A significant improvement can be achieved by applying the {Pad\'e } approximation directly to the time expansion of $c(t)$. That is, we form the rational function \begin{equation} P_{[n,m]}(t)={a_0+a_1t+\cdots+a_nt^n\over b_0+b_1t+\cdots+b_mt^m}, \end{equation} where $b_0=1$ and the rest of the $n+m+1$ coefficients are chosen such that the expansions of $P_{[n,m]}(t)$ and of $c(t)$ are identical up to the $n+m^{\rm th}$ order. Examination of the temporal range where the various {Pad\'e } approximants converge to the exact solution reveals the utility of this approximation. For the so-called ``diagonal approximants'' $P_{[n,n]}$, one finds that the approximant and the exact solution agree to within 0.001 up to $\tau_0=3.1$, 3.8, 4.5, 5.3 for $n=9$, 10, 11, 12, respectively (Fig.~1(a)). Moreover, at $t=10$ the $[11,11]$ and $[12,12]$ approximants shown in Fig.~1(a) agree with the exact function up to 10\%, while for $t{\mathrel{\raise.3ex\hbox{$>$\kern-.75em\lower1ex\hbox{$\sim$}}}} 1.1$ the 25-term Taylor series diverges badly. Unlike the example above, the original function is usually unknown. However, Fig.~1(a) clearly shows that the exact solution can be well approximated over a substantial temporal range by {Pad\'e } approximants. A criterion for the time regime over which two different {Pad\'e } approximants accurately represent a function can be developed by defining a cut-off time $t_0\bigl([n,m],[n',m']\bigr)$ where the $[n,m]$ and $[n',m']$ approximants first differ by $\epsilon=0.001$. We consider the two {Pad\'e } approximants which give the largest value of $t_{\rm max}={\rm max}\ \{t_0\}$ as the best approximations to the original function. Typically, this analysis extends the convergence of the series well into the asymptotic region, and thus enables one to draw conclusions about the long time nature of the decay. Although in the example discussed above, the approximants exhibit a monotonic improvement in convergence as the order is increased, one should be cautious since monotonic convergence does not always occur. Increasing the order of the approximant eventually improves the outcome, but the existence of spurious singularities close to the physical pole in specific Pad\'e approximants may lead to a transitory decrease in the quality of the results. In our application of the {Pad\'e } method, we have adopted the following procedure to estimate the exponent governing the long-time decay. We define a quality of fit function $g(\alpha)$ by \begin{equation} g(\alpha)={\langle t^{-\alpha}| P_{[n,m]}(t)\rangle\over {\Bigl( \langle t^{-\alpha} |t^{-\alpha}\rangle \langle P_{[n,m]}(t)| P_{[n,m]}(t)\rangle \Bigr)^{1/2}}}, \end{equation} with the inner product $\langle {f_1}| {f_2}\rangle \equiv\int_{t_1}^{t_2} f_1(t)f_2(t)dt$. This inner product measures how ``close'' a given {Pad\'e } approximant $P_{[n,m]}(t)$ is to $t^{-\alpha}$ in the time interval $[t_1,t_2]$, and for a perfect fit one has $g(\alpha)=1$. The optimal exponent $\alpha_{\rm opt}$ is given by the value that maximizes $g(\alpha)$. There are several constraints concerning the integration interval $[t_1,t_2]$. First, the short-time regime has to be excluded since the concentration does not follow a power law. Second, the cutoff $t_2$ should be as large as possible in order to minimize the effects of the short-time corrections to the asymptotic decay, but small enough so that the {Pad\'e } approximant is still accurate. Hence, we use $t_2=\!t_{rm max}$ with the definition given above and define $t_1\!=\!t_{rm max}-1$ in order to get a uniform integration length. Finally, the two approximants corresponding to the largest value of $t_{\rm max}$ will yield two different exponents and we average these two values to quote the estimated exponent $\alpha_{\rm est}$. An averaging procedure is chosen here since it is not possible to determine {\it a priori} which approximant is closer to the original function. We now illustrate this fitting procedure by applying it to the known series for the reaction $A+A\to 0$. By comparing all possible pairs of approximants given the $25$-term time series expansion we find that $P_{[11,11]}$ and $P_{[12,12]}$ yield the largest value of $t_{\rm max}=4.3$. Since the resulting values of $g(\alpha)$ are very close to unity, and weakly varying on $\alpha$, it is useful to locate the minimum of $f\equiv 1-g$ {\it vs.}\ $\alpha$ (see Fig.~1(b)). Following this procedure, we find $\alpha_{\rm opt}=0.512$ and $\alpha_{\rm opt}=0.508$ for the [11,11] and [12,12] approximants, respectively, with the integration range $[3.3,4.3]$. Thus, we quote the estimated value $\alpha_{\rm est}=0.51$ for the decay exponent. This differs by 2\% from the exact value $\alpha=1/2$. The dependence of the estimated exponents on the fitting procedure can be examined by varying the parameters used. The tolerance parameter $\epsilon$ was varied by $10\%$ and the above analysis showed a $0.01\%$ change in the resulting $\alpha_{\rm est}$, due to a corresponding change in the integration limits. Moreover, the length of the integration interval was varied by $10\%$, with the upper limit kept fixed, and in this case the estimated exponent agreed with the original estimate to within $0.02\%$. These tests demonstrate the weak dependence of the analysis outcome on the fitting parameters. To demonstrate the rate by which the estimated exponents approach the ultimate value, we performed the aforementioned analysis with a varying order of the time series expansion. The resulting estimated exponents were $\alpha_{\rm est}=0.517$, $0.507$ and $0.505$ for a Taylor series of the 20, 30 and $40$ order, respectively. We conclude that convergence to limiting value does indeed occur, although the rate of convergence is quite slow. \section{Application to reaction models} In this section we apply the power series expansion method to two reaction processes where only Monte Carlo simulation results and scaling arguments for the asymptotic behavior are presently available. The first is a ballistic annihilation process with asynchronous particle motion, and the second is a diffusive three species annihilation process. \subsection{Asynchronous Ballistic Annihilation} Consider a ballistic annihilation process in one dimension in which each particle moves ballistically with velocity $v_0$ or $-v_0$ with equal probabilities. Either from the exact solution [8] or from a qualitative random-walk argument, the concentration decays asymptotically as $c\sim t^{-1/2}$. However, if particles move one at a time, {\it i.e.}, asynchronous dynamics, there is diffusion in addition to the primary ballistic motion. This diffusive transport component permits two particles with the same velocity to annihilate. We now examine the asymptotic behavior of this reaction process by the series expansion approach. For this asynchronous ballistic reaction process, the rate operator can be written as \begin{equation} {\cal L}={\cal L}_{AO}+{\cal L}_{OB}+{\cal L}_{AB}+{{\cal L}_{AA} \over 2}+ {{\cal L}_{BB} \over 2}, \end{equation} where the last two terms account for the self-reaction of particles with the same velocity. We consider an initial state where each particle can assume a velocity $\pm v_0$ with equal probability. The initial state function then equals the sum of all $N$-particle configurations of velocities divided by $2^{N}$ (see Sec.~II). With our method, we obtain series for the concentration to $16^{\rm th}$ order (see Table). The two {Pad\'e } approximants that agree over the largest temporal range are $[7,7]$ and $[7,8]$ and the upper cutoff time is $t_{\rm max}=2.5$ (Fig.~2). We do not consider the next diagonal approximant $[8,8]$ since this approximant leads to a reduced temporal range of internal consistency. Thus, we perform the integration over the interval [1.5,2.5]. By maximizing $g(\alpha)$, the values of the optimal exponents are found to be $\alpha_{\rm opt}=0.716$ and 0.717 respectively. We conclude an estimated value of $\alpha_{\rm est}=0.72$, which yields a 4\% deviation when compared with the value $\alpha=3/4$ suggested by scaling arguments and by Monte-Carlo simulations. For reference, we outline the scaling argument which suggests that the decay exponent $\alpha$ should equal 3/4 [9]. In the asynchronous ballistic reaction, particles with velocity $+v_0$ perform a random walk with the bias in the $+x$ direction and similarly for particles with velocity $-v_0$. Hence, diffusion is introduced in addition to the ballistic motion. We denote the diffusion coefficient by $D$. This representation is useful since exact solutions are known for the two limiting cases, $D=0$ and $v_0=0$, and these can be incorporated to produce the general case result. We first write an ansatz for the concentration for $D,v_0\ne 0$ and match it to the known limiting solutions at the appropriate crossover times. The concentration in the general case must have the form \begin{equation} c(t)=c_0^{\beta} \left({1\over v_0t}\right)^{\gamma} \left({1 \over \sqrt{Dt}}\right)^{\delta}, \end{equation} that is, $c(t)$ is a geometric average of the fundamental concentrations $c_0$, $1/v_0t$ and $1/\sqrt{Dt}$. To ensure the proper dimensions of $c(t)$, the exponents $\beta$, $\gamma$ and $\delta$ must obey the constraint $\beta+\gamma+\delta=1$. In the case of a small drift [3], the exact asymptotic form of the concentration for $v_0=0$, $c_D(t)\sim(Dt)^{-1/2}$ holds for $\sqrt{Dt}>v_0t$ and matches Eq.~(20) at the crossover time $t_1=D/v^2$. Imposing $c_D(t_1)=c(t_1)$ and equating corresponding powers of the various coefficients yields equation $\gamma+\delta=1$. In the same way, the exact drift-limited solution [8] $c_v(t)\sim (c_0/v_0t)^{1/2}$ should coincide with Eq.~(20) at $t_2=1/Dc_0^2$. This condition leads to $\gamma=1/2$. Consequently, $c(t)$ has the general form \begin{equation} c(t)\sim \left({1\over Dv_0^2t^3}\right)^{1/4}. \end{equation} This suggested asymptotic decay was confirmed by simulation studies. Interestingly, the decay of the concentration in the general case is faster than that of the two limiting cases of $v_0=0$ and $D=0$. Moreover, although one can claim that the diffusive nature of the particles does not play a role for sufficiently large times, this scaling argument suggests that this is not the case. It is seen that the series expansion does give useful information about the long time behavior, despite the relatively short time domain where the concentration function is obtained. \subsection{Three species annihilation} In the $n-$species annihilation process, $n$ distinct species move diffusively and whenever two different species land on the same site an annihilation event occurs [10]. In one dimension the concentration was predicted to obey a power law decay $c\sim t^{-\alpha(n)}$ with the value $\alpha(n)=(2n-3)/4(n-1)$. This prediction appears to be confirmed by simulation studies. In the limiting case $n\to\infty$, one recovers the known $\alpha=1/2$ decay, since this case corresponds to single species annihilation. To develop a series expansion for this model we again define the rate operator and an appropriate initial state. We consider the initial state of all possible $N$ particle configurations divided by $3^N$. The rate operator now includes diffusion for each species and reaction for each pair of non-identical species, \begin{eqnarray} {{\cal L}}&=&{\cal L}_{AO}+{\cal L}_{OA}+{\cal L}_{BO}+{\cal L}_{OB}+ {\cal L}_{CO}+{\cal L}_{OC}+\nonumber\\&& {\cal L}_{AB}+{\cal L}_{AC}+{\cal L}_{BC} \end{eqnarray} with the definitions of Eqs.~(3) and (4). In this case we obtain the concentration to $13^{\rm th}$ order as listed in the table. Since the expansion is obtained to a relatively small order, we have relaxed the tolerance $\epsilon$ to $0.005$, so that the integration interval will exclude the range $t<1$. The ``best'' approximants here are $[6,6]$ and $[6,7]$ and the corresponding cut-off time is $t_{\rm max}=2.2$ (see Fig.~3). The optimal exponents are found by maximizing the fit function $g(\alpha)$ to be $\alpha_{\rm opt}=0.372$ and 0.380 respectively. Hence, we estimate the exponent by the average value $\alpha_{\rm est}=0.376$. This estimate is surprisingly close to the suggested $3/8$ value and has a 2\% deviation when compared with the MC simulation result of $0.37$ [10]. Note that the reduction of $\epsilon$ increases the discrepancy between the two optimal exponents. \subsection{Conclusions} We have presented a method for calculating a power series expansion in time of basic physical quantities, such as the concentration, in simple chemical reaction models. The series expansion approach has the advantage of being exact and being applicable to general classes of non-equilibrium models. A particular reaction process can be defined by a rate operator, and for simple reaction processes, this operator is specified by its action on only a small number of nearest neighbor pair configurations. We have also suggested an analysis technique, based on a direct application of the {Pad\'e } method to obtain the concentration over an extended temporal range. The asymptotic behavior of the series is determined by matching the power series to the two ``closest'' {Pad\'e } approximants. We have illustrated our general approach with two specific reaction models in section IV. Reasonable estimates for the asymptotic behavior were extracted from the time expansion. On the other hand, our analysis method was not entirely robust. For some processes, the exponent estimates converged slowly towards their long time limiting behavior. For example, a restricted asynchronous ballistic annihilation model, where self annihilation is excluded, yields poor estimates for the expected decay exponent given the $16^{\rm th}$ order expansion. Therefore, we conclude that series expansion is not guaranteed to be helpful. Theoretically, if one obtains enough terms the long time behavior will eventually emerge, but this may be beyond the reach of available computing resources. Although we applied the series expansion to spatially homogeneous systems, the method can be easily applied to study heterogeneous problems, such as the probability distribution of diffusing particles in the vicinity of a trap, and the nature of the reaction front in models of epidemic spread and biological waves, for example. \vskip .2in We thank S.~Redner for numerous useful discussions and for careful reading of the manuscript. We gratefully acknowledge ARO grant \#DAAH04-93-G-0021 for partial support of this research. Acknowledgement is also made to the Donors of The Petroleum Research Fund, administered by the American Chemical Society, for partial support of this research. \begin{thebibliography}{99} \bibitem{} M.~Doi, {\sl J. Phys. A} {\bf 9}, 1465; 1479 (1976). \bibitem{} P.~Grassberger and M.~Scheunert, {\sl Fortschr.\ Phys.} {\bf 28}, 547 (1980). \bibitem{} A.~A.~Lushnikov, {\sl Sov.\ Phys.\ JETP} {\bf 64} 811 (1986). \bibitem{} B.~Friedman, G.~Levine and B.~O'Shaughnessy, {\sl Phys. Rev. A} {\bf 46}, 7343 (1992). \bibitem{} R.~Dickman and I.~Jensen, {\sl Phys. Rev. Lett.} {\bf 67}, 2391 (1991); I.~Jensen and R.~Dickman, {\sl J. Stat. Phys.} {\bf 71}, 89 (1993). \bibitem{} S.~Song and D.~Poland, {\sl J. Phys. A} {\bf 25} 3913 (1992). \bibitem{} J.~L.~Spouge {\sl Phys. Rev. Lett.} {\bf 60} 871 (1988). \bibitem{} Y.~Elskens and H.~L.~Frisch, {\sl Phys. Rev. A} {\bf 31}, 3812 (1985). \bibitem{} E.~Ben-Naim, S.~Redner and F.~Leyvraz {\sl Phys. Rev. Lett.} {\bf 70}, 1890 (1993). \bibitem{} D.~Ben-Avraham and S.~Redner, {\sl Phys. Rev. A} {\bf 34}, 501 (1986). \end{thebibliography} \section*{Table Caption} {\bf Table 1} Expansion coefficients in the time power series for the concentration for the two representative reaction processes discussed in the text: (a) Ballistic two-species annihilation with self-reaction and (b) Three-species diffusive annihilation $A+B\to 0$, $A+C\to 0$ and $B+c\to 0$. \vskip 3pt \smallskip \hrule width 8.9cm\vskip 2pt\hrule width 8.9cm\vskip2pt \tabskip=1em plus2em minus.5em\halign to \hsize{\hfil#&\hfil#&\hfil#&\qquad#&\qquad#&\qquad#\cr &\hfill Ballistic\hfill&\hfill Three Species\hfill\cr $n\ $&\hfill$c_n$\hfill&\hfill$c_n$\hfill&\cr 0 & $0.500000000000000$& $0.333333333333333$&&&\cr 1 &$-1.000000000000000$&$-0.888888888888888$&&&\cr 2 & $1.500000000000000$& $2.074074074074074$&&&\cr 3 &$-1.666666666666667$&$-3.621399176954732$&&&\cr 4 & $1.416666666666667$& $4.949245541838134$&&&\cr 5 &$-0.935416666666667$&$-5.462094192958391$&&&\cr 6 & $0.468402777777778$& $4.928191840674693$&&&\cr 7 &$-0.152715773809524$&$-3.595704294012478$&&&\cr 8 &$-0.001765581044385$& $2.000766523547165$&&&\cr 9 & $0.048626063154228$&$-0.657130132013536$&&&\cr 10 &$-0.045346447200260$&$-0.157411205671068$&&&\cr 11 & $0.028467120232814$& $0.434194697238727$&&&\cr 12 &$-0.013827197289845$&$-0.337141906840114$&&&\cr 13 & $0.005071234556065$& $0.074818523203021$&&&\cr 14 &$-0.001031397263083$&&&&\cr 15 &$-0.003347397123149$&&&&\cr 16 & $0.000536792117891$&&&&\cr} \vskip 3pt\hrule width 8.9cm\vskip 2pt\hrule width 8.9cm \section*{Figure Caption} {\bf Fig.1} Approximating the time dependence of the concentration from its power series in the case of single species annihilation. (a) Shown are the exact solution (solid), the 25-term Taylor series (dashed) and the $[11,11]$ (dotted) and $[12,12]$ (dashed) Pad\'e approximants. (b) The fit function $1-g(\alpha)$ for the $[11,11]$ (solid) and $[12,12]$ (dashed) approximants. {\bf Fig.2} (a) Time dependence of the concentration for the ballistic annihilation model. Shown are the $[7,7]$ (solid), $[7,8]$ (dotted) and $[8,8]$ (dashed) Pad\'e approximants. (b) The fit function $1-g(\alpha)$ for the $[7,7]$ (solid) and $[7,8]$ (dashed) approximants. {\bf Fig.3} (a) Time dependence of the concentration in the case of three-species annihilation model. Shown are the $[6,6]$ (solid) and $[6,7]$ (dashed) Pad\'e approximants. (b) The fit function $1-g(\alpha)$ for the $[6,6]$ and $[6,7]$ approximants. \end{multicols} \end{document}