Characterization of the phase transitions in the ZiffGulariBarshad model with spatial diffusion under two different prescriptions
Abstract
In this work, we modify the ZiffGurariBarshad (ZGB) model by including the spatial diffusion of oxigen atoms and of carbon monoxide molecules, both adsorbed on the lattice, in order to study its effect on the continuous and discontinuous phase transitions of the model through time dependent Monte Carlo simulations. We use an optimization method based on a concept known as coefficient of determination to construct color maps and obtain the phase transitions when the parameters of the model vary, i.e., the rates of adsorption of carbon monoxide and of diffusion. This method was proposed recently to locate the nonequilibrium secondorder phase transitions and has been successfully used in both systems: with defined Hamiltonian and with absorbing states. We also use two different approaches: in the first one, we consider or the adsorption process XOR (Exclusive OR) the diffusion process in each trial (never both in the same trial); in the second one, we consider the adsorption process OR the diffusion process in each trial (same trial can include both). Both approaches present two interesting results: the continuous phase transition point of the original ZGB model disappears when the diffusion rate increases. However, the discontinuous phase transition is not destroyed even when the diffusion rate approaches to 1, indicating that there exists a line of discontinuous phase transition points.
I Introduction
In a pioneering work back in 1986, Ziff, Gulari, and Barshad proposed a model that mimics reactions on a catalytic surface ziff1986 () to describe some nonequilibrium aspects of the production of carbon dioxide molecules through the reaction of carbon monoxide molecules with oxygen atoms. This stochastic model, also known as ZGB model, has attracted much interest due to its simplicity and rich phase diagram with continuous and discontinuous phase transitions which separate absorbing phases from a reactive steady state phase ziff1986 (); meakin1987 (); dickman1986 (); fischer1989 ().
In addition to the scientific interest, the reason for the increasingly study of catalytic reaction models on surfaces lies in the possible applications to the technology bond1987 (); zhdanov1994 (); marro1999 (). In this respect, the ZGB model has been vastly studied nowadays and is considered a prototype for the study on catalytic surfaces due to its inherent reaction processes that also take place in industry which deals with oxidation of CO on transitionmetal catalysts.
Nevertheless, some aspects of the catalytic reaction can not be explained by this model. While some experimental works on platinum confirm the existence of discontinuous irreversible phase transitions (IPT) in the catalytic oxidation of CO golchet1978 (); matsushima1979 (); ehsasi1989 (); christmann1991 (); block1993 (), there is no experimental evidence of the continuous IPT. From the theoretical point of view, the continuous phase transition has been studied by several authors and the results support that this transition belongs to the directed percolation (DP) universality class janssen1989 (); grinstein1989 ().
Modified versions of the ZGB model has been proposed in order to obtain systems with actual catalytic processes and/or to eliminate the continuous phase transition of the original model. Some modifications include CO desorption fischer1989 (); dumont1990 (); jensen1990a (); albano1992 (); tome1993 (); kaukonen1989 (); buendia2009 (), diffusion ehsasi1989 (); fischer1989 (); kaukonen1989 (); jensen1990a (); grandi2002 (), impurities hoenicke2000 (); buendia2012 (); buendia2013 (); hoenicke2014 (); buendia2015 (), attractive and repulsive interactions between the adsorbed molecules satulovsky1992 (), surfaces of different geometries meakin1987 (); albano1990 () and with hard oxygen boundary conditions brosilow1993 (), etc. In addition, it has been studied through several techniques, such as simulations, meanfield theories, series analysis, etc marro1999 ().
In this paper, we perform shorttime Monte Carlo (MC) simulations in order to explore the behavior of the phase transitions of the ZGB model when both atoms and molecules adsorbed on the catalytic surface are able to move, i.e., when diffusion is allowed. By using timedependent MC simulations and an optimization method known as coefficient of determination trivedi2002 (), we construct maps and obtain the transitions when the parameters of the model varies. As shown in the next section, the only two parameters considered in our study are the rates of adsorption of molecules and the rates of diffusion. This technique has been used in the study of reversible systems roberto2012 (); roberto2013a (); roberto2013b (); roberto2014 () and was considered recently in the study of an epidemic model to determine its critical immunization probability roberto2015 (), and to obtain the continuous transition point and the upper spinodal point of the ZGB model fernandes2016 (). In addition, two different ways of performing the trials are proposed.
Our work is organized as follows: In the next section, we present the model with diffusion and show the two manners of moving the atoms/molecules considered in this paper. In Sec. III we describe some finitesize scaling relations in nonequilibrium systems with absorbing states and show the power laws obtained in the criticality as well as the optimization method know as coefficient of determination. Our main numerical results are presented in Sec. IV. Finally, in Sec. V we summarize and conclude our work.
Ii The Model
In this work, we study the dynamical behavior of a modified version of the ZGB model, first devised in 1986 by R.M. Ziff, E. Gulari, and Y. Barshad ziff1986 (), in order to describe the production of carbon dioxide molecules through the reaction of carbon monoxide molecules and oxygen atoms on a catalytic surface. In other words, the ZGB model is a dimermonomer model which simulates the catalysis between the carbon monoxide molecule and the oxygen atom ziff1986 (); evans1991b (). In our approach, in addition to the well known reactions presented in the original model,
(1)  
(2)  
(3) 
there are other two reactions that are related to the diffusion of molecules or atoms adsorbed on the surface:
(4)  
(5) 
In these reactions, stands for vacant sites on the surface, and and refer, respectively, to the gas and adsorbed phases of the atoms/molecules. The indexes and stand for nearest neighbor sites and are explicitly presented in Eqs. (4) and (5) in order to represent the diffusion of the atoms/molecules if there exists some vacant site in the vicinity of them. As shown in the above equations, and molecules exist only in the gas phase while atoms are present only in the adsorbed phase and molecules are able to be in both phases.
In order to simulate such a model, one can think of the catalytic surface as a regular square lattice with its sites occupied by molecules or by atoms or even be empty . The simulation is carried out as follows ziff1986 (); albano1996 (); marro1999 (): By following Eq. (1), the molecule in the gas phase is chosen to impinge on the surface with a rate . The molecule strikes the lattice in a site previously chosen at random. If the site is vacant , then the molecule is adsorbed on it. Otherwise, if the site is occupied by a molecule or by an atom, the trial ends, the molecule returns to the gas phase, and a new molecule is chosen. On the other hand, from Eq. (2), if an molecule in the gas phase is chosen to hit the lattice, it does so at a rate . In this case, a nearestneighbor pair of sites is chosen at random, and, if both sites are vacant , the molecule dissociates into a pair of atoms and are adsorbed on the chosen lattice sites. However, if one or both sites are occupied, the trial ends, the molecule returns to the gas phase, and a new molecule is chosen. As can be seen, these rates are relative ones, , and the model, as presented in Ref. ziff1986 (), has a single free parameter: . Eq. (3) is related to the reaction between the atom and the molecule, both adsorbed on the lattice. Immediately after each adsorption event [Eqs. (1) and (2)], the nearestneighbor sites of the adsorbed molecule/atom are checked. If an pair is found, a molecule is formed and quits the lattice, leaving two empty sites on it. However, if there is the possibility of formation of two or more pairs, a pair is chosen at random to quit the lattice. The last two equations represent the diffusion process proposed in this work: A site is chosen at random and, if it is occupied by an atom or by a molecule, Eqs. (4) and (5) respectively, a nearestneighbor site is also chosen at random. If is vacant, the atom/molecule moves to it with a rate leaving the site empty. It is worth to mention that if there is diffusion, Eq. (3) must be used again to verify if there is formation of molecule or not.
The ZGB model possesses two phase transitions. The first one is a continuous phase transition and occurs at the critical point voigt1997 (); fernandes2016 (). The second transition is discontinuous and occurs at ziff1992 (). For the surface becomes irreversibly poisoned (saturated) by atoms (poisoned state) and for the surface becomes irreversibly poisoned by molecules (poisoned state). The poisoned state, or absorbing phase, represents states in which, once reached, the systems become trapped and can not escape anymore. On the other hand, for there exists a reactive steady state with sustainable production of molecules. So, both and are irreversible phase transition (IPT) points between the reactive and poisoned states.
In this study, we are wondering if the diffusion of the atoms/molecules influences the first and second order phase transitions of the ZGB model. In order to answer this question, we separate the MC simulations into two different moments: the first moment is related to the diffusion process given by Eqs. (4) and (5), and the second one deals with the catalytic reaction given by Eqs. (1), (2), and (3). In addition, we proposed two different ways of performing the simulations:

