\documentstyle[aps,multicol,epsf]{revtex}
\begin{document}
\title{Inhomogeneous Two-Species Annihilation in the Steady State} 
\author{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 steady-state geometrical properties of the reaction
interface in the two-species annihilation process, $A+B\rightarrow 0$,
when a flux $j$ of $A$ and $B$ particles is injected at opposite
extremities of a finite domain.  By balancing the input flux with the
number of reactions, we determine that the width $w$ of the reaction
zone scales as $j^{-1/3}$ in the large flux limit, and that the
concentration in this zone is proportional to $j^{2/3}$.  This same
behavior is deduced from the solution to the reaction-diffusion
equation.  In the low flux limit, the concentration is almost
independent of position and is proportional to $\sqrt{j}$.  In the
latter case, the local reaction rate reaches maximum at the edges of the
system rather than at the midpoint.  When the two species approach at a
finite velocity, there exists a critical velocity, above which the
reactants essentially pass through each other.  Results similar to those
in one dimension are found in two- and three-dimensional radial
geometries.  Finally, we apply the quasistatic approximation to our
steady-state solution to recover the known time dependence for the
reaction zone width for the case of initially separated components with
no external input.
\end{abstract}
%\noindent INSPEC Numbers: 0250, 0560, 8220.
\begin{multicols}{2} 
A fundamental ingredient that controls the kinetics of the two-species
annihilation process, $A+B\to 0$, is the rate of reaction in the
interfacial region between a domain of $A$'s and a domain of $B$'s.  A
natural way to study this interface is to prepare a system with
initially separated components and then monitor the reaction rate as a
function of time, as first studied by G\'alfi and R\'acz [1]. This
type of geometry is also relatively more amenable to experimental
investigation [2] than the more extensively studied case of a
homogeneous, but random initial condition. In a mean-field
description, the kinetics is described by the reaction-diffusion
equations 
\begin{eqnarray}
{\partial c_A\over\partial t}&=&D_A\nabla^2 c_A -k c_A c_B,\nonumber\\ 
{\partial c_B\over\partial t}&=&D_B\nabla^2 c_B -k c_A c_B,
\end{eqnarray} 
where $D_i$ is the diffusion coefficient of species $i$ and $k$ is the
reaction constant.  For the case of initially separated components, the
initial conditions are $c_A(x,t=0)=c_0$ for $x>0$ and $c_A(x,t=0)=0$ for
$x<0$, and conversely for the concentration of $B$'s.  On the basis of
scaling arguments, and through solutions of the reaction-diffusion
equation, G\'alfi and R\'acz${}^1$ found that the width $w$ of the
reaction zone increases with time $t$ as $t^{1/6}$, and that the
reaction rate vanishes as $t^{-2/3}$.  Numerical simulations appear to
confirm this result in two dimensions [3,4] but different exponents
appear to hold in one dimension [5,6]. Related geometrical
properties of the reaction interface have also been studied for the
random initial condition [7].

We discuss here new results for the complementary situation in which
particles are confined to a finite $d$-dimensional bar with equal fluxes
of $A$'s and $B$'s injected at opposite ends of the system.  Our goal is
to understand the geometrical properties of the ensuing steady state.
Notice that in the steady state, we can redefine the concentrations by
$c_A\Rightarrow D_Bc_A$ and $c_B\Rightarrow D_Ac_B$ to map the problem
to the case of equal diffusion coefficients for the two species.  Thus
without loss of generality, we posit equal diffusion coefficients for
the two species.  This system is described by the steady state equations
\begin{eqnarray}
D\nabla^2 c_A&=&kc_Ac_B,\nonumber\\
D\nabla^2 c_B&=&kc_Ac_B,
\end{eqnarray}
with the boundary conditions 
\begin{eqnarray}
Dc'_A\big\vert_{x=L}&=&j,\quad Dc'_B\big\vert_{x=L}=0, \nonumber \\
Dc'_A\big\vert_{x=-L}&=&0, \quad Dc'_B\big\vert_{x=-L}=-j.
\end{eqnarray}
The reflecting boundary condition imposed on the $A$'s at $x=-L$, and on
the $B$'s at $x=L$ ensures that the reactants remain within $[-L,L]$.
We emphasize that these reaction-diffusion equations provide a
mean-field description of the reaction interface in which all spatial
fluctuations and the role of the spatial dimensionality is ignored.
Within this approximation, we shall determine the extent of the reaction
zone and the spatial distribution of the two species.

A rough estimate for these two quantities can be obtained from simple
heuristic arguments (Fig.~1).  Consider first the situation of a large
input flux, so that the concentration profile is linear near the domain
boundaries with the magnitude of the slope proportional to $j/D$.  If we
define the reaction zone as the region for which the concentrations of
both species are non-negligible, then the concentration in the reaction
zone should be of the order $jw/D$, where $w$ is the reaction zone
width.  Consequently, the number of annihilation events per unit time is
of order $kc_Ac_Bw$, which is obtained by integrating the reaction term
over the reaction zone.  This number should equal
the flux of particles $j$ entering the domain.  Therefore, balancing
these two rates gives
\begin{equation}
kc_Ac_Bw\sim k\Bigl({wj\over D}\Bigr)^2w\sim j,
\end{equation} 
or 
\begin{equation}
w_0\sim\Bigl({D^2\over jk}\Bigr)^{1/3}.
\end{equation}

Thus we conclude that the width of the reaction zone is proportional to
$j^{-1/3}$, and that the typical concentration in the reaction zone is
\begin{equation}
c_A, c_B\sim jw/D\sim\left({j^2\over{Dk}}\right)^{1/3}.
\end{equation}  
These results apply as long as the width $w$ is much less than (and
independent of) $L$.  From Eq.~(5), this corresponds to the flux $j$
being greater than a threshold value $j_0$, which scales as
$j_0=D^2/kL^3$.  When $j<j_0$, the rate balance argument fails,
indicative of a different scaling behavior for this case.  We will
outline an alternative approach for this limiting situation based on the
approximate solution to the reaction-diffusion equations (see below).

For a quantitative analysis, it is convenient to consider the
reaction-diffusion equations for the difference, $c_-=c_A-c_B$, and the
sum, $c_+=c_A+c_B$, from which the behavior of $c_A$ and $c_B$ can
easily be inferred.  From Eq.~(2), $c_-$ obeys $c_-(x)''=0$, with the
boundary conditions $Dc_-(x=\pm L)'=j$.  The corresponding solution is
$c_-(x)=j(x-x_0)/D$, with $x_0$ determined by the constraint that the
number of $A$'s equals the number of $B$'s in the domain, since the two
species are injected and annihilate at the same rate.  Thus, imposing
$\int_{-L}^{+L}c_-(x)dx=0$, gives $x_0=0$.  This condition also
specifies the location of the center of the reaction zone.

Using Eq.~(2), the sum of the concentrations obeys the equation
$Dc_+(x)''=2kc_Ac_B$, which can be rewritten in a closed form
by using $4c_Ac_B=c_+^2-c_-^2$ to yield,
\begin{equation}
c_+(x)''=\,{k\over 2D}\,\biggl[\,c_+^2-\Bigl({jx\over D}\Bigr)^2\,
\biggr],
\end{equation}
with the boundary conditions $Dc_+(x=\pm L)'=\pm j$.  Since $c_+(x)$ is
an even function of $x$, it is sufficient to consider $x>0$ only.
Notice also that $c_+\cong c_-$ for positive $x$ outside the reaction
zone, since the concentration of the minority species $B$ is negligible
in this region.  (In fact, the ansatz $c_+=|c_-|$ satisfies the Eq.~(7)
everywhere, except at the origin.)

For large flux, $j>j_0$, we attempt a power-series solution for $c_+(x)$
in the reaction zone, namely $c_+(x)=c_0+c_2x^2+c_4x^4+\cdots\ $.
Substituting a three-terms series truncation for $c_+(x)$ in Eq.~(7), we
find
\begin{equation}
c_+(x)\,=\,{w_0j\over D}\,\left[\,\left({4\over5}\right)^{1/3}+
\left({1\over10}\right)^{2/3}\Bigl({x\over
w_0}\Bigr)^2-\left({1\over 40}\right)\Bigl({x\over w_0}\Bigr)^4\,\right],
\end{equation}
for $|x|<w$, 
with $w_0=(D^2/jk)^{1/3}$.  In general, the coefficient of $x^{2n}$ in
the series is proportional to $j/(Dw_0^{2n-1})\sim j^{2(n+1)/3}$, with
the numerical prefactor rapidly converging to a limiting value as more
terms in the series approximation fo $c_+(x)$ are retained.

We may now determine the reaction zone width by equating $c_+$ to $c_-$
at $x=w$.  Up to the order given in Eq.~(8), this prescription yields
$w\cong 1.179w_0$, in agreement with our previous flux balance argument.
From Eq.~(8), we also find that the concentration in the reaction zone
is of the order of $c_0\propto (j^2/Dk)^{1/3}$.  Outside the reaction
zone, $c_+$ matches smoothly with $|c_-|$, as illustrated in Fig.~1.
It is worth noting that the above results for $c_0$ and $w$ can be
inferred by determining the scaled variables that bring Eq.~(7) into
dimensionless form.

It is instructive to determine the concentration of the
minority species outside the reaction zone by direct means.  For the
concentration of $B$'s at large $x$, we substitute $c_A=c_B+(jx/D)$ in
Eq.~(2) to give $Dc_B(x)''=kc_B(c_B+jx/D)$.  In the domain $x\gg w$,
$c_B(x)\ll jx/D$ and therefore the differential equation reduces to the
Airy equation,
\begin{equation}
c_B(z)''=z\,c_B(z),\qquad z=x/w_0.
\end{equation} 
The limiting case $x\gg w$ corresponds to $z\gg 1$, for which the
appropriate (decaying) solution is
\begin{equation}
c_B(x)\sim \,{w_0j\over D}\,\Bigl({w_0\over x}\Bigr)^{1/4}
\exp\big(-2(x/w_0)^{3/2}/3\big).
\end{equation}
In this expression, the prefactor has been chosen to match the small-$x$
power series representation for $c_+(x)$ when $x\cong w_0$.

Thus in the high flux limit, the concentration of $B$'s is approximately
$-jx/D$ for $x<-w$ and exponentially small for $x>w$ (Eq.~(10)).  In the
reaction zone itself, $|x|<w$, the concentration of $B$'s can be deduced
by writing $c_B$ as $(c_+-c_-)/2$ and using the above expressions for
$c_+$ and $c_-$ in the reaction zone.  This yields $c_B\propto j^{2/3}$.
The concentration of $A$'s is merely the mirror image of $c_B(x)$ about
$x_0=0$.

In the low flux limit, $j<j_0$, the approach presented above fails
to satisfy the boundary conditions.  However, for this case, we exploit
the fact that the concentration must be slowly varying in $x$.
Therefore we postulate that $c_+(x)=c_0+c_2x^2$ over the entire domain
$|x|\leq L$.  With this assumption, the boundary conditions on $c_+$ at
$x=\pm L$ now yields $c_2=j/2DL$.  By substituting this form into
Eq.~(7), we then find that the constant $c_0$ is given by
$c_0=\sqrt{2j/kL}$.  With these two leading terms, the remaining terms
in an infinite power series representation of the solution can be
evaluated, and it is straightforward to show that these higher-order
terms are negligible.  Thus, the approximate
form for $c_+(x)$ for $j\ll j_0$ is
\begin{equation} 
c_+(x)\cong\, \sqrt{{2j\over kL}}\, +\, \biggl({j\over 2DL}\biggr)\,
x^2,
\end{equation}
from which the concentration of each species can be reconstructed by
$c_{A,B}=\bigl(c_+(x)\pm jx/D\bigr)/2$ (Fig.~2).  Notice that the
magnitude of the variation in concentration $(\sim j)$ is small compared
to the value of the concentration itself $\big(\sim
\sqrt {j}\,\big)$.  More precisely, the ratio
$\big(c_+(L)-c_+(0)\big)/c_+(0)$ is proportional to $(j/j_0)^{1/2}$.  

A quantity which characterizes the spatial extent of the reaction is the
local reaction rate $R(x)\equiv
kc_A(x)c_B(x)=k\bigl(c_+^2-c_-^2\bigr)/4$.  Using the
previously-obtained expressions for $c_+$ and $c_-$, we find
\begin{equation}
R(x)\propto\cases
{(jw)^2-{\rm const.}\,\times(jx)^2, & $j>j_0$;\cr
         \cr
          j +{\rm const.}\,\times j^{3/2}x^2, & $j<j_0$;\cr}
\end{equation}
with both constants positive, so that the reaction rate exhibits a
unimodal to bimodal transition at $j\cong j_0$ (Fig.~3).  For large
flux, the reaction rate is sharply peaked around $x=0$, with a width
$w_0$, but in the small flux limit, particles are more likely to react
near the boundary of the domain rather than in the center.

In the limit of vanishing flux, we can, in principle, find the exact
expression for the reaction rate.  This limit is defined by the average
time between the injection of an $AB$ pair, $t_j=1/j$, being large
compared to the average time to diffuse across the domain, $t_D\propto
L^2/D$, so that the system is occupied by at most two particles.  To
solve for the kinetics of the two-particle system, we define $x_1=L-x_A$
and $x_2=L+x_B$, where $x_A(x_B)$ is the position of particle $A(B)$.
Both variables satisfy $0\leq x_1,x_2\leq 2L$, and a reaction occurs
whenever $x_1+x_2=2L$ (Fig.~4).  By this formulation, the interacting
system is mapped onto a two-dimensional random walk in the
first-quadrant triangle $x_1+x_2\leq 2L$, with reflecting boundary
conditions for $x_1,x_2=0$, and absorbing boundary conditions for
$x_1+x_2=2L$.  To specify the initial conditions, we assume that the
input flux at one boundary is uncorrelated with that at the other
boundary.  Therefore, the first particle attains its asymptotic uniform
distribution in the domain before the second particle is injected.  For
the equivalent two-dimensional problem, this leads to the initial
condition $p(x_1,x_2,t\!=\!0)=\big(\delta (x_1)+\delta (x_2)\big)/4L$.
The reaction rate corresponds to the total flux at the absorbing
boundary $x_1+x_2=2L$.  One way to find this flux is by 
exact enumeration of the probability distribution for the aforementioned
initial distribution of random walks.  This method confirms that the
reaction rate is indeed bimodal (Fig.~4).

An interesting and natural generalization is the situation where the two
species have a superimposed drift toward each other with velocity $v$,
as would be the case if the two species where oppositely charged and
placed in an external electric field.  In the high flux limit, $j>j_0$,
the density profile exhibits a boundary layer of order $D/v$ from the
reaction center, as long as $v$ is greater than $D/L$.  Outside the
boundary layer the concentration of majority species is $j/v$ and that
of the minority species is exponentially small.  Inside the boundary
layer the drift is negligible and the concentration profile merely
reduces to that of the zero drift case in the high flux limit.  This
approach is valid as long as the boundary layer $D/v$ is larger
than the reaction zone width $(D^2/jk)^{1/3}$, or $v\!<\!v_0$, with the
critical velocity given by
\begin{equation}v_0=(Djk)^{1/3}.\end{equation}  

In the high velocity limit, we find the amusing feature that the density
profile of each species is highly concentrated near the boundary
opposite the entrance point.  This result can be understood by first
rewriting the reaction-convection-diffusion equation in dimensionless
form.  Introducing $c_{A,B}=(v^2/Dk)\tilde c_{A,B}$ and $x=D{\tilde x}/v$
yields the scaled equations,
\begin{eqnarray}
{\tilde c''_A}+{\tilde c'_A}&=&{\tilde c_A}{\tilde c_B},\nonumber \\
{\tilde c''_B}-{\tilde c'_B}&=&{\tilde c_A}{\tilde c_B},
\end{eqnarray}
with the boundary conditions 
\begin{eqnarray}
&{\tilde c'_A}+{\tilde c_A}\big\vert_{{\tilde x}={\tilde L}}=
\epsilon,\qquad
&{\tilde c'_B}-{\tilde c_B}\big\vert_{{\tilde x}={\tilde L}}=0,\nonumber\\
&{\tilde c'_A}+{\tilde c_A}\big\vert_{{\tilde x}=-{\tilde L}}=0,
&{\tilde c'_B}-{\tilde c_B}\big\vert_{{\tilde x}=-{\tilde L}}=-\epsilon.
\end{eqnarray}
The scaled flux $\epsilon=\tilde j=(v_0/v)^3$ is vanishingly small when
$v>v_0$.  To leading order, then, the input is zero, and this implies
that the reaction term must vanish, as can be seen by integrating
Eq.~(14) over the length of the box and using the boundary conditions.
This fact suggests a perturbative approach, namely ${\tilde c_A}=\tilde
c_{A0}+\epsilon\tilde c_{A1}+\cdots\ \ $ and ${\tilde c_B}=\tilde
c_{B0}+\epsilon\tilde c_{B1}+\cdots\ \ $.  Solving for ${\tilde c_A}$ gives
${\tilde c_A}=c_0\exp(-\tilde x -{\tilde L})$ in the vicinity of $-\!{\tilde L}$ and ${\tilde c_A}$
vanishing elsewhere, and conversely for the $B$'s.  Thus the particle
are confined to a boundary layer whose width is of order unity in scaled
units where $\tilde L=Lv/D\gg 1$.  This zero-order solution can be used
in Eq.~(14) to obtain the corrections to this leading behavior.
Transforming back to the initial variables, we find, to leading order,
that the density is proportional to $v^2/Dk$ within a boundary layer of
order $D/v$, and that corrections are of the order $j/v$.  This gives a
reaction rate which is strongly peaked near the boundaries of the box
(Fig.~5).  In the low flux limit, we find that this same behavior still
applies as long as $v>D/L$.

We can also determine the nature of the reaction zone in a
higher-dimensional analog of the finite linear domain.  Consider the
radially symmetric situation in which the reaction takes place between
two concentric hyperspherical shells of different radii, with $A$'s
injected at the outer radius, and $B$'s injected at the inner radius.
The concentration depends on $r$ only, and the corresponding boundary
conditions are $Dr_A^{d-1}c_A(r_A)'=j$, $Dr_B^{d-1}c_A(r_B)'=0$, and
$Dr_B^{d-1}c_B(r_B)'=-j$, $Dr_A^{d-1}c_B(r_A)'=0$.  The equation for
the concentration difference is
\begin{equation}
{1\over r^{d-1}}\,{\partial \over \partial r}\biggl
(r^{d-1}\,{\partial c_-(r)\over \partial r}\biggr)=0,
\end{equation}
which, together with the boundary conditions, gives
\begin{equation}
c_-(r)=\cases{
{j\over D}\,\log(r/r_0) &
 with $r_0=\exp\,\left(\,{a^2\log a-b^2\log b\over a^2-b^2}
\,-\,{1\over 2}\,\right)$\quad $d=2$\cr
\cr
{j\over D}\,\left(\,{1\over r_0}-{1\over r}\,\right)&
with $r_0={2(a^3-b^3)\over 3(a^2-b^2)}$ \quad $d=3$\cr}
\end{equation}

As in the one dimensional case, the position of the reaction zone center
is determined by requiring that $\int_{r_B}^{r_A} c_-(r)\,r^{d-1}dr=0$.
When both radii are large with their difference remaining finite, then
$r_0\cong(a+b)/2$.  However, if the inner radius is much smaller than
the outer radius, then $r_0\cong a/\sqrt{e}$ for $d=2$ and $r_0\cong
{2a/3}$ for $d=3$.  The procedure for finding the individual
concentrations of each species closely follows that used for the one
dimensional case.  In the high flux limit, $c_+$ is simply equal to
$|c_-|$ outside the reaction zone, and $c_+$ has a power series
representation (which is not necessarily even about $r_0$) within the
reaction zone.  Following the same method of analysis as in one
dimension, we find the same scaling behavior, namely, $w\sim j^{-1/3}$
and $c(r=r_0)\sim j^{2/3}$.  Similarly, in the low-flux limit, we find
that the concentration scales as $c\sim \sqrt{j}$, and that the
magnitude of the spatial variation of $c$ is of order $j$.

By applying the quasi-static approximation to our description of the
steady-state, we can also determine the time dependence of the width of
the reaction zone in the G\'alfi-R\'acz problem, where the two components
are initially separated and with no external source of particles.  The
basis of this approximation is to neglect the time derivative in the
diffusion equation and shift this time dependence to a moving boundary
condition [8]. For the G\'alfi and R\'acz problem, the growing
depletion layer is the source of the moving boundary which produces the
requisite time-dependent flux.  For the initial conditions
$c_A(x,t\!=\!0)=c_0H(x)$ and $c_B(x,t\!=\!0)=c_0H(-x)$, with $H(x)$ the
Heaviside step function, and taking $D_A=D_B=D$ in Eq.~(1), then
$c_-(x,t)$ satisfies $\dot c_-=Dc_-''$, with the corresponding initial
condition $c_-(x,t\!=\!0)=c_0\bigl(2H(x)-1\bigr)$.  The solution of this
equation is
\begin{equation} 
c_-(x)=c_0\,{\rm erf }\biggl({x\over 2\sqrt{Dt}}\biggr),
\end{equation}
with the error function defined by ${\rm erf}(y)=
{2\over{\sqrt{\pi}}}\int_0^y \exp(-u^2)du$.  This solution leads to a
depletion zone of the order of $\sqrt{Dt}$.  The (time-dependent) flux
into the reaction zone is thus given by $j(t)\sim Dc_-(x\!=\!0,t)'\sim
Dc_0/\sqrt{Dt}$.  We now use this flux in the steady-state expression
for $w$ (Eq.~(5)).  This yields $w\sim (D^3t/k^2c_0^2)^{1/6}$, in
agreement with the G\'alfi and R\'acz result.  {\it A posteriori}, using
the high-flux limit is appropriate, since the flux $j\sim t^{-1/2}$ is
much greater than $j_0\sim L^{-3}\sim t^{-3/2}$.

In conclusion, we have investigated the geometrical properties of the
reaction zone for two-species annihilation, when particles of each
species are injected at a fixed rate $j$ from opposite edges of the
system.  Our solution to the reaction-diffusion equation provides a
mean-field description of a one-dimensional reaction zone geometry.  In
the limit of large flux, the reaction zone width is proportional
to $j^{-1/3}$, and the concentration within this zone is of order
$j^{2/3}$.  In the opposite limit, the concentration assumes a nearly
constant value of order $j^{1/2}$ throughout the system, in which the
maximal spatial variation is of order $j$.  The former case leads to a
sharply localized reaction rate near the domain center, while in the
latter case, the reaction rate is maximal near the extremities of the
domain.  Our approach is easily extended to a variety of potentially
relevant situations.  Particularly noteworthy is the case of a finite
superimposed drift of the reactants towards each other.  When the drift
velocity exceeds a threshold value, there is a strong tendency for the
reactants to pass through each other, a potentially destabilizing
mechanism of the reaction interface.

We thank M. J. Stephen for helpful discussions.  We also gratefully
acknowledge grants from the ARO and NSF for partial support of this
research.

\section*{References}

\begin{thebibliography}{99}

\bibitem{} L.~G\'alfi and Z.~R\'acz, {\sl Phys. Rev. A} {\bf 38}, 3151 (1988).
\bibitem{} Y.-E. L. Koo, L. Li, and R. Kopelman 
{\sl Mol. Liq. Crit.} {\bf 183}, 187 (1990);
Y.-E. L. Koo and R. Kopelman, {\sl J. Phys. A} {\bf 65} 893 (1991).
\bibitem{} Z.~Jiang and C.~Ebner, {\sl Phys. Rev. A} {\bf 42}, 7483 (1990).
\bibitem{} H. Taitelbaum, S. Havlin, J. E. Kiefer, B. Trus, and G. H. Weiss,
{\sl J. Stat. Phys.} {\bf 65}, 873 (1991).
\bibitem{} S. Cornell, M. Droz, and B. Chopard, {\sl Phys. Rev. A} 
{\bf 44}, 4826 (1991).
\bibitem{} M. Araujo, S. Havlin, H. Larralde, and H. E. Stanley, preprint.
\bibitem{} F. Leyvraz and S. Redner, {\sl Phys. Rev. Lett.} 
{\bf 66}, 2168 (1991)
; S. Redner and F. Leyvraz, {\sl J. Phys. A.} {\bf 65} 1043 (1991).
\bibitem{} S.~Redner and D.~ben-Avraham, {\sl J. Phys. A.} {\bf 23} L1169 
(1990).
\end{thebibliography} 

\section*{Figure Captions}

\noindent{\bf Fig. 1.} Reaction geometry in the steady state.  $A$
particles are injected at the right edge and $B$'s are injected at the
left edge, both at the same rate $j$.  Shown are the spatial
distribution of $A$'s and $B$'s (dashed), and the sum of their
concentrations $c_+(x)$ (solid) in the limit of large flux, $j=100$,
based on the numerical solution to the steady-state reaction-diffusion
equation for a system with $L=k=D=1$.

\noindent{\bf Fig. 2.} Concentration of $B$'s in the low flux limit
for the case $j=0.01$.  The order of magnitude and the variation in
this concentration over the extent of the domain are indicated.

\noindent{\bf Fig. 3.} Spatial variation of the local reaction rate
$R(x)=kc_A(x)c_B(x)$ in the case of high flux, $j=100$ (solid), and
low flux, $j=0.01$ (dashed).  The latter data has been multiplied by a
factor of 10,000.

\noindent{\bf Fig. 4.} Equivalence between two annihilating random
walks in a finite interval of length $2L$, and a single
two-dimensional random walk in a triangular domain with reflecting
boundary conditions for $x_1=0$, $x_2=0$, and absorbing boundary
conditions for $x_1+x_2=2L$.  Also shown is a plot of the flux
reaching this boundary, or equivalently, the local reaction rate
$R(x)$, as a function of position along the boundary.  This data is
based on exact enumeration inside a triangle of base 40 after 2000
time steps.

\noindent{\bf Fig. 5.} The density profile in the high velocity
limit. Shown are the densities as evaluated by exact enumeration with
the following parameters $L=30$, $v=0.2$, $j=10^{-4}$ and $D=k=1$.
\end{multicols}

\end{document}