Potts-model critical manifolds revisited

The computation of critical thresholds in percolation and the Potts model has long been the domain of Monte Carlo methods [2–4] or transfer matrix techniques [5, 6] of similar accuracy. Recently, we developed a radically different technique called the method of critical polynomials [7–13] which gives us the ability to compute percolation and Potts-model thresholds to precisions that are orders of magnitude greater than that obtained with traditional tools. The idea is based on the conjecture that the roots of a particular graph polynomial, to be defined below, provide estimates for critical points that become more accurate by increasing the size of the graph on which it is computed. This conjecture is essentially equivalent [12] with the universality conjecture for crossing probabilities [14], so it is probably not in serious doubt (at least for percolation). However, computing the polynomial on large graphs is a significant challenge for which, at present, the transfer matrix appears to be the best tool.

A recent advance in this area was made in [1], which allowed the polynomial to be calculated on graphs of hundreds of edges in some cases, producing estimates of Potts and percolation thresholds far exceeding that possible with traditional techniques, such as Monte Carlo. It also allows us to study the $v\lt -1$ regime of the Potts model on various lattices, something that was very difficult previously. Here we implement the algorithm of [1] in parallel on the Cab and Vulcan supercomputers at Lawrence Livermore National Laboratory, improving the resolution of the phase diagrams and the accuracy of the percolation thresholds found in that work. Thus, this paper should be considered an addendum to [1] and, although we will give a brief description of the critical polynomial method, the reader interested in learning the details of the transfer matrix algorithm is referred there.

We consider here the Potts model [15] in the Fortuin–Kasteleyn cluster representation [16]. On a graph $G=(V,E)$ composed of vertices V and edges E, the partition function is given by

$$\begin{eqnarray}&&Z=\displaystyle \sum _{A\subseteq E}{v}^{| A| }{q}^{k(A)},\end{eqnarray} \tag{ 1 }$$

where the sum is over all the subsets, A, of E, $| A|$ is the number of edges in A and k(A) is the number of connected components including isolated vertices. Each term in the sum can be thought of as a configuration in which each edge contributes a factor of v if it is present in the configuration (or 'open') and a factor of one otherwise. The sum (1) restricted to configurations consistent with the occurence of a particular event is the weight of that event. On a finite graph, weights are polynomials in q and v.

To compute the critical polynomial on a given lattice, we must first choose a basis, B. This is a subgraph of the lattice along with a specification of how copies of B are arranged to form the infinite lattice. This is defined by identifications between the vertices on the boundary, which we call terminals. Once B is chosen, we compute the weights of two events; 2D and 0D. To understand what these are, consider an arbitrary configuration on B. If, after B has been tiled to form the infinite lattice with each copy of B containing the same configuration, every copy of B is connected by open edges, we say that the event 2D occurs. Likewise, if we can traverse an infinite number of copies of B but some remain unconnected, then we have the event 1D. Finally, if one cannot get anywhere on open edges, we have 0D. The weights of these events are denoted ${P}_{2{\rm{D}}}(q,v)$, ${P}_{1{\rm{D}}}(q,v)$ and ${P}_{0{\rm{D}}}(q,v)$. These three possibilities are shown in figure 1. The critical polynomial is then given by

$$\begin{eqnarray}&&{P}_{B}(q,v)={P}_{2{\rm{D}}}(q,v)-{{qP}}_{0{\rm{D}}}(q,v).\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

The root of this polynomial provides the estimates for the critical point, v_c(q). In the case of a lattice on which the critical point is known exactly [17–19], the root of the polynomial provides the exact answer. Indeed, this method was motivated by the observation that all known exact solutions are roots of polynomials, and our graph polynomial always factorises in such cases, for any size of B, shedding a factor which corresponds to the polynomial of the exact solution.

On the other hand, for unsolved problems, the roots of the critical polynomial are approximations that become more accurate with the size of B. The estimates are generally very good, even for bases of moderate size, and they get better as we increase the number of edges in B. Using an extrapolation technique allows us to locate critical points with unprecedented accuracy. Although the fact that the roots of the polynomials should approach the exact critical point as the basis approaches the infinite lattice follows from universality, as argued in [12], the rapidity of this convergence, especially when compared with standard techniques, is very surprising.