Algorithm I: Diffusion process “XOR” adsorption process – With a rate , the diffusion process is chosen (first moment). On the other hand, the adsorption process occurs with a rate , and, in this case, the process is governed by the parameter as well as by the rules defined previously (second moment). In summary, both the diffusion and adsorption processes are allowed but not in the same trial.

Algorithm II: Diffusion process “OR” adsorption process – The diffusion process of the atoms/molecules occurs with a rate and then the adsorption process can occur with a rate . In summary, both process can occur in the same trial.
Iii Time dependent Monte Carlo simulations
In this work, we carry out shorttime MC simulations to study the behavior of the ZGB model with diffusion of the atoms/molecules adsorbed on the lattice. In order to perform the numerical simulations, we take into consideration that, for systems belonging to the DP universality class, the finitesize scaling near criticality can be described by:
(6) 
where is the density of molecules, stands for the average on different evolutions of the system, is the dimension of the system (for the ZGB model, ), is its linear size of a regular square lattice, and is the time. The indexes and are dynamic critical exponents, and , , and are static ones. Here, denotes the distance of a point to the critical point, , which governs the algebraic behaviors of the two independent correlation lengths: the spatial one which behaves as and the temporal one, .
The density of molecules adsorbed on the lattice is given by
and when the sites are occupied by molecules, otherwise, it is equal to zero.
As shown in Ref. fernandes2016 (), the dynamic and static critical exponents of the model can be obtained by using the Eq. (6) and performing timedependent MC simulations with two different initial conditions : (i) the lattice is completely empty, i.e. there exist only vacant sites () and (ii) the lattice is completely filled with atoms but a random site which remains empty (). When considering the first condition, one obtains
(7) 
where , and the second condition produces
(8) 
The numerical computing of the exponents , , and , follows a straightforward procedure: First, by considering these two different initial conditions, it is possible to obtain the exponent in a very simple way, leading to the following power law . In addition to the analysis of these power laws in the study of the original ZGB model fernandes2016 (), this idea has been applied successfully in a large number of spin systems: for example, the Ising model, the and Potts models Silva20021 (), Heisenberg model Fernandes2006c () and even for models based on the generalized Tsallis statistics Silva2012Tsallis (). It was also introduced recently in systems without defined Hamiltonian, as can be seen in Ref. Roberto2004 (). Second, in order to compute , another power law is obtained when taking into account the following derivative Grassberger1996 (): which is numerically represented by , where is a tiny perturbation needed to move the system slightly off the criticality, yelding a power law decay that only depends on , i.e, .
However to determine these exponents by using this tool of power laws, the critical parameters must be accurately determined. Can this task also be determined via timedependent MC simulations? In this paper we are concerned with the influence of the diffusion on the existence of the continuous and discontinuous irreversible phase transitions of the ZGB model. So, instead of focusing on the estimate of the critical exponents, we decided to take into account a very simple and efficient optimization method to obtain an overview of the phase transitions when the rate of adsorption of molecules on the lattice and the rate of diffusion vary from zero to one, i.e., and .
This method considers the robustness of the power laws, and it is based on a quantity known as coefficient of determination. It is given by
(9) 
where is obtained for each pair (), is the number of MC steps, and are, respectively, the intercept and the slope of a linear function, and . Here differently from which denotes an average over the MC steps, denotes an average over the different runs of simulation, i.e., an average over the different time evolutions considering different set of random numbers.
This coefficient ranges from 0 to 1 and measures the ratio: (expected variation)/(total variation), so that the bigger the , the better the linear fit of the data in logscale. When the system is out of criticality, there is no power law and . However, at criticality ( and ) or close to it, it is expected that possesses a power law behavior given, for instance, by Eq. (7) and approaches to 1.
Therefore, we perform nonequilibrium MC simulations for each pair and , with and , in order to obtain color maps for the coefficient of determination as function of each pair . Such maps are able to show the existence (or not) of phase transitions when diffusion of atoms and molecules are allowed. The pair of critical values corresponds to . A similar procedure has been used by the authors to the study the effects of mobility on diluted lattices in an epidemic model dasilva2015 (). In that work, they showed that mobility influences the immunization rates in such a way that the greater the mobility, the bigger the immunization rate.
It is important to mention that this technique can also be extended to study the weak firstorder phase transitions fernandes2016 (), since these transitions possess long correlation lengths and small discontinuities and therefore behave similarly to secondorder phase transitions Schulke2001 (); Albano2001 (). It has been conjectured that near a weak firstorder transition there exist two pseudocritical points: one point is just below (lower) the firstorder phase transition point, and the other one is just above (upper) it. These pseudocritical points are known as spinodal points Schulke2001 ().
Iv Results
In this work, we perform nonequilibrium MC simulations in order to obtain color maps for the phase transitions of the ZGB model with diffusion of atoms and molecules adsorbed on the lattice. The results of our simulations are obtained by considering the coefficient of determination given by Eq. (9). As stressed at the end of Sec. II, we proposed two different approaches for each trial: Algorithm I, which considers that only one process (adsorption or diffusion) can occur in a given trial, and the Algorithm II in which both adsorption and diffusion processes are able to happen simultaneously.
In order to have a clue of the phase transitions of the model without spending too much computational time, we first performed rough simulations to be able to scan the entire space of parameters and , i.e., and , and construct complete color maps for the model. With these results in hand, we focused our attention in some regions of the color map where the transitions occur and refined our simulations, improving the accuracy of the results. In this study, we consider and the parameters used in each stage are resumed in Table 1.
Stage  window  

