Critical properties of the contact process with quenched dilution

TL;DR

This study investigates the critical behavior of the contact process on a diluted lattice, revealing its universality class and phase diagram through Monte Carlo simulations and mean-field theory.

Contribution

It introduces a novel simulation approach combining percolation cluster generation with contact process dynamics on diluted lattices.

Findings

Critical exponents align with the random transverse-field Ising model

Phase diagram characterized by mean-field theory considering site correlations

Activated scaling relation confirms dynamic critical behavior

Abstract

We have studied the critical properties of the contact process on a square lattice with quenched site dilution by Monte Carlo simulations. This was achieved by generating in advance the percolating cluster, through the use of an appropriate epidemic model, and then by the simulation of the contact process on the top of the percolating cluster. The dynamic critical exponents were calculated by assuming an activated scaling relation and the static exponents by the usual power law behavior. Our results are in agreement with the prediction that the quenched diluted contact process belongs to the universality class of the random transverse-field Ising model. We have also analyzed the model and determined the phase diagram by the use of a mean-field theory that takes into account the correlation between neighboring sites.

Tables

Table 1. Table 1: Critical exponents obtained by numerical simulations on a square lattice at p = 0.8 𝑝 0.8 p=0.8 and λ c = 2.1075 subscript 𝜆 𝑐 2.1075 \lambda_{c}=2.1075 (second column). The third and fourth columns show results for the quenched diluted contact process [ 10 , 11 ] in d = 2 𝑑 2 d=2 whereas the last two columns show results for the random transverse-field Ising model [ 16 , 20 ] also in d = 2 𝑑 2 d=2 .

ref.	this work	[10]	[11]	[16]	[20]
$\theta$	$0.145\left(8\right)$		$0.15\left(3\right)$
$\psi$	$0.55\left(3\right)$	$0.48\left(7\right)$	$0.51\left(6\right)$	$0.42\left(6\right)$	$0.48\left(2\right)$
$\delta$	$1.88\left(11\right)$		$1.9\left(2\right)$
$\varphi$	$1.87\left(10\right)$			$2.5\left(4\right)$
$\beta /{\nu }_{⟂}$	$0.964\left(2\right)$	$0.95\left(2\right)$	$0.96\left(2\right)$	$1.0\left(1\right)$	$0.982\left(15\right)$
$\beta$	$1.11\left(6\right)$		$1.15\left(9\right)$
${\nu }_{⟂}$	$1.15\left(6\right)$		$1.20\left(15\right)$	$1.07\left(15\right)$	$1.24\left(2\right)$
$d_F$	$1.036\left(2\right)$			$1.0\left(1\right)$	$1.018\left(15\right)$

Equations

$$\mathrm{S}+\mathrm{E}\to \mathrm{U}+\mathrm{E},\mathrm{rate}a,$$

$$\mathrm{S}+\mathrm{E}\to \mathrm{E}+\mathrm{E},\mathrm{rate}b.$$

$$\frac{d}{dt}{P}_S=-(a+b)P_{ES},$$

$$\frac{d}{dt}{P}_E=b{P}_{ES},$$

$$\frac{d}{dt}{P}_{ES}=-\frac{a+b}{k}{P}_{ES}-(a+b)\mu P_{ESE}+b\mu P_{ESS},$$

$$\frac{d}{dt}{P}_{US}=-(a+b)\mu P_{ESU}+a\mu P_{ESS},$$

$$\frac{d}{dt}{P}_{EU}=\frac{a}{k}{P}_{ES}+a\mu P_{ESE}+b\mu P_{ESU},$$

$$\frac{dx}{dt}=-(a+b)v,$$

$$\frac{dy}{dt}=bv,$$

$$\frac{dv}{dt}=-\frac{a+b}{k}v-(a+b)\mu \frac{v^2}{x}+b\mu \frac{v(x-v-u)}{x},$$

$$\frac{du}{dt}=-(a+b)\mu \frac{vu}{x}+a\mu \frac{v(x-v-u)}{x},$$

$$\frac{dw}{dt}=\frac{a}{k}v+a\mu \frac{v^2}{x}+b\mu \frac{vu}{x},$$

$$x=s^k,$$

$$y=p(1-s^k),$$

