Phase diagram of the chromatic polynomial on a torus

We study the zero-temperature partition function of the Potts antiferromagnet (i.e., the chromatic polynomial) on a torus using a transfer-matrix approach. We consider square- and triangular-lattice strips with fixed width L, arbitrary length N, and fully periodic boundary conditions. On the mathematical side, we obtain exact expressions for the chromatic polynomial of widths L=5,6,7 for the square and triangular lattices. On the physical side, we obtain the exact ``phase diagrams'' for these strips of width L and infinite length, and from these results we extract useful information about the infinite-volume phase diagram of this model: in particular, the number and position of the different phases.


Introduction
The two-dimensional (2D) q-state Potts model [1] is one of the most studied models in Statistical Mechanics. Despite many efforts over more than 50 years, its exact free energy and phase diagram are still unknown. The ferromagnetic regime of the Potts model is the best understood case: exact (albeit not always rigorous) results have been obtained for the ferromagnetic-paramagnetic phase transition temperature T c (q) (at least for several regular lattices), the order of the transition (continuous for 0 ≤ q ≤ 4), the phase diagram, and the characterization in terms of conformal field theory (CFT) of the universality classes. (See e.g., Ref. [2].) The antiferromagnetic regime is less understood, partly because universality is not expected to hold in general (in contrast with the ferromagnetic regime). We know the exact free energy along some curves of the phase diagram (q, T ) (where T is the temperature) for some regular lattices [2]; but in this regime, this (partial) solubility of the model does not imply criticality (as it is been observed for the ferromagnetic regime with 0 ≤ q ≤ 4).
The standard q-state Potts model at temperature T can be defined on any undirected finite graph G = (V, E) as follows: on every vertex of the graph x ∈ V , we place a spin σ x taking an integer value in the set {1, 2, . . . , q}, where q ∈ N. The spins interact via the Hamiltonian where the sum is over all edges e = x, y ∈ E joining nearest-neighbor sites, J is the coupling constant, and δ is the Kronecker delta. The ferromagnetic (resp. antiferromagnetic) regime corresponds to J > 0 (resp. J < 0). The partition function of the Potts model on a graph G at inverse temperature β = 1/T ≥ 0 can then be written as Potts (1.2) A very useful representation of the Potts model is due to Fortuin and Kasteleyn [3]. They showed that the partition function (1.2) can be rewritten as where v = e Jβ − 1, the sum is over the 2 |E| spanning subgraphs (V, E ′ ) of G, and k(E ′ ) is the number of connected components (including isolated vertices) of the spanning subgraph (V, E ′ ). As a result, the partition function (1.3) is a polynomial in both variables q and v. Thus, we can promote (1.3) to the definition of the model, which allows us to consider q as an arbitrary complex number. We shall here focus on the zero-temperature antiferromagnetic case J = −∞ (i.e., v = −1), whose restriction to q ∈ N can be interpreted as a coloring problem. In this case, P G (q) = Z G (q, −1) is known as the chromatic polynomial of G. Its evaluation for a general graph is a hard problem. More precisely, Oxley and Welsh have proven that the computation of the coefficients a k of the chromatic polynomial P G (q) = |V | k=1 a k q k of a general graph G (including bipartite graphs) is #-P hard [4].
In this paper, we will consider only strip graphs of certain regular lattices (i.e., square and triangular). For a family of strip graphs {G n } (or, more generally, recursive families of graphs [5]), the chromatic polynomial P Gn (q) of any member of the family G n can be computed from a transfer matrix T and certain boundary condition vectors u and v: (1.4) Thus, the computation time grows as a polynomial in n for fixed strip width L. However, it still grows exponentially in the strip width (due to the L dependence of dim T). It is therefore of interest to devise new and efficient algorithms to be able to handle as large widths L as possible.
Remark. It is important to stress that the Potts spin model has a probabilistic interpretation (i.e., has non-negative Boltzmann weights) only when q is a positive integer and v ≥ −1 (i.e., positive temperature). The Fortuin-Kasteleyn random-cluster model (1.3), which extends the Potts model to non-integer q, has non-negative weights only when q, v ≥ 0. In all other cases, the model belongs to the "unphysical" regime (i.e., some weights are negative or complex), and some familiar properties of statistical mechanics need not hold. For example, even for integer q ≥ 2 and real v in the antiferromagnetic range −1 ≤ v < 0, where the spin representation exists and has non-negative weights, the dominant eigenvalue λ ⋆ of the transfer matrix in the cluster representation need not be simple (i.e., the eigenvector may not be unique); and even if simple, it may not play any role in determining Z Gn (q, v) because the corresponding amplitude may vanish. Both these behaviors are of course impossible for any transfer matrix with positive weights, by virtue of the Perron-Frobenius theorem. It is important to note that the dominant eigenvalues in the cluster and spin representations need not be equal, because the former may have a vanishing amplitude (see below). In any case, the Potts model can thus make probabilistic sense in the cluster representation at parameter values where it fails to make probabilistic sense in the spin representation, and vice versa. It is worth mentioning that for q = 4 cos 2 (π/n) with integer n ≥ 3 and planar graphs G, the Potts-model partition function Z G (q, v) admits a third representation, in terms of a restricted solid-on-solid (RSOS) model [6,7] (See also Ref. [8] for an recent study with a more comprehensive list of references). Even though the zero-temperature limit of the Potts antiferromagnet may not have a probabilistic interpretation in the Fortuin-Kasteleyn representation [3] for non-integer q, it can nevertheless being studied with the standard tools of Statistical Mechanics and CFT with appropriate modifications. In particular, we expect that critical "unphysical regions" can be described with non-unitary field theories. (See [9,10] for a more detailed discussion on these points.) It is well established that antiferromagnetic models in general, and the chromatic polynomial in particular, are very sensitive to the choice of boundary conditions. Indeed, different choices may lead to quite different thermodynamic limits. In our earlier publications [11][12][13] we studied free and cylindrical boundary conditions: free boundary conditions in the longitudinal direction and free or periodic boundary conditions in the transverse direction, respectively. More recently [9], we considered cyclic boundary conditions: namely, free boundary conditions in the transverse direction, and periodic boundary conditions in the longitudinal direction. In particular, we observed that in the latter case, the phase diagram of the triangular-lattice chromatic polynomial was more involved than for free or cylindrical boundary conditions. In this paper we shall deal with fully periodic boundary conditions (i.e., periodic in both directions). Periodic boundary conditions in the longitudinal direction (i.e., toroidal and cyclic) can be easily implemented in the spin representation as P Gn = tr T n . In the Fortuin-Kasteleyn representation, this implementation is less obvious as shown in Ref. [9]. For the triangular-lattice model with periodic boundary conditions on the transverse direction we should also use some additional tricks [11,13,14].
The results existing in the literature about exact chromatic polynomials for strips of the square and triangular lattices with toroidal boundary conditions are limited to L ≤ 4. The earlier results were obtained by recursive use of the contraction-deletion theorem. These works include the pioneering paper by Biggs, Damerell, and Sands [15] (where the solution for a square-lattice strip of width L = 2 was presented), and those of Chang and Shrock [16][17][18]. In Ref. [17] the latter authors gave the exact solution for the full Potts-model partition function (i.e., with arbitrary v) for square-lattice strips of widths L = 2, 3. In Ref. [18] the solution for the triangular-lattice chromatic polynomial for width L = 3 was presented, while in Ref. [16] they gave the solution for both lattice strips of width L = 4. Finally, in 2006 Chang and Shrock exhibited a transfer-matrix formalism to compute the full Pottsmodel partition function for regular lattice strips with toroidal boundary conditions [19]. 1 In particular, they gave the exact partition function for square-lattice strips of widths L ≤ 4, and triangular-lattice strips of widths L = 2, 3. 2 In Ref. [19], they also obtained structural properties of the partition function and the chromatic polynomial. Additional structural properties were discussed by one of us [20].
We have several motivations for this work. As explained above, the exact solution of the 2D q = 2 state Potts model at general v is still unknown, even for the square lattice (and a fortiori for more complicated regular lattices), and even in the thermodynamic limit (and a fortiori on general finite lattices). Therefore, exact solutions of this model on strips of finite fixed width and arbitrary length are valuable in their own right (even in the particular case v = −1). In addition, the chromatic polynomial P G (q) is an object of great interest to mathematicians (see e.g., Ref. [21] for a recent survey). As the computation of the chromatic polynomial is #-P hard [4], new and efficient methods for computing this object on large graphs are also of great interest to this community.
The second motivation of this work is to deepen our physical understanding of the critical points of the q-state Potts model in the antiferromagnetic −1 ≤ v < 0 and unphysical regimes v < −1. Indeed, we can extract this information from our finite-width strips using CFT and finite-size scaling (FSS) [22][23][24]. In particular, we will consider the following way to attain the thermodynamic limit. From the chromatic polynomial of a strip of size L × N we can compute its chromatic zeros in the complex q-plane. From these zeros it is very hard to obtain infinite-volume quantities, as they have strong FSS corrections. Thus, we consider the infinite-length limit N → ∞. Then, the chromatic zeros accumulate along certain limiting curves B L and around isolated limiting points. This phenomenon is a consequence of the Beraha-Kahane-Weiss theorem [25]. Our methods allow us to obtain the exact values of the relevant physical quantities in this limit. Their FSS corrections are found to be smaller than when both L and N are finite. The extrapolation to the true infinite-volume limit L → ∞, using standard FSS techniques, therefore becomes more precise by employing this order of limits.
Remark. Toroidal boundary conditions are rather special compared to the other boundary conditions studied in previous papers: namely, free, cylindrical, and cyclic [9,[11][12][13]. Strip graphs with the latter boundary conditions are planar; thus, the four-color theorem [26] applies. However, strip graphs with the former boundary conditions can be regarded as graphs embedded on a torus. In this case, we have only a seven-color theorem. This follows from the upper bound obtained by Heawood [27] in 1890 on the chromatic number of a graph G embedded on an orientable surface of genus g ≥ 1 χ(G) ≤ H(g) = 1 2 7 + 1 + 48g (1.5) combined with the fact that K 7 (with χ(K 7 ) = 7) can be embedded on a torus, as shown in Figure 1. Ringel and Youngs [28] finally proved that the maximum chromatic number of any graph G on an orientable surface of genus g ≥ 1 is given by H(g). This graphtheoretic difference between graphs embedded on a surface or on a torus may have physical implications that are worth studying: in particular, is there any additional structure (e.g., phase-transition lines) in the interval q ∈ (4, 7) for the chromatic polynomial of a strip graph embedded on a torus?
The main part of our physical (as opposed to mathematical) understanding of the mechanism for generating chromatic zeros stems from Saleur's analysis of the so-called Berker-Kadanoff (BK) phase [29]. This is a massless phase with algebraic decay of correlations. For the square lattice, this phase is bounded by the curves v ± (q) = −2 ± 4 − q (1.6) and all this region is attracted (at fixed q) by the curve It is important to stress four points here: First, the BK phase corresponds to generic values of q ∈ (0, 4). More precisely, for the Beraha numbers B p the BK phase does not exists due to the vanishing of the amplitude of some of the leading eigenvalues and/or cancellations among them [29,30]. Second, the chromatic-polynomial subspace v = −1 intersects the BK phase for q ∈ [0, 3]. Thus, because of the attractive nature of the BK critical curve (1.7), we can study its conformal properties by considering the chromatic-polynomial case in the window q ∈ [0, 3]. Third, the upper boundary of the BK phase v + (q) (1.6) can be identified as the antiferromagnetic critical curve [31]. Fourth, the exact free energy of the chromatic-polynomial line is not known. However, from our previous studies [9,11,12], we can draw a qualitative picture of the chromatic-polynomial phase diagram for real q: the system is disordered in the region q ∈ (−∞, 0) ∪ (3, ∞), while it is critical in the interval q ∈ (0, 3). The nature of the system in this critical interval depends on the boundary conditions: for free and cylindrical boundary conditions there is a single phase [11,12]; but for cyclic boundary conditions, we find two different phases: the intervals q ∈ (0, 2), and q ∈ (2, 3) corresponding to two distinct values of a topological order parameter (See Section 2.2) [9]. If we define the value q c as the largest value of q such that the system is disordered for all q > q c and the system is critical at q = q c [32], then for the square lattice we have that q c (sq) = 3 [33]. The situation for the triangular lattice is slightly different. This model is solvable on the curve [34] v 3 + 3v 2 = q (1.9) The lower branch of the curve (1.9) can be identified with the lower bound of the BK phase v − (q), while the middle branch of (1.9) is expected to play the same role as the curve (1.7) [i.e., it is the renormalization-group (RG) attractor governing the BK phase]. Thus, for the triangular lattice at least one curve is missing (if we assume a RG flow similar to that of the square lattice): the upper bound of the BK phase [i.e., the analogue of the curve v + (q) (1.7)], and the antiferromagnetic critical curve. In the limit q → 0, we have found [10] that there is a single phase-transition line separating the (critical) BK phase and the (non-critical) high-temperature phase. This line is interpreted as the antiferromagnetic critical curve for this model, and its slope is given approximately by dv dq q=0 = −0.1753 ± 0.0002 (1.10) The position of these two curves is an open problem (Some preliminary numerical results about the antiferromagnetic critical curve were published in [10, Figure 2]. The full curve will be published elsewhere [35]). Theoretical [30] and numerical [10] studies show that the BK phase does exist in the triangular-lattice model, and the results suggest that universality holds for this phase (i.e., the critical properties of the square-and triangular-lattice BK phase coincide). Even though some important curves of the phase diagram of the triangular-lattice chromatic polynomial are not known, Baxter [14] found its exact free energy, albeit with peculiar boundary conditions that amount to not giving a weight q to clusters with non-trivial homotopy. We comment further on that solution in Sections 5-6 (see also Ref. [13]). The resulting phase diagram consists of three different phases in the complex q-plane: the disordered phase contained the interval q ∈ (−∞, 0) ∪ (4, ∞), and there are two critical phases. One contains the interval q ∈ (0, q 2 (tri)), where q 2 (tri) is given by q 2 (tri) = q B ≃ 3.8196717312 (1.11) and the third one contains the interval q ∈ (q 2 (tri), 4). From this phase diagram, we obtain that q c (tri) = 4. Our previous studies [9,13], gave a similar qualitative picture of that phase diagram. In these papers, we have defined q 0 (L) as the largest value of q at which the limiting curve B L crosses the real q-axis. We expect that the limit q 0 = lim L→∞ q 0 (L) exists, and it is equal to q 2 (tri). 3 In Ref. [13] we reanalyzed Baxter's solution and found that if his three analytical expressions for the free energy are correct and exhaustive, the phase transition should not occur at q 2 (tri) = q B , but rather at Unfortunately, our numerical data were not conclusive enough to judge which is the correct value of the transition point. For free and cylindrical boundary conditions [13], the corresponding phase diagram contains only the above three phases. For cyclic boundary conditions [9], we found a richer phase structure within the interval q ∈ (0, 4). In particular, in the interval q ∈ (0, q 2 (tri)), there should be phase transitions at the Beraha numbers . If q 2 (tri) = q B , then the last such transition occurs at q = B 14 ≈ 3.8019377358; but if q 2 (tri) = q G , then the last transition is at q = B 12 = q G . Finally, the chromatic-polynomial line v = −1 has a non-void intersection with the BK phase for the triangular lattice. In particular, we expect that the intersection should include the whole interval q ∈ (0, q 2 (tri)). Thus, the critical properties extracted from this interval are related directly to those of the BK phase. In the interval q ∈ (q 2 (tri), 4), the universality class should be different from that of the BK phase. For other boundary conditions, we have been unable to extract any meaningful physical information for this interval; so this question remains open. In the ferromagnetic regime, it is well-known that results obtained from toroidal boundary conditions are closer to the true thermodynamic limit that for the other boundary conditions. Although the argument leading to this conclusion may not apply for antiferromagnetic or "unphysical" models, it opens an opportunity to try to answer the above-posed questions.
This paper is organized as follows: in Section 2 we show in detail how to obtain the chromatic polynomial of a strip with toroidal boundary conditions, and some structural properties: dimensionality of the transfer matrix and properties of the amplitudes. Sections 3 and 4 contain the detailed description (when possible) of the chromatic polynomial for squareand triangular-lattice strips, respectively. In these two sections we also describe the limiting curve for each regular lattice. The less-mathematically inclined readers can skip most of these sections. In Section 5 we discuss the finite-width results and extract the physically interesting quantities: phase diagram of the chromatic polynomial for the square and triangular lattices, and critical properties of the different phases. Finally, in Section 6, we draw our conclusions. In Appendix A we prove two useful combinatorial identities needed in Section 2.2.

Chromatic polynomials with toroidal boundary conditions
We want to build the transfer matrix for square-and a triangular-lattice strips of width m and arbitrary length n with toroidal boundary conditions. The goal is to write the transfer matrix T(m) as a product of two matrices where H (resp. V) represents the horizontal (resp. vertical) bonds of the lattice, and to write its chromatic polynomial in the form for some left and right vectors u and v id . The subindex P denotes the periodic boundary condition in the corresponding direction.

The transfer-matrix formalism for toroidal boundary conditions
The main problem in writing a transfer matrix in the Fortuin-Kasteleyn representation is the non-local factor q k(E ′ ) in (1.3) [11]. For free boundary conditions in the longitudinal direction (i.e., free and cylindrical), this can be efficiently handled by keeping track of the connectivity state v P of the top layer. This connectivity is related to a certain partition P of the single-layer vertex set {1, 2, . . . , m} [11,13].
For periodic boundary conditions in the longitudinal direction (i.e., cyclic and toroidal) we must keep track of the connectivity state of both the top and bottom rows: each connectivity state v P is now given by a certain partition P of the vertex set of the top row {1, 2, . . . , m} and the vertex set of the bottom row {1 ′ , 2 ′ , . . . , m ′ }. In Figure 2 we show an example of a square-lattice strip of size 4 × 2. As for cyclic boundary conditions [9], initially the top and bottom rows are identical. Then, we enlarge the strip by adding new layers via the action of the transfer-matrix operator T(m) exclusively on the top-row vertices (while the bottom-row sites remain unchanged). At the end, when we have built a strip of n + 1 rows, we identify the top and bottom rows.
Let us now explain in more detail how to obtain the chromatic polynomial for a squarelattice strip of width m, length n, and toroidal boundary conditions.
We characterize the state of the top and bottom rows by a connectivity state v P , which is associated to a partition P = {P 1 , P 2 , . . . , P k } of the set {1, 2, . . . , m, 1 ′ , 2 ′ , . . . , m ′ }. We then introduce the operators acting on the connectivity space {v P }. Here I is the identity, D x is the detach operator that detaches site x from the block it currently belongs to, and J x,y is the join operator that amalgamates the blocks containing sites x and y, if they were not already in the same block.
(Further details can be found in Ref. [11].) The matrices V and H are defined in terms of the above operators as follows: Notice that periodic boundary conditions have been implemented by adding the operator Q m,1 in (2.4a). We can similarly define a matrix H ′ sq acting only on the bottom row.
Then, the chromatic polynomial for a square-lattice strip of size m × n with toroidal boundary conditions can be written as in (2.1)/(2.2), where the two vectors are those of cyclic boundary conditions [9]. The vector v id denotes the partition (i.e., we start with the top and bottom rows identified). The left vector u T acts on a connectivity state v P as where |P ′ | denotes the number of blocks in the connectivity state v ′ P obtained from v P by Thus, u T acts on v P by identifying the top and bottom rows, and then by assigning a factor of q to each block in the resulting partition. The definition of the transfer matrix for a triangular-lattice strip of width m, length n, and toroidal boundary conditions is more involved because of the periodic boundary conditions in the transverse direction. As explained in Refs. [11,13], in this case one has to consider a triangular lattice of width m + 1 with free boundary conditions in the transverse direction, and then identify the columns 1 and m + 1. (See Figure 3 for an example of size 4 × 2.) Then, the matrices H and V take the form where the operators Q i,i+1 in (2.8b) represent the oblique bonds e = (i, j + 1), (i + 1, j) . With this definitions, the chromatic polynomial can be written as in (2.1)/(2.2) with the same right and left vectors as for the square-lattice case (2.5)/(2.6). A further simplification can be obtained if we take into account the fact that v = −1. Then, H is a projector (i.e., H 2 = H), and we can use instead of the transfer matrix T and the basis vectors v P , the modified transfer matrix T ′ and the modified basis vectors w P = H(m) · v P (2.10) Note that w P = 0 if P has any pair of nearest-neighbor sites of the top row in the same block. Thus, the space {w P } has lower dimension than the original connectivity space {v P }.
Notice that we start with the state w id = H · v id , and that at the end we identify the top and bottom rows. Then, it is very useful to work directly with the modified connectivity basis This choice implies that w P = 0 also if P contains any pair of nearest-neighbor sites of the bottom row in the same block. For simplicity, we will drop hereafter the prime in the modified transfer matrix (2.9) and the hat in the modified basis vectors (2.11). Then, the chromatic polynomial can be written as P m P ×n P (q) = u T · T(m) n · w id (2.12) Finally, to lighten the notation, we will write the connectivity states v P using Kronecker delta-functions. For instance, v {{1,2},{1 ′ ,3 ′ },{3},{2 ′ }} will be written as δ 1,2 δ 1 ′ ,3 ′ . The identity connectivity state (2.5) will be denoted simply as 1.

Structural properties of the chromatic polynomial
Our first goal is to compute the number C m,m of non-crossing connectivities in the interior of an annulus with m points on the inner rim and m points on the outer rim (See Figure 4). This problem can be reduced to a simpler one by introducing the quantities d(m, ℓ). Given m points on a circle and a non-crossing connectivity P exterior to that circle, we can compute the number of ways d(m, ℓ, P) we can connect that particular connectivity P to infinity by means of ℓ mutually non-connected paths (bridges). Then, the d(m, ℓ) are the sum over all The values of d(m, ℓ) are given by the following formula [19] where C m is the Catalan number which gives the number of non-crossing connectivities with m points on a circle (or on a line). The number of connectivities C m,m can be obtained easily from the d(m, ℓ) (2.14) by using the ℓ paths as bridges and by matching two such bridge connectivities with the same value of 0 ≤ ℓ ≤ m If we insert (2.14) for the d(m, ℓ), then we obtain a closed expression for C m,m : where we have used the identity (A.1a) proven in Appendix A. The asymptotic growth of C m,m as m → ∞ is exponentially fast with a large constant: In practice, we build the C m,m non-crossing connectivities in a recursive way. We start with all the C m−1,m−1 non-crossing connectivities on an annulus with m−1 points on its inner and outer rims, and then we add a new site on the inner rim. The number of non-crossing connectivities for this annulus is obtained again by matching the ℓ-bridge connectivities Then, we add a new site on the outer rim, and we obtain the C m,m connectivities (2.17). A closed formula for C m,m−1 (2.19) can be obtained if we use the combinatorial identity (A.1b) proven in Appendix A. After some algebra, we obtain valid for all m ≥ 2. This formula also grows exponentially fast in m in the large-m limit: In Table 1 we show the values of C m,m and C m,m−1 for 1 ≤ m ≤ 7. The condition v = −1 implies that connectivities P containing nearest-neighbor sites in any sub-block do not contribute. The number of non-crossing non-nearest-neighbor connectivities for an annulus with m points on its inner and outer rims will be denoted by D m,m . This number gives therefore the dimension of the modified transfer matrix (2.9) for a squareor triangular-lattice strip of width m and toroidal boundary conditions. These numbers are displayed for 1 ≤ m ≤ 7 in Table 1.
As in previous works [9,[11][12][13], we can reduce this dimension by considering the equivalence classes modulo symmetries of non-crossing non-nearest-neighbor connectivity states. For the triangular (resp. square) lattice, we should consider equivalence classes modulo rotations (resp. rotations and reflections with respect to any diameter of the strip). The number of generated equivalence classes modulo rotations (resp. reflections and rotations) of noncrossing non-nearest-neighbor connectivity states for a strip of width m will be denoted by TriTorus ′ (m) (resp. SqTorus ′ (m)). The adjective "generated" means that we only take into account those classes of connectivity states that can be produced by applying the various operators P i and Q i,j on the right vector w id . This distinction is relevant for the square lattice, as there are perfectly legal non-crossing non-nearest-neighbor partitions that cannot be generated when we start from the right vector w id . (See e.g., the discussion on this issue in Ref. [20, end of Section 2.1].) The simplest example corresponds to m = 2 and the connectivity state δ 2,1 ′ δ 1,2 ′ . In general, out of the m possible connectivity states with ℓ = m bridges, only the state v id (2.5) can be realized on a square-lattice strip. The other m − 1 non-crossing connectivity states cannot be realized because of the lattice structure. However, on the triangular lattice, which has a larger set of "vertical" edges, one can produce those m − 1 connectivity states without violating the non-crossing condition.
Finally, as for cyclic boundary conditions, for a given width m we find many repeated eigenvalues. Thus, we define TriTorus(m) (resp. SqTorus(m)) as the number of distinct eigenvalues (with non-zero amplitude) for a triangular-lattice (resp. square-lattice) strip of width m with fully periodic boundary conditions. In Table 1 we have displayed the quantities TriTorus ′ , SqTorus ′ , TriTorus and SqTorus for m ≤ 7 (when available).
The symbolic computation of the transfer matrix and the u T and w id vector has been carried out using a mathematica script for the smallest widths 2 ≤ m ≤ 4, and a C program for the larger widths m ≤ 7. Indeed, for m ≤ 4, both methods perfectly agree. We refer to Refs. [9,12,13] for details about the symbolic computation. For larger widths m ≤ 15, we have computed numerically the transfer matrix using another C program.
2) Even though the final dimensions of the transfer matrices are TriTorus ′ (m) or SqTorus ′ (m), the actual computations are carried out in the larger space of non-crossing connectivities (of dimension C m,m ). This fact limits the maximum width we are able to handle.
3) We have cross-checked our results using the trivial identity P m P ×n P (q) = P n P ×m P (q) (2.22) for all 2 ≤ m, n ≤ 7.