Using a transfer matrix algorithm, one of us [1] computed the critical polynomials in (q, v) of all the Archimedean lattices for square bases of size 5 × 5, and computed individual percolation and Potts thresholds for 7 × 7. An example of such a basis is shown in figure 2 and, as explained in [1], all the Archimedean lattices can be handled with such bases. The accuracy of these estimates is orders of magnitude greater than that obtained with standard techniques, such as Monte Carlo. In this paper, we implement a parallel version of the algorithm in [1] to push the basis size to n = 8 for percolation thresholds and to n = 6 for phase diagrams of all the Archimedean lattices. Many interesting features of these phase diagrams become apparent when we examine them for various n and it is surely safe to say that at this point our ability to numerically compute these plots far outstrips our theoretical understanding of them.

Zoom In Zoom Out Reset image size

The most efficient way to compute the critical polynomial is to use a transfer matrix algorithm [12, 13]. Our earlier transfer matrix computations of critical polynomials would compute the weights of the 0D, 2D and 1D configurations and we could then set ${P}_{2{\rm{D}}}={P}_{0{\rm{D}}}$. However, this is wasteful because P_1D is never used for anything. The trick then would seem to be efficiently generating the weights while discarding the 1D configurations on the fly. This problem is partially solved by Jacobsen's algorithm [1] by recasting the problem in terms of a loop model on square bases that are periodic in the x-direction (figure 2). At the same time the number of terminals is effectively reduced, leading to further significant gains in performance.

In figure 2 each grey square can contain an arbitrary arrangement of vertices and edges, provided that the interaction with the remainder of the basis takes place only through the four vertices at the corners. The transfer matrix contains an operator ${\check{{\mathsf{R}}}}_{i}$ materialising the propagation through each grey square at position x = i, and ${\check{{\mathsf{R}}}}_{i}$ factorises as a product of simpler operators, ${{\mathsf{H}}}_{i}$ and ${{\mathsf{V}}}_{i}$, responsible for the addition of one 'horizontal' or 'vertical' edge to the lattice. Finally, ${{\mathsf{H}}}_{i}$ and ${{\mathsf{V}}}_{i}$ are simply expressed in terms of elementary operators ${{\mathsf{E}}}_{i}$ which are the generators of a periodic Temperley–Lieb algebra [20]. The relevant factorised expressions for ${\check{{\mathsf{R}}}}_{i}$ corresponding to any Archimedean lattice have been given in [1] and below we shall just use those expressions, with however one significant improvement for the cross lattice (see section 7.5).

This setup just outlined allows us to discard 1D configurations that wrap through the periodic direction. We still end up computing some 1D weights but the size of the transfer matrix is greatly reduced over the more naïve approach in [12, 13]. Again, we refer the reader to [1] for a detailed description of this improved transfer matrix computation.

The eleven Archimedean lattices are shown in figure 3. By definition, an Archimedean lattice is such that each vertex is surrounded by the same types of faces, appearing in the same cyclic order.

Zoom In Zoom Out Reset image size

The Potts model is exactly solvable on the square, triangular and hexagonal lattices. We henceforth focus on the remaining eight cases which are unsolved.

The algorithm used in this work is a parallelised version of that described in detail in [1]. At its heart, the computation is the multiplication of a vector of states by a matrix, and the parallel strategy is simply to distribute the vector among the processors such that each contains an equal number of states. Then if the result of an elementary ${{\mathsf{E}}}_{i}$ operation produces a state that belongs on a different processor, that state's weight is sent to its appropriate destination. Unfortunately, there is no way to partition the states such that no communication is needed between the processors, so the goal of a good parallelisation scheme is to minimise the amount of inter-processor communication.

To achieve this goal we modify slightly the ordering of the states from the one given in [1]. Here we take care to ensure that when implementing the bijection between connectivity states and integers, the reduced state describing the bottom row is the most slowly varying⁶ . This ensures that blocks of states with the same bottom-row connectivity will tend to be placed on the same processor. Because the action of the transfer matrix is implemented with ${{\mathsf{H}}}_{i}$ and ${{\mathsf{V}}}_{i}$ operators acting on the top row, the bottom row is often unchanged upon performing these operations. Even so, the connectivity on the bottom row might still change if two strings are attached⁷ , but this is rare enough that ordering the states in the way we describe significantly cuts down on the inter-processor communication. Beyond this, we take no other care in distributing states across tasks. If the system has N total states, and we use p processors, p is chosen as a divisor of N and the first processor gets the first N/p states, the second the next N/p etc.

