Critical decay exponent of the pair contact process with diffusion
Abstract
We investigate the onedimensional pair contact process with diffusion (PCPD) by extensive Monte Carlo simulations, mainly focusing on the critical density decay exponent . To obtain an accurate estimate of , we first find the strength of corrections to scaling using the recently introduced method [S.C. Park. J. Korean Phys. Soc. 62, 469 (2013)]. For small diffusion rate (), the leading correctionstoscaling term is found to be , whereas for large diffusion rate () it is found to be . After finding the strength of corrections to scaling, effective exponents are systematically analyzed to conclude that the value of critical decay exponent is irrespective of . This value should be compared with the critical decay exponent of the directed percolation, 0.1595. In addition, we study two types of crossover. At , the phase boundary is discontinuous and the crossover from the pair contact process to the PCPD is found to be described by the crossover exponent . We claim that the discontinuity of the phase boundary cannot be consistent with the theoretical argument supporting the hypothesis that the PCPD should belong to the DP. At , the crossover from the mean field PCPD to the PCPD is described by which is argued to be exact.
pacs:
05.70.Ln, 05.70.Jk, 64.60.HtUTF8mj
I Introduction
The pair contact process with diffusion (PCPD) is an interacting particle system with diffusion, pair annihilation (), and branching by pairs (). The PCPD was introduced in 1982 by Grassberger Grassberger (1982), but had remained unnoticed in the statistical physics community for about 15 years. Since Howard and Täuber Howard and Täuber (1997) (re)introduced the socalled ‘bosonic’ PCPD in 1997, the PCPD has been captivating statistical physicists and the effort to understand the critical behavior of the PCPD has continued until now Ódor (2000); Carlon et al. (2001); Hinrichsen (2001); Park and Kim (2002); Barkema and Carlon (2003); Kockelkoren and Chaté (2003); Janssen et al. (2004); Noh and Park (2004); Park and Park (2005a, b); de Oliveira and Dickman (2006); Park and Park (2006); Kwon and Kim (2007); Smallenburg and Barkema (2008); Park and Park (2009); Schram and Barkema (2012); Gredat et al. (2014).
It was accepted almost without question that the upper critical dimension of the PCPD is 2 and the PCPD does not belong to the directed percolation (DP) universality class in higher dimensions Ódor et al. (2002). On the other hand, the one dimensional PCPD has gained its notorious fame from the beginning because of its strong corrections to scaling. Consequently, numerical studies have reported scattered values of critical exponents (see Table I of Ref. Smallenburg and Barkema (2008) for a summary of reported values of critical exponents) and, in turn, as many hypotheses concerning the critical behavior were suggested as the number of research groups involved (for a review of the early various scenarios, see Ref. Henkel and Hinrichsen (2004)). It still remains an open question whether the one dimensional PCPD belongs to the DP class (DP hypothesis) or forms a different universality class just like the higher dimensional PCPD (nonDP hypothesis). Since we are mainly interested in the PCPD in one dimension, by the PCPD in the following we will exclusively mean the onedimensional PCPD, unless dimensions are explicitly stated.
To support the DP hypothesis, Hinrichsen Hinrichsen (2006) provided a theoretical argument as to why the PCPD should belong to the DP class. This argument begins with the numerical observation that the dynamic exponent of the PCPD is smaller than 2, which means critical clusters spread superdiffusively. Since isolated particles can spread at most diffusively with dynamic exponent 2, diffusive motion of isolated particles can be regarded as frozen in comparison to the critical spreading and, in turn, the effectively frozen isolated particles can at best play the role of isolated particles of the pair contact process (PCP) without diffusion which is known to belong to the DP class. Accordingly, the PCPD should belong to the DP class. Although this argument is not unquestionable (see Ref. Park and Park (2008) for a critique), it is quite persuasive once the dynamic exponent of the PCPD is accepted to be smaller than 2 as numerical studies suggest. We will discuss this argument at the end of Sec. IV.1. Barkema and his colleagues have reported numerical results to support the DP hypothesis Barkema and Carlon (2003); Smallenburg and Barkema (2008); Schram and Barkema (2012).
At the same time, numerical studies supporting the nonDP hypothesis are also available. The critical behavior of the driven pair contact process with diffusion Park and Park (2005a, b) seems to suggest that the PCPD cannot be described by a field theory with a single component field, which makes the PCPD not satisfy the prerequisite of the DP conjecture Janssen (1981); Grassberger (1982). Besides, the existence of nontrivial crossover from the PCPD to the DP Park and Park (2006, 2009) was invoked to support the nonDP nature of the PCPD.
Since the controversy arises from numerical difficulty of finding accurate value of critical exponents due to strong corrections to scaling, it is necessary to tame the corrections to scaling at our disposal. To this end, this paper exploits a systematic method suggested in Ref. Park (2013) to find corrections to scaling without prior knowledge of leading asymptotic behavior. Once the strength of corrections to scaling is determined, the effective exponent can be systematically analyzed to get an accurate estimate of the critical exponent. In this paper, we find the critical decay exponent by extensive Monte Carlo simulations, using the method briefly mentioned above.
The paper is organized as follows: Section II consists of two parts. To be selfcontained, Sec. II.1 introduces the dynamic rules of the PCPD and describes the expected behavior of the order parameter in each phase. Also an algorithm to simulate the stochastic dynamics is detailed with comparison to previous studies. In Sec. II.2, a method to estimate the leading correctionstoscaling term is explained. Numerical estimate of the critical decay exponent is presented in Sec. III, using the method explained in Sec. II.2. Section IV studies crossover behavior from two extreme points of the model, (without diffusion) and with a suitable time rescale (mean field), to the PCPD. Section V summarizes the work.
Ii Model and Method
ii.1 Model
The PCPD is defined on a onedimensional lattice of size with periodic boundary conditions. Each site is either occupied by a particle () or empty () and every site can accommodate at most one particle. The dynamics of the PCPD are defined by the following transition events,
(1) 
where and . For bookkeeping purposes, we introduce a stochastic process at every site which takes 1 (0) if site is occupied (vacant) at time . We define the particle density, , the pair density, , and the triplet density, , as
(2) 
where means the average over ensembles. The evolution equation of is
(3) 
When (critical point), both and are expected to show asymptotic powerlaw behavior with corrections to scaling as
(4)  
(5) 
where , , , are constants and and contain higher order terms which decay faster than and , respectively. In one dimension, it is believed that equals , whereas the mean field theory assumes ; see Sec. IV. Also it is believed that , numerical evidence of which will be provided in Sec. III. On this account, we will drop the primes in the symbols of exponents for in what follows and we will refer to and as the critical decay exponent and the leading correctionstoscaling exponent (LCSE), respectively. In the active phase ( within our model definition), both and approach certain nonzero values exponentially as and in the absorbing phase (), and decrease to zero faster than for nonzero .
There are two important limiting cases. When , this model corresponds to the PCP Jensen (1993) which has infinitely many absorbing states and belongs to the DP class. Meanwhile, taking limit with kept finite, the (site) mean field theory becomes exact. Hence there are two crossover behaviors at (from the PCP to the PCPD) and at (from the mean field PCPD to the PCPD). In Sec. IV, we will study these two kinds of crossover behavior.
To simulate the model, we employ the following algorithm. At first, we introduce
(6) 
which makes and interpreted as probability. Now assume that there are particles at time . We randomly choose one particle among particles with equal probability and choose one of two nearest neighbors of the chosen particle with equal probability. Let us assume that the chosen particle is located at site and the selected neighbor site is (). If , the particle at site moves to site with probability , but with probability , nothing happens. In the case , two particles at sites and are removed with probability (pair annihilation), the site becomes occupied with probability (branching), or with probability nothing happens. If is already 1 in the branching attempt, no change in the configuration happens. After the above attempt, time increases by .
Notice that the PCPD studied in Refs. Park and Park (2005a, 2006) corresponds to the case with up to a timerescale factor (time of the PCPD in Refs. Park and Park (2005a, 2006) corresponds to of the case with in this paper). Also note that the simulation algorithm employed in Ref. Ódor (2003) is slightly different from that used in this paper. But if we set where is the diffusion parameter used for simulations in Ref. Ódor (2003) and if we rescale time appropriately, the simulation results in Ref. Ódor (2003) can be directly compared to those obtained by the algorithm in the above. For example, the reported critical point of the case with in Ref. Ódor (2003) is consistent with that in Ref. Park and Park (2006) which is .
ii.2 Corrections to scaling
A systematic way to find the critical decay exponent simultaneously together with the critical point is to investigate the behavior of the effective exponent defined as ()
(7) 
which, by Eq. (4), is expected to behave at the critical point as
(8) 
In the time regime where higher order terms are negligibly small but the leading correction term is not negligible, a plot of the effective exponent against with the correct value of should be a straight line if the system is at the critical point. Meanwhile, it veers up (down) if the system is in the active (absorbing) phase as gets smaller. From the expected behavior in each phase, the critical exponent and the critical point can be found simultaneously by investigating the behavior of effective exponents near criticality, once the LCSE is known. Hence the information about the LCSE is indispensable in order to estimate the critical decay exponent and the critical point accurately via the effective exponents.
A systematic method to find corrections to scaling without knowing was recently suggested Park (2013). The idea is that at the critical point the double ratio of at three different time points should behave asymptotically as
(9) 
where is a (fixed) constant. Although is not necessarily the same as in Eq. (8), we will use the same value for and in this paper and we will drop the subscript in in the following. We introduce a correctionstoscaling function (CTSF) as the logarithm of the left hand side of Eq. (9),
(10) 
which behaves at criticality as
(11) 
The behavior of at offcriticality is also of interest. If the system is slightly away from the critical point with , is indistinguishable from Eq. (11) up to , where is the critical exponent describing the divergence of correlation time. Since can be understood as a curvature of the curve as a function of , that is, , gets larger (smaller) with if the system is in the active (absorbing) phase. Hence, if the coefficient of in Eq. (11) is positive, a typical behavior of around the critical point looks like Fig. 1. These curves are obtained from simulations of the PCPD with . To be specific, the system size is and the number of independent runs are 9800, 7200, and 1600 for (active), 0.258 112 (critical), and 0.258 114 (absorbing), respectively. Due to the double derivativelike nature of , the curves obtained from numerical simulations can be very noisy as seen in Fig. 1 unless the number of independent simulation runs is very large. In this respect, alone may not be a good measure to find the critical point and the analysis of the effective exponents which are less noisy than should be accompanied.
Although the behavior of looks qualitatively similar to around the critical point, in the active phase actually should approach zero as because in this limit saturates to a finite number with zero curvature. We will soon see such a long time behavior in the active phase from an exactly solvable model. In the absorbing phase, the behavior of in the limit of infinite time depends on the asymptotic behavior of . If decreases exponentially, so does . On the other hand, if decreases as a powerlaw in the absorbing phase like the PCPD, should also approach zero as . Even though of the PCPD should approach zero in all phases as , the deviation of at some point from the critical is conspicuous as Fig. 1 illustrates. Thus, such infinite time behavior does not limit the practical usefulness.
The sign of is not necessarily positive and an example of the case with a negative can be found by the following equation
(12) 
with initial condition . The solution is
(13) 
Since the coefficient of the leading correctionstoscaling term at the critical point () is negative, it is appropriate to draw vs on a double logarithmic scale as in Fig. 2 which makes the curve in the active (absorbing) phase veer down (up) contrary to Fig. 1. The inset of Fig. 2 shows that in the active phase approach 0 as as argued before. Since decreases exponentially in the absorbing phase of this example, also decreases exponentially to .
Iii Critical Decay Exponent
This section investigates the critical decay exponent of the PCPD for various ’s via the analysis of the effective exponents along with the corresponding CTSFs. For all numerical analyses in this section, the system size is and the initial density is 1. Since and are equally important, the CTSFs for both and are studied and will be denoted by and , respectively.
At first, we will present the results for the case with . As we will show later, the critical point is found to be , where the number in parentheses indicates the uncertainty of the last digit. Figure 3 shows doublelogarithmic plots of and vs at for . These data are obtained from 3000 independent runs at the designated value of . It seems that shows a powerlaw behavior as while has no symptom of powerlaw decay up to . On the other hand, and become almost overlapped after , which implies that not only the LCSE but also the coefficients of the leading correctiontoscaling terms of and are identical. Since exhibits a more or less clean powerlaw behavior and eventually follows , we estimate to be 0.138 from . Note that this estimate is comparable to that in Ref. Schram and Barkema (2012).
It is worthwhile to discuss the implication of the difference between the two CTSFs, , which shows a clean powerlaw behavior with albeit . By definition of the CTSF, can be regarded as the CTSF for . Hence if one analyzes rather than and separately with the assumption of in Eqs. (4) and (5), one may wrongfully conclude that the corrections to scaling of the PCPD are weaker than the actual strength . In fact, the cancellation of the leading correctionstoscaling term in was already anticipated in Ref. Schram and Barkema (2012) and we confirmed it through the direct numerical analysis.
Actually, the cancellation of the leading correctionstoscaling term in the ratio of two quantities is not unusual. An immediate example arises when we analyze Eq. (3). Since , and at the critical point should be
(14) 
where and are constants satisfying , contains all terms decaying faster than but slower than , and should be strictly smaller than . Otherwise, the leading power of the left hand side of Eq. (3) cannot be the same as that of the right hand side of the equation. Hence, at the critical point should have the form
(15) 
Unlike the cancellation of the leading correctionstoscaling terms in , however, we could not find any theoretical reason why and should have exactly the same leading correctionstoscaling term. This can be a theoretical challenge of the PCPD.
To check the consistency of Eq. (11), we analyze ’s for various values of ( and ). We first fit for using a fitting function with two fitting parameters and . From the fitting, we estimate and the coefficient of the leading correctionstoscaling term to be . The straight line lying on for (top curve) in Fig. 4 is the result of this fitting. Then, we compare for and with the estimated values of and to ’s for corresponding ’s, which shows an excellent coincidence.
We think this coincidence provides a numerical evidence for the absence of logarithmic corrections in the leading correctionstoscaling term. When we derive Eq. (11), we tacitly assumed that the leading term has no logarithmic corrections. If there happens to be such corrections, the above procedure should exhibit a systematic deviation for different values of . Also a nice powerlaw behavior of for about four decades, as seen in Fig. 4 suggests that logarithmic corrections, even if exist, are negligible. Furthermore, the clean powerlaw behavior of also provides an indirect evidence against logarithmic corrections.
Finding the LCSE to be 0.138, we now analyze the effective exponent . Since shows a cleaner powerlaw behavior than , we analyze calculated from . In Fig. 5, obtained using is drawn as a function of with near criticality. From this plot, it is clear that the system with (0.111 159) is in the active (absorbing) phase and the critical point should be . Although the number of independent runs for both offcritical simulations is only 200, the effective exponents give a clear illustration of the expected behavior in each phase. We fit for using a linear function to obtain that . Since a fitting result of varies from 0.17 (for ) to 0.176 (for ) with , we conclude that .
Figure 5 shows that at criticality becomes an almost straight line in the region . This behavior is indeed consistent with the analysis of . We observed in Fig. 4 that for exhibits a nice powerlaw behavior from . Recall that is calculated using at , , and . Thus, this numerical observation implies that the leading correctionstoscaling term becomes dominant from . Thus for should be almost straight from , as seen in Fig. 5.
Since the estimated critical exponent 0.173 is close to the DP exponent 0.1595, it would be an interesting practice which value of can predict the DP critical exponent. By trial and error, we found that gives the DP critical exponent; see the inset of Fig. 5. Although 0.105 is different from the estimated 0.138 by 25%, it is indeed hard to exclude the possibility of the DP critical scaling. Thus, the analysis of the PCPD with alone may not be enough to conclude that the PCPD does not belong to the DP class. To make the estimate 0.173 more convincing, we analyze other cases with different values of .
Before delving into the case with different diffusion probability, we will show that the estimated critical point is rather insensitive to the estimates of and . If is small as in the present case, we can approximate as long as is not so large. Thus, the effective exponent at criticality should be insensitive to . As can be seen in Fig. 6, at is more or less insensitive to , which strongly supports that the density at exhibits the critical scaling up to the simulated time. Meanwhile, if the system is in the active (absorbing) phase and , at given should increase (decrease) significantly as gets smaller. The inset of Fig. 6 depicts at and for . By comparing this figure with with Fig. 5, it is easily recognized that at offcriticality is significantly affected by the change of . Although we plotted vs in Fig 6 for convenience, different choice of does not affect the observed behavior of the effective exponents under the change of . Also note that does not play any role in the above discussion. Hence, we conclude that the estimated critical point is accurate irrespective of whether we are using the right value of and .
Now, we will analyze the case of . Figure 7(a) depicts the CTSFs as functions of on a double logarithmic scale at which will turn out to be the critical point. The data are collected from 7000 independent simulation runs and is used. Unlike the case with , shows a powerlaw decay from about . in the short time regime decays faster than but after , and are almost indistinguishable. Since shows a more stable powerlaw behavior than , we estimate to be 0.15 which is comparable to the above estimate. One can also see that decays faster than .
In Fig. 7(b), corresponding to for is drawn against . The effective exponent for (0.152 476) results from 2400 (2500) independent simulation runs. This figure shows that the critical point is and a linear fit of for gives which is consistent with the estimate for the case of . Also note that is almost straight in the region , which is consistent with the behavior of .
We would like to emphasize that the estimate of the critical point is rather insensitive to the accuracy of and , so the accuracy of is less questionable than the exponents. Using the data of the density at the critical point, we also estimate, by trial and error, the value of which gives the DP exponent. At this time, the desired value of becomes which is quite different from the estimated . Furthermore, 0.077 differs by 25% from the case of . That is, for our numerical data to be consistent with the DP hypothesis, the LCSE has to vary continuously with significantly.
Since the case of was also studied in Ref. Schram and Barkema (2012) which supports the DP hypothesis, it is worth while to compare our results with those in Ref. Schram and Barkema (2012). First, the critical point in this paper is more accurately estimated than in Ref. Schram and Barkema (2012); see Table I of Ref. Schram and Barkema (2012). Second, a value close to the DP exponent was obtained from the system at which is actually in the active phase. Because the density in the active phase decays slower than at the critical point, it is not surprising that the estimated value of in Ref. Schram and Barkema (2012) is smaller than 0.173. Interestingly, however, a similar estimate of was attained when the system at was analyzed; see the sixth row of Table 1 in Ref. Schram and Barkema (2012).
The results up to now seem to suggest that the LCSE obtained from the behavior of CTSFs is about 0.15 for the PCPD in general. Quite intriguingly, however, the PCPD with has relatively weak corrections to scaling. In Fig. 8(a), the CTSFs for at the critical point of the case with are depicted as functions of on a doublelogarithmic scale. Note that in this figure is the middle curve in Fig. 1. Unlike the previous two cases, decays as rather than . Note that the powerlaw behavior of is observed from which means the LCSE in becomes dominant from .
Using this LCSE, we depict calculated from with as a function of in Fig. 8(b) which again shows a typical behavior of near criticality. From the behavior of we estimate . A linear fit of for suggests , which is again consistent with the estimates from the cases of different ’s. Since the LCSE is dominant from , the effective exponent should be a straight line for , as seen in Fig. 8 (b).
We also investigated which value of can give the DP exponent for the high diffusion case. Unlike the low diffusion cases, however, the DP exponent was hardly observed by varying , which seems to imply that the critical behavior of the PCPD with cannot be consistent with the DP hypothesis.
Since needs not be universal, appearing dependence of on is not contradictory to our common sense formed by the renormalization group (RG) theory. Still, what kind of mathematical structure is behind the change of with is an interesting question.
We think there are two possible scenarios. The first obvious scenario is that is a continuous function of . Without a fixed line in the RG sense, however, it is hard to expect such a continuously varying exponent albeit nonuniversal, so this scenario does not look plausible. The second one is that the correction term is actually present even in the case with but the coefficient is very small.
The implication of the second scenario can be clearly stated by the following toy equation
(16) 
In this toy example, the leading correctionstoscaling term becomes dominant only when and before this time plays the role of the leading correctionstoscaling term. Thus, is well described by and the corresponding effective exponent drawn as a function of rather than the true asymptotic behavior gives the correct value of 0.15; see Fig. 9. We think the second scenario is plausible, but more investigation is necessary to fully understand how the mathematical structure of the corrections to scaling changes with . This also can be a theoretical challenge of the PCPD. In any case, is a useful tool to find the correct value of the critical exponent, although it may not predict the true LCSE as in the toy example.
To conclude this section, we found that the critical decay exponent of the PCPD is robust against with value . Although this value differs from the DP value only by 8%, the consistent estimate for a wide range of supports the nonDP hypothesis. In particular, the DP hypothesis is not consistent with the case with up to the simulation time in this work (). Since the correctionstoscaling for the high diffusion case is much weaker than those for the low diffusion case, it seems promising to find other exponents accurately by investigating the PCPD for large .
0  0.077 0905(5) 