Amplitudes for toroidal boundary conditions
The transfer matrix for a strip with toroidal boundary conditions has a block structure. This property also holds for cyclic boundary conditions [9], and in both cases the origin is the same. Given a certain connectivity P of the top and bottom rows (see e.g., Figure 4), we can regard it as a bottom-row connectivity P bot , a top-row connectivity P top joined by ℓ bridges. Notice that each top-row block can be connected to a bottom-row block by at most one bridge, and vice versa. Thus, the number of bridges is an integer 0 ≤ ℓ ≤ m.
As explained in Ref. [9], the full transfer matrix has a lower-triangular block form (if we order the connectivity states in decreasing order of bridges).
The reason is the application of the H and V operators cannot increase the number of bridges of a given connectivity state. Furthermore, as the transfer matrix does not act on the bottom-row connectivity and they act on the space of non-crossing connectivities, each diagonal block T ℓ,ℓ has also a diagonal-block form where each sub-block T ℓ,ℓ,j is characterized by a certain bottom-row connectivity P bot and a position of the ℓ bridges. Its dimension is given by the number of top-row connectivities that are compatible with P bot , and the number ℓ and relative positions of the bridges. The number of blocks N ℓ (0 ≤ ℓ ≤ m) for square-and triangular-lattice strips of width m are displayed in Tables 2 and 3, respectively. In particular, as for cyclic boundary conditions [9], this structure means that the characteristic polynomial of the full transfer matrix T can be factorized as follows In practice, for each strip width 2 ≤ m ≤ 7, we have performed the following procedure: 1. Using a C program, we compute the full transfer matrix T(m) and the left u and right vectors w id . These objects are obtained in the basis of non-crossing non-nearestneighbor classes of connectivities that are invariant under the right symmetries (depending on the lattice structure). The dimension of this transfer matrix is SqTorus ′ (m) or TriTorus ′ (m). With these objects, we can compute the chromatic polynomial of any strip of finite length n using Eq. (2.12).
2. Given the full transfer matrix T(m), we compute the the SqTorus(m) or TriTorus(m) distinct eigenvalues λ i (q). To do this, we split the transfer matrix into blocks T ℓ,ℓ,j , each block characterized by a bottom-row connectivity and the number ℓ and position of the bridges (modulo the corresponding symmetries). By diagonalizing these blocks we obtain the whole set of distinct eigenvalues λ i (q) and their multiplicities k i .
3. The final goal is to express the chromatic polynomial in the form We have used this procedure to compute the transfer matrix up to width m = 7 for both the square and the triangular lattices. The next cases were unmanageable with our current computer resources. 6 The number of distinct eigenvalues found in each block T ℓ,ℓ can be found in Tables 4 and 5 for the square and triangular lattices, respectively.
In Ref. [9], we found several properties of the eigenvalues λ i and their corresponding amplitudes α i for cyclic boundary conditions. In particular, that all eigenvalues within a block T ℓ,ℓ have the same amplitude α (ℓ) (see below), and that any two distinct blocks T ℓ,ℓ and T ℓ ′ ,ℓ ′ (with ℓ = ℓ ′ ) have no common eigenvalues. These properties were essentially due to the presence of a quantum-group symmetry in the Potts model with cyclic boundary conditions. These properties were also used to simplify the algorithm to compute symbolically the transfer-matrix eigenvalues λ i . Unfortunately, this quantum-group symmetry is broken when toroidal boundary conditions are considered, and none of these properties hold in our case.
Chang and Shrock [19] have found that the chromatic polynomial for toroidal boundary conditions (2.26) has a structure similar (but not identical) to that for cyclic boundary conditions [9,29,36,37]. For a square-or triangular-lattice strip of width m, length n and toroidal boundary conditions, the chromatic polynomial can be written as j are polynomials in q of order ℓ, and in contrast with cyclic boundary conditions, are not the same within a given ℓ sector. The eigenvalues λ ℓ,j (m) are algebraic functions of q and, in contrast with cyclic boundary conditions, a given eigenvalue can appear in more than one ℓ sector.
Chang and Shrock [19] found a family of basic amplitudes b (ℓ) given by and for ℓ ≥ 2, by conditions, a given eigenvalue for toroidal boundary conditions can appear in several sectors (each of them characterized by a different value of ℓ. See text.) Once the Ansatz for the amplitudes is fixed, we solve the linear equations for the numerical coefficients in α i (q) using as many chromatic polynomials P mP×nP (q) as needed. After finding the solution, we check that it can reproduce chromatic polynomials not used in the determination of the amplitudes.
where the coefficients α (ℓ) are the amplitudes obtained for cyclic boundary conditions [29] Finally, Chang and Shrock argue that the coefficients b (ℓ) j cannot be equal to the above b (ℓ) . For m = 2, 3 they find that the coefficients b (2) and b (3) split into two coefficients (up to a positive integer factor): In our numerical study, we have found that this also true for m = 4, 5. The basic amplitudes are given by Note that the expression for b In general, an eigenvalue λ i which only appear in a single T ℓ,ℓ block has an amplitude given by a positive integer multiple of either (2.28)/(2.31) or (2.32)/(2.33). All the amplitudes b (ℓ) j and b (ℓ) (2.31) are simply polynomials in q of order ℓ. When an eigenvalue λ i appears in more than one block T ℓ,ℓ , then its amplitude is just a certain linear combination (with integer coefficients) of the corresponding amplitudes given above. Thus, the order of this polynomial amplitude is just the largest value of ℓ involved.
Remarks. 1. After we empirically obtained the amplitudes (2.33), one of us obtained the general formula for the different amplitudes b (ℓ) j [20]. Indeed, our results fully agree with the general formula. The number of distinct amplitudes b (ℓ) j for a given ℓ was shown to be the number of integer divisors of ℓ [20].
2. As we shall discuss in the next Section, the sub-block T L,L of the square-lattice strip of width L contains the eigenvalue λ = (−1) L with amplitude b (L) (2.31). Indeed, each b (ℓ) can be written as a linear combination of the amplitudes b This fact was first noted by Chang and Shrock [19]; and has been proven in Ref. [20, c.f., Eq. (4.7)] by showing how many times each b (ℓ) j appears in a given b (ℓ) . 3. In our previous study of the chromatic polynomial for cyclic boundary conditions [9], we could only obtain the exact chromatic polynomial up to widths m = 6. The reason was that the computation of the characteristic polynomial of a symbolic matrix of dimension ∼ > 100 was unfeasible with mathematica's standard built-in functions (based essentially on Gaussian elimination). In this paper, we have been able to overcome this problem by using a different division-free algorithm, which computes the characteristic polynomial of a matrix with elements belonging to a ring (not to a field). We have used a mathematica implementation [38] of the Samuelson-Berkowitz-Abdeljaoued algorithm [39]. 8 With this algorithm we have been able to deal with symbolic matrices of dimensions up to 308. Unfortunately, for L = 7 we need to compute the symbolic characteristic polynomial of much larger matrices (of dimensions 483 and 532). These computations are beyond our current computer facilities, even with the help of the above-described improved algorithm. 9 Further details will be published elsewhere.