$$v=0,$$

$$u=q{s}^{k-1}(1-s^{k-1}),$$

$$w=pq(1-s^{k-1}),$$

$$p{s}^{k-1}-s+q=0,$$

$$\mathrm{E}+\mathrm{I}\to \mathrm{I}+\mathrm{I},\mathrm{rate}c,$$

$$\mathrm{I}\to \mathrm{E},\mathrm{rate}r,$$

$$\frac{d}{dt}{P}_I=-r{P}_I+c{P}_{IE},$$

$$\frac{d}{dt}{P}_{IE}=-r{P}_{IE}+r{P}_{II}-\frac{c}{k}{P}_{IE}-c\mu P_{IEI}+c\mu P_{IEE},$$

$$\frac{d}{dt}{P}_{IU}=-r{P}_{IU}+c\mu P_{IEU},$$

$$\frac{dz}{dt}=-rz+cg,$$

$$\frac{dg}{dt}=-rg+r(z-g-h)-\frac{c}{k}g-c\mu \frac{g^2}{y}+c\mu \frac{g}{y}(y-g-w),$$

$$\frac{dh}{dt}=-rh+c\mu \frac{gw}{y},$$

$$\alpha =\frac{r}{c}=\frac{k-1}{k}\left(1-\frac{q(1-s)}{p(1-s^k)}\right),$$

$$\xi \sim (\ln t/t_0{)}^{1/\psi },$$

$$N_I\sim (\ln t/t_0{)}^{\theta },$$

$$P\sim (\ln t/t_0{)}^{-\delta }.$$

Full text

1 Introduction

In many experiments on condensed matter, quenched disorder may be present either because it is an unavoidable feature of the sample or because disorder is deliberated introduced in the sample. In either case, if we wish to describe the properties of these systems by statistical mechanical models, quenched disorder should be taken into account in these models. In some cases the quenched disorder is irrelevant in the sense that it does not change the critical behavior. In other cases, the quenched disorder is a relevant feature that changes the critical behavior of the pure system. According to a criterion due to Harris [1], a spatially quenched disorder will be irrelevant with respect to the critical properties if the inequality $d\nu_{\perp}>2$ is obeyed for the pure system, where $\nu_{\perp}$ is the spatial correlation length exponent and $d$ is the dimension of the system. For models belonging to the directed percolation universality class, such as the contact process [2, 3, 4, 5], this inequality is not fulfilled for $d<4$. We should thus expect a change in the critical properties of the contact process with quenched disorder, as is the case of the quenched diluted contact process, which is the object of our study here. Numerical simulations of the quenched diluted contact process in two dimensional lattices [6, 7, 8, 9, 10, 11] indeed confirm the change in the critical properties.

A remarkable critical behavior of the quenched diluted contact process is the slow activated dynamics, of the logarithmic type, instead of the usual power law type. This result was advanced by Hooyberghs et al. [12, 13] by mapping the evolution operator of the stochastic process describing the quenched diluted contact process into a random quantum spin-1/2 operator, and by the use of a renormalization group approach. This critical behavior places the quenched diluted contact process into the universality class of the random transverse-field Ising model [14, 15, 16, 17, 18, 19, 20, 21]. The slow activated dynamics of the quenched diluted contact process has been confirmed by numerical simulations in two dimensions [8, 9, 10, 11].

Here we study a quenched diluted contact process in which the quenched dilution is obtained by the removal of a fraction of sites of the lattice. The remaining sites form then clusters of site percolation. We aim to study the critical properties of the contact process with quenched dilution by a method in which the percolation clusters are understood as related to the stationary state of stochastic models for the spreading of disease [22, 23, 24]. The use of an epidemic process turns out to be a procedure to create percolation clusters as efficient as the ordinary method of simply creating random vacancies and then using a clustering algorithm to find the percolating cluster. A straightforward numerical approach to the diluted contact process is to consider all the remaining sites of a lattice after a certain fraction of them has been removed [6, 7, 10]. Other methods such as ours consider instead just the sites of the percolating cluster [9, 11]. In this case the total computer time should include the time it takes to generate the percolating cluster. However, this time is very short, representing in our approach less than 1% of the total computer time.

