- Email: [email protected]
Monte Carlo study of the critical behavior of a diffusive epidemic process
Physica A 295 (2001) 49–52
www.elsevier.com/locate/physa
Monte Carlo study of the critical behavior of a di$usive epidemic process U.L. Fulco, D.N. Messias, M.L. Lyra ∗ Departamento de F sica, Universidade Federal de Alagoas, 57072-970 Macei o - AL, Brazil
Abstract We perform Monte Carlo simulations of a symmetric 1D di$usive epidemic propagation process. The system presents a continuous phase transition between a stationary vacuum state and a steady reactive state. Renormalization group calculations have predicted that the dynamical critical exponent z = 2 and the correlation exponent = 2. However, recent numerical results have reported = 2:21. In this work we employ a 6nite size scaling analysis of order parameter and relaxation time data at the vicinity of the critical point and show that, once corrections to scaling are taken into account, the numerical data are consistent with the RG c 2001 Elsevier Science B.V. All rights results. Further, we report a precise estimate for =. reserved. PACS: 75.70.Ak; 75.40.Mg; 75.40.Cx Keywords: Phase transitions; Di$usive epidemic propagation; Finite size scaling
Recently, the propagation of an epidemic process in a population of >uctuating density has been studied using the Wilson renormalization group method [1]. In the model healthy and sick individuals (A and B particles) di$use independently with di$usion constants DA and DB . Upon contact, sick individuals may infect healthy ones at a rate k1 . They also can spontaneously recover at a rate k2 . Therefore, a competition between the contamination process (creation of B particles) and the recovery process (annihilation of B particles) takes place. For low concentrations of the average total density , the stationary state is characterized by a global extinction of the epidemics characteristic of stochastic systems with absorbing states [2– 4]. Above a critical density c there is a stable steady-state regime with a >uctuating 6nite density of sick ∗
Corresponding author. Fax: +55-82-214-1645. E-mail address: [email protected] (M.L. Lyra).
c 2001 Elsevier Science B.V. All rights reserved. 0378-4371/01/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 1 ) 0 0 0 5 0 - 4
50
U.L. Fulco et al. / Physica A 295 (2001) 49–52
individuals. Near c the system exhibits a phase transition with the average density of sick individuals B acting as the order parameter [1,5]. This model can be formulated as a reaction-di$usion-decay process A + B → 2B, B → A. The d=1 version of the symmetric DA =DB case has been recently investigated in a Monte Carlo simulation by de Freitas et al. [6]. They reported = = 0:197(2). By 6tting extrapolated order parameter data away from the critical point, they found = 0:435(10) and = 2:21(5). Theoretically the critical properties of the stationary state was 6rst studied by Kree et al. in the context of a population density in a polluted environment [7]. Taking into account the total density >uctuations and using 6eld-theoretic renormalization group techniques, Kree et al. [7] computed the critical exponents to be = − =8 in 6rst order in = 4 − d, = 2=d and z = 2 in all orders in . The value of = 2 found on the numerical simulation in 1D, indicates that higher order symmetry-breaking terms in the functional action seem to be relevant at low dimensions or, alternatively, that corrections to scaling are needed for the system sizes simulated. In this work, we perform a Monte Carlo simulation of the above di$usion-reactiondecay process in d = 1 and DA = DB . We will use a 6nite size scaling analysis of order parameter and relaxation time data to precisely locate the critical point and to directly compute the critical exponents =, z and . In our simulations, we consider each site as being a locus where an arbitrary number of particles can be located. The contamination process takes place only when A and B particles are found on the same site. The whole di$usion-reaction-decay process is done in three stages. First, we consider only the di$usive motion. Particles A and B are chosen to move with probabilities pA =pB = 12 , respectively. Once a particle is chosen to move, it is displaced to one of the two neighboring sites with equal probability. In the next stage, we consider only the contamination process. Each A particle that is in the same site of at least one B particle, is transformed in a B particle with probability k1 . In the last stage each B particle can be replaced by an A particle with probability k2 . Our simulations are performed using k1 = k2 = 12 . The time unit will be considered as the time needed to perform the above three stages over all particles. To estimate the critical density c we assume that the order parameter B satis6es the scaling relation B (; L) = L−= f [(L1= ( − c )] :
(1)
The preceding relation imply that the set of auxiliary functions g(L; L ; ) = ln[B (L; )=B (L ; )]=ln(L=L )
(2)
intersect in a single point (c ; =). In Fig. 1 we plot the average relative density of sick individuals B = at the stationary regime versus the total density for a system of linear size N = 400. In the inset we employ the above phenomenological renormalization to 6nd c = 4:24(1).
U.L. Fulco et al. / Physica A 295 (2001) 49–52
51
Fig. 1. The stationary value of B = as a function of the total density for pA = pB = 0:5. The results were obtained from 3200 realizations over a chain with L = 400 sites. The stationary state is assumed to be achieved after 12800 runs over the entire chain. For each realization, 400 microscopic con6gurations of the stationary regime were used to account for temporal >uctuations. In the inset we show the function g(L; L ; ) for several pairs (L; L ) The common point determines c = 4:24(1).
At the vicinity of the inactive-active transition the typical relaxation time of the sick individuals concentration is expected to scale as = Lz [( − c )L1= ] ;
(3)
so that at the critical point the relaxation time scales as ˙ Lz . There are several distinct forms of measuring typical relaxation times in a simulated dynamical process. Particularly, whenever the order parameter vanishes assymptotically, we can measure ∞ t (t)dt B = 0∞ : (4) B (t)dt 0 In Fig. 2 we compute the exponent z by plotting ln() at criticality for several values of ln L. The data are very well 6tted by = L2 (1 + a=L). These are consistent with the RG prediction and indicate that a small correction to scaling is needed to properly 6t these 6nite size data. In the inset we also report the size dependence of ≡ B =|c ˙ L−= and 1= (d )=d|c ˙ L1= . The data are consistent with = 2 if a small correction to scaling is allowed and with = = 0:226(15). In summary, the present results corroborate the renormalization group prediction that the correlation exponent = 2 and relaxation exponent z = 2. We have found that small corrections to scaling are needed to properly extract the correct exponents from Monte Carlo data from the 6nite size systems simulated. It would be valuable to have an analytical study of the above system employed directly in 1D to verify if the order parameter exponent can be correctly derived without the need of higher order terms in the action functional.
52
U.L. Fulco et al. / Physica A 295 (2001) 49–52
Fig. 2. ln() versus ln L at the critical point. The line is a 6t considering a 6rst order correction to scaling. It corroborates the renormalization group prediction that z = 2. The inset shows ln and ln(1= (d =d)) versus ln L at the critical point where ≡ B =. The lines are 6ts considering a 6rst order correction to scaling. The results are consistent with = 2 and = = 0:226(15).
The authors would like to thank L. Lucena and G. Viswanathan for fruitful discussions. This work was partially supported by the Brazilian research agencies CNPq and CAPES and by the Alagoas State agency FAPEAL. References [1] F. van Wijland, K. Oerding, H.J. Hilhorst, Physica A 251 (1998) 179. [2] N.G. Van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1992. [3] J. Marro, R. Dickman, Nonequilibrium Phase Transitions in Lattice Models, Cambridge University Press, Cambridge, 1999, and references therein. [4] R. Dickman, J.K.L. da Silva, Phys. Rev. E 58 (1998) 4266. [5] K. Oerding, F. van Wijland, J.P. Leroy, H.J. Hilhorst, J. Stat. Phys. 99 (2000) 1365. [6] J.E. de Freitas, L.S. Lucena, L.R. da Silva, H.J. Hilhorst, Phys. Rev. E 61 (2000) 6330. [7] R. Kree, B. Schaub, B. Schmittmann, Phys. Rev. A 39 (1989) 2214.
www.elsevier.com/locate/physa
Monte Carlo study of the critical behavior of a di$usive epidemic process U.L. Fulco, D.N. Messias, M.L. Lyra ∗ Departamento de F sica, Universidade Federal de Alagoas, 57072-970 Macei o - AL, Brazil
Abstract We perform Monte Carlo simulations of a symmetric 1D di$usive epidemic propagation process. The system presents a continuous phase transition between a stationary vacuum state and a steady reactive state. Renormalization group calculations have predicted that the dynamical critical exponent z = 2 and the correlation exponent = 2. However, recent numerical results have reported = 2:21. In this work we employ a 6nite size scaling analysis of order parameter and relaxation time data at the vicinity of the critical point and show that, once corrections to scaling are taken into account, the numerical data are consistent with the RG c 2001 Elsevier Science B.V. All rights results. Further, we report a precise estimate for =. reserved. PACS: 75.70.Ak; 75.40.Mg; 75.40.Cx Keywords: Phase transitions; Di$usive epidemic propagation; Finite size scaling
Recently, the propagation of an epidemic process in a population of >uctuating density has been studied using the Wilson renormalization group method [1]. In the model healthy and sick individuals (A and B particles) di$use independently with di$usion constants DA and DB . Upon contact, sick individuals may infect healthy ones at a rate k1 . They also can spontaneously recover at a rate k2 . Therefore, a competition between the contamination process (creation of B particles) and the recovery process (annihilation of B particles) takes place. For low concentrations of the average total density , the stationary state is characterized by a global extinction of the epidemics characteristic of stochastic systems with absorbing states [2– 4]. Above a critical density c there is a stable steady-state regime with a >uctuating 6nite density of sick ∗
Corresponding author. Fax: +55-82-214-1645. E-mail address: [email protected] (M.L. Lyra).
c 2001 Elsevier Science B.V. All rights reserved. 0378-4371/01/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 1 ) 0 0 0 5 0 - 4
50
U.L. Fulco et al. / Physica A 295 (2001) 49–52
individuals. Near c the system exhibits a phase transition with the average density of sick individuals B acting as the order parameter [1,5]. This model can be formulated as a reaction-di$usion-decay process A + B → 2B, B → A. The d=1 version of the symmetric DA =DB case has been recently investigated in a Monte Carlo simulation by de Freitas et al. [6]. They reported = = 0:197(2). By 6tting extrapolated order parameter data away from the critical point, they found = 0:435(10) and = 2:21(5). Theoretically the critical properties of the stationary state was 6rst studied by Kree et al. in the context of a population density in a polluted environment [7]. Taking into account the total density >uctuations and using 6eld-theoretic renormalization group techniques, Kree et al. [7] computed the critical exponents to be = − =8 in 6rst order in = 4 − d, = 2=d and z = 2 in all orders in . The value of = 2 found on the numerical simulation in 1D, indicates that higher order symmetry-breaking terms in the functional action seem to be relevant at low dimensions or, alternatively, that corrections to scaling are needed for the system sizes simulated. In this work, we perform a Monte Carlo simulation of the above di$usion-reactiondecay process in d = 1 and DA = DB . We will use a 6nite size scaling analysis of order parameter and relaxation time data to precisely locate the critical point and to directly compute the critical exponents =, z and . In our simulations, we consider each site as being a locus where an arbitrary number of particles can be located. The contamination process takes place only when A and B particles are found on the same site. The whole di$usion-reaction-decay process is done in three stages. First, we consider only the di$usive motion. Particles A and B are chosen to move with probabilities pA =pB = 12 , respectively. Once a particle is chosen to move, it is displaced to one of the two neighboring sites with equal probability. In the next stage, we consider only the contamination process. Each A particle that is in the same site of at least one B particle, is transformed in a B particle with probability k1 . In the last stage each B particle can be replaced by an A particle with probability k2 . Our simulations are performed using k1 = k2 = 12 . The time unit will be considered as the time needed to perform the above three stages over all particles. To estimate the critical density c we assume that the order parameter B satis6es the scaling relation B (; L) = L−= f [(L1= ( − c )] :
(1)
The preceding relation imply that the set of auxiliary functions g(L; L ; ) = ln[B (L; )=B (L ; )]=ln(L=L )
(2)
intersect in a single point (c ; =). In Fig. 1 we plot the average relative density of sick individuals B = at the stationary regime versus the total density for a system of linear size N = 400. In the inset we employ the above phenomenological renormalization to 6nd c = 4:24(1).
U.L. Fulco et al. / Physica A 295 (2001) 49–52
51
Fig. 1. The stationary value of B = as a function of the total density for pA = pB = 0:5. The results were obtained from 3200 realizations over a chain with L = 400 sites. The stationary state is assumed to be achieved after 12800 runs over the entire chain. For each realization, 400 microscopic con6gurations of the stationary regime were used to account for temporal >uctuations. In the inset we show the function g(L; L ; ) for several pairs (L; L ) The common point determines c = 4:24(1).
At the vicinity of the inactive-active transition the typical relaxation time of the sick individuals concentration is expected to scale as = Lz [( − c )L1= ] ;
(3)
so that at the critical point the relaxation time scales as ˙ Lz . There are several distinct forms of measuring typical relaxation times in a simulated dynamical process. Particularly, whenever the order parameter vanishes assymptotically, we can measure ∞ t (t)dt B = 0∞ : (4) B (t)dt 0 In Fig. 2 we compute the exponent z by plotting ln() at criticality for several values of ln L. The data are very well 6tted by = L2 (1 + a=L). These are consistent with the RG prediction and indicate that a small correction to scaling is needed to properly 6t these 6nite size data. In the inset we also report the size dependence of ≡ B =|c ˙ L−= and 1= (d )=d|c ˙ L1= . The data are consistent with = 2 if a small correction to scaling is allowed and with = = 0:226(15). In summary, the present results corroborate the renormalization group prediction that the correlation exponent = 2 and relaxation exponent z = 2. We have found that small corrections to scaling are needed to properly extract the correct exponents from Monte Carlo data from the 6nite size systems simulated. It would be valuable to have an analytical study of the above system employed directly in 1D to verify if the order parameter exponent can be correctly derived without the need of higher order terms in the action functional.
52
U.L. Fulco et al. / Physica A 295 (2001) 49–52
Fig. 2. ln() versus ln L at the critical point. The line is a 6t considering a 6rst order correction to scaling. It corroborates the renormalization group prediction that z = 2. The inset shows ln and ln(1= (d =d)) versus ln L at the critical point where ≡ B =. The lines are 6ts considering a 6rst order correction to scaling. The results are consistent with = 2 and = = 0:226(15).
The authors would like to thank L. Lucena and G. Viswanathan for fruitful discussions. This work was partially supported by the Brazilian research agencies CNPq and CAPES and by the Alagoas State agency FAPEAL. References [1] F. van Wijland, K. Oerding, H.J. Hilhorst, Physica A 251 (1998) 179. [2] N.G. Van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1992. [3] J. Marro, R. Dickman, Nonequilibrium Phase Transitions in Lattice Models, Cambridge University Press, Cambridge, 1999, and references therein. [4] R. Dickman, J.K.L. da Silva, Phys. Rev. E 58 (1998) 4266. [5] K. Oerding, F. van Wijland, J.P. Leroy, H.J. Hilhorst, J. Stat. Phys. 99 (2000) 1365. [6] J.E. de Freitas, L.S. Lucena, L.R. da Silva, H.J. Hilhorst, Phys. Rev. E 61 (2000) 6330. [7] R. Kree, B. Schaub, B. Schmittmann, Phys. Rev. A 39 (1989) 2214.