In the serial code, n = 5 was the maximum size accessible for computing full Potts-model polynomials in (q, v). In [1], the transfer matrix was also used to compute individual thresholds for q = 1, 2 and 3 using a Newton–Raphson technique and n = 7 was the limit for computing these. Using the parallel algorithm we push these limits to n = 6 for the full polynomials and n = 8 for q = 1 (percolation) critical points, giving us the clearest look yet at the phase diagrams of all the Archimedean lattices and earning typically an extra digit of accuracy in the percolation thresholds.

To compute the full polynomials, whose coefficients are very large integers, we perform the calculation modulo a set of prime numbers. The results for each prime are then combined into the final answer using the Chinese remainder theorem. However, to compute just the percolation thresholds for n = 8, we use the real numbers provided by the CLN arbitrary precision library for C++ [22]. CLN objects are equipped with functions for writing numbers to and from strings, and to communicate between processors we send the string and convert it to a CLN number at the destination. This adds somewhat to the usual communication cost.

Although the transfer matrix has the same dimension for a given n, independent of the lattice, there is significant variation in the numbers of edges and vertices in each problem. Therefore, computing the full (q, v) polynomial requires more resources, in terms of both time and memory, for some Archimedean lattices than others, with the $(3,{12}^{2})$ and cross lattices being the most difficult. The number of processors we chose for each problem is governed both by memory and speed considerations.

We performed most runs on the Cab supercomputer at Lawrence Livermore National Laboratory. Each processor is 2.6 GHz with 2 gigabytes of memory, and we use a maximum of 4092 of these for each computation. Vulcan has 1.6 GHz processors with 1 gigabyte of RAM each and is designed for massively parallel runs. As such, we use it only to compute the phase diagram of the cross lattice, the most difficult problem we address, for which we use 30690 processors. For the n = 8 percolation thresholds we used 4004 processors in all cases, with the time the runs took varying between 6.5 to about 16 h depending on the lattice and the quality of the initial guess. For the full phase diagrams, the number of processors and computation time varied considerably and these are given in section 7. All told, we used about 10 million CPU h for the entire project.

In [1] it was shown how each of the Archimedean lattices can be represented in the form shown in figure 2, enabling a transfer matrix calculation of the graph polynomial ${P}_{B}(q,v)$ by specifying a corresponding $\check{{\mathsf{R}}}$-matrix. We use here the same bases, except for the cross lattice for which we have now found a better representation containing twice as many vertices and edges for any given n.

Table 1.  Nomenclature of the Archimedean lattices and their duals (Laves lattices). The labels (d)–(k) refer to the eight lattices (see figure 3) on which the q-state Potts model is unsolved. The notation is that of Grünbaum and Shephard [21]. The right part of the table gives the number of vertices and edges for square bases of size n × n grey squares. In addition we state any parity constraints on n. Note that the cross lattice comprises an improvement over ([1] table 3).

	Lattice	Notation	Vertices	Edges	Parity of n
(d)	Kagome	$(3,6,3,6)$	$3{n}^{2}$	$6{n}^{2}$	Any
(e)	Four-eight	$(4,{8}^{2})$	$4{n}^{2}$	$6{n}^{2}$	Any
(f)	Frieze	$({3}^{3},{4}^{2})$	$2{n}^{2}$	$5{n}^{2}$	Any
(g)	Three-twelve	$(3,{12}^{2})$	$6{n}^{2}$	$9{n}^{2}$	Any
(h)	Cross	$(4,6,12)$	$6{n}^{2}$	$9{n}^{2}$	Even
(i)	Snub square	$({3}^{2},4,3,4)$	$2{n}^{2}$	$5{n}^{2}$	Even
(j)	Snub hexagonal	$({3}^{4},6)$	$3{n}^{3}$	$\frac{15}{2}{n}^{2}$	Even
(k)	Ruby	$(3,4,6,4)$	$3{n}^{2}$	$6{n}^{2}$	Even

