2770
J. Phys. Chem. 1995, 99, 2770-2776
Nonlinear Chemical Dynamics in Low-Dimensional Lattices and Fractal Sets A. Tretyakov,? A. Provata, and G. Nicolis* Center for Nonlinear Phenomena and Complex Systems, UnivertitC Libre de Bruxelles, Campus Plaine, CP 231, 1050 Brussels, Belgium Received: August 1, 1994; In Final Form: November 22, 1994@
+
We report an extensive analytical and numerical study of the trimolecular svstem A 2X t 3X on lattices with r&ricted geometry, including fractal sets. It is shown that the asymptotic value of (A)/(@ attained starting from a lattice completely covered by X is largely determined by the mean coordination number, while fractal dimension seems to play little role. If a small external flux of A is allowed, the behavior is no longer dominated by coordination number, but by the existence of stable configurations of A-sites.
1. Introduction One of the main preoccupations of Stuart Rice has been to elucidate the connection between microscopic dynamics at the molecular level and collective behavior at the macroscopic level. As a tribute to his brilliance and to the impact of his work on so many generations of physical chemists over more than 30 years, we present, in this essay, some new results on the transition from microscopic to macroscopic behavior in connection with the limits of validity of mean-field description in nonlinear chemical dynamics. The traditional phenomenological study of chemical reactions is based on the rate equations describing the evolution of the concentrations {Xi(r,t)} of the various species in space and in time. 1-3
+
axi(r,t)/at= Div2Xi(r,t) V,({X,>>
(1)
where vi stands for the overall effect of the chemical reactions and Di is the diffusion coefficient of the species i . To be more specific, consider the chemical reaction model:
A
+ 2X
3X
(2)
According to the above scheme the rate equations for the species A and X are
+
aA(r,t)/at= DAV2A(r,t)- kA(r,t)X2(r,t) kX3(r,t) (3)
where DA and DX are the diffusion coefficients for the species A and X, respectively, and k is the kinetic constant (hereafter we shall assume that the forward and backward kinetic constants are equal). By solving eqs 3 and 4, we can get the spatiotemporal behavior of the two chemical species. In the above two equations space does not seem to play any important role except in the diffusion term. Consequently, species with locally high concentrations will tend to diffuse toward regions of low concentration. Recently, with the expansion of the field of fractals there has been considerable interest in phenomena that take place in fractal sets or in sets of low dimen~ionality.~,~ From studies of random walks and diffusion it is now known that low dimenPermanent address: Graduate School of Information Science, Tohoku University, Aramaki, Aoba-ku, Sendai 980, Japan. Abstract published in Advance ACSAbstracts, February 1, 1995. @
0022-365419512099-2770$09.00/0
sionality can have nontrivial effects. In particular, in a mixture of chemically reacting particles the space-induced limitations of diffusion might not allow a representative particle to travel far enough to interact with all the other particles in the system, as is tacitly assumed in eq 1. In previous work by J. W. Turner and two of the present authors the role of spatial dimensionality in the behavior of chemically reacting systems has been analyzed for two representative nonlinear reaction models:6 the trimolecular model, eq 2, and the bimolecular model A X 2X. We briefly summarize the main conclusions: (i) In a one-dimensional regular lattice with nearest neighbor interactions, starting from a configuration containing only X particles, the trimolecular reaction system stabilizes in a nonequilibrium, locally frozen asymptotic state in which the ratio r of A to X particles is a constant number, r = 0.38, quite different from the mean-field ratio IMF = 1. (ii) The same trimolecular system reaches the mean-field limit, as far as equilibrium properties are concerned, in a twodimensional square lattice with nearest neighbor interactions. (iii) The bimolecular reaction system A X 2X agrees, as far as equilibrium properties are concerned, with mean-field predictions in both one and two space dimensions. In the present paper the earlier work on the status of mean field in nonlinear chemical dynamics is extended to the case where the reaction is occumng on a variety of sets with restricted geometry, including fractal sets spanning a wide range of (noninteger) dimensionalities up to D = 2. Furthermore, other types of regular two-dimensional lattices than the square lattice are considered. Our principal goal is to identify some key features behind the failure of mean-field behavior and, in particular, to assess the role of the mean coordination number. In section 2 a series of numerical simulations of the equilibrium properties of the trimolecular system on a wide variety of lattices, including some typical fractal sets, are reported. This study brings out the determining influence of the mean coordinationnumber as well as a surprising correlation with simple mean-field analysis. A more detailed analytical study, devoted to this correlation is carried out in section 3. In section 4 the behavior of systems showing similar behavior under the conditions of sections 2 and 3 is studied in the presence of an external flux. It is shown that under this constraint of opening to the external environment the behaviors can become quite different; in other words, the mean coordination number is no longer the only parameter governing the validity of the mean-field behavior. Furthermore, it is shown that in the presence of a first-order structural phase transition,
+
+
0 1995 American Chemical Society
Low-Dimensional Lattices and Fractal Sets
J. Phys. Chem., Vol. 99, No. 9, 1995 2771
TABLE 1: Asymptotic (A)/(X)Ratios for the Trimolecular System in a Variety of Regular, Fractal, and Random Lattices coordination fractal geometry number dimension (A)/@) square
hexagonal modified square percolation clusters 0.4 0.5 0.592 75 0.7 0.75 0.8 0.9 Sierpinski carpets 110011101 110011011 111010011 101110111 110101111 111001111 111010111 random Sierpinski carpets 1 2 3 4 4 (4 x 4)
Sierpinski gasket 1 - dline comb strip of width 2 Cantor sets 111111000 110111011 110110011 random Cantor sets 1 2 3 4
4 3 3 1.6 2 2.4 2.8 3 3.2 3.6 1.6 2.2 2.2 2.4 2.4 2.6 2.6 3.3 2.7 2.1 1.6 2.7 4 2 2 3 1.7 1.33 1.2 1.7
1.5
exponential distribution of intervals power law distribution of intervals every site connected to every other 3 4 5
1.2 0.9 1.o 1.25 1.5 1.o 1.25 1.5 2 3 4
2 2 2 1.90
1.63 1.63 1.63 1.77 1.77 1.77 1.77 1.89 1.77 1.63 1.46 1.79 1.58 1 0.82 0.89 0.89 0.95 0.89 0.82 0.73
1 0.92 0.90 0.26 0.38 0.52 0.70 0.80 0.87 0.97 0.27 0.47 0.45
0.58 0.57 0.65 0.64 0.93 0.69 0.44 0.27 0.68 1 0.38 0.33 0.81 0.27 0.16 0.07 0.29 0.21 0.14
0.08 0.21 0.26 0.30 0.28
0.30 0.32 0.33 0.57 0.73
controlled by the concentrations of reactants, the open system can exhibit regular time oscillations similar to those reported in experimentson heterogeneous catalysis. The main conclusion and suggestions for further study are compiled in section 5 .
2. Long-Time Behavior of the Trimolecular Model: Monte Carlo Simulations In an attempt to capture the dependence of the asymptotic ratio r E (A)/(% on the structure of the underlying lattice, Monte Carlo simulations have been performed on a variety of lattices, both random and deterministic. The simulation algorithm was as follows. A site is chosen at random. If it has two neighboring X-sites, the state of the site is changed from A to X or from X to A depending on the current state of the site. At the beginning the lattice is filled completely by X particles. We use lattices containing more than 1000 sites, and usually it was sufficient to take 100 steps per site to reach a state which could be considered as asymptotic, although in some cases longer runs had to be performed. The results of the simulations are summarized in Table 1 and Figure 1, together with some results that could be obtained analytically, and will be treated in detail in the next section. In the rest of this section we describe the various lattices used in the Monte Carlo simulations and comment on the r values obtained in each case.
1 ,
0.6
I
I
I
1
1
1
0.5 0.4
0.3 0.2 0.1
0.5
1
I
I
1
I
I
I
1
1.5
2
2.5
3
3.5
4
Mean coordination number
Figure 1. Asymptotic (A)/(X)ratios for the trimolecular system in a variety of regular, fractal, and random lattices versus the mean coordination number.
Figure 2. Modified square lattice with coordination number 3. The regular periodic lattices considered include square lattice, hexagonal lattice, and a modification of the square lattice obtained by adding some extra sites and having coordination number 3, the same as the hexagonal lattice (see Figure 2). For the hexagonal lattice it turned out that r = 0.92, definitely different from the r = 1 of the square lattice. On the other hand, the modified square lattice gave the value r = 0.90, quite close to that of the hexagonal lattice. These results demonstrate that for the trimolecular model non-mean-field behavior is possible on regular two-dimensional lattices, and behavior on a square lattice cannot always be equated with behavior in two dimensions, as was tacitly assumed in ref 6. We next discuss percolation clusters constructed by removing sites from a square lattice. Since the fraction of open sites is connected to the coordination number as p = coord.number/4 the coordination number can be changed continuously from 0 to 4. It should be noted that the percolation system at the percolation threshold (p = 0.592 75) has fractal structure with fractal dimension D = 1.90.7 On the other hand, as seen from Figure 1, the (A)/(% ratio for percolation clusters changes continuously with coordination number. No particular behavior can be detected at the percolation threshold, which indicates that long-range interactions are not important in determining the r value.
To span the range of fractal dimensions between one and two, we used a family of Sierpinski carpets prepared by cutting a square into nine equal-size smaller squares, removing some of them, and repeating this procedure on each of the smaller squares. In the context of a lattice, this process was conducted as follows: take a 3 x 3 matrix filled by ones and zeros, like 1 1 1 1 0 1 1 1 0 each 1 corresponding to a site that can be occupied by A or X.
2772 J. Phys. Chem., Vol. 99, No. 9, I995
Tretyakov et al.
t-4
t-4 Figure 3. Sierpinski carpet, iteration 2.
Figure 4. Sierpinski gasket, iteration 3.
Replace each 1 by the whole matrix
number, one might expect for the Sierpinski gasket (D= 1.58) the same (A)/(X)ratio as for the square lattice. This is confirmed by our simulations. As an example of fractal sets with dimension less than 1, we used a family of Cantor-like sets. Deterministic Cantor-like sets could be treated exactly and will be discussed in the next section; random ones were constructed similarly to random Sierpinski carpets and are summarized in Table 1. One should note that, when plotted against the coordination number (Figure l), the data for random Cantor sets, random Sierpinski carpets, and percolation clusters seem to form a smooth continuous line in the range of mean coordination numbers spanning the interval from 1 to 3. In summary one can conclude that the order of magnitude of the asymptotic (A)/(* ratio is largely determined by the mean coordination number, while fractal dimension seems to bear little relevance to the problem. Results for various structures seem to group around the results obtained for the percolation clusters, which are special in the sense that the mean coordination number can be controlled. To gain further insight on this result, it is instructive to explore the following simplistic approximation for the percolation system. Let us note that the concentrations of sites with coordination numbers 0 and 1 that cannot change state according to the reaction rules and, consequently, are covered by X, are p(1 - p)4 and 4p2(1 - P ) ~ respectively. , Assuming that the rest of the sites can take A or X states with equal probability, we have
1 1 1 1 0 1 1 1 0 and replace each 0 by a matrix of zeros of the same size
0 0 0 0 0 0 0 0 0 This process can be repeated on the resulting 9 x 9 matrix and so forth. After the required number of iterations has been taken, connectivity is determined on the basis of the nearest neighborhood. Henceforth for simplicity we denote generating matrices as binary numbers obtained by putting the raws of the matrix next to each other; for example, 1 1 1 1 0 1 1 1 0 will be written as 111101110. The classical Sierpinski carpet (Figure 3), obtained by repeatedly removing the middle square, is thus encoded as 111101111. The results, summarized in Table 1, show that carpets with the same fractal dimension may lead to different r values if the mean coordination number is different, which is hardly surprising in view of the already mentioned results for hexagonal and square lattices. Random Sierpinski carpets were constructed in a similar way, except that each time the generating matrix is applied, its elements are mixed in a random manner, while the number of zeros is preserved. The random fractals constructed in this way are obviously characterized by the dimensions of the generating matrix and the number of zeros in it. Besides random carpets generated by 3 x 3 matrices, we considered one case of a random fractal obtained using a 4 x 4 generating matrix. We notice that as the fractal dimension of random fractals approaches 1, the asymptotic r value does not converge to the one for a one-dimensional regular lattice ( r = (A)/(X)= [5 (5)'12]/[5 -t (5)'12] = 0.38, see ref 6), but definitely falls below it (Table 1). On the other hand, percolation clusters with mean coordination number 2 lead to a (A)/(X)ratio similar to the one for the one-dimensional regular lattice, although the details of the structure are quite different. As one more example, the random carpet with coordination number 1.6 gives an r value similar to the deterministic carpet 110011101 having a similar coordination number, although fractal dimensions are different. We conclude that coordination number is a more appropriate parameter than fractal dimension. One further illustration of that point are the results of our simulations on the Sierpinski gasket (Figure 4):anticipating that the asymptotic ratio is largely determined by coordination
r=
1 - q4 - 4pq3 1 q4 4pq3
+ +
(5)
where q e 1 - p . As one can see from Figure 1, the shape of the curve obtained by simulation is reproduced quite well, although eq 5 gives a somewhat higher value of r. 3. Long-Time Behavior: Analytic Results
In the previous section we saw that the (A)/(X)ratio achieved from the all-X-covered state for a vast variety of systems is largely determined by the mean coordination number. Here we present some exactly solvable cases in an attempt to see whether the specificity of structure may lead to deviations from this overall behavior. We start with the particularly simple case of a lattice where every site is connected to every other site. In the asymptotic state all allowed configurations (in which every A-site has at least two X neighbors) have the same statistical weight, the only nonallowed configurationsbeing one configuration with all sites covered by A and n configurations containing only one X-site. Counting all conceivable configurations, one can write
Low-Dimensional Lattices and Fractal Sets
-n _- ( X ) ( 2 " - n - l ) + l 2
xn+Ox 1
2"
J. Phys. Chem., Vol. 99, No. 9, I995 2773
(6)
where n is the size of the system and n/2 gives the mean number of X-sites in a system filled completely at random. Equation 6 together with the obvious relation ( X ) (A) = n gives for the (A)/(@ ratio
+
(7) As expected, in the large n limit r attains a mean-field value of r = 1 (in fact, solutions for this kind of an ideally connected lattice are often equated with mean-field solutions), while for n = 3 (coordination number 2) and n = 4 (coordination number 3), respectively, we have r = 0.33 and r = 0.57. The value of r for n = 4 noticeably falls out of the range of values obtained for other systems with coordination number 3 (see Table 1 and Figure 1). System with n = 5 (coordination number 4) gives r = 0.73, far from the mean-field value obtained for the square lattice and for the Sierpinski gasket. This example shows that finite systems and, possibly, infinite systems formed as collections of disconnected finite elements with narrow distribution of sizes are likely to demonstrate a pathological behavior. Now let us consider one more case of a system with coordination number 2, this time infinite. Consider the structure given in Figure 5, to which we refer henceforth as a comb. The state of sites on the dead ends cannot be changed; starting from an X-covered state, they remain X, while for sites on the line all states are allowed with the exception of the statistically insignificant all-A-covered state. The (A)/(X)ratio, obtained as (1/2)/(1/2 1) = 1/3 = 0.33, is the same as in the case of a system of three sites with every site connected to every other and somewhat different from the case of a regular onedimensional lattice, for which r = 0.38. Thus, using data obtained by exact analytical methods, we demonstrate that two infinite systems may have different r values, even if the coordination number is the same.
+
We now turn to the system depicted in Figure 6, an infinite strip of width 2. This system has coordination number 3, the same as in the hexagonal lattice. Let us note that there are four ways a configuration, connected ergodically to the all-X-filled configuration, can be realized on a strip of length i. Either a strip made up by the first i - 1 layers taken on its own is ergodically connected to a strip of length i - 1 completely covered by X. Or the same is true for a strip composed by the first i - 2 layers while the (i - 1)th layer is covered by A. Or the first i - 1 layers are covered by intervals, ergodically connected to all-X-covered intervals, alterating from one side of the strip to another, with the other side of the stip covered by A, and leaving exactly one A-covered layer between them. Or the same is true for the first i - 2 layers with the (i - 1)th layer covered by A. Notice that to produce a configuration ergodically connected to the all-X configuration on the whole i-layer strip, the ith layer should be occupied by at least one X. Counting all possibilities leads to the following system of equations, in which the mean number of X-covered sites in a strip of length i is written as X::
X: = [(3/2
+ X;-,)4n;-, + (4/3 + X!-,)6n;-, + (4/3 + X:-,)3n:-, + (4/3 + X;-,)3n;-,]/n; n: = 4n:-, + 6n:-, + 3n:-, + 3n:-,
(8) (9)
Figure 5. Comblike structure with mean coordination number 2.
B
Figure 6. Strip of width 2. All sites have coordination number 3.
Here Xi' is the mean number of X-sites for configurations with all-X ergodically connected intervals alterating from one side of a strip to another, and n: are numbers of configurations. Solving this system gives the ratio (A)/(X)= (2i - X:)/Xf = 0.81 in the limit of large i, which is somewhat lower than the simulation results for the hexagonal lattice and close to the percolation clusters with p = 0.75 and, correspondingly, mean coordination number 3 (see Table 1). The examples of the hexagonal lattice, the modified square lattice, percolation clusters with p = 0.75, and the infinite strip of width 2, all having the same coordination number, show that while the mean coordination number determines the order of magnitude of r, the actual values of this ratio are distributed over some interval according to the peculiarities of individual lattices. To clarify how the (A)/(X)ratio is formed on fractals, we solve some cases of fractal sets. We start with Cantor set-like fractals, obtained by cutting an interval into nine parts, removing some of them, and repeating this process with each of the obtained intervals. We denote this process using a binary number, like 110110011, with zeros corresponding to removed intervals. In the context of a lattice each 1 corresponds to a site that can be occupied by A or X, and in the iteration process each 1 is replaced by the whole 110110011, while 0 is replaced by a combination of nine zeros. From the point of view of the (A)/(X) ratio, since solution for a single interval is known,6 the problem is solved once the distribution of intervals is known. It should be noted that, in general, systems with fractal dimension less than 1, as our Cantor-like sets, can always be reduced to a collection of finite sets. In the case of the rule given by 110110011, as one can easily see, the system is formed by intervals of sizes 2 and 4. The numbers of intervals at iteration step i are given by
2774 J. Phys. Chem., Vol. 99, No. 9, 1995
T -b Figure 7. Sierpinski carpet generated by matrix 110011101, iteration 3.
with ni = 3 and n; = 0. Solving these equations for large i and averaging the mean numbers of A and X over the resulting distribution of intervals using the relation (L- 1)/2
k(kPk-‘) !FO
(4= (L-1)12
c (:-“-‘I
,
O=L-
(4
(16)
k=O
for an interval of length L,6 a formula obtained by counting all configurations ergodically compatible to the configuration in which all sites of the interval are covered by X, one obtains r = 0.071. In the case of a Cantor-like set generated by 110111011 one has, on the other hand,
and the ratio r = 0.156. Any deterministic Cantor-like set can be solved using a similar approach. As can be seen in Table 1, the resulting r values do not correlate with dimension. That is not very surprising since the distribution of empty spaces, crucial for determining the fractal dimension, does not play much of a role in determining the (A)/(@ ratio. One might also notice from Table 1 and Figure 1 that deterministicCantor sets, unlike ranodm Cantor sets, don’t fall on the overall r-coordination number curve very well, because the distribution of the finite sets comprising them is rather narrow. We finish our discussion of deterministic fractals with the Sierpinski carpet corresponding to the matrix 110011101 (see the previous section for an explanation of how the carpet is formed). This case is relatively easy to solve as it becomes apparent from Figure 7, depicting the second iteration of the set. We see that the set is formed by segments of 1 - d lines. The numbers of lines of length k, nf, can be obtained by solving the following system:
where di is the size of the longest segment at iteration i. Although this case looks quite peculiar, the value of r = 0.266 obtained fits well with the values for systems with a similar coordination number. To complete our survey of results that can be obtained using analytical methods, we should mention that it is possible to construct systems with different coordination numbers and (A)/ ( X ) ratios by collecting intervals according to an appropriate distribution. The (A)/(X)ratio for any interval distribution can be obtained by averaging eq 16. Some results for exponential and power law distributions are given in Table 1 and Figure 1, showing once again that although the mean coordination number determines the order of magnitude of the (A)/(X)ratio surprisingly well, generally speaking, systems with different geometry will have different ratios even if the coordination number is the same.
4. Structural Stability with Respect to Small Perturbation of the Reaction Scheme In this section we elaborate further on the deviation from mean-field behavior in systems with the same coordination number by analyzing their relative responses when the reaction scheme is slightly modified to account for inflow of chemical from the outside world. Most of the systems both in nature and in chemical technology are open systems exchanging material with the environment. It is by now well established that when operating under nonequilibrium conditions, open systems can show a variety of behaviors, including multistability, oscillations, and c h a ~ s . ~ Our . ~ J first ~ objective here will be to see the effect of this mode of operation on mean-field behavior. For this purpose the following modified trimolecular scheme will be introduced:
k2
x-s*
+x,
(29)
Here S* denotes an empty but active adsorption site on the lattice; Ag and X, denote molecules of A and X species in a fluid (e.g. gaseous) phase in contact with the lattice; A and X denote the same species in the adsorbed phase. It is understood that the adsorption of A is fast and that the desorption of X-the “final product” of the reaction-is slow. Under these conditions the rate of mass exchange with the outside world, occurring in steps 27 and 29, is determined by the rate of step 29: once an X molecule leaves the lattice in step 29, the resulting empty adsorption site S* is immediately occupied by A via step 27. We simulate this process as follows: reactions A -I- 2X 3X and 3X A 2X are carried out, as previously, at a rate k = 1. Additionally, at a rate k2 we take a site at random, and if it happens to be X, we change it to A. For a system with every site connected to every other site the transition matrix with flux of A taken into account can be constructed quite straightforwardly. Distinguishing only con-
-
+
-
J. Phys. Chem., Vol. 99,No. 9,1995 2775
Low-Dimensional Lattices and Fractal Sets
0 O. 4e 55 t
0.251
I I
Time ”.-
0
