Dynamic response of an irreversible reaction system to a pulsed perturbation

. The dynamic response of the ZGB surface reaction lattice gas model, for the catalyzed reaction A + (1 = 2)B 2 → AB, is studied by means of Monte-Carlo simulations in the neighborhood of its second order irreversible phase transition (IPT). It is found that shortly after driving a stationary con(cid:12)guration into the absorbing state, the relaxation of the system can be well described by a stretched exponential behavior. The dependence of the relaxation characteristic time and the induced changes on the coverage of the reactants, on both, the intensity and the period of the pulsed perturbation, are systematically investigated. The obtained insights can straightforwardly be extended to a wide variety of irreversible systems exhibiting second order IPT’s, such as directed percolation, forest (cid:12)re models, the contact process, branching annihilating walkers, catalyzed reactions, etc.


Introduction
The dynamic response to an external perturbation of model systems close to reversible phase transitions, has recently been studied extensively [1][2][3][4][5][6][7].The archetype model for such studies is the well-known Ising magnet.The dynamical response of Ising systems to an oscillatory magnetic field has led to the discovery of interesting phenomena such as dynamic hysteresis [1] and a fluctuation-induced symmetry breaking transition [6].Also the anomalous behavior of the pulsed susceptibility has been reported.Very recently, the dynamic response of an Ising system to a pulsed magnetic field has been studied by means of Monte-Carlo simulations and solving numerically the mean-field equation of motion [7].Also, the nonexponential relaxation of the fully frustrated Ising model has been reported very recently [5].Within this context, one should also mention interesting experimental observations of the dynamic scaling behavior of the magnetic hysteresis in thin ferromagnetic Fe/Au(001) films [8].
So far, the study of the dynamic response of an intrinsically irreversible system close to irreversible phase transitions (IPT's) has not been addressed yet.IPT's between a reactive (or active) regime and an inactive (or absorbing) state have been reported in a wide variety of model systems, e.g.directed percolation [9], the contact process [10], branching annihilating walkers [11], forest fire models [12], models of the dynamic evolution of living individuals [13] a e-mail: ealbano@isis.unlp.edu.ar and more significantly in irreversible reaction systems [14] (for reviews see e.g.[15,16]).
One of the better known irreversible reaction system is the dimer-monomer ZGB model, as earlier proposed by Ziff, Gulari and Barshad [14].The ZGB model deals with a catalyzed lattice-gas monomer-dimer reaction of the type A + (1/2)B 2 → AB.Reactant's A and B 2 are adsorbed on the surface of the catalyst; i.e. according to the Langmuir-Hinshelwood mechanism; with probabilities Y A and Y B , respectively.Taking Y A + Y B = 1, the ZGB model has a single parameter; namely Y A ; which plays the role of an external pressure due to the reactants in the gas phase in contact with the surface.For Y A → 1 (Y A → 0) the surface of the catalyst becomes inactive due to complete saturation with A-(B-)species, respectively.However, between these two inactive states there is a reactive regime as it is shown in Figure 1.More precisely, for Y A ≤ Y A1 ∼ = 0.3906 the final state of the system is a catalyst's surface fully covered (or poisoned) with B-species.Similarly, for Y A ≥ Y A2 ∼ = 0.5256 the surface becomes poisoned by A-species.IPT's between the reactive regime and the B-(A-)poisoned states are of second (first) order, respectively, as it can be observed in Figure 1.For more details on the ZGB model see e.g.[14,15].Within this context, the aim of this work is to investigate the dynamic response of the ZGB model, close to the second order IPT, to a pulsed perturbation of the external pressure.The understanding of this physical situation may straightforwardly be generalized to other systems exhibiting second order IPT's such as forest fire models, the contact process, branching annihilating walkers, and many others irreversible reaction systems [15].Also, the irreversible second order transitions exhibited by these systems have extensively been studied using a variety of theoretical approaches [9,[17][18][19][20].The main finding is that such kind of IPT belongs to the universality class of directed percolation [17] and therefore it is very well understood from the theoretical point of view [9].Furthermore, the approach proposed in this work is the natural previous step for further investigations of the dynamic response of irreversible systems to more complex perturbations.It should be noticed that this field remains unexplored, so we have restricted ourselves to present and discuss numerical results which are aimed to gain the first insight on the field and to stimulate further theoretical work.Also, this kind of study has motivations in the field of applied catalysis since, in some cases, oscillatory variations of the external pressure may enhance the rate of production of catalyzed reactions [21,22].

Description of the ZGB model and the simulation technique
In order to simulate the ZGB model we have used a square lattice of side L = 512 with periodic boundary conditions to represent the catalytic surface.The Monte-Carlo algorithm for the simulation is as follows: (i) A or B 2 molecules are selected randomly with relative probabilities Y A and Y B = 1−Y A , respectively.These probabilities are the relative impingement rates of both species, which are proportional to their partial pressures.Due to the normalization, Y A + Y B = 1, we used Y A as a parameter.If the selected species is A, one surface site is selected at random, and if that site is vacant, A is adsorbed on it.Otherwise, if that site is occupied, the trial ends and a new molecule is selected.If the selected species is B 2 , a pair of nearest neighbor sites is selected at random and the dimer becomes adsorbed only if they are both vacant.(ii) After each adsorption event, the nearest neighbors sites of the added molecule are examined in order to account for the reaction.If more than one [B(a),A(a)] pair is identified, a single one is selected at random and removed from the surface.
The Monte-Carlo time step (MCTs) (t) involves L 2 trials, so each site of the lattice is visited once, on average, during each time unit.Simulations are started with empty lattices taking a fixed working value (Y Aw ) for the parameter of the model.Y Aw is always selected within the reactive regime of the system.The first 10 4 time steps are disregarded in order to allow the establishment of a stationary state and then the perturbation is switched on, that is rectangular pulse of intensity ∆p is applied during a period of time τ p .So, during the operation of the perturbation the actual value of the tunable parameter is 3906 the system is driven to the B-poisoned state (see Fig. 1) and therefore one may expect the surface to become saturated by B-species if the period of the pulse is long enough.
Just after the application of the pulse and subsequently, several quantities of interest; e.g. the coverages with A-and B-species given by θ A and θ B , and the rate of AB production R AB , respectively; are recorded as a function of time.All these quantities are averaged over 10 3 different realizations.pressure is Y Aw = 0.4600 and the parameters characterizing the pulse are τ p = 20 MCTs and ∆p = 0.10.So, in this example the system is shortly driven into the B-poisoned state with Y * A = 0.3600.Consequently the coverage with B-species increases in contrast to the behavior of both the coverage with A-species and the rate of reaction.After suppression of the pulse the system quickly relaxes towards the stationary state at Y Aw .Unless otherwise stated, we will focus our attention to the relaxation properties of the B-species because, on the one hand, all measured quantities relax showing the same behavior and, on the other hand, θ B exhibits the best signal/noise ratio.

Results and discussion
Figures 3a and 3b show two typical relaxation curves of θ B .The first one (Fig. 3a), is taken for Y Aw = 0.4600, i.e. well inside the reactive regime, while Y * A = 0.3600 is rather close to the critical edge given by Y A1 ∼ = 0.3906.In the second example (Fig. 3b, we have taken Y Aw = 0.4175 close to the critical edge so the system has been driven deeply into the poisoned state with Y * A = 0.3175.It is found that the relaxation of the coverage can be very well fitted by a stretched exponential curve of the form where θ Bmax is the maximum value of the coverage with B-species which is just achieved when the peak is suppressed, θ ∞ is the coverage after proper relaxation; i.The relaxation time in the case of Figure 3b is longer that for the example of Figure 3a.This is due to the fact that in the former the system has been driven deeper into the poisoned state.However, this preliminary result has to be taken with caution because in the limit Y * A → 0 a different behavior is observed, as it will be discussed below. It should be noted that the conditions under which stretched exponential behavior occurs in physical systems are controversial and of current interest.In many cases it is very difficult to distinguish an stretched exponential from simple exponential or even power law behavior (for an excellent discussion on this issue see e.g.[23]).
In the examples of Figure 3 the stretching exponents b are close to unity, causing further difficulties in the identification of the actual behavior.So, in order to support the validity of our claim, let us define the function X(t) = [θ B (t)−θ ∞ ]/θ Bmax (see Eq. ( 1)).So, in the asymptotic regime (t → ∞), at least one of the following behavior may hold: (i) if log[X(t)/X(t+1)] vs. t decays and goes to zero one has an stretched exponential.(ii) However, if the same plot approaches a positive decay constant one has a simple exponential.Finally, (iii) if a log-log plot of [X(t)/X(t+1)] vs. t/(t+1) shows a lineal decay with slope α < 0 the behavior corresponds to a power law (exponent α) multiplied by a simple exponential decay.Figure 4 shows a plot of log[X(t)/X(t + 1)] vs. t for the relaxation of θ B after the application of a pulse with τ p = 15 MCTs, ∆p = 0.10 and Y Aw = 0.3925.In this example the data is averaged over 10 4 different configurations in order to improve the signal/noise ratio, since X(t)/X(t + 1) is quite sensitive to small fluctuations of the stochastic system.The observed behavior of the data (full dots) strongly supports the stretched exponential decay which can be followed as far as t = 150 MCTs > τ p = 14 MCTs.The horizontal straight line (triangles) corresponds to the expected simple exponential behavior which is far from the obtained data.The inset shows that a log-log plot of [X(t)/X(t+1)] vs. t/(t + 1) exhibits pronounced curvature, so behavior (iii) can also be disregarded.For the example of   5) shows the behavior that would be obtained for a simple exponential decay.The insert shows a log-log plot of [X(t)/X(t + 1)] vs. t/(t + 1).More details in the text.
In the following paragraphs we would like to present a systematic numerical study of the dependence of the relevant quantities which characterize the relaxation process; namely θ Bmax and τ B ; on the parameters of the applied pulse; namely ∆p and τ p .
Let us first discuss results obtained applying different pulses to stationary configurations generated within the reactive regime with Y Aw = 0.3925; i.e. rather close to the second order critical threshold.For this value of the parameter the average coverage with B-species is close to θ ∞ ∼ = 0.9002 ± 0.0005.Figures 5a and 5b show plots of the dependence of τ B and ∆θ Bmax on ∆p, respectively.For peaks of very small amplitude, more exactly for ∆p < 0.002, one has that Y * A = Y Aw − ∆p lies within the reactive regime and the relaxation time is expected to remain almost unchanged, as observed in Figure 5a.However, τ B also remains constant even after driving the system deep in the absorbing state increasing the peak intensity up to ∆p < 0.1.This observation is in contrast with the monotonic increment of ∆θ Bmax which can be very well fitted by a power law behavior of the form with z ∼ = 0.92 ± 0.02, for ∆p < 0.1.Within the interval 0.1 < ∆p < 0.25 the surface becomes almost saturated with B-species; i.e. θ ∞ + ∆θ Bmax ∼ 1; however full saturation is not observed because the poisoning kinetics is delayed by the requirement of two neighboring empty sites for adsorption of B 2 -species.In fact, the mechanism for the generation of such kind of pairs is provided by A-adsorption and subsequent reaction which is heavily affected by the low partial pressure of A-species.Furthermore, the mechanism is almost suppressed for Y A → 0 and one observes a slight decay in ∆θ Bmax (see Fig. 5b).Simultaneously with this decay a drastic drop of the relaxation time is observed (see Fig. 5a).This behavior can be rationalized as follows: in the limit Y A → 0 most empty sites on the surface are single sites because, on the one hand, these sites can not be occupied by the majority species in the gas phase; i.e.B 2 -dimers; and on the other hand, the generation of pairs of empty sites is quite inefficient since A-monomers is the minority species in the gas phase.Due to the presence of these single sites ∆θ Bmax has to drop slightly.After removal of the pulse, back to Y Aw , single sites are quickly occupied by A-monomers which immediately react with neighboring B-species generating additional empty sites in a sort of chain-reaction.This process causes a drastic drop of τ B (Fig. 5a) and the relaxation of the system is faster.
Another approach to the study of the relaxation process for a fixed value of Y Aw is to apply pulses of the same intensity but different period, as it is shown in Figure 6.The dependence of both τ B (Fig. 6a) and ∆θ Bmax (Fig. 6b) on τ p can be described by a power law behavior, that is with x = 1.11 ± 0.03 and y = 0.52 ± 0.02.So, equations ( 2) and ( 4) suggest that, in the limit ∆p → 0 and τ p → 0 the dependence of ∆θ Bmax on the parameters of the pulse may be well described by the following approach Furthermore, the plateau exhibited by τ B in Figure 5a and equation (3) suggest that independent of ∆p, for ∆p → 0 and τ p → 0. The relaxation of the system for fixed parameters of the applied pulse but different values of the station- ary pressure Y Aw (within the range Y A1 < Y Aw < Y A2 ) has also been studied.Figure 7a shows plots of τ B versus ∆Y A = Y Aw − Y A1 obtained for pulses with amplitude ∆p = 0.10 and different periods.For short periods, namely τ p ≤ 100, the relaxation is slow close to criticality showing a plateau, however increasing ∆Y A and approaching Y A2 the relaxation becomes faster and almost independent of ∆Y A close to the first order irreversible phase transition.Also, as it is shown in Figure 7, finite size effects are irrelevant, at least for the side of the used lattices, namely 128 ≤ L ≤ 512.For larger periods, e.g.τ p ≥ 200 and close to the second order transition (∆Y A → 0), the system becomes irreversibly saturated with B-species.Therefore, one has an irreversible phase transition between an active regime with AB-production and an absorbing (inactive) state which is dynamically driven by the external parameter.The detailed study of this interesting phenomena is beyond the aim of the present work and will be addressed in forthcoming publications [24,25].Also, Figure 7b shows plots of ∆θ Bmax versus ∆Y A and the observed behavior is consistent with the results shown in Figure 7b.These findings can be rationalized as follows: close to the second order IPT (∆Y A → 0) the coverage with B-species is close to unity (see e.g.Fig. 1) and the pulse drives the system into the absorbing state for a certain period causing a relatively small increment of θ B (Fig. 7b) which becomes close to saturation.Consequently, when the system is driven back to the active regime with Y Aw close to Y A1 , the relaxation is very slow (Fig. 7a).In contrast, far from the IPT at Y A1 (∆Y A → Y A2 − Y A1 ), the coverage with b-species decreases (Fig. 1) and greater changes of θ B are possible, even if the system is not driven into the absorbing state as it happens for ∆Y A ≥ 0.10.In this case the relaxation becomes faster because Y Aw is close to Y A2 and the higher pressure of A-species allows a quick evolution into the stationary state.

Conclusions
A numerical study on the dynamic response of the ZGB surface reaction lattice gas model is presented.The reaction system, under stationary conditions, is shortly drive into the absorbing state and subsequently, the relaxation process back to the initial state is studied by means of Monte-Carlo simulations.After the application of short pulses, the relaxation process can be well described by a stretched exponential behavior.The dependence of the relaxation characteristic time and the induced changes on the coverage of the reactants, on both, the intensity and the period of the pulsed perturbation, are systematically investigated.A novel irreversible transition driven by a pulsed variation of the external parameter is identified.The obtained results can straightforwardly be extended to a wide variety of irreversible systems exhibiting second order IPT's, such as directed percolation, forest fire models, the contact process, branching annihilating walkers, catalyzed reactions, etc.So, we expects that these numerical results may stimulate theoretical studies on the dynamic response of irreversible systems.In contrast to the case of reversible systems, this field remains unexplored.

Fig. 1 .
Fig. 1.Plots of the rate of AB production (RAB) and the surface coverage with A (θA) and B (θB) species versus YA, for the ZGB model.

Figure 2
Figure2shows a typical response of the reaction system to a pulsed external perturbation.The initial starting
e. the same coverage as in the stationary state at Y Aw ; τ B is the characteristic relaxation time and b is the stretching exponent of the stretched exponential curve.For the case of Figure 3a we have obtained τ B ∼ = 15.23±0.02< τ p = 20 and b ∼ = 0.963 ± 0.003.Similarly, for the case of Figure 3b one has τ B ∼ = 35.93±0.05> τ p = 20 and b ∼ = 0.862±0.002.
Figure 4 we obtain b ' 0.65, so as in the examples shown in Figure 3, the stretched exponent in close to unity.Increasing the period of the pulse b also increases and it becomes almost impossible to distinguish the type of decay, e.g. for τ p > 45 MCTs.

Fig. 4 .
Fig. 4. (•) Plot of log [X(t)/X(t + 1)] vs. t (for the definition of X(t) see the text).Data taken for the relaxation of the system with τp = 15 MCTs, YAw = 0.3925 and ∆p = 0.10.The horizontal line(5) shows the behavior that would be obtained for a simple exponential decay.The insert shows a log-log plot of [X(t)/X(t + 1)] vs. t/(t + 1).More details in the text.

Fig. 5 .
Fig. 5. (a) and (b) show plots of the dependence of τB and ∆θBmax on ∆p obtained for YAw = 0.3925, respectively.In (b) the slope of the straight line is z = 0.92 (for ∆p < 0.1), and has been drawn for comparison.