Print version ISSN 0103-9733
Braz. J. Phys. vol.33 no.3 São Paulo Sept. 2003
Corrections to finite size scaling in percolation
P.M.C. de OliveiraI; R.A. NóbregaI; D. StaufferI, II
IInstituto de Física, Universidade Federal Fluminense, av. Litorânea s/n, 24210-340, Niterói, Brasil
IIInstitute for Theoretical Physics, Cologne University, D-50923 Köln, Euroland
A 1 = L-expansion for percolation problems is proposed, where L is the lattice finite length. The square lattice with 27 different sizes L = 18; 22;... 1594 is considered. Certain spanning probabilities were determined by Monte Carlo simulations, as continuous functions of the site occupation probability p. We estimate the critical threshold pc by applying the quoted expansion to these data. Also, the universal spanning probability at pc for an annulus with aspect ratio r = 1=2 is estimated as C = 0.876657(45).
In reference  a square lattice is viewed as a torus, thus without frontiers, where sites are randomly occupied with probability p. Wrapping percolation probabilities can be defined within this geometry, counting configurations which wrap along the horizontal and/or/xor vertical directions. The probability Rh, for instance, counts all configurations wrapping along the horizontal direction, no matter which occurs vertically. As a function of p, Rh corresponds to a plot like figure 1 (to be quoted later) for a finite lattice. In the thermodynamic limit L ® ¥, this plot approaches a step function, Rh = 0 below, and Rh = 1 above the critical threshold pc. The point pL() exemplified on the plot serves as an estimator for pc: by measuring a sequence of such values, for larger and larger sizes L1, L2 ... LN, one can extrapolate this sequence for L ® ¥. Reference  presents the figure pc = 0.59274621(13), the most accurate available today, obtained from more than 7 × 109 Monte Carlo samples. Starting from an empty lattice, filled site by site at random, each sample provides a single number to the statistics, namely the precise number n of occupied sites for which the proper wrap (horizontal for Rh) appears for the first time. This computational strategy of filling up the lattice site by site and storing data onto bell-shaped n-histograms is equivalent to that of , and similar to early works . However, the multi-step strategy introduced in  involves many other components (see ). The authors of  assume a pL() pc ~ L(2+1/n) dependence, where n = 4/3 is the correlation length critical exponent. As 2 + 1/n = 2.75 is a sufficiently high value, data up to only L = 128 were needed. Also, the authors adopt the value * = 0.521058290, corresponding to that particular wrapping geometry and exactly known  as the limiting value for Rh at pc, called Pinson's number in . Outside this value or within other geometries, the not-so-high exponent 1 + 1/n or even the smaller 1/n were observed .
Here, we propose the mathematical form
for estimators pL obtained from quantities like Rh, where the cutoff M is conveniently chosen according to the numerical accuracy available.
We apply this formula to the L × L square lattice, with a set of different increasing lengths such that the numbers of sites grow by a factor of . Table I shows an example for = 0.9, for which the traditional Chi-square fitting  gives pc = 0.59274675(88), in agreement with  although within a larger error bar. We adopted M = 4 in equation (1), compatible with our smallest lattice size L = 18, since 184.75 = 1 × 106 still falls inside our numerical accuracy, whereas the next term 185.75 = 6 × 10-8 would be outside. The quality of this fit can be appreciated by the so-called goodness-of-fit Q , a quantity between 0 and 1. The fit is considered believable  for values of Q > 0.1. In our example, we get Q = 0.857 for table I. All our data to be discussed hereafter, for many other values of between 0.5 and 0.99, present the same degree of accuracy, giving credit to our proposal, equation (1). In spite of these accurate results, one cannot rule out some possible higher terms deviating from (1), for example that proposed in .
II Measured quantity
First, let's explain the spanning probability we adopted within the torus, instead of the wrapping probability . We consider two parallel lines distant L/2 from each other, on the L × L square lattice. For each sample - again obtained by filling up the initially empty lattice, site by site at random - we count the precise number n of occupied sites for which these lines become connected for the first time, no matter which occurs around the other direction. This approach has a big advantage over the wrapping probability around the whole torus: From the same sample we can count n just L times instead of only once! The parallel lines can be numbered (i, i + L/2) for i = 1, 2, ...L/2 along the horizontal direction, with the same procedure repeated vertically. Thus, the statistics is multiplied by a factor of L. In table I, for instance, the sampling counting 109 for L = 18 corresponds to an n-histogram with 18 ×109 accumulated units (the total area below the bell-shaped curve). In the same table, for L = 1594, the much smaller sampling counting 4 ×106 corresponds indeed to almost the same statistics, i.e. 6 ×109 accumulated units below the curve. This trick allows us to test a wide range of lattice sizes, and verify the validity of our proposal (1). The further computational time one needs in order to implement this trick is negligible: we simply keep in memory the top and bottom (left and right) extreme lines for each already formed cluster of neighbouring occupied sites. Thus, for each new included site, only the last updated cluster must be verified.
The spanning probability function is obtained by superimposing a lot of step functions, one for each counted n, and dividing the result by the total number (L × sampling counts). An average is then performed, weighted by C(L2, n) pn(1 p) where C is a combinatorial factor, yielding the p-continuous curves shown in figure 1.
The 3-digits error bars shown in table I were obtained by dividing the whole set of data into, say, S = 10 sub-sets, independently calculating pL() for each sub-set. The error bars are the standard deviation of this distribution divided by . This last division is based on the supposition that the whole data is normally distributed. In order to verify the validity of this approach, we repeated the same procedure with S = 20. Indeed, the error bars are approximately the same, independent of S. For safety, we adopted always the largest between both error bars so obtained. For intermediate lattice sizes (from L = 74 up to 402), we used S = 10 and 5, instead of 20 and 10.
We also simulated larger L × L lattices for = 0.5, up to L = 24000, with free instead of periodic boundary conditions, within a poorer statistics. The results are also compatible with equation (1).
III Data Analysis
By fitting data with equation (1), we have a set of M + 2 parameters to be determined, namely pc, A0(), A1(), ... AM(). The error bars for these quantities come from the Chi-square fitting procedure , as a consequence of the primary error bars directly measured for the crude data (that of table I, for instance). Thus, the accuracy obtained for each parameter, in particular the interesting one pc, involves a series of accumulated errors. Instead of taking care of pc, let's turn our attention to A0() for a while. Fig. 2 shows a plot of this value as a function of .
One can see that A0 vanishes for the particular value * = 0.984786(11). For the geometry we adopted, this is the equivalent of the above quoted Pinson's number [5,6], i.e. the universal value of our spanning probability at pc in the thermodynamic limit. Recently, Cardy  tried to calculate this value for an annulus, i.e. an L1 × L2 square lattice where periodic boundary condition is adopted only along the direction of L1. The quoted universal value corresponds to the critical spanning probability between the two free-boundary lines, and depends only on the ratio r = L2/L1. Unfortunately, within his theoretical approach, Cardy was forced to leave out spanning configurations which also wrap along the direction of L1, the periodic boundary.
Our geometry with the complete L × L torus divided in two halves by the lines i and i + L/2 corresponds to the r = 1/2 Cardy's geometry counted twice, in parallel. Thus, our * can be related to Cardy's number C by 1 * = (1 C)2, i.e. C = 0.876657(45). Compared with Cardy's exact value C Cwrap = 0.7569977963 &91;10&93; we found the contribution Cwrap = 0.119659(45) for the spanning configurations which also wrap around the periodic direction.
The critical occupation probability pc can be obtained from our expansion (1) by fixing any value for . The best choice is = *, for which the leading finite-size term in (1) becomes L11/n instead of L1/n. In this case, we get pc = 0.59274621(33), where the error bar is estimated by the same above-quoted procedure of dividing the whole data set into S = 10 sub-sets. For some unknown particular reason, the wrapping probabilities adopted in  works better yet, the leading finite-size term in (1) being L21/n. In  a numerical evidence for that behaviour is given, by plotting L(pc) * against L (assuming some previously determined value for pc) and verifying a power-law dependence with an exponent very close to 2. Indeed, for our spanning probability, the same plot gives an exponent very close to 1.
We propose the finite-size expansion (1) for spanning probabilities in percolation, where pL() is defined in Fig. 1.
This approach can be used to calculate various critical quantities. We applied it to a particular geometry within the site percolation problem on a L × L square lattice, considered as a torus. Taking two parallel lines distant L/2 from each other, our spanning probability counts configurations for which these two lines are connected. The universal value * = 0.984786(11) and the critical occupation pc = 0.59274621(33) are obtained.
As a by-product, we propose the universal value C = 0.876657(45) for the critical spanning probability within a L × L/2 annulus, i.e. a square lattice with periodic (free) boundaries along the direction of L (L/2). This probability corresponds to all configurations for which the frontiers separated by L/2 are connected. Cardy  determined the exact figure C Cwrap = 0.7569977963, where Cwrap corresponds to spanning configurations which also wrap along the direction of L, discounted in his approach.
The whole computer time for obtaining these data was 20 thousand hours, on a dozen computers, typically powered by an Athlon 1GHz processor.
Six months ago, the first author presented a preliminary version of this work in Campos do Jordão, São Paulo, Brazil, in honour of the 60th birthday of professor Silvio Salinas. Now, some other articles are presented in honour of the 60th birthday of the last author. In contrast, the younger authors dedicated this work to the senior professor Silvio Salinas.
We are indebted to Carlos Tomei, Antonio Brito Serbeto, John Cardy and Robert Ziff for helpful discussions, and Brazilian agencies CNPq and FAPERJ for support.
 C. Domb, E. Stoll, and T. Schneider, Contemp. Phys. 577 (1980); [ Links ]N. Jan, T. Lookman, and D. Stauffer, J. Phys. A16, L117 (1983); [ Links ]P.M.C. de Oliveira, S. Moss de Oliveira, and S.L.A. de Queiroz, Physica A175, 345 (1991). [ Links ]
 R.M. Ziff and M.E.J. Newman, Phys. Rev. E66, 016129 (2002). [ Links ]
 H.T. Pinson, J. Stat. Phys. 75, 1167 (1994). [ Links ]
 R.M. Ziff, C.D. Lorenz, and P. Kleban, Physica A266, 17 (1999). [ Links ]
 R.M. Ziff, Phys. Rev. Lett. 69, 2670 (1992). [ Links ]
 W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Numerical Recipes in C, Cambridge University Press (1988), section 14.3. [ Links ]
 J.-P. Hovi and A. Aharony, Phys. Rev. E53, 235 (1996). [ Links ]
 J. Cardy, J. Phys. A35, L565 (2002). [ Links ]
Received on 2 May, 2003