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