For q = 1, the root ${v}_{{\rm{c}}}$ of ${P}_{B}(1,v)$ that approximates the percolation threshold ${p}_{{\rm{c}}}$, via the relation $p=v/(1+v)$, was obtained in [1] for sizes $n=1,2,\ldots ,7$. For each of the eight unsolved Archimedean lattices we have extended this computation to n = 8; see table 2.

Table 2.  Bond percolation thresholds ${p}_{{\rm{c}}}(n)$ for n × n bases of size n = 8 on the non-solvable Archimedean lattices.

Lattice	${p}_{{\rm{c}}}$
Kagome	0.52440499951414108517718152880414539618186347478997
Four-eight	0.67680315546810921539888642691310355045728163134828
Frieze	0.41964027675236652348940526360240730997388408890937
Three-twelve	0.74042079889005115125415547646085000986933464791292
Cross	0.69373316681403389479582584979836234112904294069266
Snub square	0.41413788425992466879112356169148391265046003220755
Snub hexagonal	0.43432836999043911236626317783310197323151394980874
Ruby	0.52483149832115089795511006350993499311132979786833

For each lattice the behaviour of ${p}_{{\rm{c}}}(n)$ appears to have the form

$$\begin{eqnarray}&&{p}_{{\rm{c}}}(n)={p}_{{\rm{c}}}+{{An}}^{-w}.\end{eqnarray} \tag{ 3 }$$

Approximate values of the exponent w can be obtained from nonlinear three-point fits to the data for each lattice. These estimates for w are shown in table 3. The most reliable values are obtained for the lattice for which there is no parity constraint on n (see table 1) and in those cases we provide an error bar in table 3 which is obtained by comparing the three-point fits for $n\in \{6,7,8\}$ and for $n\in \{5,6,7\}$. In the cases with parity constraint, we just give the value for the fit with $n\in \{4,6,8\}$ to one decimal place.

Table 3.  Estimates for the finite-size scaling exponent w of (3), and final results for the percolation thresholds ${p}_{{\rm{c}}}$, for bond percolation on the non-solvable Archimedean lattices.

Lattice	w	${p}_{{\rm{c}}}$
Kagome	6.346 (2)	0.52440499917 (1)
Four-eight	4.20 (8)	0.676803125 (2)
Frieze	3.5	0.41964037 (2)
Three-twelve	6.392 (1)	0.7404207988509 (4)
Cross	4.7	0.69373295 (9)
Snub square	4.1	0.414137858 (1)
Snub hexagonal	5.1	0.43432830 (2)
Ruby	4.9	0.52483144 (2)

It appears that the lattices with three-fold rotational symmetry (kagome, three-twelve, cross, snub hexagonal, and ruby) might have a common value of w that we can estimate as⁸

$$\begin{eqnarray}&&{w}_{3}=6.2\pm 0.2\end{eqnarray} \tag{ 4 }$$

independently of the lattice. The remaining lattices which have four-fold symmetry (four-eight, snub square) or two-fold symmetry (frieze) appear to share another value of w that we estimate as

$$\begin{eqnarray}&&{w}_{4}=4.2\pm 0.2.\end{eqnarray} \tag{ 5 }$$

To extrapolate ${p}_{{\rm{c}}}(n)$ to the thermodynamic limit, $n\to \infty$, we use Monroe's implementation [23] of the Bulirsch–Stoer (BS) extrapolation scheme [24]. This assumes that the data has the form (3) and we must provide the corresponding value of w. It goes without saying that this procedure would gain substantially in precision and reliability if there was an analytic prediction for w. In the absence of any such result, we perform for each lattice the extrapolations using both the w of table 3 and the estimates w₃ or w₄ (depending on the symmetry class⁹ ). The final results and error bars on ${p}_{{\rm{c}}}$, obtained by comparing the last few BS extrapolants for both values of w, are displayed in the right column of table 3.

We should stress that this extrapolation procedure is somewhat more cautious than the one used in [1] and accordingly produces more conservative error bars for series of the same length. The fact that we have now one more data point for each lattice obviously makes the thresholds of table 3 more reliable and/or more precise than those in [1].

