Velocity of front propagation in 1-dimensional autocatalytic reactions

 

C. Warren[‡], E. Somfai[*], and L.M. Sander[†]
Dept. of Physics, The University of Michigan,
Ann Arbor, MI 48109-1120, USA

 

Received 1 December 1999

 

 

We study a discrete model of the irreversible autocatalytic reaction A + B ® 2A in one dimension. Looking at the dynamics of propagation, we find that in the low-concentration limit the average velocity of propagation approaches v = q/2, where q is the concentration, and, in the high concentration limit, we nd the velocity approaches v = 1 - e^-q/2.

 

 

I Introduction

In this paper we study a discrete model of the irreversible autocatalytic reaction A + B ® 2A in one dimension. We analyze the velocity of reaction propagation as a function of the total concentration, q. This reaction is a simple example of a chemical reaction model, but we can also think of it as a representation of the spread of an infection. The A particles are infected, and we allow them to diffuse, i.e. perform a random walk, and when a 'sick' particle encounters a 'healthy' one (a B particle) then the B is instantly infected. We assume that the A's and the B's see each other only by infection. Otherwise they are independent random walkers which do not interact. The quantity of interest is the speed of the spread of the infection. If we start with one A on the extreme left of the system, we are asking for the time dependence of the location of the rightmost A.

Models of this type have evoked a good deal of interest [1-4] because of several unexpected, fascinating features of the process. As we will see, in the limit of small q the front propogation is dominated by fluctuation effects. What attracted particular attention [1] was the realization that, as a result, continuum modeling breaks down completely for this reaction.

Specifically, we expect that the mean concentration of A particles would be described by the Fisher-Kolmogorov equation [5], as we can see by writing a conventional reaction-diffusion equation:

Here k is a rate constant, D the diffusion constant, and we have used the fact that a+b = q. The last equation becomes the Fisher-Kolmogorov equation ¶_t u = ¶_xx u +u(1 - u) after changing variables. Now from the standard theory [6], the velocity of the front approaches v_c = 2 independent of initial conditions.

However, simulations [1] showed that the front velocity was linear in q for small q, and thus very much smaller than expected. This is understandable from the work of Brunet and Derrida [2] who pointed out that discreteness has an anomalously large effect on systems which obey the Fisher-Kolmogorov equation in the continuum limit. In references [2, 3] models were introduced which interpolated between the results of [1] and the Fisher-Kolmogorov equation via a very slow crossover. The models involved a very large density of particles with a small reaction rate. They found that the velocity depression was given by v ~ v_c - K / ln²(q) where K is a constant.

A question remains, however: What is the mechanism for the small velocity for q ® 0 in the original A + B ® 2A process? In a trivial model of independent random walkers, it is startling that there is any interesting dynamics. The total density at any point is clearly given by a Poisson distribution. It turns out however, that the front is not a typical point, and this is the key to the unexpected behavior.

We find that the small velocity is a giant fluctuation effect which goes qualitatively as follows. Suppose we assume that the velocity is a monotonically increasing function of q, and, for small q recall that there are large fluctuations in the local density. The front will move quickly through high density regions, and get stuck in the low density ones. Thus, on average, the motion will be dominated by configurations where the front is behind a gap in the distribution of B's. The front motion will be random, and not advance, as long as the rightmost A cannot convert a B. Our simulations and analysis in the next section support this picture.

Also in the spirit of studying the model in its own right we consider the case of large q for a parallel updating version of the process. Clearly the velocity must approach one lattice constant per unit time since there will always be a conversion at each step. However, the nature of this approach turns out to be tricky, and we are able to give only a partial analysis.

 

II Model and simulation results

Consider a 1-D lattice of length L [7] populated with random walkers randomly distributed with concentration q. The leftmost particle is of type A, and all of the other particles are of type B. The particles make random steps with parallel updating, i.e., all the particles move simultaneously. Any number of particles are allowed to occupy a site. If a B particle encounters or passes an A particle, it becomes an A particle. The rightmost A particle defines the propagation front, and we are interested in the velocity of this front.

As we pointed out above, the particles follow a simple random walk and thus are Poisson distributed. However, taking into account particle types and following the front, the distribution of particles at the front or near the front is not so distributed. Fig. (1) shows the average density of particles from simulation for various sites around the front for q = 0.2. Ignoring the fit for the moment, the density (conditioned on there being a front at i = 0), is depleted to the right of the front. At higher concentrations, Fig. (2) shows the probability distribution of particles at various sites around the front in a simulation with average concentration q = 2; Fig. (3) is for q = 4, and Fig. (4) for q = 8.

 

Figure 1. Mean number of particles for sites about the front from a simulation with an average of 0.2 particles/site (in the rest frame of the front).

 

 