1 – Rough simulations  200  1000  
2 – Refined simulations  500  5000 
Following, we present our main results obtained from the Algorithm I and Algorithm II.
iv.1 Algorithm I  Diffusion process “XOR” adsorption process
In Fig. 1 (a) we show the color map which represents the possible regions of phase transitions of the ZGB model with diffusion.
As can be seen, there are two regions that deserve more attention whereas their points are candidates to the continuous and discontinuous phase transition points, i.e., . The first one is around the continuous phase transition of the original ZGB model ( and ), and the second one is close to (the discontinuous phase transition point of the original ZGB model) and ranges from 0 to 1. At the top of this figure (), there is a third region: a yellow horizontal strip that ranges from . This region means that, in the most of cases, only the diffusion is allowed since and, although , the points in this region are not phase transition points as in the first two regions.
These “fake” phase transition points present exponent as suggested by Fig. 1 (b), which shows the values of exponent for the same ranges of and . Differently from the “fake” points of the horizontal strip, the discontinuous phase transition (candidate) points along the vertical strip in Fig. 1 (a) present as again can be observed in Fig. 1 (b). Such behavior includes the known , which reinforces our confidence that diffusion transforms this point into an entire line of discontinuous phase transition points.
Moreover, the positivity for the exponents of these discontinuous phase transition points seems to be a good indicator to separate them from the continuous phase transition points since the last ones present negative exponents , that can be verified by taking the elliptical yellow region around (continuous phase transition points) in Fig. 1 (a) and their corresponding exponents in Fig. 1 (b).
So, by looking into the Figs. 1 (a) and (b), we are able to conclude that the diffusion eliminates the continuous phase transition and moves the discontinuous phase transition toward higher values of . With these results in hand, we simulate the following two regions (refined simulations):