Square-lattice numerical results
In this section we will analyze the results for the square-lattice strips of widths 2 ≤ L ≤ 7. We have checked our results using (in addition to the already known cases L ≤ 4 [16,17,19]) the identity (2.22).
For each L, we give a detailed (when possible) description of the eigenvalues λ i and amplitudes α i . We have included the already known cases L ≤ 4 for completeness and, more importantly, to show in detail how our method works. In particular, for L = 2, 3 we give the explicit expression of the transfer matrix and the left and right vectors. The exact expressions for the new cases L = 5, 6 are very lengthy, so we refrain from writing down all the needed formulae. Instead, we have included a mathematica file named transfer6 sqT.m available as part of the on-line version of this paper in the cond-mat archive at arXiv.org. For L = 7, even though we have obtained the symbolic form of the transfer matrix and the left and right vectors, we have been unable to compute some of the eigenvalues (because of the large dimensionality of some sub-blocks), and thus, the expression of the amplitudes.
Finally, for each strip, we compute the corresponding limiting curve and isolated limiting points (where the chromatic zeros accumulate in the infinite-length limit), and analyze their main properties.

L = 2
The chromatic polynomial for this strip is exactly the same as the chromatic polynomial for a square-lattice strip of with L = 2 with cyclic boundary conditions [9]: in the basis where we show by vertical and horizontal lines the lower-triangular block structure of this matrix. The right w id and left u vectors are given by As in [9], we find that the three diagonal blocks for ℓ = 0, 1, 2 bridges give the following eigenvalues and amplitudes: • ℓ = 0: The eigenvalue is λ 0 = q 2 − 3q + 3 with amplitude b (0) = 1.
• ℓ = 2: The eigenvalue is Each diagonal block T ℓ,ℓ is characterized by the bottom-row connectivity P bot = 1.
In Figure 5(a) we have shown the chromatic zeros for the strips of lengths n = 10 and n = 20. We also show the limiting curve B 2 where the chromatic zeros accumulate in the infinite-length limit.
The limiting curve B 2 splits the complex q-plane into four regions. The regions inside the outer parts of B 2 are characterized by a dominant eigenvalue belonging to the ℓ = 1 sector, while in the outer region the eigenvalue λ 0 dominates.
The curve B 2 crosses the real q axis at two points: q = 0 and q = 2. This latter point is in fact a multiple point. Finally, there is a single isolated limiting point at q = 1.

L = 3
The chromatic polynomial for this strip was computed by Chang and Shrock [17] (in fact, they computed the full Potts-model partition function [19]).
• ℓ = 2: The four eigenvalues are The amplitudes are given respectively by b 2 , and b In Figure 5(b) we have shown the chromatic zeros for the strips of lengths n = 15 and n = 30, as well as the limiting curve B 3 where the chromatic zeros accumulate in the infinitelength limit. This curve splits the complex q-plane into three regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The region containing the real segment (2, 3) is dominated by the ℓ = 2 sector. The rest of the complex q-plane is dominated by the ℓ = 0 eigenvalue λ 0 (3.4c).
The curve B 3 crosses the real q axis at three points: q = 0, q = 2, and q = 3. We also find a pair of complex conjugate T points at q ≈ 2.4692972375 ± 1.4013927796 i. Finally, there is a single isolated limiting point at q = 1.

