A Disproof of Tsallis' Bond Percolation Threshold Conjecture for the Kagome Lattice

In 1982, Tsallis derived a formula which proposed an exact value of 0.522372078. .. for the bond percolation threshold of the kagome lattice. We use the substitution method, which is based on stochastic ordering, to compare the probability distribution of connections in the homogeneous bond percolation model on the kagome lattice to those of an exactly-solved inhomogeneous bond percolation model on the martini lattice. The bounds obtained are 0.522394 < p c (kagome) < 0.526750, where the lower bound shows that the value conjectured by Tsallis is incorrect.


Introduction
Since the origins of percolation theory in the 1950's, rigorous determination of the exact percolation thresholds of lattices has been a major challenge.Early 1980's results [12,13,29,30] determined the exact bond percolation thresholds of the square, triangular, hexagonal and bow-tie lattices and the site percolation thresholds of the triangular and kagome lattices, using graph duality properties and the star-triangle transformation.No further progress on this problem was made until 2006, when a generalized star-triangle transformation and the concept of self-dual planar hypergraphs led to a method for solving for bond percolation thresholds for an infinite class of two-dimensional lattices [25,38,44,46].However, this method does not produce solutions for several common models, the most well-studied being the kagome lattice bond model and the square and hexagonal lattice site models.For lattices for which there is no exact solution, knowledge regarding the percolation threshold is gained via simulation estimates and mathematically rigorous bounds.
Figure 1: Induced subgraphs of the kagome lattice (left) and dice lattice (right).
The most studied unsolved two-dimensional bond percolation model is the kagome lattice model, for which there is an extensive literature of estimates and bounds.Before the 1990s, simulations, approximations, heuristic derivations, and renormalization approaches [3,4,19,20,21,24,28,39,41,42,43] produced a collection of estimated values ranging from 0.449 to 0.526.Since the late 1990s, as more efficient simulation algorithms were developed and computing power increased, the bond percolation threshold estimates [5,9,10,11,22,23,26,27,40,45,47] agree on the first four decimal places, at 0.5244, with the highest-precision simulation values clustering near 0.524405.Particularly notable is a mathematically rigorous 99.9999% confidence interval [0.52415, 0.52465] by Riordan and Walters [23].Mathematically rigorous bounds for the percolation threshold have also improved over the years, as indicated in the following  [18] This article is motivated by the disparity between the simulation estimates and the conjectured exact value by Tsallis [28], the root of p 3 − p 2 − p + 1 − 2 sin(π/18) = 0.For example, a recent estimate is 0.524404999134 by [9], while the numerical value of the Tsallis conjecture is 0.522372078 . . .We obtain mathematically rigorous bounds of 0.522394 < p c (kagome) < 0.526750, the electronic journal of combinatorics 22(2) (2015), #P2.52 which disprove the Tsallis conjecture by showing that in fact his value is strictly smaller than the kagome lattice bond percolation threshold.Note that the bounds are improvements of [18], and are approximately centered on the recent simulation estimates.Note also that, previous to the Tsallis conjecture, Wu [39] conjectured a percolation threshold of 0.524430 for the kagome lattice bond model.The discrepancy between Wu's conjecture and the consensus estimates is much smaller than for the Tsallis conjecture, and our bounds are not sufficient to disprove Wu's conjecture.
Since the dual graph of the kagome lattice is the dice lattice, the bounds imply the following bounds for the dice lattice bond percolation threshold: 0.473250 < p c (dice) < 0.477606.
Our bounds are obtained by the substitution method, using stochastic ordering to compare an unsolved percolation model to an exactly-solved model.In addition to the kagome lattice bond model, it has been applied to several other percolation models [17,18,33,34,35,36].The initial motivation for the development of the substitution method was to understand the approach of Ottavi, who proposed bounds of 0.522372 p c 0.528924 for the kagome lattice bond percolation threshold.A discussion of the substitution method, including the flaws in Ottavi's approach, can be found in [1,Ch. 6].
We compare the kagome lattice bond model to a two-parameter inhomogeneous model on the martini lattice, which is exactly solved by the generalized star-triangle transformation approach.The generalized star-triangle transformation showed how to generate a large collection of new lattices for which the percolation threshold can be exactly solved.The martini lattice has similarities to the kagome lattice which make it well-suited to be used as the reference lattice for calculating bounds for the kagome lattice bond percolation threshold.
The paper is organized as follows: Section 2 describes the martini lattice, the solution for the critical surface of the two-parameter inhomogeneous model on it, and the substitution regions of the martini and kagome lattices used.Computational issues are discussed in section 3. The apparently intractably large computation can be reduced substantially by using planarity and symmetry.A two-stage algorithm is used to determine the bounds by first generating a candidate value that is expected to be close to the bound, which greatly reduces the set of possibilities that need to be checked in the second stage.Finally, section 4 summarizes our conclusions and their consequences, and discusses further investigations in progress.