Region 1: and , and

Region 2: and .
Figure 2 (a) shows the region 1 and 2 (b) shows the region 2, using the parameters presented in the second line of Table 1.
From Fig. 2 (a), one can see that is higher for smaller values of meaning that the continuous phase transition of the original ZGB model is destroyed for . Details about elliptical region near of expected critical point and can be observed. On the other hand, the discontinuous phase transition remains even for , as can be observed in Fig. 2 (b).
In Fig. 3 we show the results from a simple procedure in which for each value of , we span all different values of , considering the region 1: 3 (a) and the region 2: 3 (b), and we plot only for the points corresponding to the highest value of found by the procedure. We can observe that Fig. 3 (a) indicates that only for values really low of (yellow points) we have candidates to continuous transition points, suggesting that the diffusion does not change the continuous transition point of the original ZGB model. On the other hand, Fig. 3 (b) shows that a notorious line of discontinuous transition points are created.
iv.2 Algorithm II  Diffusion process “OR” adsorption process
Finally, we performed simulations following the Algorithm II in which both adsorption and diffusion processes can occur in the same trial, i.e., diffusion does not exclude the possibility of adsorption. Thus, we carried out similar simulations and obtained plots as those shown previously in the Figs. 1, 2, and 3. The results obtained by considering the Algorithm II are collected and presented in Fig. 4.
In Fig. 4, the plots (a) and (b) correspond to the plots (a) and (b) of Fig. 1 respectively, the plots (c) and (d) correspond to the plots (a) and (b) of Fig. 2 respectively, and finally, the plots (e) and (f) correspond to the plots (a) and (b) of Fig. 3.
Exactly as shown for the Algorithm I, the continuous phase transition is destroyed for and the discontinuous phase transition is maintained even for large values of . In fact, it is important to notice that the discontinuous phase transition points remain practically the same for as can be observed in Fig. 4 (a) which shows the coefficient of determination . This behavior is corroborated by the value of the exponent according to Fig. 4 (b). There is a small variation as suggested by the refined simulations of the region of discontinuous phase transition considered in Fig. 4 (d) but very different from Fig. 2 (b) which shows a strong curvature to the right side. We observe that the “fake” right discontinuous points disappear since diffusion and adsorption can occur simultaneously. However, regardless of the algorithm, Figs. 2 (b) and 4 (d) indicate that there exists a line of discontinuous phase transition for the ZGB model with diffusion. In addition, Figs. 2 (a) and 4 (c) are very similar and show that the diffusion eliminates the continuous phase transition of the original ZGB model.
The plot 4 (e) shows that highest values of the coefficient of determination are obtained only at the vicinity of the critical point of the original ZGB model, i.e., the diffusion destroys the continuous phase transition, as adressed before. However, it is important to notice that power laws consistent with continuous phase transition points are obtained for . And finally, 4 (f) shows that highest coefficient of determination suggests that the discontinuous phase transition points present small variations, , for which is close of the discontinuous phase transition point of the original ZGB model, and .
V Conclusions
In this work, we studied the effects of the spatial diffusion on the phase transitions of the ZiffGulariBarshad (ZGB) model by using an alternative method that optimizes the coefficient of determination in order to localize the critical parameters of the continuous phase transition point and to obtain an estimate of the upper spinodal point (one of the pseudo critical points of the model) of the weak discontinuous phase transition point.
Our results show that critical behavior is strongly changed by introducing the diffusion of the atoms and molecules adsorbed on the lattice, leading to the extinction of the criticality even for small values of . On the other hand, we have an extension of discontinuous phase transition points for all values of the diffusion rates indicating that there exists a line of discontinuous transition points. For the Algorithm I, this line has a small curvature to the right and “fake” discontinuous transition points appear when . For the Algorithm II these points are eliminated whereas both diffusion and adsorption processes are able to occur simultaneously. In addition, all points on the line of discontinuous phase transition obtained with the Algorithm II are close to the transition point of the original ZGB model.
Acknowledgments
This research work was in part supported financially by CNPq (National Council for Scientific and Technological Development).
References
 (1) R.M. Ziff, E. Gulari, and Y. Barshad, Phys. Rev. Lett. 56, 2553 (1986).
 (2) P. Meakin and D.J. Scalapino, J. Chem. Phys. 87, 731 (1987).
 (3) R. Dickman, Phys. Rev. A 34, 4246 (1986).
 (4) P. Fischer and U.M. Titulaer, Surf. Sci. 221, 409 (1989).
 (5) G.C. Bond, Catalysis: Principles and Applications (Clarendon, Oxford, 1987).
 (6) V.P.Z. Zhdanov and B. Kazemo, Surf. Sci. Rep. 20, 113 (1994).
 (7) J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, U.K., 1999).
 (8) A. Golchet and J.M. White, J. Catal. 53, 266 (1978).
 (9) T. Matsushima, H. Hashimoto, and I. Toyoshima, J. Catal. 58, 303 (1979)
 (10) M. Ehsasi, M. Matloch, J.H. Block, K. Christmann, F.S. Rys, and W. Hirschwald, J. Chem. Phys. 91, 4949 (1989).
 (11) K. Christmann, Introduction to Surface Physical Chemistry (Steinkopff Verlag, Darmstadt, 1991), pp. 1274.
 (12) J.H. Block, M. Ehsasi, and V. Gorodetskii, Prog. Surf. Sci. 42, 143 (1993).
 (13) H.K. Janssen, B. Schaub, and B. Schmittmann, Z. Phys. B: Condens. Matter 73, 539 (1989).
 (14) G. Grinstein, Z.W. Lai, and D.A. Browne, Phys. Rev. A 40, 4820 (1989).
 (15) M. Dumont, P. Dufour, B. Sente, and R. Dagonnier, J. Catal. 122, 95 (1990).
 (16) E.V. Albano, Appl. Phys. A 54, 2159 (1992).
 (17) T. Tomé and R. Dickman, Phys. Rev. E 47, 948 (1993).
 (18) H.P. Kaukonen and R.M. Nieminen, J. Chem. Phys. 91, 4380 (1989).
 (19) I. Jensen and H. Fogedby, Phys. Rev. A 42, 1969 (1990).
 (20) G.M. Buendia , E. Machado, and P.A. Rikvold, J. Chem. Phys. 131, 184704 (2009).
 (21) B.C.S. Grandi and W. Figueiredo, Phys. Rev. E 65, 036135 (2002).
 (22) G.L. Hoenicke and W. Figueiredo, Phys. Rev. E 62, 6216 (2000).
 (23) G. M. Buendía and P.A. Rikvold, Phys. Rev. E 85, 031143 (2012).
 (24) G. M. Buendía and P.A. Rikvold, Phys. Rev. E 88, 012132 (2013).
 (25) G.M. Buendía, P.A. Rikvold, Phys. A. 424, 217 (2015).
 (26) G.L. Hoenicke, M.F. de Andrade, and W. Figueiredo, J. Chem. Phys. 141, 074709 (2014).
 (27) J. Satulovsky and E.V. Albano, J. Chem. Phys. 97, 9440 (1992).
 (28) E.V. Albano, Surf. Sci. 235, 351 (1990).
 (29) B.J. Brosilow, E. Gulari, and R.M. Ziff, J. Chem. Phys. 98, 674 (1993).
 (30) K. S. Trivedi, Probability and Statistics with Realiability, Queuing, and Computer Science and Applications, 2nd ed. (John Wiley and Sons, Chichester, 2002).
 (31) R. da Silva, J.R. Drugowich de Felício, A.S. Martinez, Phys. Rev. E. 85, 066707 (2012).
 (32) R. da Silva, N. Alves Jr., J.R. Drugowich de Felicio, Phys. Rev. E 87, 012131 (2013).
 (33) R. da Silva, H.A. Fernandes, J.R. Drugowich de Felício, W. Figueiredo, Comput. Phys. Commun. 184, 2371 (2013).
 (34) R. da Silva, H.A. Fernandes, J.R. Drugowich de Felício, Phys. Rev. E 90, 042101 (2014).
 (35) R. da Silva, H.A. Fernandes, J. Stat. Mech. P06011 (2015)
 (36) H.A. Fernandes, R. da Silva, E.D. Santos, P.F. Gomes, and E. Arashiro, Phys. Rev. E 94, 022129 (2016).
 (37) W. Evans and M. S. Miesch, Phys. Rev. Lett. 66, 833 (1991).
 (38) E.V. Albano, Chem. Rev. 3, 389 (1996).
 (39) C.A. Voigt, R.M. Ziff, Phys. Rev. E 56 R6241 (1997)
 (40) R.M. Ziff and B.J. Brosilow, Phys. Rev. A, 46 4630 (1992).
 (41) R. da Silva, N.A. Alves and J.R. Drugowich de Felício, Phys. Lett. A 298, 325 (2002).
 (42) H.A. Fernandes, Roberto da Silva, and J.R. Drugowich de Felicio, J. Stat. Mech.: Theor. Exp., P10002 (2006).
 (43) R. da Silva, J. R. Drugowich de Felicio, A. S. Martinez, Phys. Rev. E, Statistical, 85, 066707 (2012).
 (44) R. da Silva, R. Dickman, and J.R. Drugowich de Felício, Phys. Rev. E 70, 067701 (2004)
 (45) P. Grassberger and Y. Zhang, Physica A 224, 169 (1996)
 (46) R. da Silva and H.A. Fernandes, J. Stat. Mech. P06011 (2015).
 (47) L. Schulke and B. Zheng, Phys. Rev. E 62, 74827485 (2000)
 (48) E. Albano, Physics Letters A 288, 73–78 (2001)