The full polynomial ${P}_{B}(q,v)$ was obtained in symbolic form for $n=1,2,\ldots ,5$ in [1]. We have extended these computations to n = 6. As in our preceding work [1, 12, 13] the polynomials are available in electronic form as supplementary material to this paper (available at stacks.iop.org/jpa/49/125003/mmedia)¹⁰ .

The zero set ${P}_{B}(q,v)=0$ in the real (q, v) plane was shown previously [1, 11–13] to provide precise information about the phase diagram of the corresponding Potts models. In the next sections we show these diagrams for the unsolved Archimedean lattices.

Actually, even for solved problems, the set of solvable curves does not give access to all relevant information in the antiferromagnetic region $v\lt 0$. It was argued in [11] that the roots of the critical polynomial will produce extra curves inside the so-called Berker–Kadanoff phase [25] which will form, in the thermodynamic limit, vertical rays in the (q, v) plane at the Beraha numbers

$$\begin{eqnarray}&&q={B}_{k}\equiv 4{\mathrm{cos}}^{2}\displaystyle \frac{\pi }{k}\end{eqnarray} \tag{ 6 }$$

for certain integer k. On the basis of conformal field theory results and cognate studies of partition function zeros [26, 27], as well as actual evidence from critical polynomials, it was conjectured in [13] that such vertical rays will exist for any even $k=4,6,8,\ldots$, as far as allowed by the extent of the Berker–Kadanoff phase.

For the case of the square lattice the Berker–Kadanoff phase extends all the way to $(q,v)=(4,-2)$, and explicit calculations of ${P}_{B}(q,v)$ with n up to five [1] fully confirmed the validity of the conjecture. Accordingly, we see no reason to compute this polynomial also for n = 6. However, vertical rays in a Berker–Kadanoff-like phase appear to be a generic feature, present in all the Archimedean lattices, and the number and position of such rays permit us to shed light on the extent of the Berker–Kadanoff phase for each lattice. The n = 6 critical polynomials discussed below will extend the results of [1] in this respect.

We should stress that the physics of the antiferromagnetic regime is very difficult to access using standard numerical techniques, in particular for the non-probabilistic regime $v\lt -1$, but critical polynomials capture this regime very well.

7.1. Kagome lattice $(3,6,3,6)$

To get the n = 6 critical polynomial, we needed to perform the calculation with 13 different primes, with each prime taking about 5 h on 3861 processors. The phase diagram for this lattice is shown in figure 4 along with those from the smaller bases already reported in [1]. The n = 6 curve is very similar to the union of the n = 3, 4 and 5 curves, probably indicating that our n = 6 result is indeed an accurate picture of the true phase diagram. The only vertical rays to appear are at q = 2 and 3 and the n = 6 polynomial does not give a hint that any further rays are going to appear. Also, the Berker–Kadanoff phase appears as though it may not extend to $(q,v)=(4,-2)$ as it does for most Archimedean lattices. A curious persistent feature is the small bump that rises to the left of the q = 2 ray. As for most of the features in this and the subsequent phase diagrams, there is not yet any analytical explanation for this.

Zoom In Zoom Out Reset image size

Close scrutiny of figure 4 reveals a number of parity effects, confirming and extending observations made in [1, 11, 13]. For instance, the lower boundary of the Berker–Kadanoff phase—the curve emanating from $(q,v)=(0,-3)$ and going towards $\approx (3,-2)$—appears only for even n, and a small part of the upper boundary is visible only for odd n. Maybe the most significant feature of the n = 6 plot is that it confirms the existence of a bubble to the right of q = 3, previously seen only for n = 3. We conjecture that it will appear when n is any multiple of 3.

7.2. Four-eight lattice $(4,{8}^{2})$

Although the unit cell of the $(4,{8}^{2})$ lattice has the same number of edges as the kagome lattice, it has more vertices and therefore the n = 6 polynomial is slightly more difficult to compute due to its higher order in q (144 versus 108). To do the computation with one prime took around 7 h on 4092 processors and we used 14 primes. The phase diagram arising from the n = 6 polynomial is shown in figure 5 along with those for n up to five reported previously in [1]. Unlike in the kagome case, we have vertical rays here at several Beraha numbers (${B}_{4},{B}_{6},\ldots ,{B}_{12}$) and it certainly seems plausible that more would appear if we could compute higher order polynomials. In fact, this phase diagram is very similar to that of the square lattice, which is perhaps not surprising. Here, though, we have some additional features such as the finger-like structures that emanate from the lower left of the diagram. The line that exits through the lower right (for even n only) is similarly not present in the square lattice.