The stochastic model we use to generate clusters of site percolation is defined as follows [23, 24]. Each site of a regular lattice is occupied by an individual that can be in one of three states: susceptible, exposed or immune. A susceptible individual, in the presence of an exposed individual, becomes exposed with certain probability $p$ and immune with the complementary probability. The exposed and immune individuals remain forever in these states. Starting with a single exposed individual in a lattice full of susceptible individuals, a cluster of exposed individual is generated such that at the stationary state it is exactly mapped into a cluster of site percolation [23, 24], with $p$ being identified as the probability of a site occupation.

Once a cluster of site percolation is generated by the model of spreading of disease explained above, we simulate the contact process in the top of the percolating cluster. Only the percolating cluster is needed because a finite cluster cannot sustain an active state. That is, if we wait enough time, the absorbing state will be reached. This procedure is thus interpreted as the contact process with quenched site dilution. More details on the models will be given in the next section. In this same section we set up the evolution equations for one and two-site correlations and solve them by the use a pair mean-field approximation, which allows us to construct the phase diagram. This phase diagram shows that at the percolation threshold the critical creation rate of the quenched diluted contact process is finite.

Using the method presented above, we have obtained the critical properties and the phase diagram of the diluted contact process by numerical simulations and also by a mean-field theory. The method allowed us to obtain more accurate values for the critical exponents and thus confirming the prediction that the quenched diluted contact process belongs to the universality class of the random transverse-field Ising model. The mapping of the epidemic processes into the quenched dilution contact process, allows also to conclude that this universality class may include some models for epidemic spreading.

The contact process, and other models belonging to the universality of directed percolation, describes the transition of an active state to an absorbing state, in which a system cannot never scape. This nonequilibrium phase transition is very common in nature and may occur in various situations [25]. However, the experimental observation of the critical exponents is very difficult, as any amount of disorder should alter the critical behavior, but the exponents were eventually measured in a electrohydrodynamic convection of nematic liquid crystals [26].

2 Models and pair approximation

We begin by defining the two models by using the spreading of disease language. The two models are illustrated in figure 1. The first model (A) is the one that generates the site percolation clusters, and is thus the underlying support over which the second model (B), the contact process, is defined.

2.1 First model

Each site of a regular lattice is occupied by an individual that can be susceptible (S), immune (U) or exposed (E). The possible processes of the first model are as follows:

[TABLE]

[TABLE]

These two processes define a continuous time stochastic process whose probability distribution obeys a master equation. Instead of writing down the master equation, which gives the time evolution of the probability distribution, we write the time evolution of some marginal probability distribution such as the one-site and two-site probability distributions. Using a procedure developed earlier [24, 23] and the notations $P_{X}$, $P_{XY}$, $P_{XYZ}$ for one-site, two-site and three-site probabilities, the following time evolution equations can be derived

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $k$ is the coordination number of the regular lattice and $\mu=(k-1)/k$.

An approximate solution can be obtained by the use of the pair mean-field approach which amounts to use the approximation $P_{XYZ}=P_{XY}P_{YZ}/P_{Y}$. Using the notation $x=P_{S}$, $y=P_{E}$, $v=P_{ES}$, $u=P_{US}$, and $w=P_{EU}$, we may write

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we have taken into account that $P_{SS}=P_{S}-P_{ES}-P_{US}=x-v-u$. Equations (8)-(12) have been solved in reference [23]. At the stationary state, the solution is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $s$ is the root of the polynomial equation

[TABLE]

and $p=b/(a+b)$ and $q=1-p$. The trivial solution is $s=1$, which gives $x=1$, $y=v=u=w=0$ and corresponds to the non spreading regime. The solution $s\neq 1$ corresponds to the spreading regime and occurs only when $p>p_{c}=1/(k-1)$. We remark that the stationary solution is exactly mapped into the site percolation model with $p$ playing the role of the probability of occupancy or the fraction of occupied sites in the percolation model. The spreading regime ($s\neq 1$) corresponds to the existence of the percolating cluster. The non-spreading regime corresponds to the absence of the percolating cluster ($s=1$).

2.2 Second model

As before, an individual can be susceptible (S), immune (U) or exposed (E). In addition, an individual can also be infected (I). Thus in the second model, each site can be in one of the states: S, U, E, and I. However, the sites in states S and U remains forever in these states. The only sites that have their states modified are the E and I sites. They are modified according to the following processes

