- Email: [email protected]
Nonequilibrium critical poisoning in a single-species model
Volume 127. number
3
NONEQUILIBRIUM Ronald
PHYSICS
CRITICAL POISONING
LETTERS
A
15 February
IN A SINGLE-SPECIES
1988
MODEL
DICKMAN
Department q/Physics and Astronomy, Herbert H. Lehman College, CI’NY, Bronx, N. Y. 10468. USA
and Martin
A. BURSCHKA
Btophysics Laboratory, Rockefeller University, New’ York, NY, USA Received 28 September 1987; accepted Communicated by A.R. Bishop
for publication
7 December
1987
A simple model is introduced, which exhibits a poisoning transition similar to that observed in studies of catalytic surface reactions. Steady-state behavior and critical exponents are determined via series expansion and Monte Carlo simulations. The relation between catalytic surface models and reggeon field theory is examined.
1. Introduction In studies of nonequilibrium phase transitions [ 1,2], long-range spatial correlations have commonly been ignored. However, in several models inrecently long-wavelength troduced [3-71, fluctuations are not suppressed, and the critical exponents are not correctly predicted by a simple (mean-field) rate-equation approach. An example is the catalytic surface reaction model devised by Ziff, Gulari and Barshad (ZGB) [ 51. As the relative abundance of reactants is varied, the ZGB model exhibits phase transitions between an active steady state (with ongoing reactions), and poisoned states (surface permanently covered with a single species - an absorbing state of the Markov process). One of the transitions - a nonequilibrium critical point - is of a distinctly non-mean-field character [ 61. In an effort to understand this phenomenon, we have devised a simple model which also exhibits a critical poisoning transition. While models of the kind introduced here have not been studied in the context of surface reactions, our model turns out to be closely related to the reggeon spin model, a discrete-space version of reggeon field theory (RTF) [8-l 21, also known as the “contact process” in the 132
study of interacting particle systems [ 13,141. Earlier series analysis [ 111 and Monte Carlo [ 121 studies of the reggeon spin model have examined generalized susceptibilities [ 111 and critical dynamics [ 121. Here we emphasize properties directly relevant to catalytic surface reactions: we determine the occupation fraction 2 and the cluster-size distribution via Monte Carlo simulations, and describe a perturbation method which yields precise values for X and the critical exponent p. After defining the model and discussing its relation to the ZGB and reggeon spin models, we report our simulation and theoretical results for a one-dimensional system.
2. The model Since the primary motivation for the present study is to understand the critical behavior of the ZGB model [ 51, we begin by summarizing the essential features of this model. Carbon monoxide (CO) molecules arrive at the catalytic surface at a rate Yco, and adsorb if they encounter a vacant site. Oxygen (0,) molecules arrive at rate 1 - Yco, and adsorb if they encounter a nearest-neighborpair of vacant sites. Oxygen (0) atoms and CO molecules adsorbed at 0375-9601/88/$ (North-Holland
03.50 0 Elsevier Science Publishers Physics Publishing Division)
B.V.
Volume 127, number 3
PHYSICS LETTERS A
neighboring sites react instantly to form CO*, which quits the surface. (0 atoms are free to react as individuals, the O2 molecule in effect dissociating on adsorption.) Adsorbed species remain at their site of arrival until they react. Monte Carlo simulations [ 561 yield the following phase diagram for the ZGB model in two dimensions. For Yco close to 0, the steady state is a lattice poisoned with 0 atoms, but as YcO is increased there is a continuous transition to an active steady state, followed, for still larger Yco, by a discontinuous transition to CO-poisoning. A mean-field treatment [ 151 of the ZGB model yields a qualitatively correct phase diagram, and predictions which are in excellent agreement with simulation, in the vicinity of the discontinuous transition. Our interest here is in the continuous transition, which is not accurately described by mean-field theory. Recent simulations [ 61 of the ZGB model yield an order-parameter critical exponent j3=0.6-0.7, in contrast with the mean-field prediction /I = 1. In attempting to understand the nonequilibrium critical point of the ZGB model, it is valuable to study the simplest system which displays a transition of this kind. In the present work we simplify the ZGB model in two respects. We relegate the minority species (CO molecules in the ZGB model), to the background, and consider only a single adsorbing species, which requires only a single vacant site for adsorption, Transition rules for the Markov process are as follows. The probability for a vacant site to become occupied in a short time interval dt is pdt. The probability for an occupied site x to become vacant is r dt, provided that at least one of the nearest neighbors ofx is vacant. If all of the neighboring sites are occupied, desorption is forbidden. The parameter p is analogous to 1 - Yc, in the ZGB model, and r may be integrated as the average rate of arrival of the nonadsorbing (background) species, which is able to remove adsorbed particles from the surface. Particles interior to a cluster (i.e., with all neighbors occupied), cannot desorb since they are shielded from interaction with the background species. In treating the minority species as non-adsorbing, we neglect the possibility that there molecules may also form clusters with shielded interiors. But such clusters should be of minimal importance in the ZGB model, close to the continuous transition to O-poisoning. To highlight the fact that only one species adsorbs on
I 5 February 1988
the lattice, we shall refer to the model defined above as the “A model”. Our model is closely related to the contact process [ 13,141, an interacting particle system which has been used to model the spread of an infection. In the contact process the desorption rate is proportional to the number of vacant nearest neighbors. As will be shown below, the evolution operators for the A model and the contact process differ only by a three-site interaction. We expect the latter to be an irrelevant operator, and anticipate that the A model belongs to the same universality class as the contact process, i.e., reggeon field or directed percolation [ lo]. The ZGB and A models, and the contact process, all describe microscopically irreversible evolutions: the transition probabilities do not satisfy detailed balance. The correspondence between the contact process (in d space dimensions), and directed percolation (in d + 1 dimensions), is an explicit manifestation of the timeasymmetry. A key question regarding Markov processes with absorbing states is the existence of nontrivial (“active”, in our terminology), steady states, or, equivalently, of the ergodicity of the model in the infinitesize limit [ 141. Simulations [ $61 and mean-field theory [ 151 suggest that in the ZGB model, active steady states are possible only in two or more dimensions. The contact process is known to possess (for a certain range of parameter values), nontrivial steady states for all dimensions d>, 1 [ 13,141. We expect the same to hold for the A model, thereby affording the opportunity to study critical poisoning in a simpler context.
3. Simulation results We have performed Monte Carlo simulations of the A model in one dimension, on lattices of L = 5000-40000 sites (periodic boundary conditions), in order to determine steady-state properties as a function of the parameter A=plr. At each simulation step (one unit of time), a single trial is performed, consisting in (1) choosing the process (adsorption with probability U( 1 + A), desorption with probability l/(1 +a)); (2) choosing one of the L lattice sites at random; (3) performing the process 133
Volume 127, number 3
PHYSICS LETTERS A
at the chosen site if it is permitted by the transition rules. For I< 0.55, the simulations (from an initially vacant lattice), achieve a steady state quite rapidly (within 1O3trials per site). Deriving reliable results in the critical region (Ak0.57) is more difficult, owing to (1) long relaxation times ( w 1O4trials per site); (2) large fluctuations; (3) a tendency for small systems (L,< 104) to fluctuate into the absorbing state. Thus we were led to perform long runs (up to 10’ trials per site) on large systems (L=2 X 1O4 or 4 x 104). Cheks for finite-size effects revealed that L= 5000is sufficient for simulations with A,<0.5, L=104 forJGO.55, and L=2x lo4 for 1GO.56. For A> 0.57 we used L = 4 x 104. All results reported here were obtained using periodic boundary conditions. (A series of simulations in which sites 0 and L + 1 were kept vacant yielded results consistent with the periodic boundary simulations for I < 0.57, but showed marked boundary effects in the critical region.) Our numerical estimate for the steady state (time averaged) occupation fraction X(J) is plotted in fig. 1. x is proportional to ;I for small ,I, but increases pre-
1.0
I5 February1 988
Fig. 2. Scaling analysis of the order parameter (1-Z) age cluster size S: log S (open circles) and -log( 1-X) cles) versus -log(l,-l) (1,=0.575). The straight slopes 0.28 (order parameter data), and 0.35 (cluster
and aver(filled cirlines have size data).
x I /-.
I’
5
0 Fig. 1.Steady-state occupation fraction versus d. The points represent results of Monte Carlo simulations (Lao.5 L=5000; 1~0.55, L= 104;1~0.56, L=2x 104;130.57, L=4x 104). Solid line: perturbation theory; broken line: triplet mean-field theory. The inset shows a detailed plot of the critical region.
134
Fig. 3. Cluster-size distributions: average number density of s-mers, n,, versus cluster size, s. Open circles, 1=0.5; open squares, 6~0.57; tilled circles, Lc0.575. All data is from Monte Carlo simulations with L = 4 x 104.
Volume
127, number
3
PHYSICS
cipitously as il approaches ACx 0.57 5. The inset of fig. 1 reveals significant finite-size roundoff of the critical singularity. Logarithmic plots of the order parameter, 1--x, and the average cluster size, S, versus &-I are shown in fig. 2. To obtain a straight line in the order parameter plot, one is forced to choose &=0.575; the values 0.574 or 0.576 give noticeable curvatures. If we define critical exponents through the relations ~-XCC(L-J)~ and Sa(&-n)-Y, our simulation results yield j? r 0.28 and y x 0.35. In fig. 3 we plot n, (averaged number density of s-mers) versus s, for several values of 1. The distribution is seen to follow a power law in the critical region, a manifestation of the scale-invariance of the system at the critical point. For 1~0.575 we find n,as-‘, with T= 1.72. 4. Theory We shall now summarize our preliminary theoretical results. Mean field theory (MFT) has been found to give accurate qualitative (and, in some instances, quantitative), predictions for kinetic lattice models [ 15,161. MFT may be applied at the site, pair, triplet, etc. levels; at nth level one derives a closed set of equations of motion for the densities of n-site sequences of various kinds, by factorizing higher-order densities. In the present case, MFT yields B = 1, and II,= 2, 1, and 0.8165 at the site, pair and triplet levels, respectively, where A, is defined through x(&) = 1. The triplet MFT prediction for X(1) is plotted in fig. 1, while higher-level approximations may be expected to give better estimates for X and A,, the incorrect prediction of j? = 1 appears to be an inherent feature of MFT. To develop a more systematic approach, we turn to the master equation for the probability distribution on configuration space. This distribution and its evolution are conveniently expressed in terms of an operator formalism [ 171. (Details of this formalism and of the perturbation expansion will be furnished in a forthcoming paper [ 181.) We introduce the states 1@o,i) and )@l,i) (vacant and occupied, respectively) at site i, together with creation and annihilation operators
LETTERS
15 February1988
A
A basis set for the states of the system is comprised of the vectors
I{~,) with (Ia: The state of the system at time t is given by (3) where p( {aij, t) is the probability of finding the configuration I {ai)) at time t. Physically .admissible states must obey the positivity and normalization conditions: ( {ai} ] Y ) > 0, and ( Y” I Y ) = 1, where (Y”]=&((&]+(@l,j]). Ifwe adopt a continuous-time description, and absorb the parameter r into a resealing of time, the master equation takes the form d] Y)ldt=S]
Y) ,
(4)
where S= C Sic I
1 [Ja(l-A,)A: I
+(l-A:)Ai(l-A~1Ai_lA&IAi+,)].
(5)
S preserves the conditions required for a probability
interpretation. The corresponding evolution operator for the contact process (reggeon spin model) is given (with a suitable resealing of parameters), by SC,= 1 [A(1-A,)A: +
(1-A: MA,- ,A?..I +4+ ,A:, I )I
,
which is related to the evolution for the A model through S=S,,-
1 (l-A:)A;A,_,A:_,A,+,A:+,
.
The models differ by a three-body interaction, which one expects to be an irrelevant operator. A perturbation expansion for the steady state distribution I p) may be derived by decomposing S as S= So + V, where 135
Volume 127, number
so=
ES?= I
3
PHYSICS
1 [~(l-A;).4:+(1-.4:)A,] I
(6a)
and v= T K =-
LETTERS
A
15 February
may be determined recursively, via diagrammatic expansion [ 181. Using computer enumeration of diagrams, we have derived the following expansion for the occupation fraction X=v+;~~+;~~+4$~+9.58565
f
7 (l-A:)A,A:_,A*-,A:,,A,+,
(6b)
+22.78376
Iwo,~>=(l-~)l~o.r~+~l~~,,~ Iwl.,>=I@0,,>-ldl,,>
3
(7a)
7
(7b)
and left eigenvectors ( w0.i I = < q0.i I + (6,; I 7 ~w~.il=~~~o.rl-~~-~~~~,.il
(8a> 1
CfJb)
where E U( 1 +A). Thus Sp has a nondegenerate zero eigenvalue, with eigenvector I p ) = flj I w~,~). It is convenient to work in the basis of eigenstates of Sp, and to introduce raising and lowering operators B: lWO,r)=IVUl.i)t ~,Ivo,,)=0,
Bt+ lWl.r)=“, ~,Iv,,,)=lvo,,)
9
(9)
in terms of which
fv[B,_,
-(l-v)]&?,}
X{[(l-v)B:,, +aB,+
I -u-u)lB:,,j.
-llB;+,
(10)
From the left-most factor B: it is evident that ( p” 1VI Y ) = 0 for any I Y ) . (In fact, any normalization-conserving operator has a null projection onto ( p I . ) Let [So] - ’ be the inverse of So in the nonzero subspace. Multiplying the equation S] p) =O by [SO]-‘SO=l[SO]-’ and noting that 1To ) ( 9’ 1, we obtain (l+[SO]-‘v)(P)=]!Y).
(11)
The operator [SO] - ’ V is a sixth-order polynomial in the parameter v, while the r.h.s. of eq. (11) is independent of v. Thus successive terms in the expansion
136
v6
v’ + 5 1.26263 v8 + 125.5008 1 u9
+288.96525 Sp has eigenvalues 0 and - ( 1+A), with corresponding right eigenvectors
I988
~‘~$699.61608
+ 1710.03808 + 10501.39300
v”
vi2 +4221.27861 v’4+O(v’5)
.
u13 (12)
To estimate the critical parameters we form PadC approximants to the series for (dldv) In( 1-n). Both the [ 7,6] and [6,7] PadC approximants yield &=0.5741, the former giving /3=0.276, the latter p=O.277. We adopt the value /3=0.277 ?0.002 as our series-expansion estimate, where the (rather conservative) uncertainty reflects the lack of complete regularity in the ratios of successive coefftcients in the series for I Our value for/? is in excellent agreement with earlier studies of the reggeon spin model in 1 -d, which yielded estimates: /? = 0.277 ? 0.002 (from series analysis of generalized susceptibilities [ 111, and /3 = 0.273 + 0.006 (from Monte Carlo simulations [ 121). (The calculations of refs. [ 11,121 furnish /3 indirectly, through the scaling relation p =SV.) Our present Monte Carlo estimate, B = 0.28 is also in excellent agreement with the values quoted above. To derive a prediction for the occupation fraction, we integrate the PadC approximant to (dldv) ln( 1 -X). As is evident from fig. 1, the resulting prediction is in excellent agreement with simulation over the entire range of interest. Thus the perturbation expansion approach provides an effective method for calculating non-universal as well as critical properties, and should be useful in predicting the static properties of other nonequilibrium lattice models. We have seen that, as expected, the A model belongs to the RFT/directed percolation universality class. A more interesting question is whether the critical behavior of the original ZGB model also belongs to this class. While a definitive answer is not yet apparent, we note that for RFT in 2-d, the results of ref. [ 111 yield 8=0.586+0.015, which is marginally consistent with the Monte Carlo estimatep x 0.6-0.7 for the ZGB model [ 61. We are in the process of ex-
Volume 127, number 3
PHYSICS LETTERS A
tending our series-expansion method to treat the ZGB model in two dimensions. Also under investigation are time-dependent perturbation theory and renormalization group approaches.
Acknowledgement
We thank V. Privman for a helpful discussion. The simulations reported in this work were performed on the facilities of the CUNY University Computing Center. Graph enumerations were performed using the Cornell National Supercomputing Facility.
References [ I] G. Nicolis and I. Prigogine, Self-organization in nonequilibrium systems ( Wiley-Interscience, New York, 1977). [2] H. Haken, Synergetics (Springer, Berlin, 1983).
15 February1988
[ 31 S. Katz, J.L. Lebowitz and H. Spohn, Phys. Rev. B 28 (1983) 1655; J. Stat. Phys. 34 (1984) 497. (41 J.M. Gonzalez-Miranda, P.L. Garido, J. Marro and J.L. Lebowitz, preprint. [ 51 R.M. Ziff, E. Gulari and Y. Barshad, Phys. Rev. Lett. 56 (1986) 2553. [ 61 P. Meakin and D. Scalapino, J. Chem. Phys. 87 (1987) 73 1. (71 D. Elderfield and D.D. Vvedensky, J. Phys. A 18 (1985) 2591; D. Elderfeld and M. Wilby, J. Phys. A 20 (1987) L77. [8] H.D.I. Abarbanel, J.B. Bronzan, R.L. Sugar and A.R. White, Phys.Rep. 21 (1975) 119. [9] M. Moshe, Phys. Rep. 37 (1977) 255. [ IO] J.L. Cardy and R.L. Sugar, J. Phys. A 13 (1980) L423. [ 111 R.C. Brower, M.A. Furman and M. Moshe, Phys. Lett. B 76 (1978) 213. [ 121 P. Grassberger and A. de la Torre, Ann. Phys. I22 (1979) 373. [ 131 T.E. Harris, Ann. Probab. 2 (1974) 969. [ 141 T.M. Liggett, Interacting particle systems (Springer, New York, 1985) ch. VI. [ 151 R. Dickman, Phys. Rev. A 34 (1986) 4246. [ 161 R. Dickman, Phys. Lett. A 122 (1987) 463. [ 171 M. Doi, J. Phys. A 9 (1976) 1465, 1479. [ 181 R. Dickman and M. Burschka, in preparation.
137
3
NONEQUILIBRIUM Ronald
PHYSICS
CRITICAL POISONING
LETTERS
A
15 February
IN A SINGLE-SPECIES
1988
MODEL
DICKMAN
Department q/Physics and Astronomy, Herbert H. Lehman College, CI’NY, Bronx, N. Y. 10468. USA
and Martin
A. BURSCHKA
Btophysics Laboratory, Rockefeller University, New’ York, NY, USA Received 28 September 1987; accepted Communicated by A.R. Bishop
for publication
7 December
1987
A simple model is introduced, which exhibits a poisoning transition similar to that observed in studies of catalytic surface reactions. Steady-state behavior and critical exponents are determined via series expansion and Monte Carlo simulations. The relation between catalytic surface models and reggeon field theory is examined.
1. Introduction In studies of nonequilibrium phase transitions [ 1,2], long-range spatial correlations have commonly been ignored. However, in several models inrecently long-wavelength troduced [3-71, fluctuations are not suppressed, and the critical exponents are not correctly predicted by a simple (mean-field) rate-equation approach. An example is the catalytic surface reaction model devised by Ziff, Gulari and Barshad (ZGB) [ 51. As the relative abundance of reactants is varied, the ZGB model exhibits phase transitions between an active steady state (with ongoing reactions), and poisoned states (surface permanently covered with a single species - an absorbing state of the Markov process). One of the transitions - a nonequilibrium critical point - is of a distinctly non-mean-field character [ 61. In an effort to understand this phenomenon, we have devised a simple model which also exhibits a critical poisoning transition. While models of the kind introduced here have not been studied in the context of surface reactions, our model turns out to be closely related to the reggeon spin model, a discrete-space version of reggeon field theory (RTF) [8-l 21, also known as the “contact process” in the 132
study of interacting particle systems [ 13,141. Earlier series analysis [ 111 and Monte Carlo [ 121 studies of the reggeon spin model have examined generalized susceptibilities [ 111 and critical dynamics [ 121. Here we emphasize properties directly relevant to catalytic surface reactions: we determine the occupation fraction 2 and the cluster-size distribution via Monte Carlo simulations, and describe a perturbation method which yields precise values for X and the critical exponent p. After defining the model and discussing its relation to the ZGB and reggeon spin models, we report our simulation and theoretical results for a one-dimensional system.
2. The model Since the primary motivation for the present study is to understand the critical behavior of the ZGB model [ 51, we begin by summarizing the essential features of this model. Carbon monoxide (CO) molecules arrive at the catalytic surface at a rate Yco, and adsorb if they encounter a vacant site. Oxygen (0,) molecules arrive at rate 1 - Yco, and adsorb if they encounter a nearest-neighborpair of vacant sites. Oxygen (0) atoms and CO molecules adsorbed at 0375-9601/88/$ (North-Holland
03.50 0 Elsevier Science Publishers Physics Publishing Division)
B.V.
Volume 127, number 3
PHYSICS LETTERS A
neighboring sites react instantly to form CO*, which quits the surface. (0 atoms are free to react as individuals, the O2 molecule in effect dissociating on adsorption.) Adsorbed species remain at their site of arrival until they react. Monte Carlo simulations [ 561 yield the following phase diagram for the ZGB model in two dimensions. For Yco close to 0, the steady state is a lattice poisoned with 0 atoms, but as YcO is increased there is a continuous transition to an active steady state, followed, for still larger Yco, by a discontinuous transition to CO-poisoning. A mean-field treatment [ 151 of the ZGB model yields a qualitatively correct phase diagram, and predictions which are in excellent agreement with simulation, in the vicinity of the discontinuous transition. Our interest here is in the continuous transition, which is not accurately described by mean-field theory. Recent simulations [ 61 of the ZGB model yield an order-parameter critical exponent j3=0.6-0.7, in contrast with the mean-field prediction /I = 1. In attempting to understand the nonequilibrium critical point of the ZGB model, it is valuable to study the simplest system which displays a transition of this kind. In the present work we simplify the ZGB model in two respects. We relegate the minority species (CO molecules in the ZGB model), to the background, and consider only a single adsorbing species, which requires only a single vacant site for adsorption, Transition rules for the Markov process are as follows. The probability for a vacant site to become occupied in a short time interval dt is pdt. The probability for an occupied site x to become vacant is r dt, provided that at least one of the nearest neighbors ofx is vacant. If all of the neighboring sites are occupied, desorption is forbidden. The parameter p is analogous to 1 - Yc, in the ZGB model, and r may be integrated as the average rate of arrival of the nonadsorbing (background) species, which is able to remove adsorbed particles from the surface. Particles interior to a cluster (i.e., with all neighbors occupied), cannot desorb since they are shielded from interaction with the background species. In treating the minority species as non-adsorbing, we neglect the possibility that there molecules may also form clusters with shielded interiors. But such clusters should be of minimal importance in the ZGB model, close to the continuous transition to O-poisoning. To highlight the fact that only one species adsorbs on
I 5 February 1988
the lattice, we shall refer to the model defined above as the “A model”. Our model is closely related to the contact process [ 13,141, an interacting particle system which has been used to model the spread of an infection. In the contact process the desorption rate is proportional to the number of vacant nearest neighbors. As will be shown below, the evolution operators for the A model and the contact process differ only by a three-site interaction. We expect the latter to be an irrelevant operator, and anticipate that the A model belongs to the same universality class as the contact process, i.e., reggeon field or directed percolation [ lo]. The ZGB and A models, and the contact process, all describe microscopically irreversible evolutions: the transition probabilities do not satisfy detailed balance. The correspondence between the contact process (in d space dimensions), and directed percolation (in d + 1 dimensions), is an explicit manifestation of the timeasymmetry. A key question regarding Markov processes with absorbing states is the existence of nontrivial (“active”, in our terminology), steady states, or, equivalently, of the ergodicity of the model in the infinitesize limit [ 141. Simulations [ $61 and mean-field theory [ 151 suggest that in the ZGB model, active steady states are possible only in two or more dimensions. The contact process is known to possess (for a certain range of parameter values), nontrivial steady states for all dimensions d>, 1 [ 13,141. We expect the same to hold for the A model, thereby affording the opportunity to study critical poisoning in a simpler context.
3. Simulation results We have performed Monte Carlo simulations of the A model in one dimension, on lattices of L = 5000-40000 sites (periodic boundary conditions), in order to determine steady-state properties as a function of the parameter A=plr. At each simulation step (one unit of time), a single trial is performed, consisting in (1) choosing the process (adsorption with probability U( 1 + A), desorption with probability l/(1 +a)); (2) choosing one of the L lattice sites at random; (3) performing the process 133
Volume 127, number 3
PHYSICS LETTERS A
at the chosen site if it is permitted by the transition rules. For I< 0.55, the simulations (from an initially vacant lattice), achieve a steady state quite rapidly (within 1O3trials per site). Deriving reliable results in the critical region (Ak0.57) is more difficult, owing to (1) long relaxation times ( w 1O4trials per site); (2) large fluctuations; (3) a tendency for small systems (L,< 104) to fluctuate into the absorbing state. Thus we were led to perform long runs (up to 10’ trials per site) on large systems (L=2 X 1O4 or 4 x 104). Cheks for finite-size effects revealed that L= 5000is sufficient for simulations with A,<0.5, L=104 forJGO.55, and L=2x lo4 for 1GO.56. For A> 0.57 we used L = 4 x 104. All results reported here were obtained using periodic boundary conditions. (A series of simulations in which sites 0 and L + 1 were kept vacant yielded results consistent with the periodic boundary simulations for I < 0.57, but showed marked boundary effects in the critical region.) Our numerical estimate for the steady state (time averaged) occupation fraction X(J) is plotted in fig. 1. x is proportional to ;I for small ,I, but increases pre-
1.0
I5 February1 988
Fig. 2. Scaling analysis of the order parameter (1-Z) age cluster size S: log S (open circles) and -log( 1-X) cles) versus -log(l,-l) (1,=0.575). The straight slopes 0.28 (order parameter data), and 0.35 (cluster
and aver(filled cirlines have size data).
x I /-.
I’
5
0 Fig. 1.Steady-state occupation fraction versus d. The points represent results of Monte Carlo simulations (Lao.5 L=5000; 1~0.55, L= 104;1~0.56, L=2x 104;130.57, L=4x 104). Solid line: perturbation theory; broken line: triplet mean-field theory. The inset shows a detailed plot of the critical region.
134
Fig. 3. Cluster-size distributions: average number density of s-mers, n,, versus cluster size, s. Open circles, 1=0.5; open squares, 6~0.57; tilled circles, Lc0.575. All data is from Monte Carlo simulations with L = 4 x 104.
Volume
127, number
3
PHYSICS
cipitously as il approaches ACx 0.57 5. The inset of fig. 1 reveals significant finite-size roundoff of the critical singularity. Logarithmic plots of the order parameter, 1--x, and the average cluster size, S, versus &-I are shown in fig. 2. To obtain a straight line in the order parameter plot, one is forced to choose &=0.575; the values 0.574 or 0.576 give noticeable curvatures. If we define critical exponents through the relations ~-XCC(L-J)~ and Sa(&-n)-Y, our simulation results yield j? r 0.28 and y x 0.35. In fig. 3 we plot n, (averaged number density of s-mers) versus s, for several values of 1. The distribution is seen to follow a power law in the critical region, a manifestation of the scale-invariance of the system at the critical point. For 1~0.575 we find n,as-‘, with T= 1.72. 4. Theory We shall now summarize our preliminary theoretical results. Mean field theory (MFT) has been found to give accurate qualitative (and, in some instances, quantitative), predictions for kinetic lattice models [ 15,161. MFT may be applied at the site, pair, triplet, etc. levels; at nth level one derives a closed set of equations of motion for the densities of n-site sequences of various kinds, by factorizing higher-order densities. In the present case, MFT yields B = 1, and II,= 2, 1, and 0.8165 at the site, pair and triplet levels, respectively, where A, is defined through x(&) = 1. The triplet MFT prediction for X(1) is plotted in fig. 1, while higher-level approximations may be expected to give better estimates for X and A,, the incorrect prediction of j? = 1 appears to be an inherent feature of MFT. To develop a more systematic approach, we turn to the master equation for the probability distribution on configuration space. This distribution and its evolution are conveniently expressed in terms of an operator formalism [ 171. (Details of this formalism and of the perturbation expansion will be furnished in a forthcoming paper [ 181.) We introduce the states 1@o,i) and )@l,i) (vacant and occupied, respectively) at site i, together with creation and annihilation operators
LETTERS
15 February1988
A
A basis set for the states of the system is comprised of the vectors
I{~,) with (Ia: The state of the system at time t is given by (3) where p( {aij, t) is the probability of finding the configuration I {ai)) at time t. Physically .admissible states must obey the positivity and normalization conditions: ( {ai} ] Y ) > 0, and ( Y” I Y ) = 1, where (Y”]=&((&]+(@l,j]). Ifwe adopt a continuous-time description, and absorb the parameter r into a resealing of time, the master equation takes the form d] Y)ldt=S]
Y) ,
(4)
where S= C Sic I
1 [Ja(l-A,)A: I
+(l-A:)Ai(l-A~1Ai_lA&IAi+,)].
(5)
S preserves the conditions required for a probability
interpretation. The corresponding evolution operator for the contact process (reggeon spin model) is given (with a suitable resealing of parameters), by SC,= 1 [A(1-A,)A: +
(1-A: MA,- ,A?..I +4+ ,A:, I )I
,
which is related to the evolution for the A model through S=S,,-
1 (l-A:)A;A,_,A:_,A,+,A:+,
.
The models differ by a three-body interaction, which one expects to be an irrelevant operator. A perturbation expansion for the steady state distribution I p) may be derived by decomposing S as S= So + V, where 135
Volume 127, number
so=
ES?= I
3
PHYSICS
1 [~(l-A;).4:+(1-.4:)A,] I
(6a)
and v= T K =-
LETTERS
A
15 February
may be determined recursively, via diagrammatic expansion [ 181. Using computer enumeration of diagrams, we have derived the following expansion for the occupation fraction X=v+;~~+;~~+4$~+9.58565
f
7 (l-A:)A,A:_,A*-,A:,,A,+,
(6b)
+22.78376
Iwo,~>=(l-~)l~o.r~+~l~~,,~ Iwl.,>=I@0,,>-ldl,,>
3
(7a)
7
(7b)
and left eigenvectors ( w0.i I = < q0.i I + (6,; I 7 ~w~.il=~~~o.rl-~~-~~~~,.il
(8a> 1
CfJb)
where E U( 1 +A). Thus Sp has a nondegenerate zero eigenvalue, with eigenvector I p ) = flj I w~,~). It is convenient to work in the basis of eigenstates of Sp, and to introduce raising and lowering operators B: lWO,r)=IVUl.i)t ~,Ivo,,)=0,
Bt+ lWl.r)=“, ~,Iv,,,)=lvo,,)
9
(9)
in terms of which
fv[B,_,
-(l-v)]&?,}
X{[(l-v)B:,, +aB,+
I -u-u)lB:,,j.
-llB;+,
(10)
From the left-most factor B: it is evident that ( p” 1VI Y ) = 0 for any I Y ) . (In fact, any normalization-conserving operator has a null projection onto ( p I . ) Let [So] - ’ be the inverse of So in the nonzero subspace. Multiplying the equation S] p) =O by [SO]-‘SO=l[SO]-’ and noting that 1To ) ( 9’ 1, we obtain (l+[SO]-‘v)(P)=]!Y).
(11)
The operator [SO] - ’ V is a sixth-order polynomial in the parameter v, while the r.h.s. of eq. (11) is independent of v. Thus successive terms in the expansion
136
v6
v’ + 5 1.26263 v8 + 125.5008 1 u9
+288.96525 Sp has eigenvalues 0 and - ( 1+A), with corresponding right eigenvectors
I988
~‘~$699.61608
+ 1710.03808 + 10501.39300
v”
vi2 +4221.27861 v’4+O(v’5)
.
u13 (12)
To estimate the critical parameters we form PadC approximants to the series for (dldv) In( 1-n). Both the [ 7,6] and [6,7] PadC approximants yield &=0.5741, the former giving /3=0.276, the latter p=O.277. We adopt the value /3=0.277 ?0.002 as our series-expansion estimate, where the (rather conservative) uncertainty reflects the lack of complete regularity in the ratios of successive coefftcients in the series for I Our value for/? is in excellent agreement with earlier studies of the reggeon spin model in 1 -d, which yielded estimates: /? = 0.277 ? 0.002 (from series analysis of generalized susceptibilities [ 111, and /3 = 0.273 + 0.006 (from Monte Carlo simulations [ 121). (The calculations of refs. [ 11,121 furnish /3 indirectly, through the scaling relation p =SV.) Our present Monte Carlo estimate, B = 0.28 is also in excellent agreement with the values quoted above. To derive a prediction for the occupation fraction, we integrate the PadC approximant to (dldv) ln( 1 -X). As is evident from fig. 1, the resulting prediction is in excellent agreement with simulation over the entire range of interest. Thus the perturbation expansion approach provides an effective method for calculating non-universal as well as critical properties, and should be useful in predicting the static properties of other nonequilibrium lattice models. We have seen that, as expected, the A model belongs to the RFT/directed percolation universality class. A more interesting question is whether the critical behavior of the original ZGB model also belongs to this class. While a definitive answer is not yet apparent, we note that for RFT in 2-d, the results of ref. [ 111 yield 8=0.586+0.015, which is marginally consistent with the Monte Carlo estimatep x 0.6-0.7 for the ZGB model [ 61. We are in the process of ex-
Volume 127, number 3
PHYSICS LETTERS A
tending our series-expansion method to treat the ZGB model in two dimensions. Also under investigation are time-dependent perturbation theory and renormalization group approaches.
Acknowledgement
We thank V. Privman for a helpful discussion. The simulations reported in this work were performed on the facilities of the CUNY University Computing Center. Graph enumerations were performed using the Cornell National Supercomputing Facility.
References [ I] G. Nicolis and I. Prigogine, Self-organization in nonequilibrium systems ( Wiley-Interscience, New York, 1977). [2] H. Haken, Synergetics (Springer, Berlin, 1983).
15 February1988
[ 31 S. Katz, J.L. Lebowitz and H. Spohn, Phys. Rev. B 28 (1983) 1655; J. Stat. Phys. 34 (1984) 497. (41 J.M. Gonzalez-Miranda, P.L. Garido, J. Marro and J.L. Lebowitz, preprint. [ 51 R.M. Ziff, E. Gulari and Y. Barshad, Phys. Rev. Lett. 56 (1986) 2553. [ 61 P. Meakin and D. Scalapino, J. Chem. Phys. 87 (1987) 73 1. (71 D. Elderfield and D.D. Vvedensky, J. Phys. A 18 (1985) 2591; D. Elderfeld and M. Wilby, J. Phys. A 20 (1987) L77. [8] H.D.I. Abarbanel, J.B. Bronzan, R.L. Sugar and A.R. White, Phys.Rep. 21 (1975) 119. [9] M. Moshe, Phys. Rep. 37 (1977) 255. [ IO] J.L. Cardy and R.L. Sugar, J. Phys. A 13 (1980) L423. [ 111 R.C. Brower, M.A. Furman and M. Moshe, Phys. Lett. B 76 (1978) 213. [ 121 P. Grassberger and A. de la Torre, Ann. Phys. I22 (1979) 373. [ 131 T.E. Harris, Ann. Probab. 2 (1974) 969. [ 141 T.M. Liggett, Interacting particle systems (Springer, New York, 1985) ch. VI. [ 151 R. Dickman, Phys. Rev. A 34 (1986) 4246. [ 161 R. Dickman, Phys. Lett. A 122 (1987) 463. [ 171 M. Doi, J. Phys. A 9 (1976) 1465, 1479. [ 181 R. Dickman and M. Burschka, in preparation.
137