The percolation critical polynomial as a graph invariant
Abstract
Every lattice for which the bond percolation critical probability can be found exactly possesses a critical polynomial, with the root in providing the threshold. Recent work has demonstrated that this polynomial may be generalized through a definition that can be applied on any periodic lattice. The polynomial depends on the lattice and on its decomposition into identical finite subgraphs, but once these are specified, the polynomial is essentially unique. On lattices for which the exact percolation threshold is unknown, the polynomials provide approximations for the critical probability with the estimates appearing to converge to the exact answer with increasing subgraph size. In this paper, I show how this generalized critical polynomial can be viewed as a graph invariant, similar to the Tutte polynomial. In particular, the critical polynomial is computed on a finite graph and may be found using the recursive deletioncontraction algorithm. This allows calculation on a computer, and I present such results for the kagome lattice using subgraphs of up to bonds. For one of these, I find the prediction , which differs from the numerical value, , by only .
I Introduction
Percolation is the study of the formation of random clusters on lattices. Given an infinite lattice, we declare each bond to be open with probability , and closed with probability . The resulting random clusters grow in average size with until we reach the critical point, , above which an infinite cluster appears. The determination of is an unsolved problem except in one dimension and on twodimensional lattices that are selfdual 3uniform hypergraphs Bollobás and Riordan (2010). An example is shown in Figure 1a, where the shaded triangle, excluding corner vertices, can represent any configuration of bonds, sites and correlations. Critical thresholds on these lattices are given by the Ziff criterion Ziff (2006),
(1) 
where is the probability that all three vertices are connected through open paths within the triangle, and is the probability that none are connected. Application of (1) to find a bond threshold results in a polynomial in , with order at most equal to the number of bonds in the unit triangle. We may consider each bond of an bond triangle to have a different probability, giving a critical surface of the generic form
(2) 
where is at most first order in any of its arguments, a property referred to as multilinearity. Examples of such critical surfaces include the square lattice,
(3) 
the honeycomb lattice,
(4) 
and the triangular lattice,
(5) 
The critical probability for the lattice is then found by setting all probabilities equal. For the square lattice, for example, this gives the polynomial , and the threshold .
In this paper, I will show how such critical surfaces may be generalized to any lattice, even those which are not in the solvable class, by employing a deletioncontraction algorithm. Such surfaces may not be exact, but we will see that they give approximations that, in principle, can be made arbitrarily precise. I begin by describing the deletioncontraction property of critical surfaces derived from (1). Then I show how this can be used to define the critical surfaces and polynomials for unsolved problems. Although it appears that the definition may not result in a unique critical surface in all cases, I present an argument to show that it is in fact always welldefined. I conclude by reporting generalized critical polynomials for the kagome lattice, giving approximations to the threshold that become increasingly more precise, with the best estimate within only of the numerically determined value. Although I use the kagome lattice as an illustrative example, this procedure can be applied on any periodic lattice.
Ii Deletioncontraction
Consider the martiniA lattice, with the probability assignments shown in Figure 1b. If we set , this bond is contracted and its end vertices merged, the result being the honeycomb lattice with two bonds doubled in parallel. This doubled bond can be replaced by a single effective bond, so we have
(6) 
By setting , we delete the bond, and the result is the square lattice with bonds doubled in series. Thus, we have,
(7) 
The critical surface of the martiniA lattice is given by (1). The only way to reduce to the correct deleted and contracted surfaces and preserve the required multilinearity property is to set
(8)  
(9)  
(10) 
which has the formal appearance of an average of the two special cases. It is a basic property of critical surfaces found using (1) that they satisfy such deletioncontraction formulas. That is, if such a lattice, , has a triangular unit of bonds such that deletion and contraction of the bond results in the lattices and , then we have for the critical surface of ,
(11)  
(12) 
As another example, we find for the martiniB lattice (Figure 1c),
(13)  
(14) 
Iii Generalized critical surfaces and polynomials
The deletioncontraction formula can be used to extend the definition of the critical surface to lattices for which the threshold is not known exactly, and in these cases it will be referred to as the generalized critical surface. For the kagome lattice of Figure 2a, an unsolved problem, we may write,
(15)  
(16)  
(17) 
where and are given by (10) and (14). Generalized critical polynomials are found by setting all probabilities equal, and in this case , i.e.,
(18) 
which gives Wu’s wellknown Wu (1979); Scullard and Ziff (2006) approximation , compared to the numerical Feng et al. (2008) . Thus, a generalized critical polynomial may be defined on any lattice for which known lattices appear upon deletion and contraction of a bond, but there is no guarantee that the resulting prediction will be exact. In fact, the definition extends to any periodic lattice. If is in the solvable class, the polynomial is given by (1), and otherwise it is defined by (12), with this formula applied recursively if either or have unknown threshold. By repeated application, known lattices eventually appear and the full critical surface can be found. That it satisfies such a deletioncontraction formula makes the generalized critical surface similar to other graph invariants, such as the Tutte polynomial and its special cases (e.g., the chromatic polynomial) Bollobás (1998).
In the following, I call the subgraph of a lattice, on which different probabilities are assigned, the “base” of the process, which is then tiled to form the infinite lattice. For example, in Figures 2a, 2b, 3 and 4 are bases for the kagome lattice consisting of one, two, four, and six unit cells. I now enumerate some properties of generalized critical surfaces and polynomials:

They are unique once the base and tiling are specified Scullard (2011); Wu (2010). In particular, critical surfaces are independent of the bonds chosen to delete and contract at each step. To prove this, we procede by induction. First, we call critical surfaces for which the final result is independent of deletioncontraction bond order welldefined. Clearly, any selfdual 3uniform hypergraph has this property, as the polynomial is defined by (1); deletioncontraction must give this answer, regardless of the bond chosen, because it only imposes the correct boundary values on the critical surface. Similarly, any critical surface that is a special case of a welldefined surface is itself welldefined. Next, we show that the critical surface for a lattice, , that is not a realization of a selfdual 3uniform hypergraph, and is thus not solved by (1), is welldefined if every critical surface resulting from deleting and contracting any bond of is welldefined. The kagome lattice with the bond base is one example of such a system. Assuming that has this property, consider the generic formula (12). The question is whether this is the same as what we would obtain by choosing the bond instead,
(19) (20) i.e., whether . We can apply deletioncontraction on to the graphs , and on to to obtain , etc. We then have,
(21) (22) (23) (24) and
(25) (26) (27) (28) and clearly these two expressions are the same if , , etc. By assumption, , , and are welldefined. Now, the function arises by contracting the bonds and in , which is the same way we arrive at . Thus, these two functions must describe the same lattice with a welldefined critical surface, and so we have . By the same reasoning, , etc., and we conclude that if every critical surface arising from deleting and contracting any bond of is welldefined, then so is . Because each step of the algorithm lowers the number of bonds in the base by (at least) one, welldefined surfaces eventually appear and the general result follows by induction. A slight complication arises due to the existence of certain seemingly pathological cases. Consider, for example, setting on the square lattice. Now we have created a onedimensional system, with the trivial critical point , or . This latter form is the correct surface to use whenever the onedimensional case appears, and is indeed reflected in the formula (3). A less straightforward example appears on the square lattice when we set in (3). It is not obvious how to interpret the result of contracting this bond since it seems to collapse the lattice. However, in such cases it is better to think of the bond as not contracting to zero length but simply carrying probability and thus introducing infinitely long linear clusters on the lattice. In the square case, the system is supercritical unless we also set , which gives a critical onedimensional system. The supercritical phase of (3) is represented by , so by the above considerations we would expect the surface for to read , which is indeed what we find from (3). Even if the lattice is unknown, by keeping these “supporting” bonds uncontracted when they are used in the algorithm, we can continue until we reach cases on which we may use (1), and then the supporting bond is just given probability in the formula. Finally, consider the result of setting in the triangular formula (5), thus imposing an infinite twodimensional cluster. This system is supercritical (there is no infinite cluster at the critical point) regardless of the value of , and thus we have , which is a general rule for whenever such a situation arises. This completes the discussion of potentially problematic cases, and thus the proof that the critical surface for a general lattice and base is welldefined. When calculating the polynomials reported below, each one was computed several times using different bond paths to ensure that the code was working properly. For the largest bases of bonds, each calculation was repeated over a dozen times and the same polynomial always resulted.

If the singlecell prediction is not exact, it is generally very close, usually within of the numerically determined threshold Scullard and Ziff (2010). It is not clear why this should be so, but it has been clearly demonstrated for many different systems, including all the Archimedean lattices Scullard and Ziff (2010); Scullard (2011). In addition, proper choices of larger bases lead to predictions closer to the exact answer, with accuracy increasing with base size. It is this latter conjecture I seek to support here by computing polynomials on the kagome lattice for bases of up to bonds.