L = 4
The chromatic polynomial for this strip was computed by Chang and Shrock [16]. The full transfer matrix T(4 T ) has dimension 68. In this case, we restrict ourselves to report the eigenvalue structure: • ℓ = 0: This block has dimension five with two sub-blocks of dimensions three and two. There are three distinct eigenvalues: • ℓ = 2: This block has dimension 33 containing three sub-blocks of dimensions 8, 12, and 13. Among these 33 eigenvalues, we find 15 distinct eigenvalues. Three of them are simple: 1 , and the special one λ (2,3) = 2−q which also appears in the ℓ = 3 sector (see below). The next six eigenvalues are solutions of three second-order equations [16, Eqs. (2.14)-(2.16)], with amplitudes b 1 , and b 1 and 2b (2) 2 , respectively. • ℓ = 2, 3: The eigenvalue λ (2,3) = 2 − q appears in these two blocks with an amplitude 2 . • ℓ = 3: This block has dimension seven: in addition to λ (2,3) , we find six other eigenvalues: λ 3,1 = 5 − q with amplitude b Thus, for this strip we find 33 different eigenvalues. All of them but λ (2,3) = 2 − q, belong to a single ℓ sector.
In Figure 5(c) we have shown the chromatic zeros for the strips of lengths n = 20 and n = 40, as well as the limiting curve B 4 . This curve splits the complex q-plane into five regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The region containing the real segment (2, q 0 (4)) with q 0 (4) ≈ 2.7827657401 is dominated by the ℓ = 2 sector. The rest of the complex q-plane (i.e., the other three regions) is dominated by the ℓ = 0 sector. Notice that the special eigenvalue λ (2,3) = q − 2 is not dominant in any open set of the complex q-plane.
The curve B 4 crosses the real q axis at three points: q = 0, q = 2, and q = q 0 (4) ≈ 2.7827657401. We also find three pairs of complex conjugate T points at q ≈ 2.4184666383 ± 1.7281014476 i, q ≈ 2.6228997762 ± 1.5279208067 i, q ≈ 2.7784542034 ± 0.9708262424 i. Finally, there is a single isolated limiting point at q = 1.

L = 5
The full transfer matrix T(5 T ) has dimension 347. The relevant eigenvalues can be extracted from the diagonal blocks T ℓℓ (5 T ). We find the following eigenvalue structure: • ℓ = 0: This block has dimension six. In the basis . .} (where the dots ". . . " mean all equivalent states under rotations and reflections), it takes the form: where s k and S are defined in (3.4a)/(3.4b), and This matrix is block diagonal: the first four-dimensional block corresponds to the bottom-row connectivity P bot = δ 1 ′ ,3 ′ , while the second two-dimensional block corresponds to the trivial bottom-row connectivity P bot = 1.
In Figure 5(d) we have shown the chromatic zeros for the strips of lengths n = 25 and n = 50, and the limiting curve B 5 . This curve splits the complex q-plane into three regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The region containing the real segment (2, 3) is dominated by the ℓ = 2 sector. The rest of the complex q-plane is dominated by the ℓ = 0 sector. Notice that the special eigenvalue λ (1,2) is not dominant in any open set of the complex q-plane.
The curve B 5 crosses the real q axis at three points: q = 0, q = 2, and q = 3. We also find two pairs of complex conjugate T points at q ≈ 2.2834116867 ± 2.0404211832 i, q ≈ 2.4349888915 ± 1.9935503628 i. We also find a pair of complex-conjugate endpoints at q ≈ 2.5034648020 ± 2.0851731763 i. Finally, there is a single isolated limiting point at q = 1.

L = 6
The full transfer matrix T(6 T ) has dimension 2789. The relevant eigenvalues can be extracted from the diagonal blocks T ℓℓ (6 T ). We find the following eigenvalue structure: • ℓ = 0: This block has dimension 39, and it contains five sub-blocks of dimensions five, seven, eight (two of them), and eleven. We find eleven distinct eigenvalues. The simplest one is λ 0,1 = q 4 − 9q 3 + 30q 2 − 44q + 25 with amplitude 2b (0) = 2. Two eigenvalues are the solutions of the second-order equation The other eigenvalues come from polynomial equations of order three (with amplitude 2b (0) = 2) and five (with amplitude b (0) = 1).
• ℓ = 1: This block has dimension 494 and it contains eleven sub-blocks of dimensions ranging from 21 to 71. We find 48 different eigenvalues in this block. The simplest ones are given by λ (1,2) which also they appear in other blocks (see below). Two eigenvalues are the solutions of the second-order equation: The other eigenvalues come from polynomial equations of order ten, eleven (two of them), and twelve.
• ℓ = 4: This block has dimension 163 and it contains four sub-blocks of dimensions 26, 45, and 46 (two of them). We have found 49 distinct eigenvalues. The simplest one is 1 (2.33a). Six eigenvalues are solutions of three second-order equations The other eigenvalues come from equations of order three (four equations), four (four equations), six, and eight.
• ℓ = 5: This block has dimension 16. The simplest eigenvalues are The remaining eight eigenvalues come from two fourth-order equations. The amplitude of all these twelve eigenvalues is 2b • ℓ = 6: This one-dimensional block gives λ 6 = 1 with amplitude b (6) (2.31e).
In Figure 6(a) we have shown the chromatic zeros for the strips of lengths n = 30 and n = 60, and the limiting curve B 6 . This curve splits the complex q-plane into four regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The region containing the real segment (2, 5/2) is dominated by the ℓ = 2 sector. The two small complex conjugate regions are dominated by the ℓ = 1 sector. The rest of the complex q-plane is dominated by the ℓ = 0 sector.
The curve B 5 crosses the real q axis at three points: q = 0, q = 2, and q ≈ 2.9078023424. We also find four pairs of complex conjugate T points at q ≈ 2.1180081660 ± 2.2566257839 i, q ≈ 2.1632612447±2.2209172289 i, q ≈ 2.8895175995±1.1648906031 i, and q ≈ 2.9197649745± 0.7583704656 i. We also find a pair of complex-conjugate endpoints at q ≈ 2.0571168133 ± 2.3885607275 i. Finally, there is a single isolated limiting point at q = 1.

L = 7
The full transfer matrix T(7 T ) has dimension 20766. We have been unable to obtain the full eigenvalue structure, as some of the blocks for ℓ = 2, 3 are very large. Thus, we cannot obtain the amplitudes. We find the following (partial) eigenvalue structure: • ℓ = 0: This block has dimension 111, and it contains six sub-blocks of dimensions six and 21 (five of them). We find 21 different eigenvalues coming from polynomial equations of order six and 15. • ℓ = 4: This block has dimension 2711, and it contains eight sub-blocks of dimensions 155, 158, and 308. We find 165 different eigenvalues coming from polynomial equations of order three, four, seven, eight, eleven, 33 (two of them), and 66.
• ℓ = 5: This block has dimension 288, and it contains four sub-blocks of dimension 72. We find 72 different eigenvalues coming from polynomial equations of order four, eight, twelve, and 48.
In Figure 6(b) we have shown the chromatic zeros for the strips of lengths n = 35 and n = 70, and the limiting curve B 7 .
The curve B 7 crosses the real q axis at three points: q = 0, q = 2, and q = 3. We also find six pairs of complex conjugate T points at q ≈ 1.4876491726 ± 2.5577806355 i, q ≈ 1.6498030004 ± 2.5002880707 i, q ≈ 1.9765852319 ± 2.3664338035 i, q ≈ 2.5531414480 ± 1.8775702667 i, q ≈ 2.8395564367 ± 1.4513813377 i, and q ≈ 2.8734078918 ± 1.3454749268 i. There are two complex conjugate oval-shaped regions between the T points at q ≈ 1.488 ± 2.558 i and q ≈ 1.650 ± 2.500 i. And there is another pair of small bulb-like regions between the T points at q ≈ 2.840 ± 1.451 i and q ≈ 2.873 ± 1.345 i.
Thus, the limiting curve B 7 splits the complex q-plane into seven regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The region containing the real segment (2, 3) is dominated by the ℓ = 2 sector. The two small complex conjugate oval-shaped regions are dominated by the ℓ = 0 sector; and the two small complex conjugate bulb-like regions are dominated by the ℓ = 1 sector. The rest of the complex q-plane is dominated by the ℓ = 0 sector.
We also find a pair of complex-conjugate endpoints at q ≈ 1.6947007027±2.5327609879 i. Finally, there is a single isolated limiting point at q = 1.

Triangular-lattice numerical results
In this section we will analyze the results for the triangular-lattice strips of widths 2 ≤ L ≤ 7. We have checked our results using (in addition to the already known cases L ≤ 4 [16,18,19]) the identity (2.22).
As for the square lattice, for each L we give a detailed description of the eigenvalues and amplitudes. For the already known cases L = 2, 3, we include for the sake of clarity the full expressions for the transfer matrices and the left and right vectors. The exact expressions for the new cases L = 5, 6 can be found in the mathematica file transfer6 triT.m available as part of the on-line version of this paper in the cond-mat archive at arXiv.org. As in the square-lattice case, for L = 7 we can only present a partial description of the eigenvalue structure, as some of the sub-blocks are very large.

L = 2
The chromatic polynomial for this strip was computed by Chang and Shrock [19] (in fact, they computed the full Potts-model partition function).
• ℓ = 1, 2: The eigenvalue λ ⊙ = 0 appears in these two blocks. Its amplitude is irrelevant for the computation of both the limiting curve B 2 and the chromatic polynomials P 2 P ×n P (q).
In Figure 7(a) we have shown the chromatic zeros for the strips of lengths n = 10 and n = 20, and the limiting curve B 2 . This curve splits the complex q-plane into three regions. The outer one is dominated by the ℓ = 0 sector. The one containing the real interval q ∈ (0, 2) is dominated by the ℓ = 1 sector; and the one containing the real interval q ∈ (2, 4) is dominated by the ℓ = 2 sector.
The curve B 2 crosses the real q axis at three points: q = 0, q = 2, and q = 4. This latter point is in fact a multiple point. Finally, there are two isolated limiting points at q = 1 and q = 3.
• ℓ = 3: In addition to λ (2,3±) , we find λ 3 = −2 with amplitude b In Figure 7(b) we have shown the chromatic zeros for the strips of lengths n = 15 and n = 30, and the limiting curve B 3 . This curve splits the complex q-plane into three regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The region containing the real segment (2,3] is dominated by the ℓ = 2 sector. The rest of the complex q-plane is dominated by the ℓ = 0 eigenvalue λ 0 (4.4b).
The curve B 3 crosses the real q axis at three points: q = 0, q = 2, and q ≈ 3.7240755514. We also find a pair of complex conjugate T points at q ≈ 3.9777108575 ± 1.6284204463 i. Finally, there are two isolated limiting points at q = 1, 3.

L = 4
The chromatic polynomial for this strip was computed by Chang and Shrock [16]. The full transfer matrix T(4 T ) has dimension 99. In this case, we restrict ourselves to report the eigenvalue structure: • ℓ ≥ 0: The eigenvalue λ ⊙ = 0 is common to all blocks 0 ≤ ℓ ≤ 4. Its amplitude cannot be computed from our analysis; but it is unimportant to compute the chromatic polynomials and the limiting curve.
• ℓ = 0: This block has dimension five with two sub-blocks of dimensions three and two.
• ℓ = 1: This block has dimension 26, and it contains three sub-blocks of dimensions six, and ten (two of them). We find eight different eigenvalues, one of them being λ ⊙ = 0. The other eigenvalues come from a third-order equation [16,Eq. (4.12)] and from a fourth-order equation [16,Eq. (4.14)]. The amplitudes are in all cases b (1) = q − 1.
Thus, for this strip we find 38 distinct eigenvalues. All of them but three (λ ⊙ and λ (3,4±) ) belong to a single ℓ sector.
In Figure 7(c) we have shown the chromatic zeros for the strips of lengths n = 20 and n = 40, and the limiting curve B 4 . This curve splits the complex q-plane into three regions.
The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The region containing the real segment (2, 4) is dominated by the ℓ = 2 sector. The rest of the complex q-plane is dominated by the ℓ = 0 sector.
The curve B 4 crosses the real q axis at three points: q = 0, q = 2, and q = 4. We also find a pair of complex conjugate T points at q ≈ 3.3943280448 ± 2.1041744504 i. Finally, there are two isolated limiting points at q = 1, 3.

L = 5
The full transfer matrix T(5 T ) has dimension 653. The relevant eigenvalues can be extracted from the diagonal blocks T ℓℓ (5 T ). We find the following eigenvalue structure: • ℓ = 0: This block has dimension eight, and it contains two sub-blocks of dimensions two and six. We find six distinct eigenvalues (with amplitude b (0) = 1) coming from a fourth-order equation and the following second-order equation: • ℓ = 1: This block has dimension 125 and it contains seven sub-blocks of dimension 25. We find 24 different eigenvalues in this block. The eigenvalue λ (1,2) = −(q 2 −5q+6) also appears in the ℓ = 2 block (see below). The other eigenvalues come from polynomial equations of order three and 20. They all have the same amplitude b (1) = q − 1.
In Figure 7(d) we have shown the chromatic zeros for the strips of lengths n = 25 and n = 50, and the limiting curve B 5 . This curve splits the complex q-plane into six regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The regions containing the real segment (2, 4) is dominated by the ℓ = 2 sector. The rest of the complex q-plane (including the two complex-conjugate triangular-shaped regions) is dominated by the ℓ = 0 sector.
The curve B 5 crosses the real q axis at four points: q = 0, q = 2, q ≈ 3.8482632901 and q = 4. We also find four pairs of complex conjugate T points at q ≈ 2.7382976052 ± 2.8142366786 i, q ≈ 2.7643522733 ± 2.3252257137 i, q ≈ 3.7492171129 ± 1.9879721106 i, and q ≈ 3.9634199872±0.3751435740 i. Finally, there are two isolated limiting points at q = 1, 3.