Zoom In Zoom Out Reset image size

7.3. Frieze lattice $({3}^{3},{4}^{2})$

For the frieze lattice we used 4092 processors, and completed a computation for a single prime in about 3 h. We used 14 primes, but 13 were sufficient to get the final answer. For this problem, the lattice only fits into our basis for even n. The n = 6 phase diagram is plotted in figure 6 along with the n = 2 and 4 results from [1]. In this and subsequent cases, where only even n are possible, the n = 6 polynomial adds significant detail to the phase diagram. In particular we observe two more vertical rays at the Beraha numbers B₁₀ and B₁₂. Note that there seem to be gaps in the boundary of the Berker–Kadanoff phase for all n between q = 2 and 3. These are probably not a real feature but an effect of the basis we are using. Doing the calculation for a different basis, such as one in which we shift the identifications between the top and bottom rows, may fill in these gaps.

Zoom In Zoom Out Reset image size

7.4. Three-twelve lattice $(3,{12}^{2})$

The n = 6 polynomial for the $(3,{12}^{2})$ lattice is order 324 in v and 216 in q. This makes the computation difficult both in terms of the number of primes needed, 22, and the time it takes to do a calculation for a single prime on 4092 processors, 9 h. However, we are rewarded for this calculation with the very intricate phase diagram shown in figure 7 (along with those from previous work [1, 13]).

Zoom In Zoom Out Reset image size

Although this lattice bears some similarity to the kagome lattice, aside from the small bump (difficult to see at the scale of the plot) rising off the q = 2 ray the phase diagram is quite different. The Berker–Kadanoff phase clearly extends to $(q,v)=(4,-2)$ and has vertical rays at at least four Beraha numbers (${B}_{4},{B}_{6},{B}_{8},{B}_{10}$). But closer scrutiny of figure 7 reveals emergent features at five more Beraha numbers (${B}_{12},{B}_{14},\ldots ,{B}_{20}$) that would presumably turn into fully fledged rays if further sizes $n\gt 6$ were available. We also have the fingers that emerge from the lower left, while the region outside the Berker–Kadanoff phase for $q\gt 3$ seems to be a rather wild tangle.

7.5. Cross lattice $(4,6,12)$

For the cross lattice we have found a four-terminal representation which for a given n contains twice the number of vertices and edges than the representation used in [1]. This is accomplished by choosing the following $\check{{\mathsf{R}}}$-matrix

$$\begin{eqnarray}&&{\check{{\mathsf{R}}}}_{i}={{\mathsf{V}}}_{i+2}{{\mathsf{H}}}_{i+1}{{\mathsf{V}}}_{i}{{\mathsf{V}}}_{i+2}{{\mathsf{H}}}_{i+1},\end{eqnarray} \tag{ 7 }$$

whose action is illustrated in the grey squares of figure 8. In addition we use the trick that was dubbed 'white hexagons' in ([1] section 4.9) to insert extra structure in-between two rows of the basis. The corresponding operator reads

$$\begin{eqnarray}&&{{\mathsf{O}}}_{i}={({{\mathsf{V}}}_{i}{{\mathsf{V}}}_{i+2}{{\mathsf{H}}}_{i+1})}^{2}{{\mathsf{V}}}_{i}{{\mathsf{V}}}_{i+2}\end{eqnarray} \tag{ 8 }$$

and its action is illustrated in the pink hexagons (dubbed 'white' in [1]) of figure 8.

Zoom In Zoom Out Reset image size

Computing this polynomial required the most effort; and like the frieze lattice, the cross lattice only fits into bases with even n. In this one case we used the Vulcan supercomputer. We used 30690 processors, for which it took approximately 10 h to compute a single prime. Processors on this machine are slower than on Cab with less memory and we needed so many primarily to get the entire vector to fit in memory. The phase diagram is shown in figure 9; the curves with n = 2 and n = 4 this time do not correspond to those computed in [1] but represent the new construction (7) applied to those two cases. We can see vertical rays at the Beraha numbers, with n = 6 providing two more (B₁₀ and B₁₂) above n = 4. We also have the fingers that enter through the lower left and the line that exits through the lower right. Overall, the phase diagram is qualitatively very similar to that of the $(4,{8}^{2})$ lattice.

