# Bond percolation thresholds on Archimedean lattices from critical polynomial roots

Christian R Scullard, Jesper Lykke Jacobsen

► To cite this version:

Christian R Scullard, Jesper Lykke Jacobsen. Bond percolation thresholds on Archimedean lattices from critical polynomial roots. Physical Review Research, 2020, 2 (01), pp.012050. 10.1103/phys-revresearch.2.012050. hal-04978195

HAL Id: hal-04978195

https://hal.science/hal-04978195v1

Submitted on 5 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Bond percolation thresholds on Archimedean lattices from critical polynomial roots

Christian R. Scullard¹ and Jesper Lykke Jacobsen 2,3,4

¹Lawrence Livermore National Laboratory, Livermore, California 94550, USA

²Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, Paris, France

³Sorbonne Université, École Normale Supérieure, CNRS, Laboratoire de Physique (LPENS), 75005 Paris, France

⁴Institut de Physique Théorique, Université Paris Saclay, CEA, CNRS, 91191 Gif-sur-Yvette, France

<figure><img src="image_2.png" /></figure>

(Received 28 October 2019; accepted 6 January 2020; published 28 February 2020)

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_c = 0.524 404 999 167 448 20(1)$, whereas the best estimate using standard techniques is $p_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_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 the critical points of other models.

DOI: 10.1103/PhysRevResearch.2.012050

**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–9]. There have been great advances in the understanding of the continuum limit using conformal invariance and stochastic Loewner evolution [10–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–20] (though these are ever tightening) and numerical studies [21–23]. For example, the critical probability for bond percolation on the kagome lattice is known rigorously to satisfy [20]

$$0.522 551 < p_c < 0.526 490, \quad (1)$$

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.

a range with a width of $10^{-3}$, and the best estimate using traditional numerical techniques is $p_c = 0.524 404 99(2)$ [21], an accuracy of $10^{-8}$. It is very possible that in order to gain a complete understanding of subjects such as conformal invariance and universality, a solution to, or at least a firmer grasp of, these lattice-specific problems will be necessary. In this Rapid Communication, 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 exactly solved problems $p_c$ is the (unique) root, $p_c \in (0, 1)$, of a polynomial. For example, on the triangular lattice [Fig. 1(a)], the bond percolation threshold is given by [4]

$$p_c^3 - 3p_c + 1 = 0, \quad (2)$$

so that $p_c = 2 \sin \pi/18 \approx 0.347 296$. Similar results are obtained for the hexagonal [Fig. 1(b)] and square [Fig. 1(c)] lattices. This polynomial can be generalized unambiguously even to problems for which the exact solution is not known [24–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

$$P_B(p) \equiv Z_p(2D) - Z_p(0D). \quad (3)$$

This is clearly a polynomial in $p$ as $B$ has a finite number of edges. For reasons related to universality [27], the root of
<figure><img src="image_3.png" /><figcaption>FIG. 1. The 11 Archimedean lattices.</figcaption></figure>

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–31]. Although different

methods have been used to compute $P_B$ [25,26,32], the transfer matrix [27,29] 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 $\Lambda_{open}$ and $\Lambda_{closed}$, adjusting $p$ until

$$ \Lambda_{open} = \Lambda_{closed} \quad (4) $$

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 computer, 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 11 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 [34]. Note that the computations are actually done using $v \equiv p/(1-p)$, the bond weight in the Fortuin-Kasteleyn representation of the corresponding $Q = 1$ state Potts model [35]. Once we have the values of $\Lambda_{open}(v)$ and $\Lambda_{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 $\delta \equiv 10^{-40}$. To adjust $v$, we use Householder’s method of order 2 [36]. That is, if the equation to solve is $f(v) = 0$, then the $m$th iteration is given by

$$ v_{m+1} = v_m + 2 \frac{\left(\frac{1}{f}\right)'(v_m)}{\left(\frac{1}{f}\right)''(v_m)}. \quad (5) $$

<table><caption>TABLE I. Percolation thresholds for the nine unsolved Archimedean lattices, as compared to their previous numerical determinations. The letters refer to Fig. 1.</caption><thead><tr><th>Lattice</th><th></th><th>p<sub>c</sub> (This work)</th><th>p<sub>c</sub> (numerical)</th><th>Ref.</th></tr></thead><tbody><tr><td>Kagome</td><td>(d)</td><td>0.524 404 999 167 448 20(1)</td><td>0.524 404 99(2)</td><td>[21]</td></tr><tr><td>Four-eight</td><td>(e)</td><td>0.676 803 124 390 011 3(3)</td><td>0.676 802 32(63)</td><td>[22]</td></tr><tr><td>Frieze</td><td>(f)</td><td>0.419 640 358 863 69(2)</td><td>0.419 641 91(43)</td><td>[22]</td></tr><tr><td>Three-twelve</td><td>(g)</td><td>0.740 420 798 850 811 610(2)</td><td>0.740 421 95(80)</td><td>[22]</td></tr><tr><td>Cross</td><td>(h)</td><td>0.693 733 124 922(2)</td><td>0.693 733 83(72)</td><td>[22]</td></tr><tr><td>Snub square</td><td>(i)</td><td>0.414 137 856 591 7(1)</td><td>0.414 137 43(46)</td><td>[22]</td></tr><tr><td>Snub hexagonal</td><td>(j)</td><td>0.434 328 317 224 0(6)</td><td>0.434 328 0(5)</td><td>[33]</td></tr><tr><td>Ruby</td><td>(k)</td><td>0.524 831 461 573(1)</td><td>0.524 832 58(53)</td><td>[22]</td></tr></tbody></table>
For our purposes, $f(v) = \Lambda_{\text{open}}(v) - \Lambda_{\text{closed}}(v)$. To approximate the derivatives we use the central difference formulas

$$f'(v) \approx \frac{g_2 - g_0}{2\epsilon}, \qquad (6)$$

$$f''(v) \approx \left( \frac{g_2 - g_1}{\epsilon} - \frac{g_1 - g_0}{\epsilon} \right) \epsilon^{-1}, \qquad (7)$$

where we set $\epsilon = \sqrt{\delta}$, and

$$g_0 \equiv f(v - \epsilon), \qquad (8)$$

$$g_1 \equiv f(v), \qquad (9)$$

$$g_2 \equiv f(v + \epsilon). \qquad (10)$$

One typically needs three or four Householder iterations to get $v$ converged to the target precision $\delta$. We take care to choose the initial $v$ carefully, using extrapolations of the $v_c$ obtained for smaller $n$. Indeed, if $|v - v_c| \ll 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. 1(j)], 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 $\Lambda_{\text{open}}$ and $\Lambda_{\text{closed}}$ at $v$ and $v \pm \epsilon$, 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 [37] 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 desktop computers. 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 $\simeq 10^3$ processors. The time needed to complete a single power iteration varies with the lattice, but typically takes 3–4 h. 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 nonzero 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^*$.

<table><caption>TABLE II. Bond percolation thresholds for the kagome lattice computed on semi-infinite cylinders of width $n$.</caption><thead><tr><th>n</th><th>p<sub>c</sub></th></tr></thead><tbody><tr><td>1</td><td>0.524 429 717 521 274 793 546 879 681 534 455 071 6205</td></tr><tr><td>2</td><td>0.524 406 057 896 062 634 245 378 836 666 345 666 7920</td></tr><tr><td>3</td><td>0.524 405 092 218 718 391 406 491 710 278 995 604 5159</td></tr><tr><td>4</td><td>0.524 405 013 882 343 450 677 924 933 274 891 201 3263</td></tr><tr><td>5</td><td>0.524 405 002 666 098 533 997 468 638 043 799 737 1046</td></tr><tr><td>6</td><td>0.524 405 000 252 138 641 166 065 238 385 312 009 4089</td></tr><tr><td>7</td><td>0.524 404 999 570 802 604 857 648 689 641 603 340 3853</td></tr><tr><td>8</td><td>0.524 404 999 338 748 706 184 041 906 677 709 350 6317</td></tr><tr><td>9</td><td>0.524 404 999 247 980 209 806 701 838 958 679 653 4330</td></tr><tr><td>10</td><td>0.524 404 999 208 475 451 262 855 277 119 432 089 6813</td></tr><tr><td>11</td><td>0.524 404 999 189 755 973 511 801 309 010 812 928 2307</td></tr><tr><td>12</td><td>0.524 404 999 180 248 443 779 969 638 346 824 671 7858</td></tr><tr><td>13</td><td>0.524 404 999 175 132 845 053 820 303 018 487 609 0453</td></tr><tr><td>14</td><td>0.524 404 999 172 242 908 087 780 703 763 071 248 1530</td></tr><tr><td>15</td><td>0.524 404 999 170 540 780 670 080 646 173 449 291 2196</td></tr><tr><td>16</td><td>0.524 404 999 169 501 410 335 190 170 832 654 998 6109</td></tr><tr><td>&infin;</td><td>0.524 404 999 167 448 20 (1)</td></tr></tbody></table>

**Extrapolation.** Our extrapolation scheme for the resulting sequences is based on the empirical scaling form

$$p_c(n) = p_c + \sum_{k=1}^{\infty} \frac{A_k}{n^{\Delta_k}}, \qquad (11)$$

where the amplitudes $A_k$ and exponents $\Delta_k$ are dependent upon the lattice. There is currently little theoretical understanding of this scaling, but its form has previously been well established from fits to data [30,38]. The parameters $A_k$ and $\Delta_k$ must be determined from these fits, and the new data points added by the present work allow us to compute these for larger $k$ and thus to improve the accuracy of the extrapolations. While, given the large number of data points we have computed across a range of lattices (see Table II for the kagome thresholds for $n \le 16$), we have confidence that (11) is indeed the correct scaling behavior, it would obviously be preferable to have some means of computing the coefficients and exponents from theory. It is not yet known how to do this, but in recent work [39] Mertens and Ziff have shown that the critical polynomial can be reformulated in terms of the cluster density. Perhaps it will eventually be possible to relate the scaling theory of the latter [40] to that of the former. For now, however, we simply use (11) as an empirical law and we briefly describe below the method used to calculate the parameters $A_k$, $\Delta_k$ and the extrapolations to infinite width.

To begin, we determine $\Delta_1$ from the truncated form $p_c(n) = p_c + A_1 n^{-\Delta_1}$. To eliminate the unknowns $p_c$ and $A_1$, we form the combination $q(n) \equiv \frac{p(n)-p(n-1)}{p(n-1)-p(n-2)}$. Assuming the truncated form, we have

$$q(n) = \left(1 - \frac{2}{n}\right)^{\Delta_1} \frac{n^{\Delta_1} - (n-1)^{\Delta_1}}{(n-1)^{\Delta_1} - (n-2)^{\Delta_1}},$$

a nonlinear relation that determines $\Delta_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
<figure><img src="image_4.png" /><figcaption>FIG. 2. Effective scaling exponents $\Delta_1(n)$ and $\Delta_2(n)$, for the kagome lattice, shown against $n^{-1}$ with their extrapolations.</figcaption></figure>

unique. Figure 2(a) plots $\Delta_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 $\Delta_1(n)$], we conclude that $\lim_{n\to\infty} \Delta_1(n) = 6.00(2)$; cf. Table III. It appears safe to conjecture that this asymptotic value is $\Delta_1 = 6$ exactly for this lattice.

Next, we set $p_c(n) = \tilde{p}_c(n) + \tilde{A}_1n^{-6}$, defining a new series $\tilde{p}_c(n)$, obtained from two successive sizes, in which the leading $n^{-6}$ scaling term has been eliminated. We then repeat the above work, using now the truncated form $\tilde{p}_c(n) = p_c + A_2n^{-\Delta_2}$. The determinations of $\Delta_2(n)$ are shown against $n^{-1}$ in Fig. 2(b), and polynomial fits now lead to $\Delta_2 = 7.00(5)$, from which we conjecture that $\Delta_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 \le 16$. For the last four, the construction of their bases $B$ requires $n$ to be even [29], 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 that there exist two classes of nonsolvable 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 $\ge 6$ for the former class, and all even integers $\ge 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,38] 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 $\Delta_k$ with $k \ge 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, $\approx 3 \times 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. [29] now makes some of the error

<table><caption>TABLE III. Scaling exponents for the Archimedean lattices. Blank entries cannot be determined with sufficient precision.</caption><thead><tr><th>Lattice</th><th></th><th>Δ<sub>1</sub></th><th>Δ<sub>2</sub></th><th>Δ<sub>3</sub></th></tr></thead><tbody><tr><td>Kagome</td><td>(d)</td><td>6.00(2)</td><td>7.00(5)</td><td>8.1(2)</td></tr><tr><td>Four-eight</td><td>(e)</td><td>4.00(2)</td><td>6.0(1)</td><td></td></tr><tr><td>Frieze</td><td>(f)</td><td>4.00(2)</td><td>6.00(5)</td><td>8.1(2)</td></tr><tr><td>Three-twelve</td><td>(g)</td><td>6.04(2)</td><td>7.00(2)</td><td>8.0(1)</td></tr><tr><td>Cross</td><td>(h)</td><td>4.0(2)</td><td>6.0(5)</td><td></td></tr><tr><td>Snub square</td><td>(i)</td><td>4.000(2)</td><td>6.2(4)</td><td></td></tr><tr><td>Snub hexagonal</td><td>(j)</td><td>5.9(1)</td><td></td><td></td></tr><tr><td>Ruby</td><td>(k)</td><td>4.5(5)</td><td></td><td></td></tr></tbody></table>

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 nonsolvable Archimedean lattices. This distinction cannot be explained solely from the difference between three-fold and fourfold rotational symmetries of the lattices [31], although it certainly is remarkable that all members of the first class (kagome, three-twelve, snub hexagonal) enjoy threefold 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,41], even in the imaginary temperature regime [42]. 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 [38]. Generalizations to site percolation, coupled Potts models, or to nonplanar models are currently under investigation.

We thank Robert Ziff for sharing his numerical bond threshold of the snub hexagonal lattice reported in Table I. The work of C.R.S. was performed under the auspices of the U.S. Department of Energy at the Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 and was supported by the LLNL-LDRD Program under Project No. 19-DR-013. J.L.J. is grateful for support from the European Research Council under the Advanced Grant NuQFT. This work was granted access to the HPC resources of IDRIS under the allocation 2016-057751 attributed by GENCI (Grand Equipement National de Calcul Intensif).

[1] S. R. Broadbent and J. M. Hammersley, Proc. Cambridge Philos. Soc. 53, 629 (1957).

[2] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).

[3] M. W. Choptuik, Phys. Rev. Lett. 70, 9 (1993).

[4] M. F. Sykes and J. W. Essam, J. Math. Phys. 5, 1117 (1964).

[5] J. C. Wierman, J. Phys. A: Math. Gen. 17, 1525 (1984).

[6] C. R. Scullard, Phys. Rev. E 73, 016107 (2006).
[7] R. M. Ziff, Phys. Rev. E 73, 016134 (2006).

[8] R. M. Ziff and C. R. Scullard, J. Phys. A: Math. Gen. 39, 15083 (2006).

[9] B. Bollobás and O. Riordan, in *An Irregular Mind*, Bolyai Society Mathematical Studies Vol. 21 (Springer, Berlin, 2010), pp. 131–217.

[10] J. L. Cardy, J. Phys. A: Math. Gen. 25, L201 (1992).

[11] H. Saleur and B. Duplantier, Phys. Rev. Lett. 58, 2325 (1987).

[12] O. Schramm, Electron. Commun. Probab. 6, 115 (2001).

[13] S. Smirnov, C. R. Acad. Sci., Ser. I: Math. 333, 239 (2001).

[14] W. D. May and J. C. Wierman, Combin. Probab. Comput. 14, 549 (2005).

[15] J. C. Wierman, Combin. Probab. Comput. 12, 95 (2003).

[16] J. C. Wierman, Combin. Probab. Comput. 11, 629 (2002).

[17] J. C. Wierman, G. Yu, and T. Huang, Electron. J. Combin. 22, P2.52 (2015).

[18] W. D. May and J. C. Wierman, Combin. Probab. Comput. 16, 285 (2007).

[19] O. Riordan and M. Walters, Phys. Rev. E 76, 11110 (2007).

[20] J. C. Wierman, J. of Phys. A: Math. Theor. 50, 295001 (2017).

[21] X. Feng, Y. Deng, and H. W. J. Blöte, Phys. Rev. E 78, 031136 (2008).

[22] R. Parviainen, J. Phys. A: Math. Gen. 40, 9253 (2007).

[23] P. N. Suding and R. M. Ziff, Phys. Rev. E 60, 275 (1999).

[24] C. R. Scullard and R. M. Ziff, Phys. Rev. Lett. 100, 185701 (2008).

[25] C. R. Scullard and R. M. Ziff, J. Stat. Mech. (2010) P03021.

[26] C. R. Scullard, Phys. Rev. E 86, 041131 (2012).

[27] C. R. Scullard and J. L. Jacobsen, J. Phys. A: Math. Theor. 45, 494004 (2012).

[28] J. L. Jacobsen and C. R. Scullard, J. Phys. A: Math. Theor. 45, 494003 (2012).

[29] J. L. Jacobsen, J. Phys. A: Math. Theor. 47, 135001 (2014).

[30] J. L. Jacobsen, J. Phys. A: Math. Theor. 48, 454003 (2015).

[31] C. R. Scullard and J. L. Jacobsen, J. Phys. A: Math. Theor. 49, 125003 (2016).

[32] C. R. Scullard, J. Stat. Mech. (2011) P09022.

[33] R. M. Ziff (private communication).

[34] B. Haible, CLN – Class Library for Numbers, http://www.ginac.de/CLN/.

[35] C. M. Fortuin and P. W. Kasteleyn, Physica 57, 536 (1972).

[36] A. S. Householder, *The Numerical Treatment of a Single Non-linear Equation* (McGraw-Hill, New York, 1970).

[37] I. Jensen, J. Phys. A: Math. Gen. 36, 5731 (2003).

[38] J. L. Jacobsen, C. R. Scullard, and A. J. Guttmann, J. Phys. A: Math. Theor. 49, 494004 (2016).

[39] S. Mertens and R. M. Ziff, Phys. Rev. E 94, 062152 (2016).

[40] R. M. Ziff, S. R. Finch, and V. S. Adamchik, Phys. Rev. Lett. 79, 3447 (1997).

[41] J. L. Jacobsen and C. R. Scullard, J. Phys. A: Math. Theor. 46, 075001 (2013).

[42] H. Saleur, Nucl. Phys. B 360, 219 (1991).