[TABLE]

[TABLE]

which are the reactions of the contact process. The relation between the infection rate $\lambda$, often used in studies of the contact process, and $c$ and $r$ is given by $\lambda=c/r$. For convenience, we also make use of a parameter $\alpha$, defined by $\alpha=r/c=\lambda^{-1}$.

The initial state of the second model is chosen to be the stationary state of the first model. However, this state has no site in state I and the dynamics does not start. To start the dynamics we choose randomly one site in state E and replace it by a state I. By this procedure, a cluster of sites in state I growths in the top of the percolation cluster. It should be understood that the percolation cluster is formed by sites of type E and I. The sites of type U are at the border of the percolation cluster. The rest of the sites are in state S, and they are not connected to the sites of the percolation cluster.

The two reactions (19) and (20) show that the number of sites E and sites I is a constant implying that the sum $P_{E}+P_{I}$ is a constant. Since these two reactions do not involve the sites U and S, it follows that the number of sites U and the number of sites S are invariants, and the sum $P_{EU}+P_{IU}$ is a constant.

Again, using the procedure developed earlier [24, 23], the following time evolution equations for the one-site and two-site probabilities can be obtained

[TABLE]

[TABLE]

[TABLE]

where $\mu=(k-1)/k$. Due to the constraints stated above it is not necessary to write down the time evolution equations for the other one-site and two-site probabilities.

Using again the pair approximation and the previous notation together with the notations $z=P_{I}$, $g=P_{IE}$, $h=P_{IU}$, we may write

[TABLE]

[TABLE]

[TABLE]

where we have taken into account that $P_{II}=P_{I}-P_{IE}-P_{IU}=z-g-h$ and $P_{EE}=P_{E}-P_{EI}-P_{EU}=y-g-w$.

Equations (24)-(26) are to be solved using as initial conditions the stationary state of the first model. Since $P_{E}+P_{I}=y+z$ is invariant, it follows that $y+z=y_{0}$ where $y_{0}$ is the value of $P_{E}$ at the stationary state of the first model, given by equation (14). Analogously, $P_{EU}+P_{IU}=w+h$ is invariant implying $w+h=w_{0}$ where $w_{0}$ is the value of $P_{EU}$ at the stationary state of the first model, given by equation (17). At stationary state, equations (24)-(26) have a trivial solution $z=0$, characterizing the absorbing state, and a nontrivial solution for which $z\neq 0$, characterizing the active state. Solving for $z$, it is possible to obtain an expression for the nontrivial solution $z$. By taking the limit $z\to 0$ of the nontrivial solution we get the critical line, which is given by

[TABLE]

where $s$ is the root of the polynomial equation given by equation (18). The critical line $\alpha$ versus $q$, shown in figure 2, separates the active percolating phase from the inactive percolation phase. Notice that, when $p\to p_{c}=1/(k-1)$ we get $\alpha=2(k-1)/k^{2}=\alpha_{0}$, so that the critical line meet the vertical line $p=p_{c}$ at $\alpha=\alpha_{0}$, as shown in figure 2. It straightforward to show that the critical exponent related to the order parameter is the same as that of the pure system.

3 Scaling relations and numerical simulations

Around the critical point, the quantities that characterize the critical behavior are assumed to obey scaling relations. In the present case of the diluted contact model, for which the quenched disorder is relevant, the usual scaling relation in terms of power laws in time is replaced by power laws in the logarithm of time, called activated scaling [11, 12]. At the critical point, the space correlation length $\xi$ behaves as [12]

[TABLE]

where $\psi$ is the tunneling critical exponent [15] and $t_{0}$ is a constant. Other quantities behave similarly at the critical point, such as $N_{I}$, the number of infected individuals,

[TABLE]

and $P$, the survival probability at time $t$,

[TABLE]

From the scaling relations (28), (29), and (30) we find that

[TABLE]

[TABLE]

[TABLE]

valid at the critical point. These are useful relations because they do not depend on $t_{0}$.