The Lattices
The kagome lattice, shown in Figure 1, is one of the eleven Archimedean lattices, which are vertex-transitive graphs created by tilings of the plane by regular polygons [8,Ch. 2].It is the line graph of the hexagonal lattice, so the site percolation threshold of the kagome lattice is equal to the bond percolation threshold of the hexagonal lattice [6], which is exactly 1 − 2 sin( π 18 ) [29].For any lattice, the bond percolation threshold is less than or equal to the site percolation threshold [7], so, for a crude upper bound, the bond percolation threshold satisfies p c (kagome) 1 − 2 sin( π 18 ) < 0.652704.Its dual graph is the dice lattice, one of the 8 Laves lattices which were identified in crystallography [14,15].By Kesten's duality theorem [13] for bond percolation thresholds of dual graphs, p c (kagome) + p c (dice) = 1, so bounds for one threshold imply bounds for the other.
The substitution method compares probabilities of connections between vertices in the percolation model on an unsolved lattice with those on a solved reference lattice.Previous attempts to disprove the Tsallis conjecture, by May and Wierman [17,18], used the hexagonal lattice bond model as the reference model, since the hexagonal lattice has some structural similarity with the kagome lattice.Since then, a class of new exact solutions has been established.(See [2] and [38].)However, as pointed out by Ziff [44], the kagome lattice is not in this class, so an exact solution does not exist using current methods.The class of new solutions includes the martini lattice, shown in Figure 2, which is used as the reference lattice in this article.There are two principal reasons for using the martini lattice.First, its structure has more similarities with the kagome lattice than the hexagonal lattice.This is shown in Figure 3, which shows corresponding substitution regions from the three lattices.One may view the martini lattice as replacing three three-stars in the two-subdivided hexagonal lattice by triangles, and consider it to be "halfway between" the hexagonal lattice and kagome lattice structures.Intuitively, probabilities of connections between vertices are expected to be closer for the kagome and martini lattices than for the kagome and hexagonal lattices.
Secondly, an inhomogeneous bond percolation model on the martini lattice, with one parameter t for the triangular edges and a second parameter s for the other edges, is exactly solved.This provides an infinite collection of solved reference models for comparison to the kagome lattice bond model.Figure 4 shows the martini generator labeled with edge probability parameters, and its critical surface in [0, 1] 2 , respectively.The critical surface is the boundary in the parameter space separating regions of existence and non-existence of infinite clusters, which is given by 2s 3 t 3 − 3s 3 t 2 − 3s 2 t 3 + 3s 2 t 2 + 3s 2 t = 1 for the martini lattice model [44].Three particular points on the critical surface are labeled in Figure  Although increased computations are required, we find the optimal points on the critical surface, trading off one parameter for the other, for computing the best upper and lower bounds for the kagome lattice bond percolation threshold.
Although the martini lattice was solved using a particular generator, we are not obligated to use a union of those generators in the substitution method.However, it is important that the entire lattice can be constructed as an edge-disjoint union of copies of the generator.We use the larger substitution regions shown in Figure 3 rather than that shown in Figure 4.
These two advantages of the martini lattice allow us to obtain improved bounds with a smaller substitution region.In the terminology of May and Wierman [17,18], they used the electronic journal of combinatorics 22(2) (2015), #P2.52 a "seven tile" substitution to obtain the bounds 0.522197 < p c (kagome) < 0.526873.This article improves both bounds, by comparing a smaller "one tile" substitution region with the inhomogeneous martini lattice bond model reference.

Substitution Method Computations
Many of the computational methods are similar to those of May and Wierman [17,18].Modifications or replacements of these methods arose to deal with the fact that there is not a single reference model at criticality, but the full range of models on the critical surface for the martini lattice.The following subsections provide brief descriptions and remarks about various aspects of the computational methods.

Substitution Method
The substitution method compares two percolation models using stochastic ordering.There are several steps in implementing substitution method computations.
First, a finite connected subgraph, called the substitution region, is chosen for each model.These substitution regions must have the same number of vertices on the boundary; that is, vertices through which paths in the infinite lattice may enter or leave the region.Each lattice must be an edge-disjoint union of isomorphic copies of its substitution region.
Secondly, a probability measure on the partitions of the boundary vertices of the substitution region is computed for each model.For each of the substitution regions in Figure 3, label the boundary vertices 1 through 6 clockwise from the top.Use vertical bars to indicate separation of the boundary vertices into blocks.The set of boundary partitions is a partially ordered set with the refinement ordering, in which a partition is smaller than another if each of its blocks is contained in a single block of the other.For each parameter value p, and each partition π of the boundary vertices, P K p (π) denotes the probability that the boundary vertices in each block of π are connected by open paths, but no vertices in different blocks are connected by open paths, within the substitution region of the kagome lattice percolation model.For each pair of parameter values s, t, P M s,t (π) is defined similarly for the martini lattice substitution region.The number of boundary vertices of the substitution regions must be equal in order for the two partition probability measures to be defined on the same partially ordered set of partitions.
Thirdly, the families of probability measures P K p and P M s,t are then compared by stochastic ordering.An upset in the partial ordering is a set of partitions U such that if a ∈ U and a b then b ∈ U .A probability measure P is stochastically smaller than a probability measure for every upset U in the set of partitions.By the well-known equivalence of stochastic ordering and coupling, if (s, t) is a point on the critical surface of the martini lattice model, and P K p is stochastically larger than P M s,t , then p p c (kagome).Similarly, if P K p is stochastically smaller than P M s,t and (s, t) is a point on the critical surface of the martini lattice model, then p p c (kagome).
the electronic journal of combinatorics 22(2) (2015), #P2.52 The following subsections discuss the implementation of steps in this approach, emphasizing the differences from previous implementations.All computations discussed in the following were performed using MATLAB.

Generating Partitions
The first step in the computations is to generate the appropriate lattice of partitions of the boundary vertices.Since our substitution regions have six boundary vertices, the number of partitions is the 6-th Bell number, which is 203.However, there exist 90 different partitions which have exactly 3 blocks, and are thus mutually incomparable.Each of the 2 90 subsets of this set of partitions, except the empty set, generates an upset which would be involved in checking stochastic ordering later.
Fortunately, since both lattices are planar, the substitution regions are planar, so only non-crossing partitions may have positive probability.A non-crossing partition of [n] = {1, 2, 3, . . ., n} is a partition such that a < b < c < d ∈ [n] where a and c are in block A and b and d are in block B implies that A = B. Since the number of non-crossing partitions of six boundary vertices is the 6-th Catalan number, 132, computations are considerably reduced by restricting to non-crossing partitions.In fact, there are three non-crossing partitions which also have zero probability in both models, 14|23|56 and two rotations, so the number of partitions with nonzero probability is actually 129.Using the lattice of non-crossing partitions results in a considerable reduction of the number of upsets that must be checked.

Reduction to a Class Lattice
A further reduction of the computational burden is accomplished using symmetry of the substitution regions.Partitions which are rotations or reflections of each other, so are in the orbit of each other under the dihedral group, have equal probabilities in the percolation model.Thus, the set of 129 partitions can be reduced to equivalence classes of partitions which are related by rotation or reflection.Wierman [37] showed that, for the purpose of checking stochastic ordering, it is sufficient to consider the class lattice, which is a partially ordered set in which a class π is considered to be a refinement of a class σ if there is a partition in π which is a refinement of a partition in σ.
Both substitution regions have some rotational and reflection symmetries.The symmetry group for the kagome lattice is the dihedral group on 6 elements, but for the martini region it is the dihedral group on 3 elements.The class lattice for the martini substitution region contains 34 classes, but the class lattice for the kagome substitution region has 24.Since we must compare probability measures for the two models on a common partially ordered set, and each kagome class is a union of martini classes, we must use the class lattice with 34 classes.The partitions in each class are isomorphic both in the kagome substitution region and in the martini substitution region.

Partition Probability Functions
The substitution regions for the kagome lattice and the martini lattice both are comprised of 18 edges.Since each edge may be retained in the random subgraph or not, there are 2 18  possible configurations that may be realized.For each k, 0 k 18, the probability of a realization on the kagome substitution region with k retained edges is p k (1 − p) 18−k .For the martini substitution region, for each i, j, 0 i, j 9, any configuration with i retained edges in triangles and j retained edges in three-stars has probability t i (1−t) 9−i s j (1−s) 9−j .Each partition is the union of a set of configurations, so its probability is a sum of these terms, and is thus a polynomial.In principle, the partition probability polynomials may be computed by generating all configurations, determining which configurations comprise each partition, and adding the configuration probabilities to obtain each partition probability.However, we employed the more efficient "graph-welding" method introduced in [16].By carrying out these computations symbolically, the polynomial functions were found explicitly, in contrast to finding numerical values for each retention probability parameter value as needed.

Checking Stochastic Ordering
We follow the approach of [37] to find upper and lower bounds for the kagome lattice bond percolation threshold for a given point (s, t) on the martini lattice critical surface.With the values of s and t specified, the partition probability measure associated with the martini lattice is constant, while each nontrivial upset in the kagome probability measure has a probability function that is strictly increasing in p.Thus, for each non-trivial upset U in the class lattice, there is a unique value p U at which the two probabilities are equal.Let p max denote the maximum solution of all such upset equations, and p min denote the minimum solution.If p > p max , we have P K p (U ) P M s,t (U ) for all upsets U , so P K p st P M s,t , and thus p p c (kagome).Similarly, if p < p min , we have p p c (kagome).Between p min and p max , neither stochastic ordering holds.Therefore, the best upper and lower bounds that we can obtain from the comparison of the two models with the two substitution regions are p max and p min respectively.
However, in general, finding the largest and smallest upset probability equation solutions directly, by solving all upset equations, is impractical, since the number of upsets grows super-exponentially in the number of boundary vertices.This may be likened to trying to "find two needles in a haystack."Thus, when computing a lower bound, it is important to be able to rule out large collections of upsets that could not possibly give the minimum solution.We accomplish this using a process which is a modification of the approach of Wierman [37] and May and Wierman [17,18].

Desirable Classes
In this section, we define "desirable classes" and show that the largest lower bound p min is the solution for an upset U min which has only desirable classes as its minimal elements.(A similar result holds for the smallest upper bound.)Since we only need to solve upset the electronic journal of combinatorics 22(2) (2015), #P2.52 equations for upsets with this property, the computational burden is significantly reduced.
The desirable class approach can be adapted to other models, so the following reasoning is given more generally than for the kagome and martini lattices.We use the notation g A (p) for the probability of a set A of classes in the unsolved percolation model with parameter value p, and h A for the (constant) probability in the solved percolation model at criticality.In the context of this article, g A (p) = P K p (A), the kagome lattice class probabilities, and h A = P M s,t (A), the martini lattice probabilities for a fixed point (s, t) on the critical surface for the martini lattice.
We begin with definitions and preliminary facts about upsets.
Fact 1.If P is a finite poset, then every non-empty upset U of P has an antichain representation of the form where {x 1 , x 2 , . . ., x k } is an antichain in P and {↑ x i } = {π ∈ P : π x i } is the principal upset generated by x i .The elements {x 1 , x 2 , . . ., x k } are the minimal elements of upset U .This provides a recipe for constructing all the upsets, though it is practical only when the poset is relatively small.Fact 2. If U is an upset, and a is a minimal element of U , then U \{a} is an upset.
To apply our principal result, the interval [l low , l up ] must be chosen so that it contains the value p min .Even though p min is unknown, this may be done.An appropriate interval may be computed using smaller substitution regions, as is discussed in Remark 5 after the following theorem.It is beneficial to have a short interval, so the number of desirable classes is small, to reduce the computational effort required.
Theorem 4. Suppose that p min ∈ [l low , l up ] for some values 0 l low l up 1, with the corresponding upset denoted by U min .Then every minimal element of U min is desirable with respect to [l low , l up ].
Proof: Let A be the antichain consisting of the minimal elements of U min .We will prove by contradiction that each a ∈ A is desirable.Assume there exists an element a ∈ A which is not desirable.Then By assumption, p min ∈ [l low , l up ], so the fact that a is not desirable implies that Consider U min \{a}, which is an upset by Fact 2. By definition of p min and U min , the electronic journal of combinatorics 22(2) (2015), #P2.52 which may be decomposed into Consequently, g U min \{a} (p min ) > h U min \{a} .
Since U min \{a} is an upset, its probability is a nondecreasing function of p, so there exists p * < p min such that g U min \{a} (p * ) = h U min \{a} , which contradicts the fact that p min is the minimum of the solutions of all the upset equations.Thus, our assumption that a minimal element of U min is not desirable must be false, so we conclude that all minimal elements of U min are desirable.
Remark 5.An appropriate interval [l low , l up ] may be determined from substitution method computations from smaller regions.Notice that the 6-boundary-vertex substitution region is the edge-disjoint union of three isomorphic 4-boundary-vertex regions.Let p l and p u denote the substitution method lower and upper bounds, respectively, computed for the 4boundary-vertex regions (for the same point (s, t) on the martini lattice critical surface).By the equivalence of stochastic ordering and coupling, independent couplings on the 4-boundary vertex regions exist, which can be combined to provide a coupling on the 6-boundary-vertex region.This implies that the kagome lattice percolation model at p l is stochastically smaller than the martini lattice model on the 6-boundary-vertex region, so p l p min .Similarly, for the upper bounds, p u p max .Thus, we may choose [p l , p u ] as the interval [l low , l up ] for application of the Theorem.
Remark 6.In fact, the interval [l low , l up ] may be shortened further.Note that any upset equation solution provides a value that is at least as large as the minimum solution, and thus may be used as l up .A simple choice is to use the solution for the upset consisting of only the maximal class.For our calculations, we considered the classes in the order generated by our programs, solved the upset equations as the classes were successively added to form upsets, then used the smallest of these solutions for our initial value of l up .
To apply Theorem 4 in a straightforward manner, one would identify all desirable classes, construct all the upsets generated by sets of desirable classes, solve their upset equations, then find the smallest solution, to obtain the lower bound for the percolation threshold.However, if there are k desirable classes, the number of subsets of desirable classes is 2 k , so the computational efficiency can be improved considerably by two additional reductions in the number of classes needed to generate the equations.In calculations for this article, they reduce the typical running time for calculating a bound from more than 6 to 20 hours to 2 to 10 minutes.Remark 7. Note that the maximum class in the class poset must be in every upset, so it must be in the upset that generates the minimum upset solution.Suppose there is a desirable class C for which every class D C in the class poset is either desirable or is the electronic journal of combinatorics 22(2) (2015), #P2.52 the maximum class.By reasoning as in the proof of Theorem 4, C must be included in the upset that generates the minimum upset equation solution.Therefore, to determine the smallest solution, we need only consider the union of such classes with the upsets generated by the other desirable classes.Remark 8. Definition 3, of desirable class, and Theorem 4 were stated in the context of determining a lower bound for the percolation threshold of the unsolved model.One may reverse the inequality in Definition 3 and restate Theorem 4 in terms of downsets rather than upsets in the class lattice to obtain a companion result, since the complement of an upset is a downset and a downset is generated by its maximal elements.Most classes are desirable for the upset approach only or the downset approach only, with some desirable for both.We determine whether the upset approach or downset approach has the smallest number of desirable classes, and typically reduce the computational effort considerably by performing the calculations with that approach.
Finally, we note that the approach discussed above may be easily adapted to determine the upper bound for the percolation threshold.