Zoom In Zoom Out Reset image size

7.6. Snub square lattice $({3}^{2},4,3,4)$

Calculating the n = 6 critical polynomial for the snub square lattice required 11 primes. We used 2046 processors with which it took about 3 h for each prime. The phase diagram is shown in figure 10. Along with the usual vertical rays at the even Beraha numbers, there is an oval-shaped region straddling the one at B₈ near $v\approx -1.5$. It appears that the Berker–Kadanoff phase may not extend to $(q,v)=(4,-2)$.

Zoom In Zoom Out Reset image size

7.7. Snub hexagonal lattice $({3}^{4},6)$

To do the calculation of the snub hexagonal polynomial for a single prime required about 7 h on 4092 processors. We used 17 primes. The phase diagram is shown in figure 11. As for other lattices, we can only find the polynomial for even n, and there is a gap in the boundary of the Berker–Kadanoff phase between q = 2 and 3. This phase appears to reach $(q,v)=(4,-2)$ where there is an interesting flower-like structure in the n = 4 and 6 polynomials.

Zoom In Zoom Out Reset image size

7.8. Ruby lattice $(3,4,6,4)$

For the ruby lattice, we used 2046 processors; it took about 5 h to do the computation for a single prime, and we used 14 primes. Once again, only even n works for this problem. The phase diagrams are shown in figure 12. Our n = 6 calculation adds significant detail, such as the two rays at B₁₀ and B₁₂. The Berker–Kadanoff phase likely extends to $(q,v)=(4,-2)$.

Zoom In Zoom Out Reset image size

We have implemented a parallel version of the improved transfer matrix algorithm given in [1] for computing critical polynomials. We have computed phase diagrams for bases of size 6 × 6 and found percolation thresholds for 8 × 8 bases. Using Bulirsch–Stoer extrapolation we have found the latter to accuracies in some cases far exceeding that possible with standard techniques; see table 3. Even with large parallel calculations, as we have done, Monte Carlo and similar techniques would obtain nowhere near the accuracy found here for our best results (references to other work were given in the tables of [1]). On the other hand, there are still many non-Archimedean 2D problems for which it is not yet clear how to better Monte Carlo with critical polynomials (see, e.g., [28]).

For the phase diagrams, we are given a view of the $v\lt -1$ regime of the Potts model that is accessible only with critical polynomials, and in this work we have the most detailed picture yet for the Archimedean lattices. However, only the square lattice phase diagram is properly understood analytically [1, 11, 25], and, apart from some qualitative similarities with that case, the many interesting features in these plots are yet to be explained.

Aside from these issues, many interesting theoretical questions remain to be addressed about the method itself. For example, it would be very useful to know the exponent w in the scaling law (3) exactly. It seems clear that, unlike most other exponents that arise in these models, w is not universal; it takes a value around six for some lattices and four for others, possibly depending on the degree of symmetry. Gaining some understanding of this exponent may help to answer the broader question of why this method performs as well as it does.

Another outstanding problem concerns the relation between the phase diagrams obtained in the antiferromagnetic region from, respectively, the critical polynomial method and partition function zeros. It is conspicuous that the roots of ${P}_{B}(q,v)$ invariably accumulate at the vertical rays $q={B}_{k}$ with even k, just like partition function zeros for strip graphs with cyclic boundary conditions [26]. On the other hand, partition function zeros for strip graphs with toroidal boundary conditions appear to accumulate instead at integer q [27]. This is puzzling, since the critical polynomial is esentially a partition function on a torus, albeit with topological weights.

CRS acknowledges the helpful advice of Jim Glosli at LLNL in programming the parallel implementation. This work was partially performed under the auspices of the U S Department of Energy at the Lawrence Livermore National Laboratory under Contract No DE-AC52-07NA27344. The research of JLJ is supported by the Agence Nationale de la Recherche (grant ANR-10-BLAN-0414: DIME), the European Research Council (advanced grant NuCFT), and the Institut Universitaire de France.