Figure 2. Probability density of site i about the front from a simulation with an average of 2 particles/site (in the rest frame of the front).

 

 

Figure 3. Probability density of site i about the front from a simulation with an average of 4 particles/site.

 

 

Figure 4. Probability density of site i about the front from a simulation with an average of 8 particles/site.

 

Our simulations were performed in two different ways. For concentrations below q = 0.5, 200 walkers were simulated and lengths were adjusted accordingly, L = 200/q. Enough time steps were performed for each walker to walk on average halfway across the sample, T = L/q. To avoid initial transients calculations were also made with the first 40/q² time steps censored. For concentrations above q = 0.5, a 500-site lattice was simulated with 500 q walkers for 400 time steps, and the first 100 time steps were censored.

An earlier study [1] analyzed the velocity in an approximate fashion, using the Smoluchowski approach. In this method, centered in the rest frame of the front, B particles diffuse toward the front. The number density n follows the one-dimensional diffusion equation in the frame moving with velocity v,

where the diffusion constant D = 1. Assuming stationarity, the time derivative vanishes, and the boundary conditions n(±¥) = q, n¢(±¥) = 0 and n(0⁺) = 0 lead to:

We will see later that v µ q so that, as x ® 0⁺, n is order q². As shown by Fig. 1, this distribution agrees well with simulation. However, the analysis gives no way to find v.

 

II.1 Small q

In the low-concentration limit, q << 1, consider a region containing the front particle and the "second" particle, i.e., that nearest the front. (In case of more than one particle at i = 0 we arbitrarily declare one to be the front, and the other the second particle.) The size of the region will be ~ 1/q. Define a coordinate system in the rest frame of the front particle, whose position is defined to be at i = 0. The number density of the second particle at site i at time t is n_i(t). Note that since the front particle is treated separately and not included in n_i, the latter approaches the probability distribution for the second particle in the limit q ® 0.

The location of the second particle is the key to calculating the velocity. In fact, the only nonzero contributions to the average velocity occur when the second particle is one behind the front (i = -1) or on the front (i = 0). For all other positions, the front particle undergoes an unbiased random walk, and the front does not move on average. When i = -1, with proper renaming of particles, the front will move forward one step with probability 1/2, stay the same with probability 1/4 and move back one with probability 1/4. Thus, given i = -1, the average velocity is 1/2 - 1/4 = 1/4. When i = 0, the front will move forward one with probability 3/4 and move back one with probability 1/4, so the average velocity of the front is

From random walk dynamics, the motion of the second particle obeys a few simple rules:

· If the particle is at i = 0 at time t, n_-2(t + 1) = 1/2 and n₀(t + 1) = 1/2.

· If the particle is at i = -1 at time t, n_-1(t + 1) = 3/4 and n_-3(t + 1) = 1/4.

· Otherwise, if the particle is at i < -1 or i > 0, n_i+2(t + 1) = 1/4, n_i(t + 1) = 1/2 and n_i-2(t + 1) = 1/4. (Note that in these expressions the position of the second particle is defined with respect to the front at at time t + 1.)

For a stationary distribution, n_i(t) = n_i(t + 1) º n_i, so the rules lead to relationships between the densities:

(a) For positions away from the front, i < -2 or i > 2,

(b) For positions around the front

Equation (5) implies a linear behavior for the probabilities away from the front. We need to impose boundary conditions far from the front, because we have made the approximation of only one nonfront particle, which is only valid in the region 1/q around the front. We do this by matching to the continuum solution of Eq. (3). For (3) to order q, the density is simply q behind the front and 0 in front of the front. Thus, to first order in concentration q, the solution to the equations above is:

Thus v = q/2. Our simulation results, shown in Table 1, confirm this prediction.

 

Table 1. Simulation results for the velocity of front propagation for low concentrations.

 

Note that including the front particle, to first order, the stable total number density distribution is N_i = q for i < 0, N_i = 1 + q/2 for i = 0, and N_i = 0 for i > 0. That is we have average concentration to the left of the front and a depleted zone to the right, in agreement with simulations.

Note the bipartite nature of (5) -(6). We can consider a simpler "even-lattice" model in which only the even sites are populated and still get the same average velocity. This avoids the complicating factor of a site -1 particle passing the front particle and becoming the new front particle. For the even sites,

Maintaining the same overall density q would mean doubling the concentration at all of the sites. Looking at the form of (4) and taking into account differing number densities, we get equal contributions to the velocity from site 0 and site -1, and thus the same low velocity limit, v = q/2. However, this simplified "even-lattice" model does not approach the same limit as that of the "full-lattice" model in high concentration.

 

II.2 Large q

In the full-lattice model, for any concentration the probability of the front moving backward one site is