L = 6
The full transfer matrix T(6 T ) has dimension 5194. The relevant eigenvalues can be extracted from the diagonal blocks T ℓℓ (6 T ).
• ℓ ≥ 0: There is a single eigenvalue common to all blocks λ ⊙ = 0. As for the cases L = 2, 4, its amplitude cannot be determined; but it is irrelevant to compute both the chromatic polynomials and the limiting curve.
• ℓ = 2: This block has dimension 2712 and it contains 21 sub-blocks of dimensions ranging from 96 to 180. We find 165 different eigenvalues in this block: λ ⊙ = 0, and the solutions of polynomial equations of order three (two of them), four (two of them), ten (three of them), 28 (two of them), and 32 (two of them).
Thus, in this strip we find 476 distinct eigenvalues. All of them but five (λ ⊙ = 0 and λ (5,6,i) ) belong to a single sector.
In Figure 8(a) we have shown the chromatic zeros for the strips of lengths n = 30 and n = 60, and the limiting curve B 6 . This curve splits the complex q-plane into seven regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector. The region containing the real segment (2, 3) is dominated by the ℓ = 2 sector. The two bulb-like complex-conjugate regions appearing on the right hand-side of the limiting curve belong to the ℓ = 1 sector and they are rather special (see below). The rest of the complex q-plane (including the two complex-conjugate triangular-shaped regions) is dominated by the ℓ = 0 sector.
The curve B 6 crosses the real q axis at three points: q = 0, q = 2, and q ≈ 3.7943886378. We also find five pairs of complex conjugate T points at q ≈ 2.0352758561 ± 3.1492288108 i, q ≈ 3.3399899103 ± 2.4428522934 i, q ≈ 2.5037195029 ± 2.3940503612 i, q ≈ 4.1131475777 ± 0.7636584371 i, and q ≈ 4.1295836790 ± 0.8010545279 i. Finally, there are two isolated limiting points at q = 1, 3.
Important technical remark. As stated above, the two bulb-like complex-conjugate protruding from T points q ≈ 4.11 ± 0.76 i and q ≈ 4.13 ± 0.80 i, are rather special. Their peculiarity arises from the fact that inside these two bulb-like regions there is a pair of exactly equimodular dominant eigenvalues [with non-zero amplitudes ∼ (q − 1)]. Indeed, we do not see chromatic zeros becoming dense in the interior of these regions, but only on their boundaries.
The explanation of this phenomenon is rather simple: the two dominant eigenvalues come from a polynomial equation of order 24 belonging to the ℓ = 1 sector. A closer look at this equation reveals that its roots z k can be written as z k = z 0,k e ±πi/3 where the z 0,k are the roots of another polynomial equation of order 12 with integer coefficients (that can be easily computed). It is important to stress that the existence of this pair of equimodular eigenvalues does not mean that the chromatic zeros accumulate in this open set in the infinite-length limit. According to the Beraha-Kahane-Weiss theorem [25], the eigenvalues should satisfy a "no-degenerate-dominance" condition: there do not exist indexes k = k ′ such that λ k = ωλ k ′ with |ω| = 1, and such that the region where λ k and λ k ′ are dominant has non-empty interior. This is precisely the condition violated by the leading eigenvalues in the interior of the above bulb-like regions. Thus, in order to compute the limiting curve in these regions, we should only consider the roots z 0,k of the 12-order equation. In this way, we ensure the "no-degenerate-dominance" condition. In conclusion, the chromatic zeros cannot accumulate in the interior of these bulk-like regions; but on their boundary (where the dominant z 0,k root from the ℓ = 1 sector becomes equimodular to the leading eigenvalue from the ℓ = 0 sector).
In other strips we have already found the existence of two different eigenvalues differing by a phase [i.e., λ and λe iφ ]. However, this is the first case where one of these "degenerate" pairs is dominant on an open set of the complex q-plane. In order to correctly compute the limiting curve, one should eliminate one of eigenvalues of each "degenerate" pair to ensure the hypothesis of the Beraha-Kahane-Weiss theorem.

L = 7
The full transfer matrix T(7 T ) has dimension 40996. 10 We have been unable to obtain the full eigenvalue structure, as some of the blocks for ℓ = 2, 3 are very large. We find the following (partial) eigenvalue structure: • ℓ = 0: This block has dimension 186, and contains five sub-blocks of dimension 36, and a single block of dimension six. We find 36 distinct eigenvalues coming from two polynomial equations of order six and 30, respectively.
• ℓ = 1: This block has dimension 5488, and contains 28 sub-blocks of dimension 196. We find 196 distinct eigenvalues coming from a single polynomial equation of order 196.
• ℓ = 2: This block has dimension 20216, and contains 38 sub-blocks of dimension 532. A numerical study of these sub-blocks reveals that there are 508 distinct eigenvalues in this sector.
• ℓ = 3: The dimension of this block is 11109 and it contains 23 sub-blocks of dimensions 483. We numerically find that there are 454 distinct eigenvalues for this sector.
• ℓ = 5: The dimension of this block is 560 and it contains four sub-blocks of dimensions 140. There are 122 distinct eigenvalues coming from polynomial equations of order four, eight, 24, and 96.
In Figure 8(b) we have shown the chromatic zeros for the strips of lengths n = 35 and n = 70, and the limiting curve B 7 . This curve splits the complex q-plane into eight regions. The region containing the real segment (0, 2) is dominated by an eigenvalue from the ℓ = 1 sector; the region containing the real segment (2, 3) is dominated by the ℓ = 2 sector; and the small region containing the real segment (3.8, 4) is dominated also by the ℓ = 2 sector. There are two small lens-like complex-conjugate regions which belong to the ℓ = 1 sector. The rest of the complex q-plane (including the two complex-conjugate triangular-shaped regions) is dominated by the ℓ = 0 sector.

Square lattice limiting curves
A summary of qualitative results for the limiting curves B, as obtained in Section 3, is shown in Table 6. In particular, we have computed q 0 (L), the maximum real value in B for strip widths L ≤ 7 from the symbolic transfer matrices. In Figure 9, we display together all the computed limiting curves B L with 2 ≤ L ≤ 7. For comparison, we also include the limiting curve for a strip of width L = 10 with cylindrical boundary conditions [12]. These computations can be taken a lot further by diagonalizing the diagonal blocks T ℓℓ (L T ) of the numerical transfer matrix. Indeed it suffices to diagonalize one of its subblocks, e.g., the one with a trivial bottom-row connectivity. This amounts to working with a transfer matrix that acts on connectivities of only L (rather than 2L) points, with exactly ℓ marked clusters. For the purpose of getting only the eigenvalues, without going outside the specified sub-block, transfer matrix elements which would amount to adjoining two distinct marked clusters are set to zero.
Since further we need only the largest (or sometimes the first few largest) eigenvalues of the given sub-block, we can diagonalize using a standard iterative procedure (the socalled power method [40]) which has the immense advantage of allowing for sparse matrix factorization (i.e., the lattice edges are added one by one) combined with standard hashing techniques. This allowed us to access widths L ≤ 12. In all cases, q 0 (L) was found to correspond to the crossing (in modulus) of the leading eigenvalues of the ℓ = 0 and ℓ = 2 sectors. The precise value of q 0 (L) was then narrowed in by a Newton-Ralphson method. The final results are shown in Table 6. 11 From the curious fact that q 0 (L) = 3 exactly for L = 3, 5, 7, 9, 11 we conjecture that this is true for all square-lattice strips with odd widths: Conjecture 5.1 For square-lattice strips with fully periodic boundary conditions of widths L = 2p + 1 with p ≥ 1, we have q 0 (L) = 3.
In other words, we have found a particular sequence of strip graphs for which the value of q 0 (L) is constant and, furthermore, equal to the expected critical value q c (sq) = 3. Indeed, recall that for q = 3, the square-lattice chromatic polynomial is equal to three times the partition function of the six-vertex model with all vertex weights equal to one [33]; it is well-known that this latter model is critical (with central charge c = 1).
We next conjecture that parity effects do not affect the value of q 0 in the thermodynamic limit: Conjecture 5.2 For square-lattice strips with fully periodic boundary conditions of widths L = 2p, the limit q 0 (sq) = lim p→∞ q 0 (L) exists and is equal to q c (sq) = 3.
This conjecture is very similar to the one made in [12]; the only difference being the boundary conditions, cylindrical rather than toroidal.
As to the value of q 0 (L) for even L, the data of Table 6 is well fitted by the power-law form q 0 (L) = q c (sq) + AL −1/ν . Our estimates are q 0 (sq) = 2.999 ± 0.012 and ν = 0.49 ± 0.02. These results strongly corroborate Conjecture 5.2. A better estimate for ν can be obtained by fixing q 0 (sq) = 3 in the above Ansatz: we obtain the critical exponent ν = 0.502 ± 0.002. We conjecture that the exact value of the exponent is ν = 1/2, consistent with the transition being first-order in q.
The limiting curves B L for 2 ≤ L ≤ 7 show some regularities (see Figure 9): for real q, we have three different phases and each of them corresponds to a different ℓ sector (ℓ = 0, 1, 2). Furthermore, within each phase, the leading amplitude does not depend on L. All our empirical findings can be summarized in the following conjecture: The phase diagram of the zero-temperature square-lattice Potts antiferromagnet with toroidal boundary conditions at real q has three different phases characterized by the number of bridges ℓ as follows: Note that for all widths studied, 2 ≤ L ≤ 12, we have that 2 ≤ q 0 (L) ≤ 3.
We also find a curious behavior of one of the eigenvalues belonging to the highest sector for each strip width: The diagonal sub-block T L,L for a square-lattice strip with fully toroidal boundary conditions and width L contains the eigenvalue λ = (−1) L with amplitude b (L) .

Triangular lattice limiting curves
A summary of qualitative results for the limiting curves B is shown in Table 7. In Figure 10, we display together all the computed limiting curves B L with 2 ≤ L ≤ 7. For comparison, we also include the infinite-volume limiting curve obtained by Baxter [14]. The results in Table 7 with L ≤ 7 were obtained from the symbolic transfer matrices in Section 4, and those with L ≤ 12 were found by the numerical procedure outlined in Section 5.1. In all cases q 0 (L) is given by the crossing between the dominant eigenvalues in the ℓ = 0 and ℓ = 2 sectors.
A quick glance at Table 7 however reveals that the values of q 0 (L = 3p) are definitely not converging to 4. We therefore examine the data with L = 3p more closely. A first fit with q 0 (L) = q 0 (tri) + AL −1/ν yields q 0 (tri) ≃ 3.820 and ν ≃ 0.52 (See Table 8). This is consistent with the expected value ν = 1/2 corresponding to a first-order phase transition (in the variable q), as we have already found for the square lattice. Fixing ν = 1/2 we then arrive at the final estimate q 0 (tri) = 3.8196 ± 0.0005.

Remark.
We have performed all the fits reported in this paper using the following method: as the data is essentially exact, we have done a standard least-squares fit with as many data points as the number of unknown parameters N A in the Ansatz, so that in each fit there are no effective degrees of freedom. To estimate the error bars (due to the existence of correction-to-scaling terms not included in the Ansatz), we introduce the parameter L min such that each fit is performed with N A consecutive data points with L ≥ L min . From the variation of the estimates as a function of L min , we deduce the corresponding error bars.
Finally, from the FSS theory of first-order phase transitions [24], we expect that higher corrections-to-scaling terms should be integer powers of L −1/ν = L −2 . Thus, we have fitted our data to the improved Ansatz q 0 (L) = q 0 (tri) + AL −2 + BL −4 (5.1) The results are displayed in Table 8, and from these we conclude that q 0 (tri) = 3.8198 ± 0.0006. We now recall that the triangular-lattice chromatic polynomial was actually shown by Baxter [14] to be an integrable model. In particular, using coordinate Bethe Ansatz, this author identified three different analytic expressions for the free energy in the thermodynamic limit.
To be precise, Baxter used first the Nienhuis mapping [41] between the zero-temperature q-state Potts antiferromagnet on the triangular lattice and the critical O(n) loop model on the hexagonal lattice. That mapping includes an exact relation between q and the weight n of a (bulk) loop. Baxter then transformed the loop model into a vertex model, by the usual trick of redistributing the loop weight n as local complex vertex weights. Unfortunately, when imposing periodic transverse boundary conditions on the vertex model, the loops of non-trivial homotopy (i.e., that wrap around the periodic direction) get the incorrect weight of 2, a fact not mentioned in Baxter's paper. Actually, while giving a (boundary) weight of n ′ = n to those loops (by twisting the vertex model) would have been a nicer way to perform the analysis of the O(n) model, the subtleties of the Nienhuis mapping [41] are such that yet a different weight n ′ (q) = n would have been needed to ensure that clusters of non-trivial homotopy in the Potts antiferromagnet carry their correct weight of q.
One might think that these subtleties have no implications on the behavior in the thermodynamic limit, but the reader is reminded that we are here dealing with the antiferromagnetic regime in which boundary conditions can (and sometimes do) matter. Nevertheless, it seems reasonable to expect that this should have no consequence for the bulk free energies. Numerical analysis [13,14] of Baxter's analytically determined free energies shows that the corresponding eigenvalues become equimodular at q B ≃ 3.8196717312 . (5. 2) The value (5.2) is determined implicitly by equating two analytical expressions given in Ref. [14]. Given the excellent agreement with our above numerical determination of q 0 (tri) we can therefore conjecture that: For triangular-lattice strips with fully periodic boundary conditions of widths L = 3p, the limit q 0 (tri) = lim p→∞ q 0 (L) exists and is equal to q B of Eq. (5.2). Furthermore, the sequence q 0 (3p) is monotonically increasing in p.
The number of regions enclosed by B having a non-zero intersection with the real q-axis is 2 for L = 2, 3, 4 [cf. Figure 7 and Table 7]. However, for L = 5 a further such region appears, whose intersecting with the real q-axis is the interval (q 2 (L), q 0 (L)), where q 2 (L) ≃ 3.84826329. The value q 2 (L) corresponds to the crossing of the two largest eigenvalues in the ℓ = 2 sector and is tabulated for strip widths L ≤ 12 in Table 7. This table gives compelling evidence that q 2 (L) tends to a limit, q 2 (tri) < 4, independently of the parity of L mod 3.
As above, we first fit the data for q 2 (L) with a fixed parity, L mod 3 = 0, to the form q 2 (L) = q 2 (tri) + AL −1/ν (See Table 8). This is again favourable of ν ≈ 1/2. Repeating the fit with a fixed ν = 1/2 we arrive at the final estimate q 2 (tri) = 3.8194 ± 0.0016. Using the improved Ansatz (5.1), we get the value q 2 (tri) = 3.8199 ± 0.0017. We can therefore conjecture that: Conjecture 5.7 For triangular-lattice strips with fully periodic boundary conditions the limit q 2 (tri) = lim L→∞ q 2 (L) exists and is equal to q B of Eq. (5.2). Furthermore, the sequence q 2 (3p) is monotonically decreasing in p.
The monotonic nature of the sequences q 0 (tri) and q 2 (tri), and the existence of a common limit (as stated in Conjectures 5.6-5.7), means that this common limit is bounded by the finite-widths values: i.e., q 0 (3p) ≤ q B ≤ q 2 (3p) for all p ≥ 2. The data for L = 12 displayed in Table 7 implies that 3.81314 ∼ < q B ∼ < 3.82428 (5.3) In particular, this narrow interval excludes that this common limit would be equal to the value q G = B 12 (1.12) proposed in Ref. [9]. As for the square-lattice case, the limiting curves B L for 2 ≤ L ≤ 7 show some regularities (see Figure 10): for real q, we have up to four different phases belonging to only three different ℓ sectors (ℓ = 0, 1, 2). Again, within each phase, the leading amplitude does not depend on L. All our empirical findings can be summarized in the following conjecture: The phase diagram of the zero-temperature triangular-lattice Potts antiferromagnet with toroidal boundary conditions at real q has four different phases characterized by the number of bridges ℓ as follows: (a) q ∈ (−∞, 0) ∪ (q 0 (L), ∞) characterized by ℓ = 0 with amplitude b (0) = 1 (b) q ∈ (0, 2) characterized by ℓ = 1 with amplitude b (1) = q − 1 (c) q ∈ (2, q 2 (L)) characterized by ℓ = 2 with amplitude b (2) 1 = q(q − 3)/2 (d) q ∈ (q 2 (L), q 0 (L)) characterized by ℓ = 2 with amplitude b Note that the phase (d) only appears for finite width L whenever q 0 (L) > q 2 (L) (e.g. L = 5, 7). But, according to Conjectures 5.5-5.7, this phase does survive in the thermodynamic limit L → ∞. We believe that this phase will appear for all L ≥ 5 such that L mod 3 = 0.
For 5 ≤ L ≤ 7 we further note the appearance of two disjoint triangular-shaped complexconjugate enclosed regions, whose boundary contains three T points, one of which belongs to the branch of B that contains the point q = 2. The dominant eigenvalue inside these enclosed regions comes from the ℓ = 0 sector. Figures 7-8 provide some evidence that these regions may subsist in the thermodynamic limit.
Finally, for L = 6 we note the emergence of a pair of complex-conjugate enclosed regions analogous to those found for the square lattice. More precisely, these enclosed regions contain points for which Re q > 4, and the corresponding dominant eigenvalue comes from the ℓ = 1 sector. We find it likely that such regions will appear for all L = 3p with p ≥ 2, and will merge in the L → ∞ limit so as to form a branch that intersects the real q-axis in q = q 0 (tri) = 4.
There is a analogous conjecture for the triangular lattice to Conjecture 5.4 about one of the eigenvalues of the highest sector: The diagonal sub-block T L,L for a triangular-lattice strip with fully toroidal boundary conditions and width L contains the eigenvalue λ = (−1) L 2 with amplitude b (L) 1 .

Triangular-lattice free energy
In this section we will discuss the free energy for the triangular-lattice strips we have considered and its relation to the analytic eigenvalues found by Baxter [14].
The (real part of the) free energy per unit area f L (q) for a triangular-lattice strip of width L is given by 12 where λ ⋆ (L; q) is the dominant eigenvalue of the transfer matrix T(L), and G is the geometric factor [c.f., [9,10]] relating the free energy per unit area and per spin. With this normalization the FSS of the free energy is given by where f bulk is the bulk free energy per unit area (which is expected to be independent of boundary conditions), c(q) is the (effective) central charge, and the dots stand for higherorder corrections depending on the value of q (in particular, CFT predicts a non-universal L −4 correction term). In Figure 11, we have plotted the two most dominant free energies for each of the three first sectors ℓ = 0, 1, 2 for a triangular-lattice strip of width L = 12 around q = q B (1.11). 13 The other sectors ℓ ≥ 3 are always sub-dominant and play no role in this discussion. The relevant eigenvalue crossing occurs at q = q 0 (L = 12) ≈ 3.8131463701 < q B where the dominant eigenvalue changes from the ℓ = 2 sector (solid black line) to the ℓ = 0 sector (dotdashed red line). There is also a crossing between the two most dominant ℓ = 2 eigenvalues at q = q 2 (L = 12) ≈ 3.82428323 > q B . This crossing gives rise to a new ℓ = 2 phase (between q 0 and q 2 ) for all strips of width not a multiple of three. In both cases, the ℓ = 1 sector plays no role in the phase structure around q = q B . This computation gives numerical evidence that our scenario is correct and no new elements are at play.
Another interesting question is whether the bulk free energy f bulk (q) agrees or not with Baxter's free energies per unit area (1/G) log |g i | (i = 1, 2, 3) [14]. (See Ref. [13, Section 6.1] for a detailed description of Baxter's solution). In particular, Baxter found that the free energy is given by g 3 in the region containing the interval q ∈ (0, q B ), by g 2 in the region containing the interval q ∈ (q B , 4), and by g 1 in the rest of the complex q-plane.
In Figure 11 we observe that the agreement with Baxter's solution is very good: for q ∼ > q B the free energy is very similar to the g 2 contribution (actually it seems to be almost equal to the contribution of the ℓ = 1 sector), while for q ∼ < q B , it is similar to the contribution of g 3 .
A striking difference between our phase-transition diagrams and Baxter's solution is the transition line going through q = 2, which is absent in the latter. This line separates the phase characterized by ℓ = 1 from the one characterized by ℓ = 2 (which agrees with Baxter's solution close to q = q B ). Therefore a natural question is whether there is one eigenvalue missing in Baxter's solution and corresponding to the ℓ = 1 sector. To answer this question we have plotted in Figure 12(a) the ratio of the free energy per unit area f L (q) and Baxter's free energy (1/G) log g 3 (q) in the interval q ∈ (0, 3.5). To avoid parity effects, we have considered only lattices of widths L = 3p with p = 2, 3, 4, 5, and for each width we have considered the free-energy contribution from the three lowest sectors ℓ = 0, 1, 2 (represented by solid, dashed, and dot-dashed curves, respectively).
Note that the eigenvalue crossings occur at the points q = 0, 2. Another exact crossing, involving the second and third most dominant eigenvalues, takes place at q = 1. The fact that the loci of these crossings are subject to no finite-size corrections can presumably be explained using quantum group arguments.
The numerical agreement with g 3 is not very good as the ratio is consistently below the expected result (i.e., 1). The reason is that the FSS corrections (5.6) have not been taken into account. A better agreement can be obtained if we take into account the first correction term in (5.6). We then define the corrected free energy per unit area where c eff (q) is the effective central charge for the BK phase discussed in the next section; see (5.14). We expect that this value will be closer to the thermodynamic limit f bulk (q), as the correction terms are o(L −2 ). This plot is shown in Figure 12(b). Our first observation is that the values of the ratio are much closer to the expected value of 1 than in Figure 12(a). This shows that our assumption that CFT holds for this model is valid. Note that the eigenvalue crossings do not need to be at the expected points q = 0, 2, as the correction in (5.7) depends on both q and L. The second observation is that the agreement of the ℓ = 2 free energy with the contribution of g 3 is very good for q ∼ > 1 (less than 0.02% for L = 9, 12). For q ∼ < 1, the best agreement corresponds to the ℓ = 0 sector. However, it is difficult from this data alone to draw any definitive conclusion about the necessity of including an extra eigenvalue to the set introduced by Baxter.

Conformal properties
Following Refs. [29,30] we expect a number of the conformal properties of the chromatic polynomial to be given by the BK antiferromagnetic phase. In brief, the situation is the following. 14 The two-dimensional Potts model is exactly solvable (i.e., integrable) on the curves v = ± √ q (square lattice) v 3 + 3v 2 = q (triangular lattice) (5.8) The part of the curves with 0 ≤ q ≤ 4 and v ≥ 0 corresponds to the usual secondorder transition between a ferromagnetic and a paramagnetic phase. Its critical properties are given by the Coulomb gas (free field) construction with coupling constant g ∈ [0, 1/2]; in particular the thermal operator (corresponding to a perturbation of the critical point in the v variable) is relevant in the RG sense. Now, the analytic continuation to the region with g ∈ (1/2, 1) describes the lower (resp. middle) branch of the curves (5.8) for the square (resp. triangular) lattice. Crucially, at these critical points the thermal operator is irrelevant. Therefore, for each value of q ∈ (0, 4), these critical points will act as an RG attractor for a finite interval (v − (q), v + (q)): this is the BK phase. For the square lattice, the end-points of these intervals are given by the two mutually dual antiferromagnetic phase transition points [31] v ± (q) = −2 ± 4 − q (5.9) whereas for the triangular lattice v ± (q) are presently unknown. The upshot is that the part of the chromatic line v = −1 that intersects the BK phase will have its critical properties determined by the latter. For the square lattice this happens for Re q ∈ (0, 3), according to Eq. (5.9), in perfect agreement with Conjecture 5.2. For the triangular lattice, we obtain a consistent scenario if we suppose that the intersection with the BK phase is for Re q ∈ (0, q 2 ), with q 2 given by Conjecture 5.7.
Within the BK phase, we can parametrize q by q = 4 cos 2 π p (5.10) where p ∈ (2, ∞). The central charge is then Furthermore, the theory possesses an infinite set of bulk operators O k , with k = 1, 2, . . ., of conformal weight [42] h k (p) = Geometrically, the insertion of an O k operator at either end of the strip corresponds to a situation in which 2k cluster boundaries propagate along the length direction of the strip.
(The boundaries may of course also wind around the periodic transverse direction.) For k ≥ 2 we may equivalently think of k propagating clusters. It should be carefully noted that the case k = 1 is special: a single propagating cluster may either correspond to two or zero propagating cluster boundaries. The latter case will occur when the cluster is allowed to have cross-topology [43], i.e., to percolate across both periodic lattice directions. The corresponding bulk operator is denotedÕ 1 and has conformal weight Indeed, by duality this is nothing else than the order parameter (magnetization) operator. In the probabilistic regime (i.e., v > 0) the identity operator I is the most relevant (all other conformal weights are strictly positive). Similarly, one hash 1 < h 1 , translating the fact that a single cluster is indeed most likely to have cross-topology. In the BK phase, v < 0 and some of these statements need no longer be true, due to the presence of negative conformal weights. In fact, for any p > 2 the weight h 1 is the most negative, and so one would naively expect that the most "likely" configurations are those with two propagating cluster boundaries. This is however not true, due to a subtle effect: the amplitudes of the corresponding transfer matrix eigenvalues vanish identically on the torus [30]! So the remaining candidates for the most relevant operators areÕ 1 and O 2 . From Eqs. (5.12)-(5.13) it is easy to see that the former (resp. the latter) will dominate for 2 < p < 4 (resp. for 4 < p < ∞). This argument is of course valid in the thermodynamic limit only, but the conclusion will in fact hold true in finite size due to the underlying quantum group symmetry of the model. Accordingly we indeed observe numerically, for any width L and for both lattices, that branches of the limiting curves B cross the real q-axis at q = 0 (p = 2) and q = 2 (p = 4); the former (resp. the latter) crossing corresponds to the equimodularity of the dominant eigenvalues of the ℓ = 0 and ℓ = 1 (resp. the ℓ = 1 and ℓ = 2) sectors. However, the phase containing q = 2 + does not extend to q = 4 (p = ∞) for either lattice, since the chromatic line leaves the BK phase at q 0 = 3 (square lattice), respectively q 2 ≃ 3.8196717312 (triangular lattice).
We can further check the validity of this scenario by turning to the amplitudes of the dominant eigenvalues. For the order parameter operatorÕ 1 the amplitude is q − 1 (by a simple counting argument), and this is indeed the amplitude of the dominant eigenvalue observed numerically for p ∈ (2, 4) for both lattices.
We have checked numerically that this is true on both lattices (using strips of width L ≤ 12) to a very good precision; this check is essentially the contents of Figure 12(b). A final elusive point is the value of c in the interval q ∈ (q 2 , q 0 ) on the triangular lattice. The q 0 (tri) = 4 case is known to be equivalent to Baxter's three-colouring problem on the hexagonal lattice and hence have c = 2. Our numerical data seem to indicate that the rest of the interval (q 2 , q 0 ) may also correspond to a critical theory, the properties of which are unrelated to those of the BK phase.
We have made an attempt to obtain an expression for the central charge in the interval q ∈ (q 2 , q 0 ) by performing some fits to our free-energy data. Our first two-parameter Ansatz was given by Eq. (5.6) without any correction-to-scaling term. In Figure 13 we show our estimates as a function of q for different values of L min = 6, 9, 12. Much better estimates can be obtained by including higher-order corrections to the above Ansatz. In previous work [9,10,13], we have found that including a correction term ∼ L −4 greatly improves the estimates for the central charge. Thus, we have used the improved three-parameter Ansatz In Figure 13 we have also displayed the estimates for c(q) obtained with this Ansatz and L min = 6, 9. Their stability is superior to that of the two-parameter estimates. The solid line represents a guess for the exact value of the (effective) central charge in this regime (see (5.19) below) where p ∈ (2, ∞) is related to q by (5.10).
Let us now explain in detail how we have obtained the expression (5.16). The functional form of the Ansatz (5.16) is motivated by already known results of the central charge in other regimes [c.f., (5.11)]. Thus, our goal is to obtain the central charge c(q) as a rational function of p with integer coefficients. We started with the estimates coming from the three-parameter Ansatz (5.15). It is more useful to plot the results as a function of 1/p, so that q = 4 (p = ∞) corresponds to the origin. In Figure 14(a) we display the difference c(q = 4) − c(q) versus 1/p for the two available data sets L min = 6, 9. Both data sets fall approximately on a single curve, especially on the interval corresponding to q ∈ (q 2 , q 0 ). This empirical result implies that the largest FSS correction to our estimates for c(q) corresponds to a global shift. We then fit the data for c(q = 4) − c(q) to a polynomial Ansatz. Our best result is given by This Ansatz is motivated by several facts: by definition, the curve should go through the origin, and it is clear in Figure 14(a) that its slope at 1/p = 0 is zero. Furthermore, the addition of more terms to (5.17) only leads to less stable estimates. As our goal is to obtain a rational function of p, we computed the Padé approximant [2,2] of the above polynomial (5.17). The result is given by Finally, as we expect that the coefficients of the above rational function are integers, we systematically searched the set of closest integers which best approximates the data of In Figure 14(a) it is clear that both expressions give an excellent fit to the whole set of data points. In order to tell which integer-coefficient rational approximant was best, we plotted the difference between our estimates and the fits. In Figure 14(b) we show this difference for the above approximants (5.19). We observe that the approximant (5.19a) is slightly superior to the other one (5.19b). The final step is to fix the overall shift, i.e., the value of c(q = 4). As the triangular-lattice Potts antiferromagnet at zero temperature is critical and can be mapped onto a Gaussian model with two bosons [44,45], we expect that c(q = 4) = 2 (5.20) Thus, we arrive at the final expression (5.16). Indeed, this computation should be taken with a grain of salt: we have found a function with the expected properties (rational function of p with integer coefficients) that describes rather accurately our estimates for the central charge of our model on the interval q ∈ (q 2 , q 0 ). This does not mean that this is unique possible solution with those properties nor that the right expression for c(q) should satisfy those properties. 15 Another important point is that our conjectured expression for the central charge, although it was first devised to deal with the interval q ∈ (q 2 , q 0 ), is able to reproduce rather accurately the data outside the latter interval.

The "q = 2 branch"
At the end of Section 5.3, we have already discussed the branch of the triangular-lattice limiting curve B L going through the point q = 2, and separating the ℓ = 1 and ℓ = 2 phases. In this subsection we will study further properties of this branch in both the square and triangular limiting curves.
The first question we want to address is that of the density of chromatic zeros along this branch for finite strips of size L × N. In order to estimate this density, we first fixed a region in the complex q-plane containing the q = 2 branch. This region is In Figure 15, we show the density of the chromatic zeros ρ in the regions (5.21) for some finite strips of size L × N as a function of the aspect ratio L/N. At fixed aspect ratio, we expect the density ρ to converge to some (small) positive constant. In Figure 15, we observe that for fixed L/N, the density is decreasing in L, and it attains at L = 7 very small values ∼ < 0.02. If this monotonicity holds for larger values of L, we expect that the density of zeros along the q = 2 branch will converge to zero in the thermodynamic limit (in contrast with the behavior of the density of zeros in other parts of the limiting curve). Another way to tackle this problem is consider the parameter t introduced in Ref. [11,Section 4.1.1]. On the limiting curve B L , there are generically two equimodular leading eigenvalues λ 1 and λ 2 . Indeed, their ratio is just a phase In Ref. [11], instead of θ, we found it more convenient to work with the real parameter t, defined as In practice, we use this parameter to ensure the plots of the limiting curves B L are correct and we do not miss any small feature. Note that t = θ = 0 corresponds to two identical eigenvalues. At q = 2, we have that t = 0 for any L ≥ 3; but as we move along the q = 2 branch, t increases. Let us now study how the parameter t evolves (along this branch) as a function of |q − 2|. In Figure 16, we show (for both the square and triangular lattices) the values of t along the q = 2 branch for points in the complex q-plane with Im q = 0, 0.01, . . . , 0.10. It is clear that for fixed L, the relation is approximately linear. For the square lattice, the slope is given approximately by 1/L, while for the triangular lattice the slope is given approximately by √ 3/(2L) + 0.653/L 3 . In both cases, the discrepancies between the estimated slopes and the above results are of order O(10 −4 ). Thus, we expect that in the thermodynamic limit both leading eigenvalues λ 1 and λ 2 become equal along this branch.
Finally, we can study the behavior at q = 2 by considering the series expansion of the two equimodular leading eigenvalues λ 1 and λ 2 (corresponding to the ℓ = 1 and ℓ = 2 sectors, respectively). Indeed, such series expansions are expected to exist, as the eigenvalues λ i (q) are algebraic functions of q. On each phase, we can compute the series expansion of the leading eigenvalue at q = 2 + ǫ For L = 2, 3, this can be done symbolically by using mathematica. For 4 ≤ L ≤ 7, we can only obtain (again with mathematica) the series numerically, as the symbolic expressions for the coefficients are too lengthy. In practice, we are able to obtain the first three coefficients λ i,k (k = 0, 1, 2) in (5.24). 16 With these coefficients, we can obtain the corresponding series expansion for the free energy: For 4 ≤ L ≤ 12, we numerically obtained the contribution to the free energy f L,i (q) = (1/L) log |λ i (q)| of the two leading eigenvalues λ i close to q = 2. In particular, we obtained the values f L,i (q = 2) and f L,i (q = 2 ± 10 −5 ) for i = 1, 2. From these six values, we estimated the first three coefficients f (i) L,k (with k = 0, 1, 2) in (5.25). Indeed, we know that the free energy is continuous at q = 2, thus we expect that ∆f L (2) = f (2) L,0 − f (1) L,0 = 0, as we find in the actual computation. We are interested in computing the increment of the first derivatives of the free energy f L (q) at q = 2. In practice, we have computed the increments of the first 16 For the triangular lattice, this computation could only be performed up to width L = 6. two derivatives: 17 If E L = 0, then the transition at q = 2 is of first order. However, if E L = 0 and C L = 0, then the transition is of second order (and this argument can be generalized to higher-order phase transitions). In Figure 17 we have plotted the coefficients E L and C L versus 1/L for both the square and triangular lattices. The plots suggest that both coefficients for both lattices tend to a limit as L → ∞. If we fit the E L -data to a power law E L = E ∞ + A E L −ω E , we find that |E ∞ | ∼ < 10 −4 is indeed very small (in modulus) for L min = 10: In particular, our square-lattice data can be fitted very well with the simple Ansatz: E L (sq) = −2L −2 . These fits are depicted in the corresponding plots in Figure 17. We have chosen L min to be the largest possible value, so that there are no effective degrees of freedom in the fits. The power-law fits C L = C ∞ + A C L −ω C are not as stable as the preceding ones: we find sizeable parity effects in the estimates, so we have to split the whole data set into two smaller sets with well-defined L-parity. Thus, for each value of L min , we have fitted the data with L = L min , L min + 2, L min + 4. The estimates for even L are given for L min = 8 by: Our preferred fits are those given by (5.28), as L min is larger; these fits are depicted in the corresponding plots in Figure 17. They also show show that C L tend to some small (in modulus) number in the thermodynamic limit |C ∞ | ∼ < 5 × 10 −4 . We conjecture that for both lattices, E L , C L → 0 as L → ∞. Indeed, this does not mean that there is no transition at all at q = 2 in the thermodynamic limit, but that the transition (if any) is at least of third order.
In this section we have presented three numerical approaches to study the thermodynamic limit along the q = 2 branch. None of them is conclusive; but taken together they give strong 17 For 4 ≤ L ∼ < 7, we had two distinct ways of estimating the coefficients E L and C L , namely symbolically and numerically. We checked that the difference between those estimates was ∼ < 10 −11 for E L , and ∼ < 10 −6 for C L . support to the following scenario: even though at finite L there are two different phases coexisting along the q = 2 branch, in the thermodynamic limit there is a unique phase. This scenario is consistent with the three results discussed above: in the thermodynamic limit L → ∞: 1) The density of chromatic zeros goes to zero along this branch; 2) The ratio of eigenvalues along this branch tend to one λ 2 /λ 1 → 1; and 3) At q = 2 the increment in the first two derivatives of the free energy f L (q) (with respect to q) tends to zero. Thus, this branch seems to be somehow an artifact of the boundary conditions; in particular, of the periodic longitudinal boundary conditions, as this branch also appears in the cyclic case [9], and it is absent for free and cylindrical boundary conditions [11][12][13]. 18 This conclusion agrees with Baxter's solution [14] in the sense that only three phases survive in the phase diagram of the triangular-lattice chromatic polynomial.

Conclusions
In this paper we have extended our transfer-matrix studies [9,[11][12][13] to toroidal boundary conditions. We have obtained the exact expression of the chromatic polynomial for squareand triangular-lattice strips of fixed width 2 ≤ L ≤ 7 and arbitrary length N. The cases L = 5, 6, 7 are new in the literature (note that going to larger widths L ≥ 8 would require better computer facilities and/or improved symbolic algorithms). We have also obtained the limiting curves B L and isolated limiting points for these strips in the infinite-length limit N → ∞. The limiting curve B L for a fixed value of L provides the exact phase diagram for this semi-infinite strip.
By studying the features of these finite-L phase diagrams and using techniques from CFT and FSS, we have been able to determine some properties of the chromatic-polynomial phase diagram in the thermodynamic limit L → ∞. These findings are summarized in the following points: • Toroidal boundary conditions appear to give results qualitatively and quantitatively closer to those of the thermodynamic limit than the other standard boundary conditions (free, cylindrical, and cyclic) [9,[11][12][13]. In particular, the exact value of q c is obtained for all square lattices of (finite) odd widths, and for all triangular lattices of (finite) widths not a multiple of three.
• As for cyclic boundary conditions [9], we find that the limiting curves B L always enclose regions of the complex q-plane, thus determining true phases as in the expected infinitevolume phase diagram. Each phase is also characterized by the number of bridges ℓ 18 A similar (but not related) phenomenon appears in a recent study of the RSOS [8] representation of the Potts model when q = B p with p = B p (1.8). In the complex v-plane, the large-|v| region (which corresponds to the high-temperature regime of the model and is expected to be non-critical) displays several non-bounded outward branches separating different phases. The origin of these branches is simple: there are two almost-degenerate eigenvalues, each of them corresponding to a different topological sector. These two eigenvalues become exactly equimodular along those non-bounded outward branches, whose number increases with L. While for any finite L there are several phase-transition lines in the large-v region, in the thermodynamic limit we expect that both eigenvalues become exactly identical over all this region. Thus, we recover the standard result that there are no phase-transition lines in the large-v (high-temperature) regime of this model. associated to the leading eigenvalue; but for toroidal boundary conditions only the ℓ = 0, 1, 2 sectors contribute for both lattices. (For the triangular lattice with cyclic boundary conditions we expect to have contributions up to ℓ = 6 bridges). The other difference with the cyclic triangular strip is that there are now two different phases characterized by ℓ = 2 bridges (For cyclic boundary conditions, there is at most one phase for each value of ℓ).
• It is rather remarkable that for cyclic boundary conditions [9,29], all the Beraha numbers B p (1.8) played a special role in the BK phase, either as phase transition points (even p) or as isolated limiting points (odd p). In contrast, for toroidal boundary conditions only the finite subset of integer q play such special roles: there are phase transitions at q = 0, 2, 4 and isolated limiting points at q = 1, 3. This observation alone is a nice illustration of the crucial influence of boundary conditions on antiferromagnetic systems. Note however, that we still expect (and have to some extent numerically observed) that the full set of B p will be the loci of level crossings of sub-dominant eigenvalues.
• Seeing that the maximum chromatic number H(g) of a graph G increases with the genus g of the surface that it is embedded in according to (1.5), one might have expected that the toroidal topology studied here would have revealed further structure of the limiting curves for q ∈ (4, 7). This is not the case. The reason is that it is precisely the interplay between the properties of the BK phase and the representation theory of the quantum group U t (SU(2)) for t a root of unity that is responsible for generating the partition-function zeros. This mechanism is confined to the interval q ∈ [0, 4]. The graphs having chromatic zeros for q ∈ (4, 7) must therefore be rather exceptional and have little to do with the recursive families of graphs studied here. One might even speculate that the number of such exceptional graphs might be small in some sense, or even finite. These conclusions are expected to carry over to higher genus as well.
• For toroidal triangular strips we have found two phase transitions inside the critical region: q = 2 and q = q B (1.11).
This result is rather puzzling when compared to Baxter's exact solution [14] and the re-analysis discussed in Ref. [9] of Baxter's eigenvalues. On one hand, the value of q B coincides with that predicted by Baxter based on the result that only three eigenvalues are relevant in the complex q-plane. This conclusion stems from the Bethe-Ansatz solution of a closely related model [14]. The mapping between these two models is exact, except for the boundary conditions. Usually boundary conditions have no effect in the ferromagnetic regime; but here we are deep in the antiferromagnetic one.
On the other hand, a careful and high-precision analysis of these three eigenvalues [9] revealed that, although the outer boundary of Baxter's infinite-volume curve B ∞ was indeed correct, the inner structure was richer: there were infinitely many more branches than those included by Baxter and those indicated by the finite-strip chromatic zeros.
In particular, the point q B (1.11) does not belong to the limiting curve B ∞ , and its role is now played by q G (1.12).
• In Section 5.5 we have found some numerical evidence that the branch of the limiting curves B L going through q = 2 might disappear in the true thermodynamic limit L → ∞. This may provide an explanation of why this branch is absent from the triangular-lattice limiting curve B ∞ proposed by Baxter [14]. However, it is still not clear how to explain the existence of the q = 2 branch for finite widths from Baxter's analysis. The resolution of these discrepancies would most likely involve identifying (at least) one more analytic expression for the free energy from the Bethe Ansatz equations of Ref. [14]. Such a study might also shed some light on the possible existence of additional enclosed regions in B ∞ , such as the triangular-shaped regions discussed towards the end of Section 5.2.
• In the triangular-lattice case we obtain a novel behavior: one important feature of the infinite-volume phase diagram does depend on the parity of the width (i.e., on L mod 3). Our numerics strongly suggests that if L mod 3 = 0 (6.1) In the square-lattice case, we have that this limit for even widths is equal to q c (sq) = q 0 (L = 2p + 1) for all p ≥ 1.
• In the light of the numerical results given at the end of Section 5, we might finally speculate that the full solution of the triangular-lattice chromatic polynomial (in the sense of a Bethe Ansatz solution) should reveal not only two distinct bulk free energies for q ∈ [0, 4], but also two different sets of associated critical scaling levels throughout that interval. In other words, we suggest that there exists two "superposed" critical theories throughout the interval, one having the physics of the BK phase, and the other being some deformation of the two-boson theory [44,45] responsible for c = 2 at q = 4. This is a strong motivation for a more refined study of the Bethe equations corresponding to Baxter's g 2 solution [14]. Note that these expectations apply also to the Nienhuis loop model [41].

A Two combinatorial identities
In this section we prove the following combinatorial identities that were used in Section 2.2.  TriTorus SqTorus  1  2  5  2  2  2  2  2  2  15  52  7  5  4  4  4  3  200  780  31  13  8  11  8  4  3185  12838  361  99  68  38  33  5  52920  214956  3241  653  347  148  90  6  884268 3593700 30916  5194  2789  476  325  7 14723280 59751978 286924  40996  20766   Table 1: Dimensionality of the transfer matrix for toroidal boundary conditions. For each width m, we give the number of non-crossing connectivities C m,m (2.17), the number of non-crossing connectivities C m+1,m (2.19), the number of non-crossing non-nearest-neighbor connectivities D m,m , the number of generated connectivity classes modulo reflections (resp. reflections and translations) of non-crossing non-nearest-neighbor connectivities TriTorus ′ (resp. SqTorus ′ ), and the number of distinct eigenvalues for a triangular-lattice (resp. squarelattice) strip with toroidal boundary conditions TriTorus (resp. SqTorus). ′  2  1  2  1  4  3  1  2  4  1  8  4  5  22  33  7  1  68  5  6  67  190  72  11  1  347  6 39 494 1432 644 163 16  1  2789  7 111 2794 10224 5615 2711 288 22 1 20766 Table 2: Block structure of the transfer matrix of a toroidal square-lattice strip of width m as a function of number of bridges ℓ. For each strip width m, we quote the dimension N ℓ of the diagonal block T ℓ,ℓ (m T ) for a given number of bridges ℓ, and the total dimension of the transfer matrix SqTorus ′ (m) = dim T(m T ).  Table 3: Block structure of the transfer matrix of a toroidal triangular-lattice strip of width m as a function of number of bridges ℓ. For each strip width m, we quote the dimension N ℓ of the diagonal block T ℓ,ℓ (m T ) for a given number of bridges ℓ, and the total dimension of the transfer matrix SqTorus ′ (m) = dim T(m T ).   Table 5: Number of distinct eigenvalues for a toroidal triangular-lattice strip of width m as a function of number of bridges ℓ. For each strip width m, we quote the number of different eigenvalues L ℓ for a given number of bridges ℓ, and the total number of distinct eigenvalues TriTorus(m). Notice that the sum m ℓ=0 L ℓ (m) may be greater than TriTorus(m), as a few eigenvalues can appear on different ℓ blocks. The results marked with a † are conjectures based on the numerical evaluation of the corresponding eigenvalues.
q 0 (L) q 2 (L) Ansatz: q k (L) = q 0 (tri) + AL −1/ν L min q 0 (tri) A 1/ν q 0 (tri)  Table 8: Summary of fits of the q 0 (3p) and q 2 (3p) data (1 ≤ p ≤ 4) for the zero-temperature triangular-lattice Potts antiferromagnet. For each data set and Ansatz, we show the different estimates as a function of the parameter L min . See text for an explanation of how the fits were performed.       Figure 4: Connectivity state for a toroidal strip of width 4. (a) We show the connectivity state P = δ 1 ′ ,3 ′ δ 2,4,4 ′ , which can be seen as a bottom-row connectivity state P bot = δ 1 ′ ,3 ′ (lower solid thick line), a top-row connectivity state P top = δ 2,4 (upper solid thick line), and a bridge connecting both rows δ 4,4 ′ (dashed thick line). The thin solid lines show the top (resp. bottom) row containing the unprimed (resp. primed) sites. The thick arrow shows the transfer direction (upwards). (b) The same connectivity state depicted in such a way that the transverse periodic boundary conditions are clearer. The transfer direction is outwards (thick arrow).  : Limiting curves for square-lattice strips of width L and toroidal boundary conditions. We show the results for L = 2 (black), L = 3 (red), L = 4 (green), L = 5 (blue), L = 6 (light blue), and L = 7 (pink). We also show the limiting curve for a square-lattice strip of width L = 10 with cylindrical boundary conditions (orange) [12]. Figure 10: Limiting curves for triangular-lattice strips of width L and toroidal boundary conditions. We show the results for L = 2 (black), L = 3 (red), L = 4 (green), L = 5 (blue), L = 6 (light blue), and L = 7 (pink). We also show the infinite-width limiting curve B ∞ obtained by Baxter [14] (orange). Figure 11: Free energy per unit area f L (q) for the triangular-lattice chromatic polynomial on an infinitely long strip of width L = 12 and toroidal boundary conditions. We show the free energy f L (q) associated to a selected set of eigenvalues of the transfer matrix T(12 T ) as a function of q ∈ (3.5, 4). For each sector ℓ = 0, 1, 2, we show the dominant and the first sub-dominant eigenvalues respectively with colors red/orange for ℓ = 0, blue/cyan for ℓ = 1, and black/dark gray for ℓ = 3. Those eigenvalues that are decreasing (resp. increasing) in q are depicted with a solid (resp. dot-dashed) curve. For comparison, we also show Baxter's free energies (1/G) log g 2 and (1/G) log g 3 [14] normalized per unit area and depicted as solid pink and violet, respectively. The vertical dashed dark gray line shows Baxter's value q B (1.11).  . We show two estimates coming from the three-parameter fits (5.15) with L = 6, 9, 12 (black •) and L = 9, 12, 15 (red ). The (almost coincident) solid lines represent our two best candidates for the central charge (5.19a) (blue), and (5.19b) (red). (b) Error between our numerical estimates for the conformal charge and our conjectured expressions (5.16). We show the difference c(q) − c fit (q) between our three-parameter fits and our conjectured expressions c fit (q). Points in black (resp. red) are obtained from the fits with L = 6, 9, 12 (resp. L = 9, 12, 15). Solid [resp. empty] symbols are obtained with c fit given by (5.19a Figure 17: Coefficients of the series expansions for the increment of the first (a) and second (b) derivatives of the free energy f L (q) with respect to q evaluated at q = 2. In each plot we show the corresponding coefficient for the square (red •) and triangular (blue △) lattices as a function of 1/L. The solid lines depict the best power-law fits to the corresponding data set (5.27)/(5.28).