Study of a bimolecular annihilation process for coarsening reactants

Study of a bimolecular annihilation process for coarsening reactants A. M. Lacasta Departament de Fı ´sica Aplicada, Universitat Polite `cnica de Catalunya, Av. Dr. Gregorio Maran ˜ o ´n 50, E-08028 Barcelona, Spain F. Sague ´s Departament de Quı ´mica Fı ´sica, Universitat Barcelona, Av. Diagonal 647, E-08028 Barcelona, Spain I. M. Sokolov Theoretishe Polymerphysik, Universita ¨t Freiburg, Hermann-Herder-Strasse 3, D-79104 Freiburg i. Br., Germany J. M. Sancho Departament d’Estructura i Constituents de la Mate `ria, Universitat de Barcelona, Av. Diagonal 647, E-08028 Barcelona, Spain ~Received 21 August 1998; accepted 11 December 1998! We consider the spatio-temporal behavior of the A1B→0 bimolecular reaction in a system where both reactants tend to segregate into separated phases. Our study is based on the numerical solution of a pair of reaction–diffusion equations appropriate to capture the underlying coarsening dynamics. The interplay between reaction and coarsening leads to a complex pattern of reactants spatial distribution. At short/intermediate times two distinctive dynamical regimes are seen in the decay of overall concentration and droplet number and the behavior of droplet radii. © 1999 American Institute of Physics. @S0021-9606~99!50911-2# I. INTRODUCTION The basic bimolecular reaction A1B→0 was extensively studied during the past decade mostly due to its unusual, notoriously slow long-time fluctuation-dominated kinetics.1–4 This regime is due to the spontaneous separation of reactants and makes the A1B→0 scheme a very interesting model, e.g., for the investigations of effects of mixing procedures applied to diffusion-controlled reactions which are otherwise slow and ineffective.5,6 The reaction schemes considered in such a context are normally oversimplified because no particles interactions ~except those leading to the reaction!are considered. However, in many situations such interactions can be of great importance. In the case of screened electrostatic repulsion of particles of the same kind,7–9 the interactions lead to faster dissolution of inhomogeneities ~domains!. Contrarily, when the particles of the same kind attract, the inhomogeneities are stabilized. In this last case, corresponding to the coagulation–annihilation scheme considered in Refs. 10–12, different kinetic regimes can be observed. Thus, the overall kinetics depends strongly on the properties of the stabilized aggregates ~domains!, especially on the relation between their size and mobility. In these last situations, the aggregation process is normally considered to be irreversible, so that the domain size can decrease only if this aggregate reacts with one of the opposite kind. On the other hand, in many systems the reversibility of the aggregation process must be taken into account. Thus one can imagine a situation where large domains ~being rather immobile!dissolve ~‘‘evaporate’’!or grow ~‘‘condensate’’!rather than react, whereas their constituent particles, having higher diffusivity, react rather effectively in the bulk ~‘‘gas’’!phase. This problem is addressed in the present article where we adopt a continuous picture and describe the stabilization of domains ~‘‘droplets’’!using the models common in phase separation dynamics. Thus, the reactant’s organization is described in terms of a model for an inhomogeneous system consisting of separated droplets of the two species coexisting with a noncondensed phase where they come into contact and irreversibly annihilate. The reaction is taken into account by an additional term leading to the concentration decay. This method of description of the system exhibiting reactions and phase coarsening parallels that of Refs. 13–18, although the particular reaction and the objectives of our study are quite different. On the other hand, our study here is related to one of Ref. 19, where a phase separation modelization was used to create initial reactants distribution but was switched off after the reaction started. II. FORMULATION OF THE PROBLEM: MODEL AND SIMULATION ALGORITHM We will address our problem within a continuousmedium description, with the concentrations of reactants depending on space and time. The dynamics of the system is modeled by a pair of Cahn–Hilliard equations for each species,20–22 supplemented with a reaction term which couples them, ] ] tcA~r,t!5D¹2~2cA~122cA!~12cA!2¹2cA!2KcAcB, ~1! ] ] tcB~r,t!5D¹2~2cB~122cB!~12cB!2¹2cB!2KcAcB. Using the standard notation, cAand cBstand for the local dimensionless concentrations of the two species, DacJOURNAL OF CHEMICAL PHYSICS VOLUME 110, NUMBER 13 1 APRIL 1999 64580021-9606/99/110(13)/6458/5/$15.00 © 1999 American Institute of Physics Downloaded 08 Dec 2008 to 147.83.195.22. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp counts for the diffusion coefficient and Krepresents the constant reaction rate. All these parameters must also be understood as dimensionless. In the absence of the reaction term (K50), the two above equations are dynamically decoupled. Each species follows an independent coarsening process by growing domains rich in each component Aor Bwith cA,B.1, immersed in a background of cA,B.0. Larger droplets grow at the expense of smaller ones following a condensation–evaporation mechanism. This coarsening process is driven by the surface tension of the droplets. In the Cahn–Hilliard dimensionless picture this last quantity can be evaluated via s 5 E dx S dc~x! dx D 2 5 A 2 6,~2! where c(x) stands for the one-dimensional profile ~taken along the radial direction!of single hump-like droplet. It is well known that in this case the domain growth is characterized by a unique time-dependent length L(t) which asymptotically follows the power-law L(t);t1/3 ~Lifshitz– Slyozov behavior!.23 The introduction of the reaction term leads to the annihilation of both species in the regions where they come into contact. Thus, we are facing a competition process between the segregation dynamics that tries to form and keep growing the droplets of Aand B, and the chemical reaction that is continuously fading them. We note here that the typical fluctuation-dominated behavior stems from the disordered initial conditions. Other ~partially correlated!initial conditions could lead to a different behavior.19 Thus, to precisely characterize the simulations, we now have to describe the initial preparation of the system. The simulations start from randomly mixed configurations of droplets of reactants. To prepare this initial condition, we first integrate a single Cahn–Hilliard equation ] ] tc~r,t!5D¹2~2c~122c!~12c!2¹2c!,~3! with the initial condition prescribed by assigning to each point a concentration c(r,0)5c01 h (r) where c0is the mean concentration and h is a uniformly random number in the range @20.05,0.05#. Following Eq. ~3!, randomly located domains or droplets with c.1 form and grow in a bulk with c.0. At some time t0, the system has N0droplets of mean radius R0. The values of N0and R0depend on the mean concentration c0and time t0. Then we randomly assign the labels Aand Bamong the N0droplets and enforce if necessary the condition of equal number of droplets for both species. Obviously the droplet configuration created according to the Cahn–Hilliard dynamics as explained above leads to some dispersivity in the distribution of the number N0and of the radii R0of the droplets. To facilitate the interpretation of the different dynamical regimes of this reactive mixture, the results that are going to be presented below correspond to a distribution of equal size droplets. This naturally requires a further simple manipulation of the initial condition. However, simulations have also been conducted without this requirement, i.e., with droplets of random radii in addition to their random position. They show dynamical patterns very much like those reported here which will be commented on below. In any case, conditions to have a strictly stoichiometric mixture are always checked for any initial condition. Just before ending this section, some more technical comments are worth making concerning the numerical details employed in our simulations. We have numerically integrated Eqs. ~1!and ~3!in a square lattice L2of size L 5256 and periodic boundary conditions. We have used Euler’s method with spatial step Dx51 and time step Dt 50.02. The results have been averaged over 10 runs corresponding to different realizations of the initial condition. In all that follows we restrict ourselves to the case of very fast reaction, so we take K@D. III. NUMERICAL RESULTS Let us first consider the most important kinetic variable, i.e., the spatially averaged concentration of the reactants. Needless to say that, under stoichiometric conditions, we have ^ cA(r,t) & 5 ^ cB(r,t) & 5 r (t) for all time. The behavior of r (t) obtained numerically for a characteristic initial condition as described above is shown on a linear scale in Fig. 1. Simulations with other initial conditions ~parametrized in terms of any pair of variables chosen within the set r 0,N0, and R0) show that r (t) may slightly depend on them. In all cases considered, the temporal pattern is well represented by an initial quasilinear decay followed by a much slower reaction regime, with a well-defined crossover region in between. A useful and common strategy to interpret the different dynamic regimes in the general context of diffusion controlled annihilation processes consists in analyzing, in parallel to the behavior of the mean concentrations, the morphology and distribution of the reactants and reaction zones ~regions of high values of cAcB) in the reactive mixture. Such a distribution is shown in Fig. 2 for the same simulation conditions as Fig. 1. Representative snapshots correspond to different stages of the dynamical process studied. At very short times the reaction zones tend to concentrate in the vicinity of the droplets boundaries. This rather spurious regime is quickly followed by a fast reorganization of the reaction zones which FIG. 1. Time dependence of the concentration of the reactants, r , for parameters D50.1 and K510, corresponding to initial droplets with equal radii ~solid line!and initial droplets with distributed radii ~dotted line!. The dashed line corresponds to the pure diffusion-controlled case for the same initial condition and parameters than the solid line case. 6459J. Chem. Phys., Vol. 110, No. 13, 1 April 1999 Lacasta et al. Downloaded 08 Dec 2008 to 147.83.195.22. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp transform into a labyrinthic morphology of corrugated lines which are in close contact with a major fraction of the highly reactive droplets and separating those of unequal species. In this way, clusters of droplets of both species are progressively being formed. As time goes on, the reaction zones undergo a slow coarsening process due to the collapse and the evaporation of the droplets lying closer to the reaction zones. The remaining droplets tend more and more to concentrate far away from the reaction zones and to occupy the inner regions of the surviving clusters which become more and more rarified. The two regions in Fig. 1 can be roughly associated with this pair of distinctive regimes. Roughly speaking, the first stage represented by the quasilinear decay in the plot of r (t) corresponds to the build up of dense clusters with a high proportion of reactive droplets ~close enough to the reaction zones!. Contrarily, the second stage, with rather rarified clusters of droplets, corresponds to the slowing down of the decay rate evidenced in Fig. 1. After evaporation of almost all droplets, this last regime would eventually terminate. Figure 1 also shows the results obtained for droplets with a distribution of radii, exhibiting the same two regimes. Nevertheless, due to the dispersion of the radii, the first linear regime is shorter. In Fig. 1 we have also plotted for comparison the results obtained with the pure diffusion–reaction model,19 for the same droplets-like initial condition, showing a much faster dynamics. The final stage of the reaction would take place within a system without droplets and would correspond to one of the possible regimes considered in Ref. 19: it would either show a Zeldovich-type behavior or an exponential-like decay as a signature of the unavoidable finite-size effects. We chose the parameters of our simulations in a way to postpone the onset of these ~uninteresting!regimes. Two additional dynamical variables enable us to get a more comprehensive picture of the complex temporal behavior of such a reactive mixture. The first one corresponds to the evolution of the number of droplets N(t) and is represented in Fig. 3. Complementarily, we show in Fig. 4 a plot of the averaged droplet size R(t), calculated from the radius of gyration, RG(t), as R5 A 2RG. Again in both plots, distinctive features of the previously identified dynamical crossover are easily identified. In relation to Fig. 3, such a crossover connects a plateau of practically constant numbers of droplets with a monotonous and progressively damped decay law. The behavior is somewhat more intriguing with respect to the average radius: it first monotonically and almost linearly decreases but, after the time corresponding to the change from the plateau to a decay region in Fig. 3, starts to increase. Thus, the first stage of the reactive process is characterized by a practically constant number of continuously shrinking droplets and corresponds to the stage of active annihilation process involving a major fraction of very reactive droplets ~those relatively close to the reaction zones!. Contrarily, the constant irreversible decay in the number of droplets, with raising averaged radius, would be associated to the rather inactive stage of the annihilation process which involves progressively more and more rarified clusters of the reactive species. FIG. 2. Droplets and distribution of reaction zones depicted for the solid line case of Fig. 1, corresponding to times: ~a!t50, ~b!t51000, ~c!t 53000, and ~d!t55000. FIG. 3. Number of droplets versus t, for the same case as Fig. 2. FIG. 4. Mean radius versus t, for the same case as Fig. 2. 6460 J. Chem. Phys., Vol. 110, No. 13, 1 April 1999 Lacasta et al. Downloaded 08 Dec 2008 to 147.83.195.22. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp IV. INTERPRETATION OF THE DYNAMIC REGIMES As already mentioned, the initial kinetic regime is associated with the existence of stable reactive droplets, immersed in an almost stationary network of reaction zones. At this point it is worth emphasizing that at any time the stability of the droplets is maintained by their surface tension. The number of droplets remains almost constant in time ~up to the completion of this dynamical regime!while their radius continuously decreases. Almost each droplet is situated near some part of the reaction zone, so that this dynamical regime admits a quite straightforward analysis in terms of a diffusive controlled evaporation of a single isolated droplet in the vicinity of an adsorbing boundary. No signature of cooperative effects would be thus exhibited during this first stage. Needless to say, in our reactive system such an adsorbing boundary would correspond to the reaction zone closest to a given droplet. The problem of evaporation of a circular droplet of radius R(t) evolving under the Cahn–Hilliard dynamics @Eq. ~3!# and concentric to an adsorbing boundary located at a distance lfrom its center can be solved analytically. Use is made of a quasistatic approximation that permits us to write down a Laplacian equation for the concentration in the bulk. After invoking a sort of Gibbs–Thomson, curvaturemediated equilibrium condition at the droplet’s boundary, the steady state concentration profile c(r) is obtained. The rate of evaporation, expressed in terms of the temporal derivative of R(t), is then readily established in terms of the diffusive flux through the boundary. After computing R(t), the temporal evolution of the mean density of the isolated droplet system is then easily obtained. The whole analytical procedure is outlined in Appendix A. We simply make use here of some approximate final form of the whole dynamics, d r dt }2D s l3,~4! where the proportionality factor is a quantity which slightly increases with time @see Eq. ~A4!#. In any case, the conclusion is that a quasilinear decay of the mean density is expected due to the presence of the adsorbing boundary. Specific simulations of this behavior were conducted with different values of R0and l. Such parallel simulations enabled us to consistently check the above theoretical predictions. The quasilinear decay rate observed in our simulated reactive mixtures admits then a simple interpretation according to the result just obtained. To make this comparison a little bit more quantitative, the remaining task consists of evaluating the effective parameter lto be used in our phaseseparated reactive mixture. Such an effective parameter might be estimated from the computation of the correlation function CAB(r,t)5 ^ cA(r8,t)cB(r1r8,t) & r8/ ^ cA &^ cB & . Needless to say such a description of the whole reactive mixture cannot be maintained for all times and we can reasonably predict that it will break down as soon as larger and more rarified clusters of droplets will be progressively formed. These clusters will contain droplets situated closer to each other than to the absorbing boundaries, thus giving rise to strong cooperative effects. This leads to a completely different dynamical regime which is discussed as follows. The behavior we are facing at such an intermediate stage of reaction is of a cooperative nature involving not individual droplets but rather sets of droplets coupled through the bulk and enclosed by absorbing boundaries. In this regime, a great role is played by the peculiarities of the phase separation dynamics underlying Eq. ~3!, that did not show up at all in the earlier dynamics of the reactive process. In short, the nature of the dynamics corresponds to a competitive evaporation problem. Some results for the dynamics of a system consisting of droplets of the same type, enclosed by absorbing boundaries and corresponding to particular configurations of droplets ~for instance, disposed in arrays!have already been obtained very recently.24 The appearance of a droplet clusterization process in evaporating ~growing!droplets has been studied in Ref. 25 with a much more complicated and probably realistic modelization for its relevance in the analysis of evaporation and combustion of multicomponent ~blended!liquid fuels. At this point we recall that when entering this stage of the reaction, the system is organized in rather big and rarified clusters enclosed by the reaction zones and composed of droplets of the reactive species, most of them being situated quite far away from the boundaries of the clusters. The key idea then, is that the different droplets evaporate and in this way feed the reaction zones ~adsorbing boundaries!not independently but in a cooperative manner via exchanging matter through the bulk of the clusters. The evolution of any droplet of a particular species strongly depends on the concentration of that species in the bulk of the cluster to which the droplet belongs. For each such concentration a critical radius exists, so that the larger droplets grow while the smaller ones dissolve. The evaporation of a ~subcritical! droplet is a rather fast process, producing a considerable amount of matter diffusing in the bulk of the cluster. The matter set free during such an evaporation process diffuses not only into the direction of the reaction zones but everywhere within the cluster. This matter can be partially adsorbed by neighboring droplets. The radii of these droplets increase and in this way such droplets get transiently stabilized. The final outcome of such a mechanism is that larger and larger droplets tend to grow and become more and more stable while the clusters to which they belong become less and less populated. We note here that for the configuration of nonreactive droplets considered in Ref. 24, such a competitive mechanism is so efficient that it can lead to the fact that droplets which are further away from the boundaries evaporate earlier that those in the boundary’s vicinity. Furthermore for these particular configurations, the average radius also increases after some transient regime. Results found for this competitive dynamics would apply to the reactive mixture considered here. In this sense, preliminary results extending the situation considered in Ref. 24 to random configuration of droplets show a slower that linear decay of the concentration, likely to be diffusion controlled with a very small density in the diffusive bulk phase,26 consistent with Fig. 1. 6461J. Chem. Phys., Vol. 110, No. 13, 1 April 1999 Lacasta et al. Downloaded 08 Dec 2008 to 147.83.195.22. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp V. CONCLUSIONS We considered the bimolecular reaction A1B→0 between reactants undergoing phase separation. Initially the reactants are contained in form of equal size droplets; the overall amount of reactants correspond to stoichiometric conditions. By numerically solving the corresponding equations, we found that at short/intermediate times two dynamical regimes divided by a pronounced crossover region are clearly seen. The first regime corresponds to a reaction controlled by evaporation of single droplets towards the reaction zone; the second one corresponds to a cooperative process of redistribution of matter due to condensation within the clusters of droplets of the same kind. The mean droplet size decreases during the first regime and starts to grow with the onset of the second one. On the other hand the number of droplets during the first regime is almost constant and starts to decrease in the second. None of these regimes can be associated with standard, Ovchinnikov–Zeldovich, fluctuation-dominated kinetics. ACKNOWLEDGMENTS This research was supported by an ACCIONES INTEGRADAS-DAAD grant, by DFG, by Direccio ´n General de Investigacio ´n Cientı ´ficayTe ´ cnica ~Spain!under projects PB96-0241 and PB96-1001, and by the Comissionat per Universitats i Recerca de la Generalitat de Catalunya. We also acknowledge Fundacio ´Catalana per a la Recerca for computing support. The discussions with Professor A. Blumen are gratefully acknowledged. APPENDIX A: We study here the case of a spherical droplet of radius R(t), evolving under the dynamics of the dimensionless Eq. ~3!, but with a spherical absorbing boundary at fixed distance l. The stability of the droplet is maintained by its surface tension s . The concentration in the bulk can be evaluated analytically in the quasistatic approximation ¹2c(r)50, with the absorbing condition c(l)50 and the Gibbs– Thomson equilibrium condition at the interface c(R) 5 s /2R. For two dimensions, the concentration is given by ~A1! c~r!5 s 2R ln~r/l! ln~R/l!for R,r,l, c~r!51 for r,R. Due to the presence of the absorbing boundary, the droplet evaporates at a rate given by23 dR dt 5Ddc dr U r5R .~A2! From Eqs. ~A1!and ~A2!, the evolution of the mean density r (t)5(R(t)/l)2can be evaluated as d r dt 52 D s l2 1 Rln ~l/R!.~A3! Now, by introducing the ratio a5l/R, Eq. ~A3!can be written as d r dt 52 D s l3 a lna.~A4! During the dynamic evolution aincreases with time but, for the typical values of our simulations, the factor a/lna changes slightly. Hence, a quasilinear behavior is expected. 1A. A. Ovchinnikov and Ya. B. Zel’dovich, Chem. Phys. 28, 215 ~1978!. 2D. Toussaint and W. Wilczek, J. Chem. Phys. 78, 2642 ~1983!. 3K. Kang and S. Redner, Phys. Rev. Lett. 52, 955 ~1984!; Phys. Rev. A 32, 435 ~1985!. 4G. Zumofen, A. Blumen, and J. Klaften, J. Chem. Phys. 82, 3198 ~1985!. 5I. M. Sokolov and A. Blumen, Phys. Rev. Lett. 66, 1942 ~1991!. 6R. Reigada, F. Sague ´s, I. M. Sokolov, J. M. Sancho, and A. Blumen, Phys. Rev. Lett. 78, 741 ~1997!. 7I. M. Sokolov and A. Blumen, Europhys. Lett. 21, 885 ~1993!. 8G. Oshanin, I. M. Sokolov, P. Argyrakis, and A. Blumen, J. Chem. Phys. 105, 6304 ~1996!. 9V. N. Kuzovkov, E. A. Kotomin, and W. von Niessen, Phys. Rev. E 54, 6128 ~1996!; J. Chem. Phys. 105, 9486 ~1996!. 10I. M. Sokolov and A. Blumen, Phys. Rev. E 50, 2335 ~1994!. 11P. L. Krapivsky and E. Bennaim, Phys. Rev. E 53, 291 ~1996!. 12J. A. N. Filipe and G. I. Rodgers, Phys. Rev. E 54, 1290 ~1996!. 13S. C. Glotzer, D. Stauffer, and N. Jan, Phys. Rev. Lett. 72, 4109 ~1994!. 14S. C. Glotzer and A. Coniglio, Phys. Rev. E 50, 4241 ~1994!. 15S. C. Glotzer, E. A. Di Marzio, and M. Muthukumar, Phys. Rev. Lett. 74, 2034 ~1995!. 16S. Toxvaerd, Phys. Rev. E 53, 3710 ~1996!. 17J. J. Christen, K. Elder, and H. C. Fogedby, Phys. Rev. E 54, R2212 ~1996!. 18S. Puri and H. L. Frisch, J. Phys. A 27, 6027 ~1994!. 19J. M. Sancho, A. H. Romero, K. Lindenberg, F. Sague ´s, R. Reigada, and A. M. Lacasta, J. Phys. Chem. 100, 19066 ~1996!. 20S. J. Langer, in Solids Far From Equilibrium, edited by C. Godreche ~Cambridge University, Cambridge, England, 1992!. 21J. D. Gunton, M. San Miguel, and P. Sahni, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz ~Academic, New York, 1983!, Vol. 8. 22J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28, 258 ~1958!;ibid. 31, 688 ~1959!. 23I. M. Lifshitz and V. V. Slyozov, J. Phys. Chem. Solids 19,35~1961!. 24A. M. Lacasta, I. M. Sokolov, J. M. Sancho, and F. Sague ´s, Phys. Rev. E 57, 6198 ~1998!. 25T. Elperin and B. Krasovitov, Trans. ASME, J. Appl. Mech. 119, 288 ~1997!. 26A. M. Lacasta, F. Sague ´s, J. M. Sancho, and I. M. Sokolov, Phys. Rev. E 59, 189 ~1999!. 6462 J. Chem. Phys., Vol. 110, No. 13, 1 April 1999 Lacasta et al. Downloaded 08 Dec 2008 to 147.83.195.22. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp