High-precision percolation thresholds and Potts-model critical manifolds from graph polynomials

The graph polynomial P_B(q, v) for the Potts model depends on a finite part of the lattice, called the basis B, that generates the infinite lattice G by an appropriate set of translations, called the embedding [14]. Each regular lattice G admits an infinite number of choices for B. We shall be interested in the simplest possible family of B, called square bases in [14]. They have a chequerboard structure, shown in figure 2, consisting of alternating grey and white squares.

Zoom In Zoom Out Reset image size

The grey squares in figure 2 can contain any arrangement of lattice edges and internal vertices. Two adjacent grey squares meet at a common vertex; we shall call such shared vertices the terminals of the grey square. The lattices G which can be represented as in figure 2 are called four-terminal lattices. The Potts model on such G is not in general exactly solvable, as discussed in the Introduction. For the simplest lattices, all the grey squares contain the same arrangement of edges and internal vertices. For more complicated lattices it is necessary to let the structure of the grey squares depend on the coordinates (x, y) with some periodicity.

A rectangular (resp. square) basis is an array of n × m (resp. n × n) grey squares. We shall almost exclusively be interested in square bases. However, for a few lattices we shall need to decorate the grey squares with a periodicity that is different in the x and y directions, and in those cases we might need rectangular bases in order to respect that periodicity. The transfer matrix construction to be described below also allows for decorating (some of) the white squares by a single horizontal edge. For simplicity, we shall still say that the corresponding lattices are of the four-terminal type.

The accuracy of the critical manifold determined by P_B(q, v) = 0 increases rapidly with n, but unfortunately the same is true for the computational effort required to compute P_B(q, v). To be precise, the contraction–deletion algorithm described in [13] has time and memory requirements that grow exponentially in nm, whereas for the transfer matrix algorithm of [14] the growth is only exponential in min(n, m). The improved transfer matrix algorithm to be described here is again exponential in min(n, m) but with a smaller growth constant, enabling us to access larger sizes.

The algorithm of [13] was capable of computing P_B(q, v) for bases with up to 36 edges. Since the most interesting lattices have typically at least six edges per grey square, this means that bases of size 2 × 2 or 2 × 3 could be handled. The limit of feasibility using the transfer matrix of [14] was improved to min(n, m) = 4. In this work we further improve the transfer matrix algorithm, putting n × n square bases with n = 7 within reach. To be precise, we compute exactly the two-variable Potts polynomial P_B(q, v) up to n = 5; the exact one-variable bond percolation polynomial $P_B(q=1,v=\frac{p}{1-p})$ up to n = 6; and roots in v to 50-digit numerical precision of the equation P_B(q, v) = 0 for selected values of q up to n = 7.

The site percolation polynomials P_B(p) can sometimes be computed for even larger n, namely up to n = 11 for the square lattice, and even n = 16 for the ruby lattice (see section 7). For the site percolation problems we invariably compute the exact polynomial P_B(p), refraining from accessing any additional gain that might have been obtained by finding only the relevant root (0 < p_c < 1) with finite numerical precision.

The vertices situated at a corner of a grey square in figure 2, and not shared between two distinct grey squares, are called the terminals of the basis B. The embedding of B into G is defined by gluing distinct copies of B at the terminals⁴. This can be done in a variety of ways, but in this work we are only interested in the simplest possibility, called straight embedding in [14]. It consists of simply translating the n × m basis horizontally through multiples of ne_x, and vertically through multiples of me_y, where e_x (resp. e_y) is a unit vector in the x direction (resp. y direction).

The graph polynomial P_B(q, v) was initially defined from a deletion–contraction principle [13]. It was subsequently shown [14] that it can be equivalently written as a linear combination of conditioned partition functions, similar to (1), defined on a graph which is equal to the basis B. Consider first the partition function Z of the Potts model defined on B = (V, E), where we have imposed toroidal boundary conditions, i.e., opposite terminals are identified both horizontally and vertically. We can carry out the decomposition

$$\begin{equation} Z = Z_{\rm 2D} + Z_{\rm 1D} + Z_{\rm 0D} , \end{equation} \tag{ 5 }$$

where Z_0D is the sum over edge subsets A⊆E such that all connected components (clusters) in the subgraph G_A = (V, A) have trivial homotopy on the torus (i.e., are contractible to a point), and Z_1D is the sum over terms where there exists both a cluster and a dual cluster with non-trivial homotopy. In simpler words, Z_2D regroups the terms where there is a cluster that spans both spatial directions; the terms in Z_1D contain a cluster that spans only one, and not both, of the directions in space; and in Z_0D there are no spanning clusters.

In this set-up the graph polynomial reads [14]

$$\begin{equation} P_B(q,v) = Z_{\rm 2D} - q Z_{\rm 0D} . \end{equation} \tag{ 6 }$$

The term in Z corresponding to each A⊆E can be assigned to either Z_2D, Z_1D or Z_0D by using the Euler relation, as explained in [14]. We shall come back to this technical consideration in section 3.7.

We shall refer to (6) as the cluster representation of P_B(q, v). Below we shall use an equivalent formulation in terms of a loop model [28] defined on the medial lattice $\mathcal{M}(B)$. To this end, consider the two possible states of an edge e = (ij) ∈ E between to adjacent vertices i, j in B:

In the left picture we have e ∈ A and the edge is drawn as a thick blue line. The corresponding loops (thin red lines) reflect off the edge e. In the right picture we have e∉A and the edge is shown in dashed line style. The loops then cut through e, or equivalently, they reflect off the corresponding dual edge e*. In both cases (e ∈ A and e∉A) the loops separate the clusters in G_A = (V, A) from the dual clusters in $G_{A^*} = (V^*,A^*)$, where by definition A* consists of the edges dual to those in E∖A.

In a transfer matrix formalism the lattice is built up, row by row, starting at the bottom and moving towards the top. Each row is in turn built up, edge by edge, starting at the left and moving towards the right. One can think of this as sweeping some imaginary d − 1-dimensional surface over the lattice in a number of discrete time steps. We refer to this surface as the time slice. Note that since we are in d = 2 dimensions it is actually just a curve. More precisely, it is a horizontal line each time a row of the lattice has been completed, and a horizontal line with a kink when the row is only partially completed. At each step, the partially built lattice (the portion below the time slice) is characterized by the connectivity state of the clusters or loops intersecting the time slice, in a precise way that makes it possible to compute the partition function—or, in our case, the conditional partition functions entering in (6)—in the transfer process. The possible configurations of the edges (e ∈ A and e∉A) must be summed over, and each term in the sum induces a definite operation on the connectivity states.

A detailed description of the transfer matrix formalism in the cluster representation was given in [14]. We here need the corresponding loop representation, describing the states and operations on the thin red lines in (7). We begin with a brief outline, deferring a more precise description of a number of important points to the following subsections.

For a planar graph one may use [28] the Euler relation to rewrite the partition function (1) in the loop representation as

$$\begin{equation} Z = q^{|V|/2} \sum _{A \subseteq E} x^{|A|} n_{\rm loop}^{\ell (A)} , \end{equation} \tag{ 8 }$$

where $x = v / \sqrt{q}$, and ℓ(A) denotes the number of closed loops induced by the configuration A. The loop fugacity is $n_{\rm loop} = \sqrt{q}$. Note that because of (7) there is a local bijection between configurations of clusters and loops, so we may use the notation A⊆E to specify the configurations of loops as well.

Consider then adding a single edge to the lattice. Imagining for the moment that the transfer direction is upwards (the lattice is built from the bottom to the top), we shall refer to the horizontal edge in (7) as a 'space-like' edge. When e∉A (right picture) the two loop strands just go though and nothing happens, while for e ∈ A (left picture) the two strands are being connected and a new partially completed loop is started out. The sum over both possibilities can be described by an operator acting on adjacent loop strands at positions i and i + 1:

$$\begin{equation} {\sf H}_i = {\sf I} + x {\sf E}_i , \end{equation} \tag{ 9 }$$

where ${\sf E}_i$ are the generators of the Temperley–Lieb algebra [31] defined by

$$\begin{eqnarray} &&{\sf E}_i^2 = n_{\rm loop} {\sf E}_i , \nonumber \\ &&{\sf E}_i {\sf E}_{i \pm 1} {\sf E}_i = {\sf E}_i , \\ &&{\sf E}_i {\sf E}_j = {\sf E}_j {\sf E}_i \quad {\rm for }\quad |i-j| > 1 . \nonumber \end{eqnarray} \tag{ 10 }$$

The algebraic relations (10) can be proved graphically by gluing several diagrams of the type (7) on top of one another.

For a vertical, or 'time-like' edge, we would have the graphical correspondence

and since the situation e ∈ A now corresponds to the right picture, the corresponding operator is

$$\begin{equation} {\sf V}_i = x {\sf I} + {\sf E}_i . \end{equation} \tag{ 12 }$$

In the following subsections we explain in detail how to adapt the loop representation to the geometry of figure 2 and use it to compute the graph polynomial (6) rather than the full partition function Z. In particular, we explain how to discard the configuration contribution to Z_1D, distinguish topologically the contributions to Z_2D and Z_0D, and attribute to each diagram the correct powers of $x = v/\sqrt{q}$ and $n_{\rm loop}=\sqrt{q}$. On the basis of this, the graph polynomial (6) can the be retrieved by changing the (n_loop, x) variables back to (q, v) and providing the overall factor q^|V|/2 in (8).

The graph polynomial P_B(q, v) must be computed in the square-basis geometry (cf figure 2) which in the loop representation takes the appearance shown in figure 3. The loops live on the medial lattice $\mathcal{M}(B)$, of which a part is shown as red and blue edges in figure 3. In a nomenclature inspired by the theory of quantum integrable systems, we shall refer to the (red) horizontal edges as 'auxiliary spaces' and the (blue) vertical edges as 'quantum spaces'.

Zoom In Zoom Out Reset image size

The part of $\mathcal{M}(B)$ inside the grey squares depends of course on the lattice being studied. We denote it symbolically by the letter R. The operator constructing the corresponding part of the lattice is called the $\check{R}$-matrix and is written as $\check{\sf R}_i$. It acts on two auxiliary spaces (i, i + 1) coming from the left and two quantum spaces (i + 2, i + 3) coming from the bottom of the grey square, and produces outgoing quantum spaces (i', i' + 1) on the top and auxiliary spaces (i' + 2, i' + 3) on the right of the grey square. This labelling of spaces is shown in figure 4.

Zoom In Zoom Out Reset image size

For each of the lattices to be studied in this paper, $\check{\sf R}_i$ is a definite product of the elementary operators ${\sf H}_j$, ${\sf V}_j$ and ${\sf E}_j$ defined above, with j = i, i + 1, i + 2. This means in practice that the computation of P_B(q, v) can be adopted to any desired lattice that can be shown to have the four-terminal structure of figure 3 and for which the algebraic expression for $\check{\sf R}_i$ can be provided (see section 4).

To add one row of grey squares to the lattice, one inserts a pair of auxiliary spaces, acts with the product $\check{\sf R}_{2n-2} \cdots \check{\sf R}_4 \check{\sf R}_2 \check{\sf R}_0$, and removes the pair of auxiliary spaces. The next subsection contains the precise definitions of the connectivity states, and describes how the elementary operators act on them, and how auxiliary spaces are inserted and removed.

Note that the labelling of figure 4 implies that time flows in the north-east direction⁵ and not simply upwards as we assumed for pedagogical reasons when discussing (7) and (11). This has an incidence on the interpretation of 'horizontal' and 'vertical', since these geometrical notions have now been rotated 45° in the clockwise direction. To be precise, the epithet horizontal, or space-like (resp. vertical, or time-like) now means perpendicular (resp. parallel) to the direction of time flow. For instance, the $\check{R}$-matrix that constructs the square lattice is written as

$$\begin{equation} \check{\sf R}_i = {\sf H}_{i+1} {\sf V}_i {\sf V}_{i+2} {\sf H}_{i+1} . \end{equation} \tag{ 13 }$$

Many more examples will be given in section 4.

As already mentioned it is possible also to perform a restricted set of operations in the white squares. Since after completing a row of grey squares the auxiliary spaces are no longer at our disposal, this is essentially limited to letting the operator ${\sf H}_i$ act on the two quantum spaces within a white square in figure 3. This has the effect of adding a horizontal diagonal to the white square.

Since the terminals on the top and bottom of figure 3 must eventually be glued together, we actually need a set of two time slices. The first one is a horizontal line at the bottom of the system ($y=-\frac{1}{2}$), and the second one is a horizontal line (with a kink when a row is only partially completed) that keeps moving upwards (and the kink moving to the right, i.e., in the 'north-east direction', while completing a given row) until the lattice is completed. After completion, the two time slices are glued together in a precise way (see section 3.7) that distinguishes the contributions to Z_2D and Z_0D.

3.5.1. Description of the states

It is convenient first to describe the state space for a complete row. In this case each of the time slices is a horizontal line that intersects the 2n quantum spaces. No auxiliary spaces are involved in this case. A connectivity state describes how these 4n intersection points are pairwise connected by the loop strands.

In figure 5 we show a few examples of connectivity states. We name as an arc a loop segment that connects two points within the same time slice, and as a string a loop segment that connects points on different time slices. The number of strings is denoted s. Boundary conditions are periodic in the horizontal direction, as shown by the dashed sides of the construction boxes in figure 5.

Zoom In Zoom Out Reset image size

Suppose that the points are labelled 0, 1, 2, ..., 2n − 1 on the bottom time slice and 0', 1', 2', ..., (2n − 1)' on the upper time slice. To relate the cluster and loop representations, it is important to observe that points with an even label have a cluster on the right and a dual cluster on the left (and vice versa for the points with an odd label). The clusters (resp. dual clusters) corresponding to the connectivity states in figure 5 are shown in grey (resp. white) shading. This relation to clusters imposes important restrictions on the connectivities of loops. First, arcs connect points of opposite parities. Second, strings connect points of the same parity. In particular one may define a string to be even or odd, depending on the parity of the points that it connects. Third, s is even and there are as many even as odd strings. When s = 0, consider the configuration of arcs restricted to just one of the two time slices. We call this configuration closed if the equivalent clusters are bounded away from the other time slice (in which case one dual cluster will connect the two time slices), and open if one cluster connects the two time slices (in which case the dual clusters are bound). Then, fourth, when s = 0 a connectivity state consists of two open arc configurations, or of two closed arc configurations.

These concepts are illustrated in figure 5. Figure 5(a) shows a state with s = 2. The case of s = 0 is depicted in figure 5(b) for two open arc states, and in figure 5(c) for two closed arc states. The figures also exemplify a coding of the states which turns out to be convenient for describing the transfer matrix algorithm. On the time slice on the top of the system, the leftmost (resp. rightmost) point of an arc carries the code 1 (resp. 2). Arcs on the bottom of the system are coded similarly, but after a 180° rotation. Strings are coded as pairs of matching codes (3, 4, ...).

3.5.2. The action of Temperley–Lieb generators

3.5.3. The dimension of the transfer matrix

To determine the dimension of the transfer matrix we must count the number of states. We still suppose for the time being that the uppermost row is complete (i.e., there are no auxiliary spaces).

By cutting all the strings (if any), a connectivity state describing the full system of two time slices is transformed into a pair of reduced states each associated with one of the time slices. The reduced states consist of arcs and half-strings. When s = 0 the reduced states can be characterized as open or closed, just like the full states.

It is easy to count the number of reduced states. For s = 0, out of the 2n points there are n with code 1 and n with code 2. A moment's reflection reveals that all the ${2n \choose n}$ possible placements of these codes correspond to a valid state. For s = 2k > 0 one can similarly convince oneself that specifying the n − k points with code 1 will uniquely imply the positions of the arcs and half-strings. So there are ${2n \choose n-k}$ reduced states in general.

Conversely, a pair of reduced states with the same number of half-strings can be transformed into a full connectivity state by gluing pairs of half-strings. However, for s = 0 one obtains a valid state only by 'gluing' (or rather juxtaposing, since nothing is actually being glued!) two open or two closed half-states. Since the set of open and closed half-states are bijectively related by performing a cyclic shift, there are $\frac{1}{2} {2n \choose n}$ of each. Moreover, for s = 2k > 0 strings the gluing can be done in k inequivalent ways, since each half-string must be glued to one of the same parity and the cyclic order of strings must be respected. These observations imply that there are

$$\begin{equation} {\rm dim}(n) = \frac{1}{2} {2n \choose n}^2 + \sum _{k=1}^n k {2n \choose n-k}^2 \end{equation} \tag{ 14 }$$

connectivity states.

3.5.4. Bijection between states and integers

3.5.5. Handling auxiliary spaces

To handle a partially completed row of the lattice we need to be able to insert and remove auxiliary spaces. We also need to count the number of states in the presence of p auxiliary spaces. (The four-terminal representation requires p = 0, 1, 2 but extensions of the formalism to higher values of p may turn out to be of interest for lattices which cannot be cast in the four-terminal form.)

Suppose now that the n points in the top time slice are initially labelled from left to right:

The bottom time slice simply occupies the 2n labels following those of the top time slice, and plays no further role in the following construction. Accordingly we shall not represent it here.

Inserting the first (lower) auxiliary space amounts to adding two extra points to the left of those in the top time slice. This can be considered as two copies of the same point which are initially connected. We have now:

The 2n points of the quantum spaces have had their labels shifted by 1, in order to accommodate the labels 0 and 2n + 1 of the quantum spaces. In order to connect the latter two points, and respect the conventions for an arc on the top time slice, we attribute to them the codes c_{2n + 1} = 1 and c₀ = 2.

Similarly, the second (upper) auxiliary space is inserted by adding a pair of points with codes 1 and 2 in between those previously inserted. This looks as follows:

Note that the quantum space labels have again been shifted by 1, as have those of the first (lower) auxiliary space. This is necessary for respecting the cyclic order of the labels upon moving around the top time slice. So the complete set of four points inserted to the left of those in the top time slice have the codes c_{2n + 2} = c_{2n + 3} = 1 and c₀ = c₁ = 2. To add one complete row of grey squares to the lattice, one then applies the product of operators $\check{\sf R}_{2n-2} \cdots \check{\sf R}_4 \check{\sf R}_2 \check{\sf R}_0$. Note that before the action with the factor $\check{\sf R}_i$ the two rightmost dangling ends of the auxiliary spaces carry the labels i and i + 1, in agreement with the conventions of figure 4.

After the application of the last factor, $\check{\sf R}_{2n-2}$, the situation is as follows:

The row is basically completed, and the new quantum spaces have come out with the correct labelling 0, 1, 2, ..., 2n − 1. However, the two auxiliary spaces remain open, so the dangling ends on the left and the right will have to be glued and then removed. To this end, one applies the 'contraction' operator (discussed above when defining the Temperley–Lieb generators ${\sf E}_i$) to identify first points 2n + 1 and 2n + 2, and then 2n and 2n + 3. The four auxiliary points are then removed from the state, and we are ready to start all over and add a new row to the lattice.

The number of states dim(n, p) in the presence of p auxiliary spaces is an obvious generalization of (14). We find that

$$\begin{equation} {\rm dim}(n,p) = \frac{1}{2} {2n+2p \choose n+p} {2n \choose n} + \sum _{k=1}^n k {2n+2p \choose n+p-k} {2n \choose n-k} . \end{equation} \tag{ 15 }$$

We have tabulated these dimensions in table 2 for the values of n and p used in the computations of section 4. These dimensions should be compared with those of the transfer matrix approach of [14]. In order to compute P_B(q, v) for an n × n square basis, [14] used states which are planar partitions of the N = 4n terminal points in figure 2. The number of such states is given by the Catalan number

$$\begin{equation} {\rm Cat}(N) = \frac{1}{N+1} {2N \choose N} . \end{equation} \tag{ 16 }$$

The advantage of the present approach, as announced in the introduction, is to reduce this to dim(n, 2), thus raising the limit of practical feasibility of the computations from n_max = 4 [14] to n_max = 7. This improvement is achieved by dealing efficiently with the horizontal periodic boundary conditions, leading to the reduction from 4n to 2n of the number of terminals implied by figure 3.

Table 2. Dimension dim(n) = dim(n, 0) of the transfer matrix for a completed row, and dimensions dim(n, p) for a partially completed row with p = 1, 2 auxiliary spaces, for the sizes n used in section 4. This is compared with the dimension of the transfer matrix of [14].

n	dim(n, 0)	dim(n, 1)	dim(n, 2)	Reference [14]
1	3	10	35	14
2	36	132	490	1 430
3	500	1 900	7 245	208 012
4	7 350	28 420	109 956	35 357 670
5	111 132	433 944	1 693 692	6 564 120 420
6	1 707 552	6 708 240	26 332 020	1 289 904 147 324
7	26 501 904	104 535 288	411 945 105	263 747 951 750 360

The bijection between states and integers extends straightforwardly to the case with auxiliary spaces, the only difference being that the reduced state describing the top time slice contains 2p extra points.

We now describe some considerations on how to efficiently build the lattice by repeated applications of the transfer matrix. The general prescription for building a single row has been detailed in section 3.5.5:

• (i)  
  Insert two auxiliary spaces.
• (ii)  
  Add a row of grey squares, each corresponding to the application of an operator $\check{\sf R}_i$.
• (iii)  
  Remove the two auxiliary spaces.
• (iv)  
  (For certain lattices:) add horizontal edges in the white square by application of the ${\sf H}_i$ operators.

At the beginning of the transfer process, the top and bottom time slices coincide. Therefore the initial state consists of just a single connectivity state s₀ in which, for each i = 0, 1, ..., 2n − 1, the points i and i + 2n are connected by a string. In the conventions of section 3.5.1 this means that the corresponding coding is c_i = c_{i + 2n} = 3 + i. In other words, the initial state is therefore a unit vector where s₀ has Boltzmann weight 1, while all other states have weight 0.

For the computation of P_B(q, v), the states are represented as arrays of polynomials in the variables (n_loop, x), with coefficients that are non-negative integers. The degrees of the polynomials are |V| in the n_loop variable and |E| in the x variable, where |V| and |E| denote the numbers of vertices and edges in B. For the computation of the percolation critical polynomial P_B(1, v) it suffices obviously to employ arrays of polynomials in the single x variable. And, finally, when we only desire to find the roots of P_B(q, v) numerically, the arrays consist of real numbers to the desired numerical precision. In order to obtain the roots to (at least) 50 digits we use 100-digit real numbers and the Newton–Raphson method, where derivatives are computed from a first-order finite-difference formula with = 10⁻⁵⁰. In practice we use the CLN library [36], which efficiently handles such high-precision real numbers within our C++ implementation of the algorithm.

In all cases, the arrays have dimensions dim(n, p) given by (15), where p = 0, 1, 2 depending on the number of open auxiliary spaces.

The bijection described in section 3.5.4 is employed throughout the transfer process to translate back and forth between connectivity states (on which the action of the fundamental Temperley–Lieb operators ${\sf E}_i$ has been detailed in section 3.5.2) and integers 0, 1, ..., dim(n, p) − 1 that specify the position in the arrays of the relevant Boltzmann weights.

In the cases where the polynomials P_B(q, v) or P_B(1, v) are computed exactly, the coefficients of the polynomials are very large integers for all but the smallest values of the size n. Although the CLN library [36] offers also arbitrary-precision integers, it is more efficient to compute the result modulo a sufficient number of different primes p_i and reconstitute the exact result from the Chinese remainder theorem. The standard version of our algorithm (i.e., the one used throughout section 4) employs only additions, whereas the 'generic $\check{\sf R}$-matrix' version described in section 6.1 uses also multiplications. Since standard unsigned integers in C++ lie in the range 0, 1, ..., 2³² − 1 we therefore take p_i < 2³¹ in the former case and p_i < 2¹⁶ in the latter. The most demanding computations required of the order of 20 different primes within this scheme.

Some consideration on the application of an $\check{\sf R}_i$ operator are also in order. For most lattices—including all those of section 4—we can express $\check{\sf R}_i$ as a product of the elementary operators ${\sf H}_i$, ${\sf V}_i$ and ${\sf E}_i$. An example has been given in (13). In those cases it is usually numerically most efficient to apply each of the factors separately, i.e., to build the lattice one edge at a time. (This procedure is known as sparse-matrix factorization.) Note that the application of each of these elementary operators to a given input connectivity state produces only two output states. Expanding out the product $\check{\sf R}_i$ by hand would lead instead to up to fourteen output states (since the number of planar pairings of the eight points appearing in figure 4 is Cat(4) = 14), and with vastly more complicated coefficients. However, in some cases not all 14 output states are actually generated. This is so in particular for the kagome lattice, where only 13 states are generated. Accordingly some of the states stored in the arrays will have zero weight. For instance, for the kagome lattice we find that only 37% of the dim(n, 2) states carry non-zero weight for n = 1, but this occupation ratio increases to 74% for n = 2, 87% for n = 3 and 90% for n = 4. It is therefore hardly worth dealing with this slight waste of memory resources in this case.

However, in other situations the occupation ratio may go instead to zero for large n. This is notably the case for some of the site percolation problems of section 7. It is then more efficient to avoid storing a lot of zero coefficients in the array, and instead insert the states that are really generated (with non-zero weight) in a hash table. In particular, the bijection of section 3.5.4 becomes superfluous. We shall describe the hashing version of the algorithm in more detail in sections 6.1 and 7.

When the basis B has been completely built up, it remains to discuss how to actually compute the graph polynomial, for which we recall the definition (6):

$$\begin{eqnarray*} &&P_B(q,v) = Z_{\rm 2D} - q Z_{\rm 0D} . \end{eqnarray*}$$

The end result of the transfer process is a linear combination of connectivity states involving two time slices, such as those shown in figure 5, with coefficients that are polynomials in n_loop and x. Each state is described by its number in the canonical ordering, from which the coding (i.e., the integers c_i shown in figure 5) can be inferred from the bijection of section 3.5.4. Also, from this coding, one can rather straightforwardly construct a representation of the pairing of the 4n points (namely 2n on each time slice) induced by the loops (shown in red in figure 5) which allows, in particular, 'travel' along the loops.

The goal is now to identify, for each connectivity state, the top and bottom time slices such that the ith point on the top time slice gets glued to the ith point on the bottom time slice. Within this gluing procedure we should, on one hand, be able to distinguish which states contribute to Z_2D, Z_1D and Z_0D (those of Z_1D are then discarded, since they do not contribute to P_B(q, v)) and, on the other hand, provide some extra powers of $n_{\rm loop} = \sqrt{q}$ that have not been accounted for by the transfer process.

To count the number of loops P and, at the same time, determine whether there is a loop of non-trivial homotopy, we begin by travelling along each loop. To this end, one starts at some initial reference point and follows its arc or string to the partner point. Then one jumps to the opposite time slice (since the two are glued) and repeats the process until one comes back to the initial point. During this travel, a list of the points visited is maintained, as well as the horizontal and vertical winding numbers, w_x and w_y, incurred. Note that w_x changes as the result of an arc or string explicitly crossing the periodic boundary condition, whereas w_y changes when one jumps from one time slice to the other. If (w_x, w_y) ≠ (0, 0) the loop has non-trivial homotopy, and we are dealing with a state that contributes to Z_1D and therefore can be discarded. Note that all non-trivial loops (if any) necessarily have the same homotopy, i.e., the same values of (w_x, w_y) up to a global sign change.

Having completed the travel along the first loop, if a non-visited point still exists, this is taken as the new reference point, and we trace out the next loop. This process terminates when all points have been visited. We now know the number of loops P.

For example, the states in figures 5(b) and (c) both have a loop with (w_x, w_y) = (1, 0). This is most easily seen by considering the loop in figure 5(b) (resp. figure 5(c)) that passes through the third (resp. second) point from the left on the bottom time slice.

Suppose now that we have found that all P loops in the state have trivial homotopy. This is the case for the state in figure 5(a), which we shall henceforth use as an example. We then need to find out whether the state contributes to Z_2D or to Z_0D. To this end we use a variant of the Euler relation, i.e., we compute the quantity

$$\begin{equation} \chi = E + 2 C - V - P , \end{equation} \tag{ 17 }$$

where the quantities E, C and V will be defined below.

Corresponding to a loop configuration (red lines in figure 5) there is a corresponding cluster configuration (grey shading in figure 5). The cluster configuration can be seen as a hypergraph on the set $\mathcal{P}_1$ of 2n points (namely n on each time slice) situated to the right of the (loop) point i and to the left of point i + 1, for all even i. An area with grey shading containing d + 1 points of the hypergraph is called a hyperedge of degree d. (Note that in the definition of hyperedges we do not impose the identification of the top and bottom time slices.) Let now E be the sum of the degrees d of all hyperedges. One may think of E as the equivalent number of usual edges (not hyperedges). For instance, the state in figure 5(a) can be represented as

It has three hyperedges: e₁ and e₂ each of degree 1, and e₃ of degree 3. Therefore E = 1 + 1 + 3 = 5.

Let us provide a few details on how to construct the cluster configuration from the loop configuration. The connections between (cluster) points within a single time slice (top or bottom) can be easily inferred by travelling along the arcs on that time slice. It is more delicate to infer the connections (if any) from one time slice to the other. It is obvious that if there are s = 2k > 0 strings, there will be k such connections, and these can again be easily inferred by travelling along the strings (this is the case in figure 5(a)). However, if s = 0 the existence of connections between the clusters on the top and bottom time slices depends on whether the two reduced states are both closed or both open (see section 3.5.1). An example involving open states is provided by the following figure:

To determine in general the connections (if any) between the clusters on the top and bottom time slices we proceed as follows. Consider first the set $\mathcal{P} = \mathcal{P}_1 \cup \mathcal{P}_2$ consisting of 2n points on each time slice, which is the union of the points $\mathcal{P}_1$ on which the clusters live—i.e., those with an even loop point on the left and an odd loop point on the right, shown as solid circles in (19)—and the points $\mathcal{P}_2$ on which dual clusters live—i.e., those with an odd loop point on the left and an even loop point on the right, shown as open circles in (19). We now define a set of integer 'heights' on $\mathcal{P}$ as follows. Starting from an arbitrary initial value (chosen as 1 in (19)), and moving along the top (resp. bottom) time slice from left to right (resp. from right to left), let the height increase (resp. decrease) by one unit each time one crosses a loop opening (resp. closing), i.e., a loop point with code c_i = 1 (resp. c_i = 2). Since only height differences are defined, this procedure defines the heights only up to a global translation. The heights corresponding to the example (19) are shown next to each point in $\mathcal{P}$. It is easy to see that if, for each of the time slices taken separately, the minimum of this height profile (which is 0 in (19) for both time slices) resides at a point of $\mathcal{P}_1$ (resp. $\mathcal{P}_2$), the corresponding reduced state is open (resp. closed). In the case of a pair of open reduced states, all the points of $\mathcal{P}_1$ residing at the minimum height on the top time slice are incident on the same hyperedge as the corresponding points of minimum height on the bottom time slice. In the example (19) there are two such points on the top time slice and one on the bottom time slice, so top and bottom are connected through a hyperedge of degree 2. This concludes the construction of the cluster configuration from the loop configuration.

We finally define V and C as, respectively, the number of vertices and clusters in the hypergraph. Both of these numbers are defined after the identification of the top and bottom time slices. In particular V = n. Both of the states (18) and (19) turn out to have C = 1. We can then compute χ from (17), and we find χ = 5 + 2 − 4 − 3 = 0 for (18) and χ = 4 + 2 − 4 − 1 = 1 for (19).

In general the possible values are χ = 0, 1, 2. When χ = 0 the state belongs to the Z_0D class, and when χ = 1 or 2 it belongs to the Z_2D class. Moreover, when χ = 1 (resp. χ = 2) we can deduce that the state contains s = 2k > 0 strings (resp. s = 0 strings), but we shall not need this fact to compute (6).

Now that the nature (Z_0D, Z_1D or Z_2D) of each connectivity state has been determined, it remains only to multiply it by a certain power of n_loop that has not been accounted for by the transfer process itself. First, there is a factor of $\sqrt{q} = n_{\rm loop}$ for each vertex in the basis B, coming from the front factor of (8). The number of vertices should of course be computed up to the identification which is made by imposing the doubly periodic boundary conditions on B (i.e., gluing the left and right, and the top and bottom). Second, each state has to be multiplied by $n_{\rm loop}^P$, where we recall that P is the number of loops in the final state. Finally, the Z_2D configurations should be multiplied by a factor of $q = n_{\rm loop}^2$; this follows from the Euler relation.

A square basis for the triangular lattice is obtained by placing the following $\check{\sf R}$-matrix:

$$\begin{equation} \check{\sf R}_i = {\sf H}_{i+1} {\sf V}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf H}_{i+1} \end{equation} \tag{ 22 }$$

inside each grey square in figure 2. In addition we need horizontal diagonals on all the white squares. The resulting representation is shown in figure 6. There are four vertices and six edges per grey square.

Zoom In Zoom Out Reset image size

The critical polynomial P_B(q, v) invariably factorizes for any size n of the basis, shedding the small factor

$$\begin{equation} P_{\rm tri}(q,v) = v^3 + 3 v^2 - q . \end{equation} \tag{ 23 }$$

This is compatible with the fact [37] that the triangular-lattice Potts model is exactly solvable on the curve P_tri(q, v) = 0. In particular, for q = 1 we have the root v = −1 + 2cos (2π/9), meaning that the exact percolation threshold is

$$\begin{equation} p_{\rm c} = \frac{v}{1+v} = 2 \sin \left( \frac{\pi }{18} \right) . \end{equation} \tag{ 24 }$$

The remaining, large factor in P_B(q, v) gives additional information about the critical manifold in the regime v < 0, as we shall now see. Its roots in the real (q, v) plane are shown in figure 7. The lower boundary of the BK phase, denoted v₋(q), is the lower branch of the cubic (23). The corresponding upper boundary v₊(q) can be seen in figure 7 as the curve starting from the origin with nearly horizontal slope, passing through the top point of a series of vertical rays, and extending towards the special point (q, v) = (4, −2). Curiously, the critical polynomials miss the part of v₊(q) with 2 < q < 3. Heuristically, this is 'because' the polynomials have to trace out both v₊(q), v₋(q) and the vertical rays in a zigzag fashion that becomes increasingly complicated upon approaching q = 4. The vertical rays corresponding to k = 4, 6, 8, 10 in (21)—and to a lesser extent k = 12, 14 as well—are clearly visible from the figure. Here and in the following, we help the visual identification of such vertical rays by superimposing a number of dotted grey lines on the figures.

Zoom In Zoom Out Reset image size

In conclusion, figure 7 provides rather compelling evidence that the curves v_±(q) will merge in (q, v) = (4, −2). In particular q_c = 4 for the triangular lattice. It follows almost inevitably that the critical curve to which the BK phase is RG-attracted must be the middle branch of (23). That conclusion is backed up by a number of other studies [40–43].

We should mention that the triangular-lattice Potts model is exactly solvable on the chromatic line v = −1 [44–46]. It follows from [45, 46] that v₊(q) = −1 for some value q = 3.819 671 731... obtained by equating two infinite products. The region near (q, v) = (4, −1) exhibits some rather complicated physics and would be a suitable subject for a separate study [47].

The square lattice can obviously be obtained by removing the diagonal edges from the triangular lattice. The $\check{\sf R}$-matrix then reads

$$\begin{equation} \check{\sf R}_i = {\sf H}_{i+1} {\sf V}_i {\sf V}_{i+2} {\sf H}_{i+1} . \end{equation} \tag{ 25 }$$

The basis then looks like figure 2 with the horizontal edges removed. There are two vertices and four edges per grey square.

For any size n, the critical polynomial P_B(q, v) factorizes, shedding two small factors:

$$\begin{equation} P_{\rm sq}(q,v) = (v^2 - q)(v^2 + 4 v + q) . \end{equation} \tag{ 26 }$$

The zero set of the first factor describes the self-dual critical point of the square-lattice Potts model [5], while the zero set of the second factor yields two mutually dual antiferromagnetic critical points [38]. The model is exactly solvable on these curves [5, 38, 43, 48, 49]. In particular, we find that P_sq(1, 1) = 0, so the exact percolation threshold is

$$\begin{equation} p_{\rm c} = \textstyle\frac{1}{2} . \end{equation} \tag{ 27 }$$

The roots of P_B(q, v) in the real (q, v) plane are shown in figure 8; the curves with n ⩽ 4 were already reported in [21]. Of all the Archimedean lattices, the square lattice is the only one for which this phase diagram can be claimed to be completely understood. The boundaries of the BK phase are given by the second factor in (26), namely

$$\begin{equation} v_\pm (q) = -2 \pm \sqrt{4-q} , \end{equation} \tag{ 28 }$$

and in particular q_c = 4. In the thermodynamical limit, n → ∞, we expect an infinite set of vertical rays, corresponding to (21) with k = 4, 6, 8, .... The first few, with k = 4, 6, 8, are clearly visible in our results for n ⩽ 5, shown in figure 8, as is the precursor of the k = 10 ray (which is not yet in its correct position).

Zoom In Zoom Out Reset image size

The square and triangular lattices (and the hexagonal lattice, which is the dual of the triangular one) could have been presented in three-terminal form. This fact actually makes it possible to compute the critical manifolds, P_sq(q, v) = 0 and P_tri(q, v) = 0, exactly [11], and it is closely related to the exact solvability [5, 37, 38, 43, 48] of the models along these curves.

However, the kagome lattices—and indeed all the remaining Archimedean lattices—are not of the three-terminal type. Accordingly no exact solution is known to this day. The graph polynomial method therefore gives only approximate results, which are however very accurate, in particular in the ferromagnetic region v > 0.

The $\check{\sf R}$-matrix of the kagome lattice can be written as

$$\begin{equation} \check{\sf R}_i = {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_{i} {\sf E}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} . \end{equation} \tag{ 29 }$$

The resulting representation is shown in figure 9. There are three vertices and six edges per grey square.

Zoom In Zoom Out Reset image size

We have computed the (unique) root P_B(q, v) in the ferromagnetic regime, v > 0, for several integer values of q. The results for q = 1 are shown in table 4 to 50-digit numerical precision, in terms of the percolation probability $p = \frac{v}{1+v}$, for square bases of size 1 ⩽ n ⩽ 7. A quick glance at the table makes it obvious that these numbers converge very fast to their expected limit p_c.

The roots for 1 ⩽ n ⩽ 4 have already been reported in [21], leading the authors to propose a final estimate of p_c = 0.524 405 00(1). However, the fact that we now have three more terms in the sequence allows us to employ powerful extrapolation techniques to obtain a very accurate final value of p_c. We have chosen to apply the time-proven Bulirsch–Stoer (BS) extrapolation [50]. This algorithm requires a parameter w which can be thought of as a correction-to-scaling exponent. In an preliminary step we obtain an approximate value for w from a nonlinear fit of the data in table 4 to the form

$$\begin{equation} p_{\rm c}(n) = p_{\rm c} + A n^{-w} . \end{equation} \tag{ 30 }$$

The best results are obtained by constraining this fit to the last three available data points. We do not report the values of w found for each data set to be considered in this paper, except for the first few examples of each problem type. We shall however provide a few general remarks in the Discussion section 8. Generally speaking we find w ≈ 6 for the best behaved (bond or site) percolation problems. Obviously, the higher the value of w is, the better the precision on the final result will be.

In a second step, we insert the value found for w into the implementation of the BS algorithm described by Monroe [51]. This results in a series of approximants that are compared among themselves in order to assess a final value and error bar. We crosscheck our results by repeating the whole procedure with the last data point being eliminated, in order to ensure that the central value and error bar obtained from fits on N − 1 points are compatible with (although of course less precise than) those obtained on all N data points.

For some of the lattices for which fewer data points (i.e., sizes n) are available, some adaptations of this general procedure will be necessary. We shall return to this in the following subsections.

In table 4 and the following many tables in this paper, we compare our final result with the most accurate value known from previous numerical work. In the present case, we find w ≈ 6.36, and the relative precision on the final value of p_c is of the order of 4 × 10⁻¹¹, that is, four orders of magnitude better than the previous result [25].

We next turn to the Ising model (q = 2). In this case, all the P_B(q, v) are found to factorize into small factors. The maximum degree of the factors is d_max = 4 for n = 1, 2; d_max = 8 for n = 3, 4; and d_max = 16 for n = 5. There is precisely one of these factors, namely

$$\begin{equation} -8 - 8 v + 4 v^3 + v^4 , \end{equation} \tag{ 31 }$$

that possesses a positive root,

$$\begin{equation} v_{\rm c} = \sqrt{3 + 2 \sqrt{3}} - 1 \simeq 1.542\,459\,756\ldots . \end{equation} \tag{ 32 }$$

Moreover, (31) is a factor in P_B(q, v) for any size n. Its physical root (32) coincides with the exactly known critical point of the kagome-lattice Ising model [23, 52]. Below we shall similarly see that we recover exact results for the Ising model on any lattice.

The polynomials for the Ising model are often found to simplify under the change of variables $v = -1 + \sqrt{y}$. Recall that v = e^K − 1, but, rewriting the nearest neighbour interaction energy in Ising form, $K \delta _{\sigma _i,\sigma _j} = K_{\rm Ising} (S_i S_j + 1)$ for spins S_i = ±1, we find K = 2K_Ising, so $y = {\rm e}^{K_{\rm Ising}}$ is simply the Boltzmann factor for a pair of aligned Ising spins. In the present case (31) simplifies to

$$\begin{equation} -3 - 6t + t^2 . \end{equation} \tag{ 33 }$$

The critical points for the q = 3-state and q = 4-state Potts models are shown in tables 5–6. The exponent appearing in (30) is w ≈ 5.36 for q = 3, and w ≈ 4.80 for q = 4. In any case, using BS extrapolations as explained above, we arrive at values for v_c which are considerably more accurate than previous numerical results. Note in particular that numerical simulations of the Monte Carlo or transfer matrix type are usually particularly difficult for q = 4 because of the presence of logarithmic corrections to the scaling. By contrast, the graph polynomial method only experiences a slight decrease of w, and the precision is almost as good as for percolation. Thus, for q = 4 our final value for v_c is six orders of magnitude more precise than the previous result [24].

The phase diagram on the kagome lattice has been discussed in detail in [13], and further in [14] on the basis of the square-basis P_B(q, v) with n = 1, 2, 3, 4. In figure 10 we show again the roots of P_B(q, v) in the real (q, v) plane, but this time with the n = 5 basis included. Because of the extensive treatment of this phase diagram in [13, 14] we shall be rather brief.

Zoom In Zoom Out Reset image size

The BK phase contains vertical rays at q = B₄ = 2 and q = B₆ = 3. Unlike for the triangular and square lattices, there is no sign of the BK phase widening out towards q = 4 as n increases. Its rightmost termination might be close to the n = 3 arc extending to around q ≈ 3.2. This arc is confirmed by results from the hexagonal bases studied in [14]. On the other hand, even allowing for n mod 2 parity effects which are visible elsewhere in the phase diagram, it is curious that this arc is not confirmed by the n = 5 critical polynomial. So it might also be that the BK phase in fact terminates right at the q = 3 vertical ray. In any case, it seems certain that q_c < 4 for the kagome lattice.

The upper boundary v₊(q) of the BK phase is the nearly straight line emanating from the origin and passing through the point (q, v) = (3, −1). Indeed, the three-state zero-temperature antiferromagnet on the kagome lattice is equivalent to the corresponding 4-state model on the triangular lattice, which in turn is known to be critical with central charge c = 2 [53]. It is interesting to observe that at finite n the curve v₊(q) is approximated by various pieces from the different P_B(q, v), and that no single critical polynomial reproduces the curve completely. This is in line with observations already made in [13, 14]. In particular we see clear n mod 2 parity effects: the critical polynomials with even n (resp. odd n) are the only ones to produce the part of v₊(q) with 0 < q ≲ 1.6 and 2 < q < 3 (resp. 1.6 ≲ q < 2).

Meanwhile, the lower boundary v₋(q) of the BK phase is the nearly straight line emanating from (q, v) = (0, −3) and passing through ≈(3, −2). It is determined by the polynomials with even n. Interestingly there is a lower-lying curve, coming out from (0, −3) with infinite slope, and the space between this curve and v₋(q) is devoid of vertical rays (because it does not belong to the BK phase).

Finally, the existence of two small enclosed regions, or phases—the first a thin sliver between (2, −2) and (2, −1), and the other a triangular-shaped region above (2, −1)—is confirmed by the new n = 5 polynomial.

A four-terminal representation of the four–eight lattice is shown in figure 11. The corresponding $\check{\sf R}$-matrix reads

$$\begin{equation} \check{\sf R}_i = {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} . \end{equation} \tag{ 34 }$$

There are four vertices and six edges per grey square.

Zoom In Zoom Out Reset image size

The approximations to the bond percolation threshold p_c are given in table 7. Since also for this lattice any parity of n is possible, the BS extrapolation scheme produces very accurate values, improving considerably on the existing numerical results which are shown for comparison. We note that the exponent appearing in (30) comes out as w ≈ 4.28 in this case, so it is definitely lattice dependent.

In the Ising case (q = 2) the polynomials P_B(q, v) systematically factorize. The maximum degree of the factors is d_max = 4 for n = 1, 2; d_max = 8 for n = 3; d_max = 6 for n = 4; and d_max = 16 for n = 5. There is precisely one of these factors, namely

$$\begin{equation} -4 - 8 v - 6 v^2 + v^4 , \end{equation} \tag{ 35 }$$

that possesses a positive root,

$$\begin{equation} v_{\rm c} = \frac{1 + \sqrt{5 + 4 \sqrt{2}}}{\sqrt{2}} \simeq 3.015\,445\,388\ldots , \end{equation} \tag{ 36 }$$

and this factor occurs in P_B(q, v) for any size n. Its positive root (36) produces the exactly known critical point of the Ising model on the four–eight lattice [23, 55].

The q = 3 and q = 4 critical points v_c are shown in tables 8–9, and just like in the percolation case the BS extrapolation produces very accurate final values.

Reference [14] already contained a discussion of the phase diagram for the four–eight lattice. But to highlight the new n = 5 results, figure 12 shows again the roots of P_B(q, v) in the real (q, v) plane.

Zoom In Zoom Out Reset image size

The boundaries of the BK phase can be clearly seen. The upper boundary v₊(q) is the almost straight line emanating from the origin and extending out towards ≈(4, −2). The lower boundary starts at ≈(0, −4). The zero sets of the various critical polynomials fill in the curves v_±(q) in a zigzag fashion, while at the same time providing vertical rays at the Beraha numbers (21). The rays with k = 4, 6, 8, 10, 12 are clearly visible in the figure. It is interesting to notice the formation of narrow 'fingers' that tend to close those parts of the BK boundaries that are not provided by the principal 'zigzag' trend. For example, the lower boundary with 0 < q < 2 is produced by fingers in the critical polynomials with n ⩾ 3.

Overall, it seems likely that the BK phase for this lattice will extend all the way out to q = 4, and so we can conjecture that q_c = 4 for the four–eight lattice.

The frieze lattice is the first example of a lattice which cannot be represented in four-terminal form by using the same $\check{\sf R}$-matrix in all grey squares (x, y). Instead we have

$$\begin{equation} \check{\sf R}_i = \left\lbrace \begin{array}{@{}ll@{}}\sf H_{i+1} {\sf V}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf H}_{i+1} & {{\rm for} \quad x+y\; {\rm even}} \\ {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} & {{\rm for} \quad x+y\; {\rm odd}.} \\ \end{array} \right. \end{equation} \tag{ 37 }$$

In addition, there are horizontal diagonals on the white squares with coordinates $(x+\frac{1}{2},y+\frac{1}{2})$ for x + y even. The resulting basis, shown in figure 13, needs n to be even in order for it to produce the frieze lattice upon tiling. It has two vertices and five edges per grey square.

Zoom In Zoom Out Reset image size

The restriction that n needs to be even is somewhat problematic for our approach. We have now only three data points (n = 2, 4, 6) for the extrapolations, instead of the usual seven (n = 1, 2, ..., 7) when no parity constraint is operative. Fitting first the three data points to the form (30) we find that w ≈ 3.03. This three-point fit also provides the central value shown in table 10. Next, on the basis of other lattices for which more data points are available, we can estimate that this value of w is likely to shift by around 1/10 upon increasing the size. Performing next a two-point fit to the two largest sizes, with a fixed value of w within the range 3.03 ± 0.10, we get some alternative extrapolated values from which we can judge the size of the error bar on the central value.

Obviously this procedure leads to final results which are less precise than those of the preceding sections; but they still have a precision superior to that of existing numerical results.

The approximations for the percolation threshold p_c are shown in table 10, and those for the critical points v_c for the q = 3 and q = 4 Potts models are given in tables 11–12.

For the Ising model (q = 2) we find the now-familiar factorization of the polynomials P_B(q, v). The maximum degree of the factors is d_max = 2 for n = 2 and d_max = 12 for n = 4. But the recurrent factor that determines the ferromagnetic critical point is simply

$$\begin{equation} -1 + v \end{equation} \tag{ 38 }$$

and so

$$\begin{equation} v_{\rm c} = 1 \end{equation} \tag{ 39 }$$

in agreement with [23].

The information obtainable on the phase diagram also suffers from the parity restriction, since the polynomials P_B(q, v) are now at our disposal only for n = 2, 4. It is nevertheless clear from figure 14 that the usual features are present. The extent of the BK phase can be judged from the vertical rays at the Beraha numbers with k = 4, 6, 8. Inside the BK phase we have a critical curve coming out of the origin with a vertical tangent and extending to the point (q, v) = (4, −2). Since this is the analytical continuation of the critical curve in the ferromagnetic regime v > 0, it must be the RG attractor governing the BK phase. We should therefore have q_c = 4 for the frieze lattice.

Zoom In Zoom Out Reset image size

The curves v_±(q) bounding the BK phase are only partially represented by the roots of the polynomials P_B(q, v). Indeed the parts with 2 < q < 3 are missing altogether. We notice however that the n = 4 curve develops a small bulge near the bottom of the q = 3 vertical ray, and it is conceivable that for higher n this will develop into narrow 'fingers' that will close the curves v_±(q)—a situation that was seen clearly in figure 12 for the four–eight lattice.

With the three–twelve lattice we are back in the category of lattices that can be represented in four-terminal form for any parity of n. The appropriate $\check{\sf R}$-matrix reads

$$\begin{equation} \check{\sf R} = {\sf H}_{i+1} {\sf V}_i {\sf H}_{i+1} {\sf V}_i {\sf E}_{i+2} {\sf H}_{i+1} {\sf V}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_i . \end{equation} \tag{ 40 }$$

The corresponding basis is depicted in figure 15. It possesses six vertices and nine edges per grey square, the highest numbers that we have attained among the Archimedean lattices. Accordingly we can expect the most accurate results for the critical points of this lattice (except, obviously, for the exactly solvable three-terminal cases).

Zoom In Zoom Out Reset image size

Results for the percolation threshold p_c are given in table 13, using square bases of size n ⩽ 7. Those with n ⩽ 4 have already been reported in [21] where the final estimate p_c = 0.740 420 800(1) was proposed. The three extra data points allow us to improve considerably on this, adding four more digits to the final result reported in table 13. We have again benefited from a very high value w ≈ 6.39 of the parameter appearing in (30). This brings the relative precision to 1 × 10⁻¹², that is, more than four orders of magnitude better than the best available numerical result [22]. Note that the value of w is almost coincident with that of the kagome lattice; this is presumably linked to the fact that this lattice has the same symmetry group as the three–twelve lattice.

The graph polynomials P_B(q, v) factorize as usually in the q = 2 Ising case. The maximum degree of the factors is d_max = 4 for n = 1, 2; d_max = 8 for n = 3, 4; and d_max = 16 for n = 5. There is precisely one of these factors, namely

$$\begin{equation} -8 - 8 v - 6 v^2 - 2 v^3 + v^4 , \end{equation} \tag{ 41 }$$

that possesses a positive root,

$$\begin{equation} v_{\rm c} = \textstyle\frac{1}{2} (1 + \sqrt{3} + \sqrt{2 (6 + 5 \sqrt{3})} ) \simeq 4.073\,446\,135\ldots . \end{equation} \tag{ 42 }$$

Moreover, (41) is a factor in P_B(q, v) for any size n. Its positive root (42) furnishes the exactly known critical point [56–58].

The results for the critical point v_c of the Potts models with q = 3 and q = 4 are shown in tables 14–15. Again we obtain very high accuracies on the final estimates.

The phase diagram shown in figure 15 has a very intricate structure, containing more features than for any of the other lattices considered so far. It was discussed in detail in [14] for n ⩽ 4, but we gain extra information—and confirmation of some of the salient features—from the n = 5 critical polynomial P_B(q, v) now available.

We see of course the usual vertical rays characterizing the BK phase, visible at the Beraha numbers (21) with k = 4, 6, 8 and building up at k = 10 as well. The extent of the BK phase can be estimated from those rays. It seems clear that it extends out to the point (q, v) = (4, −2), and thus we have q_c = 4 for the three–twelve lattice.

To discuss the phase diagram (see figure 16) in some more detail, we first focus on the region 0 ⩽ q ≲ 2. We first note the formation of fingers in the n = 4, 5 curves in the range 0 ⩽ q ≲ 1.2. To the left of the lower edge of the q = 2 vertical ray one sees the formation of a triangular-shaped enclosed region with 1.3 ≲ q ≲ 2; this is visible in particular in the small wrinkle that develops in the n = 5 curve. Moreover, there is a tiny enclosed sliver with q ⩽ 2 ≲ 2.01 that is consistently visible for any parity of n. These two regions (of triangular and sliver shape) near q = 2 are reminiscent of similar features for the kagome lattice, and we believe that they subsist in the thermodynamical limit.

Zoom In Zoom Out Reset image size

For 2.8 ≲ q ⩽ 4 the situation is quite complicated. There are notably a couple of curves inside the BK phase, such that the vertical rays with k = 6, 8 are cut into three pieces. One observes just above (q, v) = (4, −2) several extra curves developing for large n which clearly tend to extend this structure to higher values of k. A similar situation happens near the lower boundary of the BK phase. Moreover, in the region 3 ≲ q ≲ 4 and −4 ≲ v ≲ −3 many curves are building up and accumulating at (q, v) = (3, −3). Clearly this latter point must have a very distinguished role in the phase diagram.

To cast the cross lattice in four-terminal form, grey squares of four different types are needed. The $\check{\sf R}$-matrix takes the form

$$\begin{equation} \check{\sf R}_i = \left\lbrace \begin{array}{@{}ll@{}}\sf H_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} {\sf V}_{i+2} & {{\rm for}\quad x \;{\rm even\; and}\; y\; {\rm even}} \\ {\sf H}_{i+1} {\sf E}_{i+2} {\sf E}_i {\sf H}_{i+1} & {{\rm for}\quad x\; {\rm odd\; and} y\; {\rm even}} \\ {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_i & {{\rm for}\quad x\; {\rm even\; and} y\; {\rm odd}} \\ {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} {\sf V}_i & {{\rm for}\quad x\; {\rm odd\; and} y\; {\rm odd}} \\ \end{array} \right. \end{equation} \tag{ 43 }$$

and the resulting basis is illustrated in figure 17. The reader is invited to check carefully from the figure that each vertex is indeed surrounded by a square, a hexagon and a dodecagon. Obviously we must require that n be even. There are on average three vertices and 9/2 edges per grey square.

Zoom In Zoom Out Reset image size

The approximations to the percolation threshold p_c are shown in table 16. Despite the parity constraint on n, we are still able to provide a final estimate that is more precise than those obtainable from numerical simulations [54].

For q = 2, the graph polynomials P_B(q, v) factorize, and the maximum degree of the factors is d_max = 8 for both n = 2 and n = 4. The unique common factor possessing a positive root is

$$\begin{equation} {-}32 - 128 v - 208 v^2 - 176 v^3 - 84 v^4 - 24 v^5 + 4 v^7 + v^8 . \end{equation} \tag{ 44 }$$

After some algebra this leads to the critical coupling

$$\begin{eqnarray} v_{\rm c} &=& \exp \left(2 {\rm arctanh } \sqrt{\frac{2}{1 + \sqrt{3} + \sqrt{-4 + 6 \sqrt{3}}}} \right) - 1 \\ &\simeq & 3.216\,563\,123\ldots \, \nonumber \end{eqnarray} \tag{ 45 }$$

in agreement with [23].

The results for the critical points v_c for the 3-state and 4-state Potts models appear in tables 17–18.

The phase diagram shown in figure 18—as witnessed by the roots of the critical polynomials P_B(q, v)—combines a number of features familiar from the four–eight lattice. There is a finger in the n = 4 curve with 0 ⩽ q ≲ 1.7 that tends to bridge the gap between the point (q, v) = (0, −4) and the lower edge of the q = 2 vertical ray. Similarly, wrinkles in both curves (n = 2, 4) tends to delimit the BK phase from above on the interval 2 < q < 3. As a result, we observe vertical rays at the Beraha numbers (21) with k = 4, 6, 8. Presumably the BK phase extends out to q_c = 4 for this lattice.

Zoom In Zoom Out Reset image size

The four-terminal representation of the snub square lattice is quite straightforward to obtain. The $\check{\sf R}$-matrix reads

$$\begin{equation} \check{\sf R} = \left\lbrace \begin{array}{@{}ll@{}}\sf H_{i+1} {\sf V}_{i+2} {\sf H}_{i+1} {\sf V}_i {\sf H}_{i+1} & {{\rm for}\quad x+y\; {\rm even}} \\ {\sf H}_{i+1} {\sf V}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf H}_{i+1} & {{\rm for}\quad x+y\; {\rm odd}} \\ \end{array} \right. \end{equation} \tag{ 46 }$$

and figure 19 contains a rendering of the corresponding basis. There are two vertices and five edges per grey square.

Zoom In Zoom Out Reset image size

The bond percolation thresholds p_c obtained from the critical polynomials are reported in table 19, and the corresponding results for the critical point v_c of the Potts model with q = 3 and q = 4 are shown in tables 20–21.

The Ising model polynomials P_B(q, v) with q = 2 factorize, the maximum degree of the factors being d_max = 6 for n = 2, and d_max = 12 for n = 4. There is a unique common factor possessing a positive root

$$\begin{equation} -4 - 8 v - 4 v^2 + 4 v^3 + 8 v^4 + 4 v^5 + v^6 . \end{equation} \tag{ 47 }$$

If we define $\omega = 37 + 27 \sqrt{2} + 3 \sqrt{315 + 222 \sqrt{2}}$, the critical coupling can be written as

$$\begin{equation} v_{\rm c} = \textstyle\frac{1}{3} ( -2 -2 w^{-1/3} + w^{1/3} ) \simeq 0.980\,730\,864\ldots \end{equation} \tag{ 48 }$$

and this coincides with the result of [23].

It remains to discuss the phase diagram, shown in figure 20. As usual, the extent of the BK phase can be seen from the vertical rays, here appearing at Beraha numbers (21) with k = 4, 6, 8. An interesting feature for this lattice is that there is a transition curve in the middle of the BK phase, reminiscent of what was found for the frieze lattice (see figure 14). This curve emanates from the origin with vertical slope and extends out to (q, v) = (4, −2). For finite n there are some gaps in this curve, here visible for n = 4 near k = 8, but these should disappear in the thermodynamical limit. The closure of the BK phase is ensured by a wrinkle developing near the lower edge of the q = 3 vertical ray in the n = 4 curve, and by fingers developing at q ≈ 4.

Zoom In Zoom Out Reset image size

The snub hexagonal lattice is a depleted version of the triangular lattice, where 1/7 of the vertices and their adjacent edges have been erased. We have constructed this lattice in two different ways.

The first construction closely parallels the one that we have used in section 4.1 for the triangular lattice. The only difference is that the ${\sf H}$ and ${\sf V}$ operators corresponding to erased edges have been replaced by just one of their two terms (the identity ${\sf Id}$ or the Temperley–Lieb generator ${\sf E}$, as the case may be). There are then $\frac{12}{7}$ vertices and $\frac{30}{7}$ edges per grey square; note that these numbers are compatible with the coordination number being 5.

This depletion representation is a bit easier to state than to actually write down. But in formal terms we arrive, after drawing things carefully, at the following $\check{\sf R}$-matrix:

$$\begin{equation} \check{\sf R} = \left\lbrace \begin{array}{@{}ll@{}}\sf H_{i+1} {\sf V}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf H}_{i+1} & {{\rm for}\quad 5 y + x = 0,2,4 \;{\rm mod}\; 7} \\ {\sf H}_{i+1} {\sf V}_i {\sf H}_{i+1} {\sf E}_{i+2} {\sf Id}_{i+1} & {{\rm for}\quad 5 y + x = 1\; {\rm mod}\; 7} \\ {\sf Id}_{i+1} {\sf E}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf H}_{i+1} & {{\rm for}\quad 5 y + x = 3\; {\rm mod}\; 7} \\ {\sf Id}_{i+1} {\sf V}_i {\sf Id}_{i+1} {\sf E}_{i+2} {\sf H}_{i+1} & {{\rm for}\quad 5 y + x = 5\; {\rm mod}\; 7} \\ {\sf H}_{i+1} {\sf E}_i {\sf Id}_{i+1} {\sf V}_{i+2} {\sf Id}_{i+1} & {{\rm for}\quad 5 y + x = 6\; {\rm mod}\; 7} \\ \end{array} \right. . \end{equation} \tag{ 49 }$$

In addition there are horizontal diagonals in the white squares at coordinates $(x+\frac{1}{2},y+\frac{1}{2})$ when 5y + x = 0, 1, 4, 5, 6 mod 7. Note that the first line in (49) corresponds to undepleted grey squares, i.e., it coincides with (22). Subsequent lines are obtained from the first one by depletion, i.e., replacing some of the ${\sf H}$ operators by the identity ${\sf Id}$, and some of the ${\sf V}$ operators by the Temperley–Lieb generator ${\sf E}$.

Note that this construction will just disconnect the erased spins from the remainder of the lattice. This means that when computing the critical polynomial from (49), the true P_B(q, v) will be multiplied by a spurious factor of q per erased spin. This will obviously not change the set of roots P_B(q, v) = 0. We adopt the convention (here and elsewhere) of dividing such spurious factors out of the polynomials that are provided in electronic form in the supplementary material available at stacks.iop.org/JPhysA/47/135001/mmedia.

Note that (49) only makes sense when n is a multiple of 7. This means that, using this construction, the only computation that we can perform in practice is to find numerically the roots in v of P_B(q, v) = 0 with n = 7. While this agrees nicely with the upper bound n_max = 7 for the feasibility of the computations, it entails two inconveniences. First, our inability to compute the full polynomial P_B(q, v) with n = 7 implies that we have no access to the phase diagram in the real (q, v) plane. Second, the fact that we get just one single estimate for v_c means that the only sensible proposal for the n → ∞ result is the n = 7 value itself. In particular, an error bar on this result can only be obtained by making the (a priori not unlikely) assumption that the distance between the n = 7 value and the would-be extrapolation is comparable to that of other lattices for which we have been able to perform the extrapolation carefully.

To elaborate on this last remark, consider for instance the n = 7 estimate for the percolation threshold shown in table 22. Let us recall that the relative deviation between the n = 7 estimate and the n → ∞ extrapolated value is 2 × 10⁻⁹ for the kagome lattice (see table 4) and 1 × 10⁻¹⁰ for the three–twelve lattice (see table 13). However, these two lattices benefited from very high values (w ≈ 6.36 and w ≈ 6.39 respectively) of the correction-to-scaling exponent in (30). For the four–eight lattice (see table 7) we found a smaller value w ≈ 4.28 and accordingly the relative precision of the n = 7 estimate comes out as 8 × 10⁻⁸. Recall also that the basis for the kagome and four–eight (resp. three–twelve) lattices has six (resp. nine) edges per grey square. From the number of edges in the basis, we thus have no a priori reason to believe that the convergence properties of the snub hexagonal lattice (whose basis has 30/7 ≃ 4.3 edges per grey square) would be significantly worse than any of those lattices. However, since the value of w for the snub hexagonal lattice is currently unknown, we could conservatively assume that the n = 7 data point has a relative precision of only 8 × 10⁻⁷, and use this assumption to provide a tentative error bar:

$$\begin{equation} p_{\rm c} = 0.434\,3281 (4) \qquad {\rm (tentative\; result.)} \end{equation} \tag{ 50 }$$

Noting that even this conservative estimate is in significant disagreement with the numerical result of Parviainen [54] (see table 22) we therefore turn to a different way of constructing the snub hexagonal lattice.

The second construction is based on the four-terminal representation shown in figure 21.⁸ It generalizes the general four-terminal set-up of figure 2 by stretching the white squares vertically, so they are now hexagons. To make these hexagons more visible in the figure we have shaded them alternately using white and pink colours, but we shall still refer to them as 'white hexagons'. The $\check{\sf R}$-matrix describing the grey squares is exactly the same as that of the snub square lattice; see equation (46). The essential characteristic of the white hexagons is that—just like for the white squares used previously—the lattice structure inside them can be constructed without the use of auxiliary spaces. The corresponding operator (the analogue of $\check{\sf R}_i$) reads

$$\begin{equation} {\sf O}_i \equiv {\sf H}_{i+1} {\sf V}_{i+2} {\sf H}_{i+1} {\sf V}_i {\sf H}_{i+1} . \end{equation} \tag{ 51 }$$

After each even row (i.e., $y = \frac{1}{2}\; {\rm mod }\; 2$) the operator ${\sf O}_i$ acts at positions i = 0 mod 4, and after each odd row (i.e., $y = \frac{3}{2}\; {\rm mod }\; 2$) it acts at i = 2 mod 4. Note that $\mathcal{O}_i$ and $\mathcal{O}_{i+4}$ commute, since each operator acts on precisely four adjacent strands, and so we do not need to specify the order of factors in these products. Moreover, since $\mathcal{O}_i$ now contains operators of the type ${\sf V}$ that propagate the system upwards, it is crucial for this representation that every other white hexagon be empty. This is indeed the case, as shown by the alternating white and pink shadings of the 'white' hexagons.

Zoom In Zoom Out Reset image size

The representation of figure 21 makes sense for even n. It contains three vertices and $\frac{15}{2}$ edges per grey square. This is a significant improvement over the first construction described in this subsection. We note that the basis with n = 2 coincides with that used in [20, figure 21(a)] in which only percolation (not the general q-state Potts model) was studied. We have checked that our result for P_B(1, v) with n = 2 is proportional to [20, equation (27)] after the usual change of variables, p = v/(1 + v).

The percolation thresholds p_c using both the first and the second constructions are shown in table 22. We have separated the results using the two different constructions by a horizontal line in this and the following two tables. Extrapolating the results with n = 2, 4, 6 leads to the estimate for the n → ∞ limit shown in table 22. This is in agreement with the tentative result (50), but is obviously more precise since we now take advantage of a well-controlled extrapolation procedure. Note that we find w ≈ 2.94, significantly lower than for the four–eight lattice, so our precautions in arriving at (50) were justified. The final result of table 22 enhances our disagreement with [54] for this lattice.

For the Ising model (q = 2) the polynomials P_B(q, v) factorize as usual. The degree of the largest factor is d_max = 8 for n = 2 and d_max = 12 for n = 4. The recurrent factor that leads to a positive real root is

$$\begin{equation} -8 - 8 v - 4 v^2 + 4 v^3 + 8 v^4 + 4 v^5 + v^6 . \end{equation} \tag{ 52 }$$

If we define $\omega = 37 + 27 \sqrt{3} + 3 \sqrt{6 (66 + 37 \sqrt{3})}$, the critical coupling can be written in the same way as (48), namely

$$\begin{equation} v_{\rm c} = \textstyle\frac{1}{3} ( -2 -2 w^{-1/3} + w^{1/3} ) \simeq 1.050\,155\,297\ldots , \end{equation} \tag{ 53 }$$

and this coincides with the result of [23]. Using the first construction, we have performed a 50-digit numerical computation to verify that this number is also a root of the n = 7 polynomial. This agreement provides compelling evidence that both constructions of P_B(q, v) are correct.

Our results for the critical point v_c in the q = 3-state and q = 4-state Potts models are shown in tables 23–24. Also in those cases we have based the n → ∞ extrapolations on the results coming from the second construction (with n = 2, 4, 6). Using the same kind of reasoning as discussed above for percolation, we observe that the n = 7 results are compatible with these extrapolated values.

The phase diagram of the snub hexagonal lattice is shown in figure 22. The extent of the BK phase can be estimated from the vertical rays at q = B_k with k = 4, 6, 8. Its upper and lower boundaries v_±(q) are interrupted by a hiatus for 2 < q < 3 with these bases, but the wrinkle developing near the bottom of the q = 3 ray in the n = 4 curve indicates that this gap may be partly filled in at larger sizes. The curve coming out of (q, v) = (0, 0) with infinite slope cuts through the BK phase. We also note a rather rich structure near q = 4 where a flower-like structure with four petals grows out of the point (q, v) = (4, −2). There are thus four (resp. eight) branches of the curve coming out that point for n = 2 (resp. n = 4), and this number might well continue to grow for larger n.

Zoom In Zoom Out Reset image size

The ruby lattice is most simply treated by considering the corresponding dual lattice. A four-terminal representation of the basis of the ruby dual lattice is shown in figure 23. The reader can check that it is indeed a quadrangulation, with each of the quadrangles being surrounded by vertices of degrees 3, 4, 6 and 4. The $\check{\sf R}$-matrix reads

$$\begin{equation} \check{\sf R} = \left\lbrace \begin{array}{@{}ll@{}}\sf V_{i} {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} & {{\rm for}\quad x+y\; {\rm even}} \\ {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_{i} {\sf H}_{i+1} {\sf V}_{i+2} & {{\rm for}\quad x+y\; {\rm odd}} \\ \end{array} \right. \end{equation} \tag{ 54 }$$

and there are horizontal diagonals on all the white squares. Clearly this representation of the ruby dual lattice requires n to be even. There are three faces (corresponding to vertices of the ruby lattice itself) and six edges per grey square.

Zoom In Zoom Out Reset image size

The graph polynomials P_B(q, v) for the ruby lattice (and their specialization P_B(p) to the percolation case) are then found from those of the dual by a simple duality transformation. Details on duality will be deferred to section 5.

Bond percolation thresholds obtained from P_B(p) are reported in table 25. Uncharacteristically for the graph polynomial approach the convergence to the thermodynamical limit is here seen to be non-monotonic. This is presumably due to the smallest result (n = 2) straying away from the general trend. Combined with the parity constraint on n, this negatively impacts the precision of the extrapolated threshold p_c for this lattice.

For the Ising case (q = 2), the graph polynomials P_B(q, v) factorize once again, and the maximum degree of the factors is d_max = 6 for n = 2, and d_max = 10 for n = 4. There is a unique common factor possessing a positive root:

$$\begin{equation} -8 - 8 v + 4 v^3 + v^4 . \end{equation} \tag{ 55 }$$

The critical coupling reads

$$\begin{eqnarray} &&v_{\rm c} = \sqrt{3 + 2 \sqrt{3}} -1 \simeq 1.542\,459\,756\ldots \end{eqnarray} \tag{ 56 }$$

in agreement with [23]. Rather curiously, this is identical to the result (32) found for the Ising model on the kagome lattice.

The critical points for the q = 3 and q = 4 models are shown in tables 26–27. As for percolation, the convergence is non-monotonic, preventing us from attaining the usual precision in the final results.

The phase diagram for the ruby lattice is shown in figure 24. The extent of the BK phase can be judged from the vertical rays at the Beraha numbers (21) with k = 4, 6, 8. The upper and lower boundaries of the BK phase are only partially brought out by the largest size n = 4. In particular, the upper boundary between q = 2 and q = 3 is still missing. The lower boundary between q = 0 and q = 2 is partially provided by the finger protruding from q = 0, and one notes the formation of a wrinkle to the right of q = 3. Interestingly there is another curve, emanating from (q, v) = (0, −4), that lies below the lower boundary of the BK phase. This curve extends towards the point (q, v) = (4, −2), and presumably the boundaries of the BK phase will also join at that point. We therefore conjecture that q_c = 4 for this lattice.

Zoom In Zoom Out Reset image size

A remarkable property of the critical polynomial P_B(q, v) is that it usually factorizes in exactly solvable cases [13, 14]. This is true in particular for the exactly solvable lattices (square, triangular and hexagonal). This is also true for the Ising model (q = 2) on any of the lattices considered. Note however that some cases also exist where P_B(q, v) fails to factorize, even though the model is known to be exactly solvable. This is so in particular for the zero-temperature antiferromagnetic 3-state model, (q, v) = (3, −1), on the kagome lattice and for the entire chromatic line, v = −1, on the triangular lattice. Indeed, Baxter has shown that the latter model is integrable [45, 46], and the former model (3-state kagome) is equivalent to a special case of the latter (4-state triangular) by means of an exact mapping [53].

Conversely, it is compelling to consider any systematic factorization of P_B(q, v) as evidence that the model may be exactly solvable (by 'systematic factorization' we mean a factorization that occurs for any value of the size n).

In this section we examine exhaustively the issue of factorization for all the Archimedean lattices. We consider the following cases: integer q = 0, 1, 2, 3, 4, the chromatic polynomial v = −1, the flow polynomial v = −q [59, 60], and the limit (q, v) → (0, 0) with fixed w = v/q, corresponding to spanning forests [61] with weight 1/w per component tree.

For q = 0, P_B(q, v) factorizes for all the lattices, producing a root v = 0. This is consistent with the observation that for all the lattices there is a branch of the critical curve going through the point (q, v) = (0, 0). This describes the problem of spanning trees, which can in turn be related to free (symplectic) fermions with central charge c = −2 [61–63]. Moreover, the triangular, kagome and three–twelve lattices have a root (q, v) = (0, −3). And the square, four–eight, cross and ruby lattices have a root (q, v) = (0, −4). By duality, the models with (q, v) = (0, v_c) are equivalent [61] to models of spanning forests on the corresponding dual lattice with weight v_c per component tree. These models can in turn be formulated as interacting fermionic theories [62, 63], and we conjecture that they are in fact exactly solvable. It follows, still by duality, that the spanning tree problem factorizes on the square lattice with w = −1/4 and on the hexagonal lattice with w = −1/3. We find moreover that spanning trees on the cross lattice factorize with w = −1/3.

The case of the snub square and snub hexagonal lattices is interesting. For q = 0, the critical polynomials of both of these lattices shed the small factor 8 + 5v + v², and we conjecture that the corresponding roots $v=(-5 \pm {\rm i}\sqrt{7})/2$ are loci of exact solvability.

The Ising case (q = 2) has been extensively discussed in the preceding sections. Cases where P_B(2, v) has a negative integer root in v occur only for v = −1 (the chromatic polynomial) and v = −2 (the flow polynomial). More precisely, v = −1 factorizes for the triangular, kagome, frieze, three–twelve, snub square, snub hexagonal and ruby lattices. And v = −2 factorizes for the hexagonal, four–eight, frieze, three–twelve, cross, snub square and snub hexagonal lattices.

For q = 3, P_B(q, v) has a root at v = −3 for the four–eight and three–twelve lattices. Note that these are three-flow problems or, equivalently, three-colouring problems of the corresponding dual lattices.

Finally, for q = 4 there is a root at v = −2 for all lattices except the kagome lattice. Motivated by this, and by the phase diagrams reported in the preceding sections, we conjecture that the BK phase extends to (q, v) = (4, −2) for all the Archimedean lattices, except the kagome lattice.

A four-terminal representation of the four–eight medial lattice is shown in figure 25. At first sight it does not appear feasible to write the corresponding $\check{\sf R}$-matrix in the usual form, namely, as a product of the single-edge operators ${\sf V}_i$, ${\sf H}_{i+1}$ and ${\sf V}_{i+2}$ (or, more generally, the Temperley–Lieb generators ${\sf E}_i$, ${\sf E}_{i+1}$ and ${\sf E}_{i+2}$) acting within the unit cell shown in figure 4. The problem is that we need an intermediate point in the spin representation, or two intermediate strands in the loop representation that would be situated between those labelled i + 1 and i + 2 in figure 4. These intermediate strands can however be eliminated once the grey square is completed, and they are not needed for connecting among themselves the grey squares of which the lattice consists.

Zoom In Zoom Out Reset image size

To avoid dealing with intermediate points, the most efficient solution is to write down directly the entire $\check{\sf R}$-matrix that propagates strands i, i + 1, i + 2, i + 3 into i', (i + 1)', (i + 2)', (i + 3)'. Recall that a single-edge operator, such as ${\sf V}_i$, consists of two terms (${\sf Id}_i$ and ${\sf E}_i$), since there are Cat(2) = 2 possible planarity-respecting pairings of the four points i, i + 1, i', (i + 1)'. Similarly, the entire $\check{\sf R}$-matrix contains in general fourteen terms, corresponding to the Cat(4) = 14 pairings of eight points that respect planarity. We have therefore written a version of the algorithm in which a lattice is specified by supplying the 14 terms of a generic $\check{\sf R}$-matrix, each of which are polynomials in n_loop and x with integer coefficients. Further remarks on this version can be found in section 3.6.

For the case at hand we remark that $\check{\sf R}_i = ({\sf B}_i)^2$, where ${\sf B}_i$ denotes the bow–tie operator, first discussed for the special case of percolation in [21, section 3.2] and subsequently generalized to the Potts model in [14, section 3.2]. Squaring the explicit expression [14, equation (27)] we therefore obtain the fourteen polynomials defining $\check{\sf R}_i$, each of which contains up to a maximum of ten monomials $x^a n_{\rm loop}^b$.

To test the general algorithm we have also investigated the case where each $\check{\sf R}$-matrix is a single bow–tie operator, $\check{\sf R}_i = {\sf B}_i$. This can be represented as in figure 26. Obviously this is just a rotated version of the four-terminal representation of the kagome lattice shown in figure 9. We have validated the 'generic $\check{\sf R}$-matrix' algorithm by verifying that in this case it gives the very same critical polynomials P_B(q, v) as were obtained in section 4.3. The choice of transfer direction made in section 4.3 is slightly more efficient for dealing with the kagome lattice, since $\check{\sf R}_i$ involves the application of six operators each with two terms (6 × 2 = 12) rather than a single operator with fourteen terms. Moreover, the approach with two terms per operator involves only very simple coefficients (1 or x).

Zoom In Zoom Out Reset image size

Returning to the four–eight medial lattice, we remark that in figure 25 there are six vertices and twelve edges per grey square (see table 28). The approximations to the bond percolation threshold p_c obtained from the unique positive root of P_B(p) are shown in table 29.

For the Ising model (q = 2) the polynomials P_B(q, v) always factorize. The maximum degree of the factors is d_max = 8 for n = 1, 2, d_max = 16 for n = 3, and d_max = 12 for n = 4. One of these factors, namely

$$\begin{equation} -64 - 128 v - 160 v^2 - 96 v^3 - 8 v^4 + 32 v^5 + 24 v^6 + 8 v^7 + v^8 , \end{equation} \tag{ 59 }$$

occurs systematically for any n. On making the change of variables $v = -1 + \sqrt{y}$ this simplifies to

$$\begin{equation} -23 - 20 y - 18 y^2 - 4 y^3 + y^4 . \end{equation} \tag{ 60 }$$

The unique positive root is

$$\begin{equation} v_{\rm c} = -1 + ( 1 + \sqrt{2} + \sqrt{10 + 8 \sqrt{2}} )^{1/2} \simeq 1.651\,582\,692\ldots . \end{equation} \tag{ 61 }$$

We expect this to be the exact critical point, although we are not aware of any exact solution of the Ising model on the four–eight medial lattice.

Moreover, the critical points for the q = 3-state and q = 4-state Potts models are given in tables 30–31.

The phase diagram for the four–eight medial lattice is rather simple; see figure 27. There is a clear vertical ray at the Beraha number (21) with k = 4; and possibly the n = 3 polynomial also indicates that a ray will emerge at k = 6, although the n = 4 result fails to confirm this. It seems likely that the BK phase will extend to the arc near q ≈ 3.3, although its upper and lower boundaries are not visible on the interval 2 < q < 3 with these critical polynomials. In any case, the absence of vertical rays for k > 6 is a clear sign that q_c < 4 for this lattice. We also note that all curves go through the point (q, v) = (0, −3) exactly.

Zoom In Zoom Out Reset image size

For the frieze medial lattice we can use the four-terminal representation depicted in figure 28. It consists of alternating rows of grey squares of the types used in the kagome and square lattices. Therefore, the $\check{\sf R}$-matrix is given by (29) on even rows, and by (25) on odd rows. It is thus convenient to use rectangular bases of size n × 2n grey squares, for any parity of n. There are 5/2 vertices and 5 edges per grey square.

Zoom In Zoom Out Reset image size

The corresponding bond percolation thresholds p_c are shown in table 32.

When q = 2 we obtain as usual a factorization of P_B(q, v). The maximum degree of the factors is d_max = 8 for n = 1, 2, d_max = 16 for n = 3, 4, and d_max = 32 for n = 5. The factor relevant for determining the critical point simplifies upon setting $v = -1 + \sqrt{y}$ and becomes

$$\begin{equation} 1 - 20 y - 10 y^2 - 4 y^3 + y^4 . \end{equation} \tag{ 62 }$$

Its physically relevant solution has an expression in terms of cube roots, which is however too lengthy to be reported here. It corresponds to v_c ≃ 1.479 990 057... in the original variable.

The critical points v_c for the q = 3-state and q = 4-state Potts models appear in tables 33–34.

The phase diagram of the frieze medial lattice, shown in figure 29, is exceedingly complicated, and arguably even more complicated than that of the three–twelve lattice (see figure 16). The upper limit v₊(q) of the BK phase contains an almost straight part from the origin to the neighbourhood of (q, v) = (3, −1). The almost straight continuation to higher q cannot be the upper limit of the BK phase, though, since it is not adjacent to the characteristic vertical rays. Rather, the continuation of v₊(q) must be provided by the n = 4 arc that bends round at (q, v) ≈ (3.5, −1.6). For n = 2 this arc is prefigured by a bubble that is not connected to the rest of the curves. It seems likely that all even n will participate in this part of the phase diagram. Note in particular that the n = 4 arc has a small wrinkle at q ≃ B₈ which would most likely turn into a vertical ray for higher (even) n. One can therefore believe that q_c > B₈ for this lattice, but it is as yet unclear whether q_c might be as large as 4. After the n = 4 arc bends round, it traces out the lower limit v₋(q) of the BK phase that continues to the point (q, v) = (0, −3) through which all curves pass exactly. The extent of the BK phase can be judged from the vertical rays at q = B_k with k = 4, 6 (and maybe 8 as just mentioned).

Zoom In Zoom Out Reset image size

Two further branches of the curve come out of (q, v) = (0, −3) and trace out a finger that extends a little further than the vertical ray at q = 2. There is a similar, broader finger coming out of (q, v) = (0, 0), whose upper side coincides with v₊(q).

Apart from these features, there is almost horizontal branch coming out of (q, v) ≈ (0, −3.7) and extending towards large q. Similarly, the bottom and top of the vertical ray at q = 3 connect to branches that extend towards large q.

A close-up of the region near (q, v) = (0, −3) is shown in figure 30. Throughout this region the curves become increasingly dense as n increases. Counting the number of fingers emanating from the v axis in the range −3.7 < v < 3 gives compelling evidence for the conjecture that the curves will, in fact, become space filling in this region when n → ∞. This is a novel feature, not seen in the phase diagrams for any of the other lattices. It is very much reminiscent of what emerged from a recent study of the phase diagram of the Potts model for a family of non-planar graphs, called the generalized Petersen graphs [66], where a number of 'critical regions' (marked by a ⋆ in [66, figure 2]) were identified throughout which the two dominant eigenvalues of the transfer matrix are exactly degenerate in norm. Such critical regions also exist in two dimensions, and in particular for the q-state Potts model on the triangular lattice close to the point (q, v) = (0, −3), as well as for q > 4 [47].

Zoom In Zoom Out Reset image size

This capability of the P_B(q, v) = 0 curves to be space filling might also offer a new interpretation of the thin fingers emanating from the v axis in many of the preceding figures (see, e.g., figure 16). Might it be that these fingers will also become more numerous and tend to fill out space for larger (i.e., not accessible in this paper) values of n? We leave this question for future investigations.

A four-terminal representation of the three–twelve medial lattice is shown in figure 31. This lattice is also known as the 2 × 2 kagome subnet [65], since it can be obtained by replacing each of the triangles of the kagome lattice by an equilateral form made of 2 × 2 = 4 triangles. In our representation there are 9 vertices and 18 edges per grey square—the highest numbers for any lattice considered in this paper. This means that the basis for the n = 7 numerical computation encompasses 18 × 7² = 882 edges, as mentioned in the abstract.

Zoom In Zoom Out Reset image size

The $\check{\sf R}$-matrix can only be computed in the sparse-matrix factorization scheme if one inserts intermediate points. In this case, three such points (or six loop strands) are needed. Just like in section 6.1, it is therefore advantageous to use the 'generic $\check{\sf R}$-matrix' version of the algorithm. The fourteen weights can be readily computed separately, from a transfer matrix that builds a single grey square, using time slices of width five points (or ten loop strands). Alternatively, this may be done by hand, defining a transfer process where time flows in the north-west (rather than the usual north-east) direction, and using the fact that $\check{\sf R}_i = \tilde{\sf S}_i {\sf S}_i$, where ${\sf S}_i$ denotes the 2 × 2 subnet operator, and $\tilde{\sf S}_i$ is its time-reflected counterpart. Either way, we obtain the fourteen polynomials defining $\check{\sf R}_i$, each of which contains up to a maximum of 30 monomials $x^a n_{\rm loop}^b$. These are obviously too lengthy to be reproduced here, but are available upon request from the author.

The bond percolation thresholds p_c are shown in table 35. As for the three–twelve lattice itself, we obtain a very substantial improvement on the precision of the existing results, here by more than four orders of magnitude.

For the Ising model (q = 2) we have the usual factorization of P_B(q, v). The maximum degree of the factors is d_max = 8 for n = 1, 2, d_max = 16 for n = 3, 4. The relevant factor for determining the critical point has degree 8 in the v variable, but changing variables through $v = -1 + \sqrt{y}$ we find the simpler polynomial

$$\begin{equation} -83 - 32 y - 6 y^2 - 8 y^3 + y^4 . \end{equation} \tag{ 63 }$$

The unique root that corresponds to a positive value of v reads

$$\begin{equation} v_{\rm c} = -1 + (2 + \sqrt{3} + \sqrt{2 (6 + 5 \sqrt{3})} )^{1/2} \simeq 2.024\,382\,957\ldots . \end{equation} \tag{ 64 }$$

Once again, we expect this to be the exact critical point, although we are not aware of any exact solution of the Ising model on the three–twelve medial lattice.

The critical points v_c for the q = 3-state and q = 4-state Potts models are given in tables 36–37.

Despite the very large size of the bases, the phase diagram of the three–twelve medial lattice (see figure 32) is not very complicated. The upper and lower boundaries v_±(q) of the BK phase are a couple of curves going out of the points (q, v) = (0, 0) and (0, −3), respectively, with finite slopes. They join via an arc at q ≈ 2.6 which is however only visible in the n = 3 result. Accordingly the only vertical ray is at q = B₄ = 2. The role of the elongated bubble at q ≈ 2.7 in the n = 3 curve is not clear and would have to be confirmed at larger sizes. An additional, lowest-lying curve goes out of (q, v) = (0, −3) vertically and continues to large q; this curve is outside the BK phase since it does not touch the vertical ray at q = 2.

Zoom In Zoom Out Reset image size

Note that the large degree of the polynomials P_B(q, v) makes the computation very memory demanding, so, departing from the general rule, we have not computed the case n = 5 for this lattice.

Figure 33 shows a four-terminal representation of the cross medial lattice. The $\check{\sf R}$-matrix is identical to that of the four–eight medial lattice (see section 6.1) for grey squares where at least one of x and y is even. When both x and y are odd the $\check{R}$-matrix is simply the identity. This corresponds formally to setting the coupling constant to infinity on two of the edges, and is represented by a coil-like symbol in the figure. We have therefore 9/2 vertices and 9 edges per grey square. Note that this representation is defined only for even n.

Zoom In Zoom Out Reset image size

The reader might want to verify the presence of hexagons and dodecagons in figure 33, besides the obvious squares. All of these polygons share each of their edges with a triangle, as they should, since the underlying (4, 6, 12)-lattice is a cubic graph.

The percolation thresholds p_c are given in table 38.

For the Ising model, the largest degree of the factors is d_max = 16 for both n = 2 and n = 4. The factor relevant for determining v_c simplifies upon setting $v = -1 + \sqrt{y}$ and becomes

$$\begin{equation} -647 - 2192 y - 2700 y^2 - 1952 y^3 - 594 y^4 - 80 y^5 - 28 y^6 + y^8 . \end{equation} \tag{ 65 }$$

The relevant root cannot be written simply, but reads numerically v_c ≃ 1.726 376 028....

The critical points for the q = 3-state and q = 4-state Potts models are displayed in tables 39–40.

The phase diagram of the cross medial lattice, shown in figure 34, is rather simple. The upper and lower boundaries v_±(q) of the BK phase go out of the points (q, v) = (0, 0) and (0, −3) with finite slope. Those curves do not continue beyond the vertical ray at q = B₄ = 2, but from the experience with other lattices we can safely assume that this is an artefact of our choice of bases. In particular, since the q = 2 ray has finite length, the BK phase must extend further to the right. Presumably it ends at the arc extending from q ≈ 3.0 to q ≈ 3.4. If so, we would expect a further vertical ray at q = B₆ = 3 to build up for larger bases (note that the lower end point of the arc in the n = 4 curve is conspicuously close to q = 3). Finally, there is another, lowest-lying curve going out of (q, v) = (0, −3) vertically which is outside the BK phase.

Zoom In Zoom Out Reset image size

A four-terminal representation of the snub square medial lattice is shown in figure 35. Its $\check{\sf R}$-matrix is a mixture of known ingredients: it is given by $\check{\sf R}_i$ of the square lattice, equation (25), when x + y is odd; by that of the kagome lattice, equation (29), when x and y are both even; and by the alternative kagome representation discussed in section 6.1 and depicted in figure 26 when x and y are both odd. This representation contains 5/2 vertices and 5 edges per grey square. Once again we must require n to be even.

Zoom In Zoom Out Reset image size

The approximate percolation thresholds p_c obtained from the positive root of P_B(p) are displayed in table 41.

The polynomials P_B(q, v) for the Ising model (q = 2) factorize, and the maximum degree of the factors is d_min = 16 for n = 2, 4. However, if we perform the change of variables $v = -1 + \sqrt{y}$, the polynomial determining y is only of degree 8:

$$\begin{equation} 1 - 72 y - 304 y^2 - 320 y^3 - 226 y^4 - 88 y^5 - 16 y^6 + y^8 . \end{equation} \tag{ 66 }$$

The largest positive root in y corresponds to the unique positive root in v, which is v_c ≃ 1.480 593 024....

Critical points of the q = 3-state and q = 4-state Potts models are given in tables 42–43.

The phase diagram of the snub square medial lattice is shown in figure 36. There are vertical rays at q = B_k with k = 4, 6, 8. The gap in the rendering of the boundaries v_±(q) of the BK phase for 2 < q < 3 is partly filled out by the bubble appearing in the n = 4 curve. The curve coming out of (q, v) = (0, 0) with infinite slope bends round and goes through the point (q, v) = ( − 3, 0). We notice just below the latter point the formation of several narrow fingers which might well tend to be space filling at larger n, in analogy with what was observed for the frieze medial lattice (see figure 30). Finally, at the bottom of the phase diagram there is an almost horizontal curve that lies below the BK phase and extends out to large q.

Zoom In Zoom Out Reset image size

Recall from section 4.9 that for the snub hexagonal lattice we have computed P_B(q, v) using two different methods. Also for the corresponding medial lattice shall we present two distinct constructions.

The first construction proceeds in analogy with that of the snub hexagonal lattice itself. Just like how the latter was obtained in section 4.9, as a depleted version of the triangular lattice, we can construct its medial by depleting the kagome (medial of the triangular) lattice. The ${\sf H}$ and ${\sf V}$ operators corresponding to erased edges are replaced by just one of their two terms (the identity ${\sf Id}$ or the Temperley–Lieb generator ${\sf E}$, as the case may be). We thus obtain 15/7 vertices and 30/7 edges per grey square.

In formal terms, the $\check{\sf R}$-matrix reads

$$\begin{equation} \check{\sf R} = \left\lbrace \begin{array}{@{}ll@{}}\sf H_{i+1} {\sf V}_{i+2} {\sf V}_{i} {\sf E}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} & {{\rm for}\quad 3 y + x = 0,3,6\; {\rm mod}\; 7} \\ {\sf Id}_{i+1} {\sf V}_{i+2} {\sf V}_{i} {\sf E}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf H}_{i+1} & {{\rm for}\quad 3 y + x = 1\;{\rm mod}\; 7} \\ {\sf E}_{i+1} {\sf V}_{i+2} {\sf E}_i {\sf E}_{i+1} & {{\rm for}\quad 3y + x = 2\;{\rm mod}\; 7} \\ {\sf E}_{i+1} {\sf E}_{i+2} {\sf V}_i {\sf E}_{i+1} & {{\rm for}\quad 3y + x = 4\;{\rm mod}\; 7} \\ {\sf H}_{i+1} {\sf V}_{i+2} {\sf V}_{i} {\sf E}_{i+1} {\sf V}_{i+2} {\sf V}_i {\sf Id}_{i+1} & {{\rm for}\quad 3 y + x = 5\;{\rm mod}\; 7.} \\ \end{array} \right. \end{equation} \tag{ 67 }$$

The first line in (67) corresponds to undepleted bow–tie patterns in the grey squares, i.e., it coincides with (29). Subsequent lines are obtained from the first one by depletion, i.e., replacing some of the ${\sf H}$ operators by the identity ${\sf Id}$, and some of the ${\sf V}$ operators by the Temperley–Lieb generator ${\sf E}$.

As with the snub hexagonal lattice itself (see section 4.9), the representation (49) only makes sense when n is a multiple of 7. So the only computation that we can perform in practice is to find numerically the roots in v of P_B(q, v) = 0 with n = 7. We therefore turn now to an alternative construction.

This second construction is based on the four-terminal representation shown in figure 37. The periodicity of the tiling is 3 horizontally and 6 vertically, so we can use it with rectangular bases of size n × 2n provided that n = 0 mod 3. Like the snub square medial lattice (see figure 35) it uses both the kagome bow–tie and its rotated counterpart. More precisely, for even y the $\check{\sf R}$-matrix is given by (29) when x + y = 0, 1 mod 3, and is that of the alternative kagome representation (see figure 26) when x + y = 2 mod 3. For odd y we have

$$\begin{equation} \check{\sf R}_i = \left\lbrace \begin{array}{@{}ll@{}}\sf E_{i+2} {\sf V}_i {\sf H}_{i+1} & {{\rm for}\quad x+y = 0\;{\rm mod}\; 3} \\ {\sf H}_{i+1} {\sf V}_{i+2} {\sf E}_i & {{\rm for}\quad x+y = 1\; {\rm mod}\; 3} \\ {{\rm Equation}\; ({29})} & {{\rm for}\quad x+y = 2\; {\rm mod}\; 3.} \\ \end{array} \right. \end{equation} \tag{ 68 }$$

Moreover, one must place a horizontal diagonal on the white squares having coordinates $(x+\frac{1}{2},y+\frac{1}{2})$ whenever x + 2⌊y/2⌋ = 2 mod 3. The reader should check that each hexagon (which contains two of the horizontal diagonals) shares each of its edges with a pentagon and each of its vertices with a triangle.

Zoom In Zoom Out Reset image size

This construction provides 5/2 vertices and 5 edges per grey square. The packing density has therefore been improved with respect to the first construction. However, due to the parity constraint on n we can now only attain n = 6, instead of n = 7.

In table 44 we show the percolation thresholds p_c using both the first and the second constructions. We have separated the results using the two different constructions by a horizontal line in this and the following two tables. Extrapolating the results with n = 3, 6 leads to the estimate for the n → ∞ limit shown in the last line of table 44. We have here supposed that the parameter w in (30) is in the same range, up to a confidence interval of ±0.2, as that found for the closely related snub square medial lattice (see section 6.5). The distance of the n = 7 result from the final value is consistent with our scaling analysis.

In the case of the Ising model (q = 2) the polynomial P_B(q, v) with n = 3 factorizes. The degrees of the largest factors are 1, 4, 6, 12 and d_max = 44. We note that the value of d_max is unusually large. The factor leading to a positive real root is the one of degree 4. It simplifies upon setting $v = -1 + \sqrt{y}$, becoming

$$\begin{equation} -3 - 6y + y^2 . \end{equation} \tag{ 69 }$$

Our conjecture for the critical point is thus

$$\begin{equation} v = -1 + \sqrt{3 + 2 \sqrt{3}} \simeq 1.542\,459\,756\ldots \end{equation} \tag{ 70 }$$

and this is seen to coincide with the result (56) for the ruby lattice. By means of a 50-digit numerical computation we have verified that the number (70) is also a root of the n = 6 polynomial, and of the n = 7 polynomial arising from the first construction. This provides compelling evidence that, on one hand, the factor (69) is recurrent for all of the critical polynomials, and, on the other hand, that our two different constructions lead to consistent results.

The results for the critical point v_c in the q = 3-state and q = 4-state Potts models are shown in tables 45–46. Also in those cases we have based the n → ∞ extrapolations on the results coming from the second construction (with n = 3, 6), and the approach is similar to that described above for q = 1. Again, the n = 7 results confirm our scaling analysis.

The phase diagram of the snub hexagonal medial lattice is shown in figure 38. Unfortunately this is based on the single size n = 3, so we cannot make very firm statements about the convergence properties. The usual vertical rays are visible at q = B_k with k = 4, 6. The rightmost arc of the bubble containing the q = 2 vertical ray may well be a precursor of the v_±(q) curves on the interval 2 < q < 3. On the other hand, the part of v₊(q) with 0 < q < 2 is clearly missing with this choice of basis. There is a narrow finger close to the point (q, v) = (0, −3) through which the curve passes twice.

Zoom In Zoom Out Reset image size

The ruby medial lattice can be represented in four-terminal form as shown in figure 39. This representation is valid only when n is a multiple of 4, and so we shall be limited to studying the case n = 4 in the following.

Zoom In Zoom Out Reset image size

As discussed in section 3.4 one can act on the two strands in a white square using the operator ${\sf H}_i$ in order to produce a horizontal diagonal. This has been done here on the white squares with coordinates $(x,y) = (\frac{3}{2},\frac{1}{2}), (\frac{3}{2},\frac{3}{2}), (\frac{7}{5},\frac{5}{2}), (\frac{7}{2},\frac{7}{2})$ mod 4. But for the present lattice we shall also need a slight variation of this trick, namely to act with the operator ${\sf E}_i$ instead. This corresponds to formally setting the coupling strength x = ∞, and renormalizing by a factor of x, meaning that the two sites sitting across the horizontal diagonal of the white square will be effectively identified, or contracted. This operation is represented by a coil-like symbol in figure 39, and we use it in the white squares with coordinates $(x,y) = (\frac{1}{2},\frac{1}{2}), (\frac{1}{2},\frac{3}{2}), (\frac{5}{2},\frac{5}{2}), (\frac{5}{2},\frac{7}{2})$ mod 4. The reader may want to verify from the figure that each triangle is surrounded by six squares (three of which share an edge with the triangle, with the other three sharing only a vertex), and similarly that each hexagon is surrounded by twelve squares (that again alternate between edge-sharing and vertex-sharing). Each hexagon has been represented as an octagon with two x = ∞ edges.

Similarly, some of the edges in the grey squares have coupling strength x = 0, i.e., they reduce to operators ${\sf Id}$ or ${\sf E}$. Thus, each hexagon comprises two spins and one dual spin that are 'free', in the sense that they do not interact with the remainder of the lattice. Since each 4 × 4 pattern contains two hexagons, this leads to an extraneous factor q⁶ by which we must divide in order to form P_B(q, v). Obviously the occurrences of x = 0 or x = ∞ do not count as edges of the lattice that we are investigating, and therefore diminish the efficiency of the representation. In the present case we have therefore a rather modest number of 3/2 vertices and 3 edges per grey square.

The $\check{\sf R}$-matrix acting on the grey squares can be formally described as

$$\begin{equation} \check{\sf R}_i = \left\lbrace \begin{array}{@{}ll@{}}\sf H_{i+1} {\sf V}_i {\sf V}_{i+2} {\sf H}_{i+1} & {\rm if }\quad (x,y)=(0,0)\; {\rm mod }\; 2 \\ {\sf H}_{i+1} {\sf V}_i {\sf V}_{i+2} {\sf H}_{i+1} & {\rm if }\quad (x,y)=(1,3) \;{\rm or }\; (3,1)\; {\rm mod }\; 4 \\ {\sf E}_i {\sf H}_{i+1} {\sf V}_{i+2} {\sf H}_{i+1} & {\rm if }\quad (x,y)=(1,0) \;{\rm or }\; (3,2)\; {\rm mod }\; 4 \\ {\sf H}_{i+1} {\sf V}_i {\sf H}_{i+1} {\sf E}_{i+2} & {\rm if }\quad (x,y)=(3,0) \;{\rm or }\; (1,2) \;{\rm mod }\; 4 \\ {\sf V}_i {\sf E}_{i+2} {\sf H}_{i+1} & {\rm if }\quad (x,y)=(0,1) \;{\rm or }\; (2,3) \;{\rm mod }\; 4 \\ {\sf H}_{i+1} {\sf E}_i {\sf V}_{i+2} & {\rm if }\quad (x,y)=(2,1) \;{\rm or }\; (0,3) \;{\rm mod }\; 4 \\ {\sf E}_i {\sf E}_{i+2} & {\rm if }\quad (x,y)=(1,1) \;{\rm or }\; (3,3) \;{\rm mod }\; 4. \\ \end{array} \right. \end{equation} \tag{ 71 }$$

The n = 4 result for the percolation threshold p_c is shown in table 47, and those for the critical points v_c of the q = 3-state and q = 4-state Potts models are reported in tables 48–49. From these single data points, which moreover do not correspond to a very large size of the basis compared to the other lattices treated in this work, it does not seem reasonable to provide a final result with error bars for the n → ∞ limit.

For the Ising model (q = 2) the factorization of P_B(q, v) with n = 4 contains factors of maximum degree d_max = 20. The factor determining v_c is of degree 16, and on setting $y = -1 + \sqrt{y}$ it becomes simpler:

$$\begin{equation} -3 - 6 y - 66 y^2 - 174 y^3 - 194 y^4 - 58 y^5 - 10 y^6 - 2 y^7 + y^8 . \end{equation} \tag{ 72 }$$

The largest real root in y determines the critical coupling as v_c ≃ 1.477 488 025.... On the basis of our experience from the Archimedean lattices we can assume this value to be the exact result.

The phase diagram in shown in figure 40. There is a clear vertical ray at the Beraha number (21) with k = 4, and another incipient ray with k = 6. Note also that the n = 4 curve passes through the point (q, v) = (0, −4) exactly. The curve emanating from that point is presumably a good approximation to the lower boundary of the BK phase. The value of q_c is however difficult to estimate from just one curve.

Zoom In Zoom Out Reset image size

We now discuss, in analogy with section 4.11, the cases of exact factorization for the seven medial lattices which are not themselves Archimedean (i.e., those discussed in section 6). Statements including the words 'all' or 'none' thus refer to those seven lattices only.

Note that for none of the medial lattices does P_B(q, v) factorize for generic values of q and v. This presumably means that the Potts model is not solvable on these lattices along any curve. Conversely, factorization does occur for isolated (integer) values of q and v, as we now discuss.

For q = 0, P_B(q, v) factorizes for all the lattices, producing a root v = 0. The resulting free fermion theories describe spanning trees [61, 62]. Moreover, the four–eight medial, frieze medial, three–twelve medial, cross medial, snub square medial and snub hexagonal medial lattices have a root (q, v) = (0, −3). And the ruby medial lattice has a root (q, v) = (0, −4). These are models of spanning forests [61] on the corresponding dual lattices, with weight v per component tree.

Unlike for the case of Archimedean lattices, we have found no medial lattices with a size-independent (finite) slope with which the curves pass through the origin (q, v) = (0, 0).

There are a couple of unusual cases for q = 0. The four–eight medial and cross medial lattices both shed the small factor 8 + 3v + v² with roots $v = (-3 \pm {\rm i}\sqrt{23})/2$. Similarly, the three–twelve medial lattice sheds the factor 6 + 3v + v² with roots $(-3 \pm {\rm i}\sqrt{15})/2$. All of these cases presumably provide loci of exact solvability.

In the Ising case (q = 2) there is a root in v = −1 (chromatic polynomial) for all the medial lattices. Unlike the case for Archimedean lattices, there are no medial lattices with a root at v = −2.

For q = 3 there is a root at v = −1 (chromatic polynomial) for the four–eight medial lattice provided that n is odd (we have checked this for n = 1, 3), while for n even (i.e., for n = 2, 4) the curves in figure 27 do not even come close to the point (q, v) = (3, −1). However, we have seen in this study that parity effects in n are extremely common, so we feel confident in conjecturing that the three-colouring problem on the four–eight medial lattice is exactly solvable¹⁰.

For q = 3 the frieze medial and snub hexagonal medial lattices both shed the small factor 3 + 3v + v² whose roots, $v = (-3 \pm {\rm i} \sqrt{3})/2$, are presumably loci of exact solvability.

Finally, for q = 4 there is a root at v = −2 for the frieze medial lattice.

Figure 41 shows a four-terminal representation for site percolation on the hexagonal lattice. There are two vertices per grey square.

Zoom In Zoom Out Reset image size

It is useful at this stage to point out the key differences with bond percolation. In bond percolation, the 'conducting units' are the edges. Since edges meet at vertices, each of the four terminals of the grey squares must be situated at the position of a vertex. On the other hand, in site percolation the 'conducting units' are the vertices. It is convenient to represent an occupied site instead as a colouring of its adjacent half-edges, so that site percolation clusters become connected components (clusters) of coloured half-edges. It follows from this picture that each of the four terminals of the grey squares must be situated at the midpoint of an edge.

The loop strands of figure 4 must turn around the clusters of coloured half-edges. It follows that each choice of occupancy of sites within a grey square induces a planar pairings of the strands i, j, k, l := i, i + 1, i + 2, i + 3 and i', j', k', l'. In the case of figure 41 the $\check{\sf R}$-matrix becomes

$$\begin{eqnarray} &&\fl \check{\sf R}_i = (ij)(kl)(i^{\prime }j^{\prime })(k^{\prime }l^{\prime }) + v \, (il)(jk)(i^{\prime }j^{\prime })(k^{\prime }l^{\prime }) \nonumber \\ &&+ v \, (ij)(kl)(i^{\prime }l^{\prime })(j^{\prime }k^{\prime }) + v^2 \, (ii^{\prime })(jk)(j^{\prime }k^{\prime })(ll^{\prime }) , \end{eqnarray} \tag{ 73 }$$

where the first term corresponds to both sites within the grey square being empty, the second and third terms correspond to one site being occupied and the other empty, and the fourth term comes from the case where both sites are occupied. This can be expressed more elegantly in terms of TL generators:

$$\begin{equation} \check{\sf R}_i = {\sf E}_{i+2} {\sf E}_i + v \, ({\sf E}_{i+2} {\sf E}_i {\sf E}_{i+1} + {\sf E}_{i+1} {\sf E}_{i+2} {\sf E}_i) + v^2 \, {\sf E}_{i+1} . \end{equation} \tag{ 74 }$$

The fact that (73) contains only four terms out of the fourteen possible means that the computation of P_B(p) can be accomplished for square bases of size up to n = 8. The corresponding thresholds p_c are shown in table 51.

The exponent appearing in (30) here comes out as w ≈ 6.35, which is the same value as was found for bond percolation on the kagome and three–twelve lattices. This seems reasonable, since all those lattices have the same threefold rotational symmetry¹². The high value of w again entails a high precision of the final value: more than two orders of magnitude better than the best numerical result [25].

Note that graph polynomials for this lattice were previously studied in [21] for hexagonal bases with up to 54 sites (compared to the square bases with up to 128 sites used here).

Site percolation on the square lattice was already discussed as an example in the introduction to section 7. The four-terminal representation is shown in figure 42. It has one vertex per grey square.

Zoom In Zoom Out Reset image size

The $\check{\sf R}$-matrix reads

$$\begin{equation} \check{\sf R}_i = (ij)(kl)(i^{\prime }j^{\prime })(k^{\prime }l^{\prime }) + v \, (ii^{\prime })(jk)(ll^{\prime })(j^{\prime }k^{\prime }) , \end{equation} \tag{ 75 }$$

as discussed previously. In terms of TL generators this reads

$$\begin{equation} \check{\sf R}_i = {\sf E}_{i+2} {\sf E}_i + v \, {\sf E}_{i+1} . \end{equation} \tag{ 76 }$$

Since there are only two terms out of fourteen possible, we have been able to obtain P_B(p) for bases of size up to n = 11. The corresponding thresholds p_c are given in table 52.

Since we have here eleven data points—the highest number for any of the problems treated in this paper—we have taken particular care to get the best possible extrapolation out of them. Applying the usual procedure we first found w ≈ 4.07 from a non-linear fit (30) to the last three data points. However, it was easily detected that this choice of w led to some unnecessary spread on the BS approximants. In fact, fitting successively w from the last three data points among the first nine, ten or all eleven points, we found w ≈ 4.139, w ≈ 4.098 and w ≈ 4.074, indicating that the true w might be slightly lower. Repeating then the BS procedure while moving w down in steps of 0.01, we have checked that the choice w = 4.03 ± 0.01 minimizes the spread of the approximants, and hence the final value given in table 52 is based on this choice. This is consistent with (and slightly more accurate than) the best numerical result [25].

The four-terminal representation for site percolation on the four–eight lattice is shown in figure 43. It can accommodate four vertices per grey square.

Zoom In Zoom Out Reset image size

The $\check{\sf R}$-matrix is now

$$\begin{eqnarray} &&\fl \check{\sf R} = (1+4v+2v^2) \, (ij)(kl)(i^{\prime }j^{\prime })(k^{\prime }l^{\prime }) + v^4 \, (ii^{\prime })(jk)(j^{\prime }k^{\prime })(ll^{\prime }) \nonumber \\ &&+ v^2 [ (ii^{\prime })(jj^{\prime })(kl)(k^{\prime }l^{\prime }) + (ij)(i^{\prime }j^{\prime })(kk^{\prime })(ll^{\prime }) \nonumber \\ &&+(il)(jk)(i^{\prime }j^{\prime })(k^{\prime }l^{\prime }) + (ij)(kl)(i^{\prime }l^{\prime })(j^{\prime }k^{\prime }) ] \nonumber \\ &&+ v^3 [ (ii^{\prime })(j^{\prime }k^{\prime })(jl^{\prime })(kl) + (ii^{\prime })(jk)(j^{\prime }l)(k^{\prime }l^{\prime }) \nonumber \\ &&+(ik^{\prime })(jk)(i^{\prime }j^{\prime })(ll^{\prime }) + (i^{\prime }k)(j^{\prime }k^{\prime })(ij)(ll^{\prime }) ] . \end{eqnarray} \tag{ 77 }$$

With ten terms out of the fourteen possible, the advantage of the hashing scheme over the complete tabulation of states is less pronounced than for the lattices treated previously. Accordingly we can treat square bases of size up to n = 7. The results for the threshold p_c are given in table 53.

We give a four-terminal representation for site percolation on the cross lattice in figure 44. It requires n to be even and hosts three vertices per grey square. The $\check{R}$-matrix is the same as for the four–eight lattice, except when x and y are both odd in which case it is replaced by the identity operator. One may check the presence of faces of degree 6 and 12.

Zoom In Zoom Out Reset image size

Table 54 provides the results for the site percolation threshold p_c.

The four-terminal representation for site percolation on the ruby lattice requires some new tricks. It is shown in figure 45, where we have supposed that n = 0 mod 4. As usual the coil-like symbol indicates an edge of infinite strength, meaning that the two vertices at its end points have been identified. The states of the two grey squares linked by a coil are now correlated: if the site in one of the squares is empty (resp. occupied) the same is true in the adjacent square, since the two sites have been identified. Identifying such conglomerates of two grey squares by the coordinates (x, y) of the leftmost one, we have the following $\check{\sf R}$-matrix:

$$\begin{equation} \check{\sf R}_i = \left\lbrace \begin{array}{@{}ll@{}}\sf E_{i+4} {\sf E}_{i+2} {\sf E}_i + v \, {\sf E}_{i+3} {\sf E}_{i+4} {\sf E}_{i+1} {\sf E}_{i+2} & {\rm if }\quad x+y=0 \;{\rm mod }\; 4 \\ {\sf E}_{i+4} {\sf E}_{i+2} {\sf E}_i + v \, {\sf E}_{i+2} {\sf E}_{i+3} {\sf E}_{i} {\sf E}_{i+1} & {\rm if }\quad x+y=3 \;{\rm mod }\; 4 \\ \end{array} \right. \end{equation} \tag{ 78 }$$

whereas the remaining grey squares (those without coils in figure 45) are described by (76), the usual $\check{\sf R}$-matrix of the square lattice.

Zoom In Zoom Out Reset image size

As a result we have an average of 3/4 sites per grey square. While this might seem mediocre, it should be remembered that the $\check{\sf R}$-matrix (78) for a conglomerate of two grey squares acts on six pairs of points, labelled i, i + 1, ..., i + 5 and i', (i + 1)', ..., (i + 5)', which may in general accommodate Cat(6) = 132 different connectivities. However, (78) contains only two terms, so the connectivity states actually generated by the transfer process are only a very small subset of the total number of states respecting planarity. Accordingly we can attain a size of n = 16, the largest for any of the problems studied in this paper. Table 55 gives the corresponding results for the threshold p_c.

Figure 46 provides a four-terminal representation of the Cairo pentagonal lattice. For the convenience of the drawing certain pairs of half-edges make an angle at the junction between neighbouring grey squares, but such a pair should of course just be considered a single edge. The reader may verify that the lattice does indeed consist of pentagons, and that going around each pentagon the degrees of the vertices are 3, 3, 4, 3 and 4 as they should be. There are on average 3/2 vertices per grey square, and we must take n even to respect the alternation of patterns.

Zoom In Zoom Out Reset image size

The $\check{\sf R}$-matrix is given by that of the square lattice, equation (75), when x + y is odd; by that of the hexagonal lattice, equation (73), when x and y are both even; and by a rotated version of the latter when x and y are both odd.

Results for the thresholds p_c are given in table 56.

A four-terminal representation of the frieze dual lattice is shown in figure 47. Its $\check{\sf R}$-matrix is that of the square lattice (equation (75)) when y is even; and that of the rotated hexagonal lattice (cf equation (73)) when y is odd. All faces are pentagons, and the reader may verify from the figure that the degrees of the surrounding vertices are (3³, 4²) as they should be. There are on average 3/2 vertices per grey square. Any parity of n defines a valid basis, provided that we use n × 2n rectangular bases in order to respect the alternation between rows.

Zoom In Zoom Out Reset image size

The percolation thresholds p_c are displayed in table 57.