At the stationary state ($t\to\infty$), the quantities that describe the critical behavior follow the usual power laws, but with exponents distinct from those of the pure system. The order parameter $\rho$, defined as the fraction of infected sites in the percolating cluster, behaves as

[TABLE]

Initially, we have simulated the first model to generate a percolating cluster. The simulation, with periodic boundary conditions, was performed as follows. At each time step, we choose at random a bond from a list of active bonds. An active bond is a pair of SE nearest neighbor sites. The site S of the chosen bond becomes E with probability $p$ and becomes U with the complementary probability $q=1-p$. The chosen bond is removed from the list and the list us updated. The time is then incremented by an amount $1/N_{a}$, where $N_{a}$ is the number of active bonds in the list. Notice that, if a site S has $n_{E}$ nearest neighbor sites in states E, then it will appear $n_{E}$ times in the list. Starting with just one E site in a lattice full of S sites, this algorithm will generate a cluster of E sites. This process stops when there is no SE bonds in the lattice. When this happens the cluster of E sites is a site percolating cluster with U sites standing in the border of the cluster, separating the E sites from the S sites.

Having generated a percolating cluster of E sites, we simulate the contact model on top of the cluster using the following algorithm. At each time step, a site is chosen at random among a list of I sites. With probability $p_{a}=\lambda/(\lambda+1)$ it becomes an E site and with the complementary probability $1-p_{a}$ a nearest neighbor site is chosen at random. If the chosen neighboring site is in state E, it becomes I, otherwise, nothing happens. The time is then incremented by an amount $1/N_{I}$, where $N_{I}$ is the number of I sites. The initial condition is formed by the cluster of E sites with one E site turned into an I site. This I site is taken as the origin.

We have performed simulations on a square lattice with $N=L^{2}$ sites with $L$ up to $L=8192$. For several values $p$ and $\lambda$, we have measured, as a function of time, the number of infected sites $N_{I}$, the survival probability $P$ and the correlation length $\xi$ defined by

[TABLE]

where the summation is over the sites occupied by an infected and $r_{i}$ is the distance from the site $i$ to the origin. Each quantity was measured by averaging over $10^{5}$ to $10^{6}$ disorder configurations, where each disorder configuration is obtained by the simulation of the first model starting with a distinct seed of random number. We have also performed simulation with smaller values of $L$. However, results coming from lattice with $L=4096$ agree, within statistical errors and up to the maximum time we have used, with those coming from $L=8192$. The statistical errors were determined by the calculation of the standard statistical deviation. The results for the three quantities $N_{I}$, $P$ and $\xi$, determined for several values of $\lambda$, are shown in figures 3, 4, and 5. The error bars in these figures are not shown, but they are less than 8%. At the critical point, $\lambda=2.1075$, they are even less reaching 1%.

Figure 3 shows the plot of the number of infected $N_{I}$ as a function of time $t$ for $p=0.8$. Fitting the expression (29) to the data points of figure 3 we estimate the critical parameter as being $\lambda=2.1075(1)$, the critical exponent as $\theta=0.13(2)$, and $\ln t_{0}=6.0(5)$. To find the exponents $\psi$ and $\delta$ and a better estimate of $\theta$ we use a procedure similar to the one used in [11] in which we first determine the quantities $\theta/\delta$, $\theta\,\psi$ and $\delta\psi$, by fitting the expressions (31), (32) and (33) to the data points. After that, we use expressions (28), (29) and (30) to find the exponents $\theta$, $\delta$ and $\psi$ by a constrained fitting, to be explained below.

From the plots of $N_{I}$ versus $P$, shown in figure 4, $N_{I}$ versus $\xi$, shown in figure 5, and $P$ versus $\xi$, we may get, respectively, $\theta/\delta$, $\theta\,\psi$ and $\delta\psi$. Since the scaling relations (31) (32) and (33) do not involve time, the estimates of these quantities are independent of the $t_{0}$, resulting in more precise values, which are found to be

[TABLE]

[TABLE]

[TABLE]

The consistency of these values can be checked by dividing equations (37) and (38). The result is $\theta\,\psi/\delta\psi=0.0758(57)$, which is in fair agreement with (36). The value of $\theta/\delta$ can be used to get the ratio $\beta/\nu_{\perp}$ between the order parameter critical exponents $\beta$ and the critical exponents $\nu_{\perp}$ related to the spatial correlation length. Using the relation $\theta/\delta=d\nu_{\perp}/\beta-2$ we find

