Bond percolation thresholds on Archimedean lattices from critical polynomial roots

We present percolation thresholds calculated numerically with the eigenvalue formulation of the method of critical polynomials; developed in the last few years, it has already proven to be orders of magnitude more accurate than traditional techniques. Here we report the result of large parallel calculations to produce what we believe may become the reference values of bond percolation thresholds on the Archimedean lattices for years to come. For example, for the kagome lattice we find $p_{\rm c}=0.524\,404\,999\,167\,448\,20 (1)$, whereas the best estimate using standard techniques is $p_{\rm c}=0.524\,404\,99(2)$. We further provide strong evidence that there are two classes of lattices: one for which the first three scaling exponents characterizing the finite-size corrections to $p_{\rm c}$ are $\Delta=6,7,8$, and another for which $\Delta=4,6,8$. We discuss the open questions related to the method, such as the full scaling law, as well as its potential for determining critical points of other models.

Introduction. Percolation is one of the simplest mathematical models with a phase transition [1]. It has served as a paradigm of such models, with basic properties that also emerge in a diverse range of systems from superconductivity [2] to black hole critical collapse [3]. In recent years the two-dimensional problem, which we focus on here, has proven to be particularly interesting, with its fascinating mix of solved and unsolved problems. Given a lattice, such as one of those shown in Fig. 1, choose each bond to be open with probability p and closed with probability 1 − p independently of all others. As p is increased, there is a critical threshold, p c , which marks a sharp transition to a regime with an infinite connected cluster. Despite the problem's apparent simplicity, the value of p c is unknown for most lattices, with exact solutions available only on a restricted class [4][5][6][7][8][9]. There have been great advances in the understanding of the continuum limit using conformal invariance and stochastic Loewner evolution [10][11][12][13], in which the details of the underlying lattice are irrelevant, but progress on unsolved lattice-dependent quantities, such as the critical probabilities of most of the Archimedean lattices, has been limited to the derivation of rigorous bounds [14][15][16][17][18][19][20] (though these are ever-tightening) and numerical studies [21][22][23]. For example, the critical probability for bond percolation on the kagome lattice is known rigorously to satisfy [20] 0.522551 < p c < 0.526490, a range with a width of 10 −3 , and the best estimate using traditional numerical techniques is p c = 0.52440499(2) [21], an accuracy of 10 −8 . It is very possible that in order to gain a complete understanding of subjects like conformal invariance and universality, a solution to, or at least a firmer grasp of, these lattice-specific problems will be necessary. In this Letter, we make progress towards this objective by pushing the precision of p c to the order of 10 −18 in the most favorable case. Critical polynomials. The method of critical polynomials originated from the observation that in all exactlysolved problems p c is the (unique) root, p c ∈ (0, 1), of a polynomial. For example on the triangular lattice (Fig.  1a), the bond percolation threshold is given by [4] so that p c = 2 sin π/18 ≈ 0.347296. Similar results are obtained for the hexagonal (Fig. 1b) and square (Fig. 1c) lattices. This polynomial can be generalized unambiguously even to problems for which the exact solution is not known [24][25][26][27]. This is done by first choosing a finite subgraph, B, called the basis, that generates the infinite lattice when copies are arranged in some periodic way. Next, assuming the percolation realization is identical on each copy of the basis, we use the label 2D for the event that there is an open cluster connecting every copy of B and 0D for the event that no infinite set of bases can be connected by open clusters. Denoting the probabilities of these events Z p (2D) and Z p (0D), the critical polynomial P B (p) is defined by This is clearly a polynomial in p as B has a finite number of edges. For reasons related to universality [27], the root of the polynomial in [0, 1] provides an estimate of the critical point that becomes more accurate as the size of B is increased. For problems with an exact solution, remarkably, the critical polynomial provides the correct answer for any choice of B-even the smallest possible- and this convenient property has been used to find some previously unknown exact solutions [25,27]. Transfer matrix. The origins, development and refinement of the method can be followed in a series of papers written over the last several years [24][25][26][27][28][29][30][31]. Although different methods have been used to compute P B [25,26,32], the transfer matrix [27,33] has proven to be by far the most efficient. It is in fact not really necessary to compute the polynomial explicitly if all one wants is the root, and in Ref. [30] it was shown by one of us that at this root, the largest eigenvalues of two transfer matrices, T open (p) and T closed (p), are equal. These two transfer matrices describe different topological sectors and build up respectively Z p (2D) and Z p (0D) on a semi-infinite cylinder of width n lattice spacings (the reader is referred to Ref. [30] for more details). The strategy is then to start with some value of p close to p c and continually compute the largest eigenvalues Λ open and Λ closed , adjusting to within some tolerance. The great advantage of this method is that it allows bases of size as large as n = 12 to be tractable on an ordinary desktop, with extrapolation used to find the estimate on the infinite lattice.
Here, we use a parallel implementation of the method, which is run on supercomputers, to reach n = 16 for bond percolation on the unsolved Archimedean lattices. These are the eleven two-dimensional lattices for which all vertices are equivalent (see Fig. 1). The square, triangular and hexagonal lattices have exact solutions, such as eq. (2), but the remaining eight are unsolved. Our results for those are shown in Table I, along with the best estimates using Monte Carlo [22] or traditional transfer matrix techniques [21].
Iterative scheme. To find the eigenvalues for a given value of p, we use the power iteration method, which involves simply multiplying an initial vector repeatedly by the transfer matrix, until it converges to the eigenvector with the largest eigenvalue. We demand convergence to 40 decimal places, which requires the use of an arbitrary precision library [35]. Note that the computations are actually done using v ≡ p/(1 − p), the bond weight in the Fortuin-Kasteleyn representation of the corresponding Q = 1 state Potts model [36]. Once we have the values of Λ open (v) and Λ closed (v), we use them to find the next value of v, and so on until the two eigenvalues are identical within our tolerance of δ ≡ 10 −40 . To adjust v, we use Householder's method of order 2 [37]. That is, if the equation to solve is f (v) = 0, then the m th iteration is given by For our purposes, To approximate the derivatives we use the central difference formulas where we set = √ δ, and One typically needs three or four Householder iterations to get v converged to the target precision δ. We take care to choose the initial v carefully, using extrapolations of the v c obtained for smaller n. Indeed, if |v − v c | 1 each second-order Householder iteration doubles the number of correct digits. Arriving at our largest size, n = 16, we are most often able to pick the initial v correct to about 10 −18 , so a single iteration gives sufficient precision for our purposes. Only in the case of the snub hexagonal lattice (Fig. 1j), we performed two iterations for n = 16 because we found our initial guess to be insufficiently accurate.
Parallelization. The parallel algorithm is used to compute Λ open and Λ closed at v and v ± , by the actual transfer matrix multiplication. It will be described more fully elsewhere, but in broad terms the goal is to distribute the components of the vector across processors in such a way that the need for inter-processor communication is minimized. Our implementation was inspired by Jensen's [38] transfer matrix enumeration of self-avoiding polygons. In that problem, Jensen identified a criterion for organizing the vector components that ensured that any transfer operation performed on data on a particular processor would not require information from any other. In our problem, as far as we know, it is not possible to replicate this exactly, but using a variant of Jensen's approach we were able to keep the need for communication to a minimum.
Up to n = 12, calculations can be done on ordinary desktops. Even then, the accuracy achieved by the method is far better than that of traditional techniques [30]; our calculations up to n = 16, should place these quantities permanently out of their reach. For the n = 16 computations, we used 10 3 processors. The time needed to complete a single power iteration varies with the lattice, but typically takes three to four hours. The number of power iterations needed to get convergence of a single eigenvalue is likewise variable and depends strongly on the initial vector used. When doing the first Householder iteration, we start with a vector with only one non-zero component. In this case, the number of power iterations needed is in the 40-60 range. However, for subsequent Householder iterations, we start with the final vector computed during the previous iteration and in this way we can reduce the number of power iterations needed to the 20-30 range as p approaches p * .
Extrapolation. Our extrapolation scheme for the resulting sequences is based on the empirical scaling form where the amplitudes A k and exponents ∆ k are dependent upon the lattice. There is currently little theoretical understanding of this scaling, and we compute the A k and ∆ k by simply fitting the actual data. Fortunately, we do now have a fair number of data points; in Table  II kagome lattice.
To begin, we determine ∆ 1 from the truncated form p c (n) = p c + A 1 n −∆1 . To eliminate the unknowns p c and A 1 , we form the combination q(n) ≡ p(n)−p(n−1) p(n−1)−p(n−2) . Assuming the truncated form, we have a non-linear relation that determines ∆ 1 (n) from three successive data points, p c (n), p c (n − 1) and p c (n − 2). Upon supposing a reasonable starting value, this determination is unique. Fig. 2a plots ∆ 1 (n) against n −1 for the kagome lattice along with a low-order polynomial fit. Trying various orders of the fit and eliminating the lowest values of n (the figure shows a fifth-order fit in n −1 to the last seven ∆ 1 (n)), we conclude that lim n→∞ ∆ 1 (n) = 6.00(2), cf. Table III. It appears safe to conjecture that this asymptotic value is ∆ 1 = 6 exactly for this lattice.
Next, we set p c (n) = p c (n) + A 1 n −6 , defining a new series p c (n), obtained from two successive sizes, in which the leading n −6 scaling term has been eliminated. We then repeat the above working, using now the truncated form p c (n) = p c + A 2 n −∆2 . The determinations of ∆ 2 (n) are shown against n −1 in Fig. 2b, and polynomial fits now lead to ∆ 2 = 7.00(5), from which we conjecture that ∆ 2 = 7 exactly. Table III compiles the scaling exponents computed in this way for the various lattices. For the first four lattices, we have data for all n ≤ 16. For the last four, the construction of their bases B requires n to be even [33], so we have only eight data points. Obviously this leads to more reliable determinations in the former cases. Taken jointly, the data of Table III suggest   two classes of non-solvable Archimedean lattices; those (kagome, three-twelve, and snub hexagonal) for which the first three exponents are 6, 7, 8, and the remaining five lattices for which they are 4, 6, 8. We find it compelling to conjecture that the complete set of exponents are all integers ≥ 6 for the former class, and all even integers ≥ 4 for the latter one.
Assuming this conjecture, we can extrapolate the p c (n) by means of the scaling form (11). This is done by identifying a stable compromise between the number of terms used in (11) and the number of low-n data points not included in the fits; see Refs. [30,39] for details. The end result is the central values and error bars given in Table I. We have checked that these are not sensitive to reasonable modifications of the values of ∆ k with k ≥ 4, and hence do not depend on the complete validity of the conjecture just made.
Discussion. Using a parallel implementation of the eigenvalue approach to critical polynomial roots, 3·10 6 CPU hours of supercomputer resources, and a comprehensive extrapolation method, we have obtained extremely accurate values of the bond percolation thresholds (see Table I) on the Archimedean lattices. In retrospect, the more naive extrapolation method used in Ref. [33] now makes some of the error bars cited there stand out as too optimistic. In the same vein, the leading scaling exponents reported in Ref. [31] can now be relegated to effective exponents for small n. However, the p c given in Ref. [30], based on the same scaling exponents as here, and also those of Ref. [31] are fully confirmed by the more precise values now obtained. This agreement-as well as other checks, including leaving out the n = 16 data point from the present analysis-lends credence to the error bars given in Table I. Our analysis for the scaling exponents reveals two different classes of non-solvable Archimedean lattices. This distinction cannot be explained solely from the difference between three-fold and four-fold rotational symmetries of the lattices [31], although it certainly is remarkable that all members of the first class (kagome, three-twelve, snub hexagonal) enjoy three-fold symmetries. The proposed scaling exponents are compatible with conformal field theory predictions [30], but it remains unclear how to derive them analytically.
In addition to the interesting open problems related to the critical polynomial method, it is yet to be determined how widely applicable it is. The definition (3) is easily generalized to the Q-state Potts model and gives excellent estimates for any Q [31,40], even in the imaginary temperature regime [41]. The eigenvalue method presented here was also adapted to compute the growth constant of self-avoiding walks and was able finally to rule out a longstanding conjecture [39]. Generalizations to site percolation, coupled Potts models, or to non-planar models, are currently under investigation.