Conversely, if the critical polynomial gives the exact threshold for a base of a single unit cell, as is the case for selfdual uniform hypergraphs, then any critical polynomial using a larger base makes the same prediction, i.e. the original polynomial always factors out. This is also a conjectured property but it would seem to be necessary for consistency, and I have found no counterexamples. While it would be ideal to have a lattice for which the exact threshold is known and is not a root of a polynomial, so that we can directly check the convergence of these polynomials, no such example is known at present and the conjecture can only be tested against numerically determined values.

In many cases, if the singlecell polynomial does not give the exact answer, then it can be shown Scullard (2011) that no critical polynomial derived from a finitesized base will give the correct percolation threshold. This is true for the kagome lattice, as discussed below.
Generalized critical polynomials have been found for all the Archimedean lattices Scullard and Ziff (2008, 2010); Scullard (2011). The polynomial (18) was originally found by Wu using his “homogeneity” Wu (1979) assumption and recently, he extended the method to include kagome subnets Wu (2010), with excellent results Ding et al. (2010). In fact, the homogeneity approximation gives a prediction for the full state Potts critical frontier, but at present it appears limited to the kagome lattice and its subnets as it relies on a transformation from the triangular lattice to these kagometype graphs. Although I focus on the kagome lattice in this paper, mostly because of the interest it has attracted over the years (e.g., Ziff and Gu (2009); Tsallis (1982); Wu (1979)), the method presented here may be applied on any periodic lattice for either site or bond percolation.
Iv Kagome lattice polynomials
In previous work, critical polynomials were found “by hand”, with the lattice and its bond unit cell probably representing the limit of what one would be inclined to do this way (see the appendix of Scullard (2011)). However, the recursive nature of the algorithm makes it an ideal problem for a computer, and here I report the results of using a program to calculate critical polynomials on the kagome lattice for bases containing , , and bonds. At each step, the program chooses a bond, and finds the two graphs that result from its deletion and contraction. It knows a small number of graphs (e.g., triangular, honeycomb, and square lattices), and it repeats the deletioncontraction algorithm recursively until it recognizes all the lattices it has found. The output is a set of function definitions, like equations (10), (14) and (17) which can be evaluated in a computer algebra package to get the critical surface and then the critical polynomial. Of course, there are many issues to overcome in programming this scheme, and a full account is given elsewhere Scullard (2012).
iv.0.1 Base of 2 unit cells
The first extension beyond a single unit cell base is that shown in Figure 2b, in which we employ a base using two unit cells which are indicated by different colors, with the tiling shown in Figure 2c. The polynomial can be written in the factored form
(29)  
We recognize the first term in brackets as the polynomial for the bond base. The second term has no real roots, and thus the prediction here is the same as the one we found previously, . There is no other way to tile this base to give a different result.
iv.0.2 Bases of 4 unit cells
Using the base consisting of unit cells shown in Figure 3a, we may tile it in two different ways as indicated in Figures 3b and 3c. Starting with Figure 3b, the critical polynomial can be written in the factored form
(33)  
Once again, the first term in brackets is just the bond polynomial (18), and the others have no roots in . The prediction is again the same as the bond estimate.
Things finally become more interesting when we use the staggered embedding in Figure 3c. In this case, the polynomial is
(34)  
and the solution on is , which differs from the numerical result by , a great improvement over the singlecell bond case. Although finding all the bases that give different polynomials is something of an art, especially as the number of bonds becomes large, there does not appear to be one that gives a different prediction using bonds.
iv.0.3 Bases of 6 unit cells
There are many options for bases of bonds, and I will discuss only three here. If we take the 6unit cell version of the base shown in Figure 3a with the embedding analogous to Figure 3b, we once again find a polynomial in which (18) appears as a factor. Thus, for this base and embedding, we get the bond estimate again. It is tempting to assume that this trend continues for larger bases and embeddings of this type.
Turning now to the base shown in Figure 4a, in which the embedding is indicated by the matching shapes on the external vertices, we get the polynomial,
(35)  
(36)  
(37)  
(38)  
(39)  
(40)  
(41)  
(42)  
(43)  
(44) 
which has solution on , differing from the numerical value by . A different base and embedding is shown in Figure 4b and has polynomial,
(45)  
(46)  
(47)  
(48)  
(49)  
(50)  
(51)  
(52)  
(53)  
(54)  
(55) 
with solution in , , slightly better than the previous estimate and within of the numerical value.
V Discussion
Clearly, we are observing the predictions converging to the exact answer with increasing base size. The results suggest that, if a generic base can be described as containing unit cells, the exact answer is approached only as both and go to since the bases appear all to provide the same incorrect prediction. From these few examples it also appears that the kagome estimates are converging from above, as they all seem to be greater than the numerical value. However, I know of no argument that guarantees this trend will continue. Also, as the polynomials respect duality, i.e., making the substitution gives a polynomial for the dual graph, the estimates for the dual of the kagome lattice, the dice lattice, would converge from below. Note also that no polynomial for the kagome lattice derived from any finitesize base will give the exact answer. This can be seen as follows. Using the base in Figure 4a, we delete many of the bonds to leave a single unit cell and a few connecting bonds, as in Figure 5a. Now, we may contract the connecting bonds to recover the kagome system in which the base is a single cell, as in Figure 5b. By uniqueness, the prediction for this case must be equation (18), which contradicts the prediction found by setting all bonds equal in the full critical surface, i.e. equation (44). Thus, if the singlecell polynomial does not provide the exact threshold, and for the kagome lattice it does not, then, although we can get arbitrarily close by using ever larger bases, no finite critical polynomial derived in this way will ever solve the problem.
There are many possible directions for further study. Aside from more firmly establishing the various conjectures, it would be interesting to try to quantify the manner of convergence to the exact solution, perhaps through a version of finitesize scaling. Another avenue is the extension to three and higher dimensions. While the polynomials are welldefined by the deletioncontraction algorithm in any dimension, preliminary results seem to indicate that they are not as successful at predicting higherD critical points. However, this will be the subject of future work.
I have presented critical polynomials for the kagome lattice, up to bases of bonds. This is the limit of the current implementation of the algorithm, as it becomes progressively more difficult to add bonds due to the exponential complexity. It is not uncommon for a large calculation to produce over a million function definitions. Nevertheless, there is room for improvement in the efficiency, as the rule used to choose the bond for deletion and contraction at each step can have a significant effect on the rate of reduction to the known simple cases. I have hardly explored this issue and presently use what amounts to a random bond selection, rejecting only those choices that lead to undue complications (such as the contraction of supporting bonds described above). Moreover, this algorithm is perfect for a parallel implementation as it would require little interprocessor communication. It remains to be determined how much extra performance can be wrought from these considerations, but the problem will hopefully be seen as an interesting computational challenge.
Acknowledgements.
I am grateful to Oliver Riordan for many valuable suggestions. I also thank Robert Ziff for the fruitful collaboration that led to this work, and Jesper Jacobsen for informative discussions. Finally, I thank an anonymous referee for providing several helpful comments. This work was performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract No. DEAC5207NA27344.References
 Bollobás and Riordan (2010) B. Bollobás and O. Riordan, in An Irregular Mind, Bolyai Society Mathematical Studies, Vol. 21 (Springer Berlin Heidelberg, 2010) pp. 131–217, arXiv:1001.4674.
 Ziff (2006) R. M. Ziff, Phys. Rev. E 73, 016134 (2006).
 Wu (1979) F. Y. Wu, J. Phys. C 12, L645 (1979).
 Scullard and Ziff (2006) C. R. Scullard and R. M. Ziff, Phys. Rev. E 73, 045102(R) (2006).
 Feng et al. (2008) X. Feng, Y. Deng, and H. W. J. Blöte, Phys. Rev. E 78, 031136 (2008).
 Bollobás (1998) B. Bollobás, Modern Graph Theory (Springer, New York, 1998).
 Scullard (2011) C. R. Scullard, J. Stat. Mech. 2011, P09022 (2011).
 Wu (2010) F. Y. Wu, Phys. Rev. E 81, 061110 (2010).
 Scullard and Ziff (2010) C. R. Scullard and R. M. Ziff, J. Stat. Mech. 2010, P03021 (2010).
 Scullard and Ziff (2008) C. R. Scullard and R. M. Ziff, Phys. Rev. Lett. , 185701 (2008).
 Ding et al. (2010) C. Ding, Z. Fu, W. Guo, and F. Y. Wu, Phys. Rev. E 81, 061111 (2010).
 Ziff and Gu (2009) R. M. Ziff and H. Gu, Phys. Rev. E 79, 020102 (2009).
 Tsallis (1982) C. Tsallis, J. Phys. C: Solid State Phys. 15, L757 (1982).
 Scullard (2012) C. R. Scullard, J. Stat. Mech. 2012, P11011 (2012).