[TABLE]

The exponents $\theta$, $\delta$ and $\psi$ are found by a procedure as follows. For each value of $t_{0}$ in the interval $5.5\leq t_{0}\leq 6.5$, we determine the exponents $\theta$, $\delta$ and $\psi$ by fitting the expressions (28), (29) and (30) to the data points. After that we choose the actual values of these exponents as the ones such that the quantities $\theta/\delta$, $\theta\,\psi$ and $\delta\psi$ are as close as possible to the values given by (36), (37), (38). This procedure leads to the following values for the exponents:

[TABLE]

[TABLE]

[TABLE]

We have also performed simulations to get the stationary properties by using systems of linear size $L=2048$. The interest quantities were obtained by the use of $10^{8}$ Monte Carlo steps after discarding $10^{7}$ Monte Carlo steps. In this case, we use as the initial state a configuration in which a fraction of sites is in the infected state. Again we determined the number of infected sites $N_{I}$ at the stationary state from which we obtained the density $\rho=N_{I}/N_{C}$ where $N_{C}$ is the number of sites of the cluster, that is, the number of I sites plus the number of E sites. Assuming the critical behavior (34), we get the value $\beta=1.11(6)$ by plotting $\rho$ as a function of $\lambda-\lambda_{c}$, as shown in figure 6 for the case of $p=0.8$. This result for $\beta$ together with the numerical value for the ratio $\beta/\nu_{\perp}$, obtained above, gives us $\nu_{\perp}=1.15(6)$.

We have also performed similar simulations for other value of $p$ and obtained the critical line, shown in figure 7. In particular, at the percolation critical point $p=p_{c}=0.59274$, we get $\lambda=3.10(1)$.

The critical exponents obtained here are shown in table 1 together with results coming from other papers on the quenched diluted contact process [11, 10] and on the random transverse-field Ising model [16, 20]. To make contact with exponents used to describe the critical behavior of the random transverse-field Ising model, we have determined from our results the fractal dimension critical exponent $d_{F}=d-\beta/\nu_{\perp}$ and the exponent $\phi$, related to the fractal dimension and the tunneling exponent by $d_{F}=\phi\psi$ [16]. We see that our results agree, within the statistical errors, to all other results cited in table 1. The results are $d_{F}=1.036(2)$ and $\phi=1.87(10)$.

4 Conclusion

We have studied the critical properties of the quenched diluted contact process through a mean-field theory and Monte Carlo simulations by using a two stage procedure. The first was the generation of the percolating cluster, obtained by the use of a stochastic lattice model whose stationary states are the clusters of percolation model. The second stage was the simulation of the contact process on the top of the percolating cluster. It should be remarked that, only the percolating clusters is necessary if we wish to study the static stationary properties because finite clusters cannot support an active state. For a finite cluster, the absorbing state will be reached if we wait enough time. As to the dynamic properties, our results show that they can also be obtained from the percolating cluster, or at least their critical properties, as can be inferred by comparing our critical exponents with other works. The present method allowed to obtain more precise critical exponents, with errors that are at most equal to 6%, confirming the prediction that the quenched diluted contact process belongs to the universality class of the random transverse-field Ising model.

The mapping of the two epidemic models into the quenched diluted contact process allows to speculate about the existence of epidemic models that are in the universality class of transverse-field Ising model. In fact, this is the case of the model illustrated in figure 8, which may be thought as a merger of the two models in figure 1. The epidemic model of figure 8 is defined on a lattice in which each site can be in one of four states: S, U, E, and I, and is composed by three catalytic reactions: S $\to$ U, S $\to$ E, E $\to$ I, and by a spontaneous reaction I $\to$ E. At the stationary states, the I and E sites form a connected cluster of sites consisting of a site percolation cluster. The E and I sites evolves then as the contact process on the top of a percolating cluster. Therefore, the model defined by rules of figure 8 is also mapped into the quenched diluted contact process and its critical properties puts the model in the universality class of the transverse-field Ising model.

Acknowledgment

We wish to acknowledge the Brazilian agency FAPESP for financial support.