Determining Bounds
The desirable class method described for checking stochastic ordering can be applied using any point (s, t) on the critical surface of the inhomogeneous martini lattice bond model.To determine the best possible bounds for these substitution regions, for each value s = .66,.67,.68, . . ., 1.00, we determined the corresponding t so that (s, t) is on the critical surface, then computed the upper and lower bounds for p c (kagome).These upper and lower bounds are graphed in Figure 5.The largest lower bound for this set of s-values was not larger than Tsallis' conjectured value.However, by iterating this discretization on successively shorter intervals, we determined the upper bound of .526750 at s = 0.804942 and the lower bound of .522394 at s = 0.806921.The results of the finest discretization are also shown in Figure 5.
The calculations were done using MATLAB on personal computers.For various points (s, t) on the critical surface of the martini lattice, the running time varied from approximately two to ten minutes.

Conclusion and Future Research
We have disproved the conjectured value of Tsallis for the bond percolation threshold of the kagome lattice, by improving both the upper and lower bounds due to May and Wierman [18], which required computations for a much larger substitution region.The improvement was possible because there is now a wider class of exactly solved bond percolation models.Among those models, the inhomogeneous bond percolation model on the martini lattice appeared to provide a better match to the connectivity of the kagome lattice than the hexagonal lattice bond model used by May and Wierman [18], and its use did in fact produce better percolation threshold bounds.
In the future, we expect to improve the methods of this article and apply them to additional unsolved percolation models.Specifically, there are two prime candidates for our next investigations: Tsallis also conjectured an exact bond percolation threshold value of p c = 0.739830 . . .for the Archimedean (3, 12 2 ) lattice, which the research literature suggests is smaller than the true value.For example, a recent high-precision estimate of Jacobsen [9] provides the value p c ≈ 0.7404207988474 (7).Our comparison with the two-parameter martini lattice for the substitution region with six boundary vertices provides the following bounds: 0.739227 < p c (3, 12 2 ) < 0.740880, which are roughly centered on Jacobsen's estimate.The lower bound does not disprove Tsallis' conjecture, but the upper bound is an improvement from the best previous upper bound, 0.741125, obtained by May and Wierman [18].
The hexagonal lattice site model is one of the most important unsolved site percolation model in two dimensions.May and Wierman [18] computed the best known bounds for it using the substitution method, comparing it to the hexagonal lattice bond model.A comparison with the inhomogeneous martini lattice model may produce improved bounds.

Figure 2 :
Figure 2: An induced subgraph of the martini lattice.

Figure 4 : 4 .
Figure 4: Left: The martini generator with two edge probability parameters.Right: The martini lattice model critical surface.

Figure 5 :
Figure 5: Graphs of bounds for various values of s: Left: Upper and lower bounds for the full range of values.Right: Lower bounds in a small interval containing the maximum.