moving forward one is

and being stationary is

where p(n_i) is the probability of site i having n_i particles.

With stationarity, p(n₁(t), ¼, n_L(t)) º p(n₁, ¼, n_L), we can use the dynamics of the model to obtain similar relations between the site averages and other site moments. In the rest frame of the front, if the front moves forward, site i at time t is "fed" by sites i and i + 2 at time t-1. If the front is stationary, site i is fed by sites i - 1 and i + 1. If the front moves backward, site i is fed by i and i - 2. The average at site i is

where f_k is the number of particles that move forward from site k, b_k is the number of particles that move backward from site k, n_k = b_k + f_k, á·|+ñ means that the average is conditioned on the front moving forward, á·|0ñ means the average is conditioned on the front remaining stationary and á·|-ñ means the average is conditioned on the front moving backward.

For example, for site i = 0,

Note that around the front, specifically for i Î {-2, -1, 0, 1, 2}, the movement of the front reveals information about the number particles at site i.

Now consider q to be large. Assume for the moment that the distribution at all sites, including i = 0, is Poisson with mean q_i, as seen in the simulations. Although this must be true far from the front, it is not clear why it is true nearby. It introduces a discrepancy of order e^-q because at site i = 0 there is a probability e^-q that there will be no particles at site i = 0. This cannot be true because i = 0 is defined as the site of the rightmost A particle.

Neglecting this inconsistency for the moment, for large q, one may use equation (13) to derive a recursion relation:

With the assumption q_i º q, this equation collapses to a relation that approaches self-consistency to order qe^-.

A similar relation may be derived for the variance of the number of particles at site i,

With the Poisson assumption of mean and variance q, this relation simplifies to an expression that is self-consistent to leading order .

Using the Poisson approximation, the velocity is

However we can handle p(n₀ = 0) in a different way by using a conditioned or truncated Poisson distribution, p_T [8]. Specifically, p(n₀ = 0) = e^-q is truncated from the distribution and distributed to the other probabilities,

where p_P(k) is the ordinary Poisson distribution. Then the velocity is

The recursion relation discrepancy for the mean is identical to that of the Poisson distribution since everything is simply divided by 1-e^-q, and the leading order of the variance is again identical, . Fig. (5) shows that this agrees with the data better than Eq. (16). We have no explanation for this.

 

Figure 5. Mean front velocity as a function of particle concentration.

 

III Summary

We have seen that the behavior of v at low concentration can be traced to the depletion zone to the right of the front. Near the front the distribution of particles is very different from the Poisson distribution, and the motion of the front is dominated by the depletion. For large q the velocity is approximated by 1-e^-q/2 and the distribution is quite close to the truncated Poisson.

In the low concentration limit, to order q, the velocity is simply v = q/2 in the truncated Poisson model, which also is correct. However, the truncated Poisson model does not satisfy the recursion relations to order q and it lacks the depletion zone ahead of the front. As shown in Fig. (5), the truncated Poisson approach gives a reasonable approximation for v for any concentration, but this should be regarded as only a convenient interpolation.

We are grateful to Martin Barlow and Igor Sokolov for useful discussions. This work is supported by DOE grant DEFG-02-95ER-45546.

 

References

[1] J. Mai, I.M. Sokolov, and A. Blumen, Phys. Rev. Lett. 77, 4462 (1996)         [ Links ]

[2] E. Brunet and B. Derrida, Phys. Rev. E 56, 2597 (1997).         [ Links ]

[3] D. A. Kessler, Z. Ner and L. M. Sander, Phys. Rev. E 58, 107 (1998).         [ Links ]

[4] J. Riordan, C. Doering, and D.ben-Avraham Phys. Rev. Lett. 75, 565 (1995).         [ Links ]

[5] R. A. Fisher, Annals of Eugenics 7, 355 (1937); A. Kolmogorov, I. Petrovsky, and N. Piscounov, Moscow Univ. Bull. Math. A 1, 1 (1937).         [ Links ]

[6] J. D. Murray, Mathematical Biology (Springer, Berlin, 1993).         [ Links ]

[7] To simulate the dynamics of a large lattice, we actually used a ring lattice of L sites to handle the far boundaries, left and right. Random A walkers that walked off the left edge were added to the right edge as B walkers, and vice versa.

[8] V.K. Rohatgi, An Introduction to Probability Theory and Mathematical Statistics (Wiley, New York, 1976).         [ Links ]

 

 

[^‡] e-mail: warrencp@umich.edu
[^*] e-mail: lsander@umich.edu
[^†] Current address: Instituut Lorentz, Leiden, NL-2333 CA, Netherlands. e-mail: ellak@lorentz.leidenuniv.nl