0.001  0.1019(1) 
0.005  0.1023(1) 
0.01  0.1028(1) 
0.02  0.1038(1) 
0.05  0.1066(1) 
0.1  0.111 158(1) 
0.133 519(3) 

0.5  0.152 4755(5) 
0.9  0.2334(1) 
0.95  0.258 112(2) 
0.99  0.2968(1) 
1 
Iv Crossover
To get a hint to the crossover behavior occurring at two limiting cases and , we found critical points for a wide range of , which are summarized in Table 1. Because the estimate of the critical points within error is relatively easy with the present computing power, we just present the resulting critical points without showing the details.
The phase boundary of the PCPD in the  plane shown in Fig. 10 has two salient features at the two end points, and . The phase boundary is discontinuous at and the phase boundary approaches the mean field critical point with infinite slope as . Each singular behavior signifies a crossover; crossover from the PCP to the PCPD at and that from the mean field PCPD to the one dimensional PCPD at . In this section, we will investigate these two kinds of crossover behavior one by one and find the corresponding crossover exponents.
iv.1 From the PCP to the PCPD
The discontinuity of the phase boundary at can be understood as follows. In this discussion, all exponents are of the DP class. Consider the system at with , where is the critical point of the PCP. If , the system is in the absorbing phase and the pair density decays exponentially if time exceeds the relaxation time . On the other hand, the particle density should approach a certain nonzero value . If , it is unlikely for isolated particles to meet each other purely by diffusion when is smaller than . Hence, effectively, the initial PCP dynamics is almost decoupled with diffusion before and only when exceeds pairs formed by diffusion can appear. Since we are considering an infinite system, the pair density is nonzero for any finite although it can be extremely small. As soon as pairs appear by diffusion, the socalled defect dynamics of the PCP begins. Since the probability that two consecutive sites are occupied by diffusion is roughly , the mean distance between two pairs formed by diffusion should be . Let be the probability that the defect dynamics continues until the cluster size becomes . should be the same order as the survival probability that the defect dynamics continues until , which for small becomes . This is because the scaling form of the survival probability is , where is a scaling function with finite value of . Then the mean distance between the starting points of two successfully spreading clusters should be . If or , a merger of two spreading clusters into a single cluster happens frequently and the system should survive indefinitely. Hence, the condition that the system dynamics continue indefinitely by any small but finite is or . Since , there should be a finite range of which satisfies the above criterion. Thus, the phase boundary should be discontinuous at . Using at the critical point of the PCP (see the inset of Fig. 11) and the DP exponents , , the validity of the above criterion gives which is comparable to the numerical result.
Due to the discontinuity, the phase boundary does not give any information about the crossover exponent . Rather, we find by data collapse using the following scaling ansatz,
(17) 
where is the steady state particle density at , is the crossover exponent, , are scaling functions, and and are the critical exponents of the DP class. We observe the best collapse when we use as shown in Fig. 11. Thus, we conclude . Note that this crossover exponent is different from that of the crossover from the DP class with infinitely many absorbing states to the DP class with finite number of absorbing states Park and Park (2007).
The existence of the nontrivial crossover behavior at has nothing to do with the change of universality classes. In fact, this crossover originates from the drastic decrease of the volume of absorbing states in the configuration space Park and Park (2007). However, the discontinuity at in the phase boundary raises a criticism on the Hinrichsen’s argument explained in Sec. I. This criticism starts from numerical observation that diffusion makes the system more active as the system at with nonzero is in the active phase. Now consider the system at and . In this case, clusters spreads even faster than the critical spreading in the long time limit. Repeating Hinrichsen’s argument, one can conclude that diffusion of isolated particles is negligible in the long time limit and, in turn, isolated particles can join a cluster by the spreading of clusters not by their own diffusion, which is the crucial feature of immobile isolated particles in the PCP. If this were the case, the critical point should approach as and the phase boundary should be continuous at just as the inhibitory route in Ref. Park and Park (2007). Obviously, this is contradictory to the numerical result. Also note that the origin of the discontinuity is the active role of isolated particles, which is completely neglected in the Hinrichsen’s argument. In other words, one cannot deduce a right conclusion by simply comparing the speed of spreading clusters with that of diffusion. Hence, we cannot neglect the effect of pure diffusion and there should be a strong correlation between diffusion and the critical cluster spreading, which can mediate the change of the universality class.
iv.2 From the mean field PCPD to the PCPD
The mean field equation for the PCPD is obtained by setting and in Eq. (3), which gives
(18) 
where . The critical point of the mean field theory is at which the density decays as
(19) 
where is the initial density (we will set ). When , is indistinguishable from and clear deviation from Eq. (19) is observable when , or . In other words, the critical density decay is observable when and the offcritical behavior dominates when . Hence, the critical exponent for the mean field theory is 2.
The relation between the mean field equation and the PCPD under limit can be understood as follows (similar argument is also found in Ref. Park (2009)). Under the limit with kept finite, any correlation generated by reaction dynamics will be removed by diffusion immediately. Since the mean field theory assumes no correlation for all time , Eq. (18) becomes exact in this limit. Even for finite , the mean field equation is an accurate approximation if the relaxation time of diffusion of a randomly chosen region is much smaller than the time between two consecutive reaction dynamics in the same region.
To find the criterion for the validity of the mean field theory for small but finite , consider the mean field dynamics at . According to the mean field solution, the mean distance at time between two nearest particles is . Now consider a region of size at time and assume that two consecutive reaction events occur at and in this region. We also assume that during , this region is not correlated with outside of this region. Since the number of particles is finite in this region, should be ; recall that is rescaled time. Since the diffusion constant in rescaled time is , the relaxation time of diffusion of this region is . If , two consecutive reaction events are uncorrelated and the mean field theory becomes accurate. On the other hand, if , two consecutive reaction events become correlated and the mean field theory fail to describe the system correctly. Hence, the crossover from the mean field PCPD to the PCPD occurs when and becomes a proper scaling parameter.
According to the above argument, the particle density in the asymptotic regime should take the scaling form
(20) 
where as above and is a scaling function. We also conjecture the scaling form of the pair density as
(21) 
where is another scaling function. From the mean field theory, and . Hence, if and , plots of against should collapse onto a single curve for sufficiently large . Furthermore, the phase boundary for should approach the mean field critical point as
(22) 
where is the critical point of the PCPD for given . Hence, the crossover exponent is .
To confirm the above argument, we simulated the PCPD for , , and at with the system size . In Fig. 12(a), we depict vs and vs on a doublelogarithmic scale. For comparison, the mean field solution is also depicted. As argued, the regime where mean field theory is accurate becomes larger as gets closer to 1. In Fig. 12(b), one can see a scaling collapse plot of vs as well as vs , which affirms that the crossover exponent is .
V Summary
To sum up, we numerically studied the critical density decay of the pair contact process with diffusion (PCPD) and estimated the critical decay exponent by investigating effective exponents after finding corrections to scaling for various diffusion strength. For small diffusion probability (), we found that the correctionstoscaling term asymptotically behaves as and for large diffusion probability () the correctionstoscaling term decays as which is weaker than that of the case with small . All the same, the analysis of the effective exponents for any with the corresponding correctionstoscaling term showed that the critical decay exponent is . Although this value is quite close to that of the directed percolation (DP) universality class which is 0.1595, the systematic deviation of for the PCPD from the DP value for any suggests that the PCPD does not belong to the DP class and forms an independent universality class.
We also studied the crossover from the pair contact process (PCP) without diffusion to the PCPD which occurs around and from the mean field theory (MFT) to the PCPD which happens around . We found that the crossover at is characterized by the discontinuity of the phase boundary and that the crossover exponent is . We showed that applying the Hinrichsen’s argument Hinrichsen (2006) which supports the DP hypothesis to this crossover leads to a contradictory conclusion to the discontinuity of the phase boundary at . The crossover from the MFT to the PCPD, occurring at , is described by the crossover exponent , which was argued to be exact.
Acknowledgements.
This work was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant No. 20110014680); and by the Catholic University of Korea, Research Fund, 2013. The author acknowledges the hospitality of the Institut für Theoretische Physik, Universität zu Köln, Germany and support under a German Research Foundation (DFG) grant within SFB 680 Molecular Basis of Evolutionary Innovations during the final stage of this work. The computation was supported by Universität zu Köln. The author also would like to thank the APCTP for its hospitality.Footnotes
 From Ref. Park and Park (2007).
 detailed analysis not shown in the paper.
 detailed analysis not shown in the paper.
 detailed analysis not shown in the paper.
 detailed analysis not shown in the paper.
 detailed analysis not shown in the paper.
 From Ref. Park and Park (2006).
 detailed analysis not shown in the paper.
 detailed analysis not shown in the paper.
References
 P. Grassberger, Z. Phys. B 47, 365 (1982).
 M. J. Howard and U. C. Täuber, J. Phys. A 30, 7721 (1997).
 G. Ódor, Phys. Rev. E 62, R3027 (2000).
 E. Carlon, M. Henkel, and U. Schollwöck, Phys. Rev. E 63, 036101 (2001).
 H. Hinrichsen, Phys. Rev. E 63, 036102 (2001).
 K. Park and I.M. Kim, Phys. Rev. E 66, 027106 (2002).
 G. T. Barkema and E. Carlon, Phys. Rev. E 68, 036113 (2003).
 J. Kockelkoren and H. Chaté, Phys. Rev. Lett. 90, 125701 (2003).
 H.K. Janssen, F. van Wijland, O. Deloubriere, and U. C. Täuber, Phys. Rev. E 70, 056114 (2004).
 J. D. Noh and H. Park, Phys. Rev. E 69, 016122 (2004).
 S.C. Park and H. Park, Phys. Rev. Lett. 94, 065701 (2005a).
 S.C. Park and H. Park, Phys. Rev. E 71, 016137 (2005b).
 M. M. de Oliveira and R. Dickman, Phys. Rev. E 74, 011124 (2006).
 S.C. Park and H. Park, Phys. Rev. E 73, 025105 (2006).
 S. Kwon and Y. Kim, Phys. Rev. E 75, 042103 (2007).
 F. Smallenburg and G. T. Barkema, Phys. Rev. E 78, 031129 (2008).
 S.C. Park and H. Park, Phys. Rev. E 79, 051130 (2009).
 R. D. Schram and G. T. Barkema, J. Stat. Mech.:Theory Exp. (2012), P03009.
 D. Gredat, H. Chaté, B. Delamotte, and I. Dornic, Phys. Rev. E 89, 010102 (2014).
 G. Ódor, M. C. Marques, and M. A. Santos, Phys. Rev. E 65, 056113 (2002).
 M. Henkel and H. Hinrichsen, J. Phys. A 37, R117 (2004).
 H. Hinrichsen, Physica A 361, 457 (2006).
 S.C. Park and H. Park, Eur. Phys. J. B 64, 415 (2008).
 H.K. Janssen, Z. Phys. B 42, 151 (1981).
 S.C. Park, J. Korean Phys. Soc. 62, 469 (2013).
 I. Jensen, Phys. Rev. Lett. 70, 1465 (1993).
 G. Ódor, Phys. Rev. E 67, 016111 (2003).
 S.C. Park and H. Park, Phys. Rev. E 76, 051123 (2007).
 S.C. Park, Phys. Rev. E 80, 061103 (2009).