Bootstrap approach to geometrical four-point functions in the two-dimensional critical $Q$-state Potts model: A study of the $s$-channel spectra

We revisit in this paper the problem of connectivity correlations in the Fortuin-Kasteleyn cluster representation of the two-dimensional $Q$-state Potts model conformal field theory. In a recent work [M. Picco, S. Ribault and R. Santachiara, SciPost Phys. 1, 009 (2016); arXiv:1607.07224], results for the four-point functions were obtained, based on the bootstrap approach, combined with simple conjectures for the spectra in the different fusion channels. In this paper, we test these conjectures using lattice algebraic considerations combined with extensive numerical studies of correlations on infinite cylinders. We find that the spectra in the scaling limit are much richer than those proposed in [arXiv:1607.07224]: they involve in particular fields with conformal weight $h_{r,s}$ where $r$ is dense on the real axis.


Introduction
In a very interesting recent paper [1], a proposal was put forward for some of the four-point correlation functions of the percolation problem in two dimensions. This proposal was part of a more general conjecture addressing various geometrical objects involving four points in the diagrammatic formulation [2] of the Qstate Potts model [3]. The case Q = 1 corresponds to percolation, and the proposal in [1] covers such objects as the probability that two of these points belong to one cluster, and the two others to another cluster. Obtaining closed-form expressions for such objects is one of the holy grails in the field. It is a far from obvious endeavour because the conformal field theory (CFT) describing percolation (and more generally geometrical features of the Q-state Potts model) is not well understood: it is non-unitary, probably involves logarithms (even for Q generic), and involves operators which are not degenerate, precluding the use of the differential equations approachà la BPZ [4].
The construction in [1] is elegant and powerful. It starts with a seemingly reasonable hypothesis for the spectrum of operators appearing in the fusion channels for the fusion of two order operators, and determines, using a clever code, the whole set of structure constants based on our knowledge of conformal blocks [5,6] together with the imposition of crossing symmetry. The results are then checked against Monte Carlo simulations, with, it is claimed, reasonably good agreement.
Although the results in [1] are appealing, they are not really consistent with what is known about the Potts model CFT and, in particular, percolation. Early work [7] has revealed indeed a much richer spectrum than the one postulated in [1], which covers only a very tiny set of the known full operator content of the theory. Of course, it could be that by some accident, the order operator in the Q-state Potts model does not couple to as many fields as one would expect, at least in the scaling limit. But it could also be that something is simply missing in the work of [1], despite the apparent numerical effectiveness of their proposal.
To investigate this question requires a long and detailed analysis, of which we present the results here. In a nutshell, we have gathered direct, in our opinion unquestionable evidence that the spectrum of the Q-state Potts model is as complex as could have been feared, that many more fields appear in the OPE of order operators in the Potts model than was conjectured in [1], and that the proposal in that paper, appealing as it may be, simply cannot be correct. It is, at best, a good numerical approximation to the true expressions for the four-point functions.
Our paper is organised as follows. In section 2 we remind the reader of basic facts and results about the Q-state Potts model and its geometrical formulation. Algebraic aspects-which constitute a crucial part of our approach-are discussed in section 3. Section 4 summarises our method of analysis, and how we extract exponents as well as amplitudes, from lattice data. Section 5 discusses our results for the spectra in the intermediate channels of four-point functions. A comparison with results in [1] is provided in section 6. In section 7, we return to the issue of divergences in the amplitudes, and re-analyse briefly our results as well as those of [1] from the point of view of degeneracies, and, potentially, logarithmic CFTs.
Since a good part of our analysis is based on extracting amplitudes from lattice data, a lot of technical aspects have to be considered both to make the program possible, and to check its validity and its limits. We have thus gathered quite a bit of material in a series of appendices. Appendix A discusses in detail many aspects of the numerical algorithms and other techniques used to obtain our results, while appendix B goes over a series of checks, including detailed comparisons, in particular, with known results for Q = 0, 2, 4.

Potts model and its correlation functions
We consider the Q-state Potts model [3] defined on a graph G = (V, E) with vertices V and edges E. There is a spin σ i = 1, 2, . . . , Q attached to each vertex i ∈ V and an interaction energy −Kδ σi,σj attached to each edge (ij) ∈ E. The partition function (in units where the inverse temperature is absorbed into the coupling constant K) reads Note that this initial formulation supposes Q to be a positive integer, Q ∈ N. This constraint can be lifted in a rewriting of Z due to Fortuin and Kasteleyn (FK) [2]. Indeed, write e Kδσ i ,σ j = 1 + vδ σi,σj with temperature parameter v = e K − 1, expand the product (ij)∈E , and perform the sum over all spins {σ} to obtain where the sum is over all 2 |E| subsets of E, and |A| denotes the number of edges in the subset. Moreover, k(A) denotes the number of connected components (also called FK clusters) in the subgraph G A = (V, A).
In the remainder of the paper we take G to be the two-dimensional square lattice. The temperature parameter will be taken at its critical value, v c = √ Q [3,8], so that the model is conformally invariant in the continuum limit. In this latter limit, we are interested in the geometry of the infinite plane, so that boundary effects are immaterial. We shall often use the trick of transforming this into the geometry of a cylinder, via an appropriate conformal mapping (details will be given below). This cylinder geometry is convenient for imposing the lattice discretisation, which is our main tool of algebraic and numerical investigations. In that case we always take the square lattice G to be axially oriented with respect to the cylinder axis, so that the row-to-row transfer matrix describes the (imaginary) time evolution of L Potts spins.

Loop model
An equivalent formulation of Z is given [9] by the loop model on the medial lattice M(G) = (V M , E M ). The vertices V M of M(G) are situated at the mid-points of the original edges E, and two vertices in V M are connected by an edge in E M whenever the former stand on edges in E that are incident on a common vertex from V . In particular, when G is a square lattice, M(G) is just another square lattice, tilted through an angle π 4 and scaled down by a factor of √ 2. There is a bijection between edge subsets A ⊆ E and completely-packed loops on M(G). The loops are defined so that they turn around the FK clusters and their internal cycles (alternatively they separate the FK clusters from their duals). One has then [9] where (A) denotes the number of loops. The loop fugacity is where the (quantum group related) parameter q will be used intensively below. The loop model will be convenient to make contact with the Temperley-Lieb (TL) algebra [10], which will be discussed below. Note also that on a square lattice, we have simply vc n = 1, so at the critical point Z depends only on (A).

Correlation functions
The Potts model allows for the definition of various correlation functions, depending on whether one uses its formulation in terms of Potts spins, FK clusters or TL loops. 1 The spin correlators are naturally defined in terms of the order parameter (or spin) operator O a (σ i ) ≡ Qδ σi,a − 1 .
More interesting and general results can however be obtained by moving to the cluster or loop formulations, in which the correlation functions acquire a geometrical content. In the same vein, the spin correlators can be analytically continued from Q ∈ N to arbitrary real values, in which case they also acquire a geometrical interpretation [12][13][14] to which we shall return in a short while. Such generalisations to Q ∈ R are not only useful, but actually indispensable in our case, since our main objective is to study correlation functions in the generic case where q is not a root of unity. Correlation functions in the loop formulation are of either electric or magnetic type, where the terminology refers to the Coulomb gas approach to CFT [7,15]. Let i 1 , i 2 , . . . , i N be a number of distinct marked vertices. Electric correlators are defined for i k ∈ V by appropriately modifying the weight of loops that contain a subset of marked vertices on their inside, and the remainder on their outside. Magnetic correlators are defined for i k ∈ V M by specifying whether given vertices belong to the same or different loops; one can also increase the set of possibilities by allowing for topological defects that insert a number of open loop strands at each marked vertex [16]. While these electromagnetic correlation functions have been intensively studied for N = 1, 2 in a variety of contexts, we wish here to recall only one recent result. Namely, the electric N = 3 correlation functions have been shown to be related, for generic values of n ∈ [0, 2], to the so-called DOZZ formula for the structure constants within Liouville field theory [14,17,18].
Our main interest here is however correlation functions defined in terms of the FK clusters. Let again i 1 , i 2 , . . . , i N ∈ V be a number of distinct marked vertices, and let P be a partition of a set of N elements. We then define where Z is given by (2), and I P (i 1 , i 2 , . . . , i N |A) is the indicator function that, ∀k, l ∈ {1, . . . , N }, vertices i k and i l belong to the same connected component in A if and only if k and l belong to the same block of the partition P. It is convenient to denote P by an ordered list of N symbols (a, b, c, . . .) so that identical symbols refer to the same block. For instance, with N = 2, P aa is the probability that vertices i 1 , i 2 belong to the same FK cluster, whereas P ab = 1 − P aa is the probability that i 1 , i 2 belong to two distinct FK clusters. In the context of four-point functions, we are therefore interested in the 15 probabilities P aaaa , P aabb , . . . , P abcd . The combinatorial properties of FK correlation functions were further discussed in [14]. It is natural to relate P P to correlation functions of the spin operator. Define where a 1 , a 2 , . . . , a N is a list of (identical or different) symbols defining a set partition P, and the expectation value · · · is defined with respect to the normalisation Z. It is straightforward to formally relate the G P to P P . Indeed, to evaluate the expectation value of a product of Kronecker deltas, we initially suppose that Q is integer, and use that spins on the same FK cluster are equal, while spins on different clusters are statistically dependent. This leads to Q-dependent relations, which can finally be extended to real values of Q by analytical continuation. For instance, with N = 2, one readily finds that G a1,a2 = (Qδ a1,a2 − 1) P aa .
In other words, the two-point function of the spin operator is proportional to the probability that the two points belong to the same FK cluster. Therefore O a (σ i ) effectively "inserts" an FK cluster at position i ∈ V and ensures its propagation until it is "taken out" by another spin operator.
Remark 1. In a recent series of works [19][20][21] we have introduced a more general class of operators O a1,a2,...,a N (σ i1 , σ i2 , . . . , σ i N ) that act on N spins according to given irreducible representations of the symmetric groups S Q and S N . These operators can enforce the propagation of more than one FK cluster, with the set of propagating clusters having specific symmetry properties. Some of the four-point functions to be considered below (namely P abab ± P abba ), with the points being considered as regrouped in two pairs, actually coincide with two-point functions of such operators, each acting on a pair of spins (N = 2).
As already stated above, for Q arbitrary, the left-hand sides of these equations are only formally defined: it is in fact the right-hand sides that give them a meaning. Note that this linear system has determinant Q 4 (Q − 1)(Q − 2) 3 (Q − 3), so it cannot be fully inverted for Q = 0, 1, 2, 3. By analogy with (8) one would expect that, in the scaling limit, the four P P of interest would be described by combinations of conformal blocks for the spin operator. In particular, the function P aaaa corresponding to the four points being in the same cluster should become, in the scaling limit, a crossing-invariant such combination. The other three would maybe not be crossing-invariant individually, but might be related with each other by crossing (or give rise, after proper combinations, to other crossing-invariant objects).
Clearly, to implement the bootstrap programme, one needs an idea of the set of conformal blocks that may appear in these geometrical correlation functions. The key question in this problem-the one that we shall pursue in the remainder of this paper-is therefore what happens in the s-channel of each of these four correlation functions, when two order operators are brought close to each other. Note that, since the conformal field theories we are dealing with are not unitary, the behavior of the G or P functions might be more complicated than in the unitary cases, and involve, in particular, logarithmic terms. Examples of such behaviours are already known for two-and three-point functions [19][20][21].
Remark 2. An important note: unless otherwise specified we will use the same notation (such as P P and G P ) for correlation functions defined on the lattice and for their scaling limits.

Lattice algebras
As mentioned above, our main exploratory tool for unravelling the structure of four-point functions is to impose a lattice discretisation and study the Potts model in the cylinder geometry. The algebraic object that propagates a row of L Potts spins axially along the cylinder axis is a linear operator called the rowto-row transfer matrix T . In this section we discuss how T can be used to build the partition function Z, and defer the more technical question about how to build the correlation functions P P to Appendix A.1. Both the algebraic definition of T and the space of states on which it acts depend subtly on the degrees of freedom defining the model. We are here interested in two different representations, viz. in terms of FK clusters and TL spins, which we now describe in turn. The key technical point is to impose a weight Q per cluster in the former case, or a weight n per loop in the latter.

FK clusters and the join-detach algebra
To build a row of an axially oriented square lattice, the transfer matrix T must first add L "horizontal" edges in some row of constant imaginary time t, and then propagate to the next row at time t + 1 by adding L "vertical" edges. It is convenient to introduce more elementary operators that add just a single edge to the lattice. Concretely, H i adds a horizontal edge between sites i and i + 1 (mod L), while V i adds a vertical edge on top of site i. We can thus write The operators H i and V i must ensure the correct building of the sum A⊆E in (2). They can be written where I denotes the identity operator, while J i and D i will be defined shortly. Each expression has two terms depending on whether the given edge e belongs to the subset A or not. In the former case, a weight v is applied. The subtle point is to obtain also the non-local weight of Q per completed cluster. To that end, T acts on states {s 1 , s 2 , . . . , s L } which are set partitions of L points describing how the sites of a row are interconnected via the parts of the FK clusters living at times prior to t. The join operator J i amalgamates the blocks of the partition corresponding to sites i and i + 1 (mod L). The detach operator D i transforms site i into a singleton, applying a weight Q if it was already a singleton beforehand. The join-detach algebra is defined by the algebraic rules emanating from these requirements: where all indices are considered modulo L. Operators associated with sites that are farther apart than in the relations given simply commute.
In two dimensions, the join-detach algebra is closely related to the Temperley-Lieb (TL) algebra [10] that we describe next. For a more general discussion, see [22]. The question of how the join-detach algebra must be adapted to accommodate the computation of correlation functions is deferred to Appendices A.1-A.2. For some applications (see Appendix A.2 in particular) we shall also need to consider the transpose of the join-detach algebra, which furnishes another geometrical representation of the TL algebra that we shall call the split-attach algebra and describe in some detail in Appendix A.8.

Loops and the Temperley-Lieb algebra
Another option is to define the transfer matrix in terms of the loops that separate the FK clusters from their duals. We recall that these loops now live on a tilted square lattice M(G). At each vertex of M(G) two pieces of loop, labelled i and i + 1 according to their horizontal position, can either bounce off a "vertical" or a "horizontal" edge of G (or its dual G * ), operations that are described in imaginary time by respectively the identity operator I and the so-called braid monoid e i : The e i generate the Temperley-Lieb (TL) algebra which has a long history [10] and is deeply associated with work on the Potts model [23,24]. We note that a horizontal cut through M(G), in between two rows of vertices, will intersect the loop pieces in N = 2L points. If we set J i = Q −1/2 e 2i and D i = Q 1/2 e 2i−1 for any i = 1, 2, . . . , L, the algebraic relations (12) become simply where n is given by (4). These are precisely the defining relation of the TL algebra.
Up to this point we have deliberately been rather loose about specifying the boundary conditions. Indeed, the TL algebra per se is associated with the Potts model on a strip-i.e., with open boundary conditions (that is, free boundary conditions for the Potts spins and reflecting boundary conditions for the loops)-and the generators e i are defined for i = 1, 2, . . . , N − 1. In the cylinder geometry-i.e., with periodic boundary conditions-a tempting possibility is to merely add a last generator "closing" the system, e N , and define the labels modulo N in the defining relations (14), so that in particular e N e 1 e N = e 1 and e 1 e N e 1 = e 1 . This natural generalisation however takes one into a sticky mathematical problem: the corresponding algebra is then seen to be infinite-dimensional, even for finite N . In a nutshell, this occurs because of through-lines or loops that can wind around the system. While what must be done with these objects is clear in the Potts model itself, this requires providing extra information that is not present in the definition of the "periodicised" Temperley-Lieb algebra. This extra information takes the mathematical form of quotients.
To define these quotients more precisely, it is useful to also introduce a translation generator u that shifts the label of the e i generators, giving rise to the following extra relations-in addition to (14a)-(14c)-with integer indices considered modulo N (that is, i ∈ Z N ): The translation operator has the diagrammatic representation The last relation (14e) is easily understood in terms of diagrams, for example for N = 4, Note also that u N is central: it commutes with all the generators e i . The resulting algebra is called the affine Temperley-Lieb algebra T a N (n). In the following we shall draw extensively on known results about its representation theory [25,26] and its relation with conformal field theory [27].

The transfer matrix sectors
While T a N (n) is infinite-dimensional, it is easy to define the finite-dimensional modules which are relevant to us. First, we fix the number of through-lines, which are the pieces of loops connecting the bottom and the top of the diagrams. The number of through-lines is denoted 2j, with j ∈ N/2-the factor of 2 comes about because one can relate each loop strand to a q-deformed representation of spin-1/2. Second, we stipulate that whenever 2j through-lines wind counterclockwise around the axis of the cylinder l times, we can unwind them at the price of a complex phase factor e 2ijlK ; similarly, for clockwise winding, the phase will be e −i2jlK [25]. This unwinding means more precisely that we equate the words in the algebra corresponding to the winding configurations with a numerical factor (the phase) times the related words without winding. This operation is known to give rise to a generically irreducible module over T a N (n), which we denote by W j,z 2 =e 2iK and call the standard module [26]. A key point is that inside the modules W j,z 2 one has the identity meaning that the central element u N of T a N (n) is replaced by the complex number z 2j . The dimensions of the standard modules W j,e 2iK are given bŷ Note that the dimensions do not depend on K, although the representations with different e iK are not isomorphic. The case j = 0 is a bit special, due to the absence of through-lines. There is no pseudomomentum, but representations are characterised by another parameter, related with the weight given to non-contractible loops (i.e., loops that close around the periodic direction). Parameterising this weight as z + z −1 , the corresponding standard module of T a N (n) is denoted W 0,z 2 and has dimension N N/2 . These modules are irreducible for generic z. As in the case j > 0, we indicate only the z 2 value, though it does not mean that the two standard modules with ±z are isomorphic. We will indicate the sign of z when it is necessary.

Potts model
The fact that we wish to apply T a N (n) to study the Potts models entails a few minor modifications of the general setup. First, since the number of sites N = 2L is even, the number of through-lines is also even, so that j ∈ N. Second, the translation operator in the Potts model shifts the L spins cyclically, meaning that the TL sites must be shifted by two units. Therefore, we are actually going to use the subalgebra in T a N generated by the e i 's and by u 2 (instead of u itself), which is why the above notation refers to z 2 rather than to z itself.
The basic ingredient is the finite-dimensional modules of T a N described above, with the loop weight n = √ Q parameterised as in (4) in order to match (3). However, to account for the particularities of the Potts model, the algebra that we are mostly interested in is a quotient of T a N , known as the Jones-Temperley-Lieb algebra JTL N (n) [29,30]. It is obtained by: (i) replacing non-contractible loops by the same weight n = √ Q as for the contractible ones; (ii) identifying diagrams connecting the same sites, even if they are non-isotopic on the cylinder; and (iii) setting u N = 1, which allows one to unwind through-lines.
Rules (i) and (ii) are only relevant in the case j = 0, where through-lines are not present. The first rule leads to z 2 = q ±2 . In this case, in fact, the affine TL module W 0,q 2 is reducible, and contains a unique simple submodule isomorphic to W 1,1 . The reason for this is that the general T a N (n) allows (in diagrammatic terms) to distinguish loop arcs that connect two given sites on the front or on the back of the cylinder, meaning that a closed loop can be given different weights depending on whether it is contractible or not. When this distinction is not needed we must identify arcs only according to which sites are being connected, as in rule (ii). Identifying non-isotopic diagrams in this way corresponds to replacing W 0,q 2 by the (unique) simple quotient W 0,q 2 /W 1,1 ≡ W 0,q 2 of dimension In technical terms, this quotient is precisely the standard module of JTL N (n) for j = 0.
Remark 3. The quotient W 0,q 2 is but one example of representations that appear more generally in T a N when q is still generic, but z takes particular values [26,28]. Indeed, the standard module W j ,z has a non-zero homomorphism to W j,z , if and only if j − j ∈ N 0 and the pairs (j , z ) and (j, z) satisfy When q is not a root of unity, there is at most one solution to (18). When there is one, the module W j,z is not irreducible, but has a unique proper irreducible submodule isomorphic to W j ,z . One can then obtain a simple module by taking the quotient of dimension The quotient W 0,q 2 appearing above is but the simplest example (with j = 0, z = q 2 and j = 1, z = 1) of this situation, and it is the only such quotient that is relevant for the Potts model at generic q.
Whenver j = 0, rule (iii) leads to a quantisation of the momentum: K = πp/M , with M a divisor of j (i.e., M |j) and with a greatest common divisor p ∧ M = 1. The modules encountered so far are thus W 0,q 2 = W 0,q −2 , and W j,e 2iπp/M for j = 0 and M |j.
On top of this, there is a small subtlety having to do with the relation between through-lines and through-clusters, by which we mean clusters that propagate along the imaginary time direction. For j ≥ 2, each of the 2j through-lines alternatingly separates a propagating FK cluster and an propagating dual cluster, implying the existence of precisely j through-clusters.
However, when we wish to impose just one through-cluster, the situation is different. Since nothing prevents this cluster from wrapping the periodic direction of the cylinder, it will in fact do so with probability one, implying that through-lines will be absent (j = 0). On the other hand, there cannot be any noncontractible loops either, since this would prevent the propagation of the through-cluster. The correct module is thus obtained by giving a vanishing weight to non-contractible loops. This is easily accomplished by setting z = ±i, leading to W 0,−1 .
The three types of modules that we have just introduced: are known to encode the full Potts model partition function on the torus [7,30,31]. Formal multiplicities for these modules are also known. For Q non-integer, they are real numbers with group-theoretical significance [20,21]. The crucial observation that will be made below is that the modules (22) are also the sufficient objects to describe the four-point correlation functions in the geometrical Potts model. An important additional fact is that actually only the modules with even values of j are necessary for the description of four-point functions. By contrast, any j ≥ 2 contributes 2 to the partition function on the torus, as has been verified in details for finite systems [32].
This last result was not totally obvious a priori. Indeed, one must in general be careful with geometrical questions in models such as the Potts model, where the set of observables is seemingly not limited, if one sways far enough from locality. It is well-known, for instance, that correlations involving several independent paths along clusters-the case of two such paths defines the celebrated backbone exponentscannot be described using T a N , and the corresponding exponents have never been identified using Coulomb gas techniques. Indeed, the numerical measurements of [33][34][35] appear to convincingly rule out any such identification for this whole class of so-called monochromatic path-crossing exponents. Similar remarks can be made about other seemingly reasonable observables, such as the shortest-path exponent [36], to mention but one example. Fortunately, then, the matters seem to be (relatively) simpler for the four-point correlation functions.

Summary of notations
For the reader's convenience we summarise here some of the notations used in this paper. They are, as far as possible, the same as those used in [27,28,[37][38][39][40][41][42].
• W j,z 2 =e 2iK -the standard modules over T a N , • W j,z 2 -the same, with P = e 2iK , • W 0,q 2 -the standard module over JTL N for j = 0, • j, e 2iK or X j,z 2 -simple modules over T a N (n).
Moreover, when discussing the transfer matrices for the Potts model correlation functions (see Appendix A, and section A.4 in particular) we shall sometimes need a lighter notation V ,k,m for the sector with propagating clusters, an integer momentum variable k = 0, 1, . . . , − 1 for the through-lines (if any), and a lattice momentum m = 0, 1, . . . , L − 1 which is the precursor of the conformal spin for a system of finite size . The notations V ,k or V with some of the variables omitted mean that these take indiscriminate values. The correspondence with the standard modules is then: • V 0 is the sector with no through-lines, and non-contractible loops have weight • V 1 is the sector with no through-lines, and non-contractible loops have weight zero: • V ,k is the sector with j = ≥ 2 pairs of through-lines and phases z 2 = e 2iπk/j : V ,k = W j,z 2 =e 2iπk/j .
4 Four-point functions in the s-channel

Generalities
We consider a general four-point function of primary operators in a CFT, which we write in the following convenient form in the plane: where we have denoted z ij ≡ z i − z j , with the exponents The antiholomorphic exponentsδ ij are obtained from the holomorphic ones δ ij by the replacement h →h, and the same convention henceforth applies to any other quantity. The parameter z denotes the anharmonic ratio z ≡ z 12 z 34 z 13 z 24 .
One finds easily that This function G(z,z) is what one usually refers to as Φ 1 (z,z)Φ 2 (0, 0)Φ 3 (∞, ∞)Φ 4 (1, 1) . It is known [43] to expand as a sum over conformal blocks where (∆,∆) are the conformal weights of the primary fields appearing in the operator product expansion relevant at small z. They define the scaling dimension (eigenvalue of the dilatation operator) ∆ +∆ and the conformal spin (eigenvalue of the rotation operator) ∆ −∆.
We shall be particularly interested in the limit z 1 → z 2 and z 3 → z 4 : this is called the s-channel (borrowing a standard terminology of particle physics due to Mandelstam). This limit corresponds to taking z → 0 in (26), so we can write the expansion One could of course similarly consider the t-channel (z → ∞) and in the u-channel (z → 1), by expanding in powers of 1 z and (z − 1) respectively. The idea of the conformal bootstrap programme is that all these expansions determine the same function G(z,z) and hence will impose constraints. The first step in any further discussion is therefore to establish the fields intervening in one of these channels, which we take here to be the s-channel.
The key question we want to address in this paper is thus to determine the set S of values of ∆,∆, which we will tackle by a combination of algebraic and numerical methods. In its crudest form, the latter involves the brute force numerical determination of a (very) large number of terms appearing in the right-hand side of (28).
Note that the determination of the set S from the knowledge of these terms will only be fully possible in "generic" cases, where none of the ∆,∆ differ by integers. Otherwise, there will be ambiguities, as a term such as ∆ + n,∆ +n (with n,n integer) may arise from a genuine primary field, or from a Virasoro descendent of some primary field with weights ∆ + p,∆ +p, with p ≤ n andp ≤n (with at least one of the inequalities being strict). The non-generic case requires Q to take particular values (with q being a root of unity); it is clearly more complicated than the generic case and will typically lead to at least some correlation functions having logarithmic behaviour. A few non-generic cases (not all of them logarithmic) will be discussed in Appendix B. But the main text is henceforth dedicated to the generic case, for which we shall determine S fully.
Our strategy is to study the expansion (28) on the cylinder, where we will be able to use, on the numerical side, transfer matrix techniques, and, on the analytic side, algebraic results. The four-point function on the cylinder follows from (23) via the conformal map Using the fact that the fields are primary, and restricting here to i = j = k = l for notational simplicity, we find where the subscript "cyl" on the left-hand side refers to the cylinder geometry, and we have set w ij ≡ w i −w j as before. The expansion variable is now Using (28) we can write this as In practice, to access the s-channel properties, we will take the points w 1 , w 2 on a given slice of imaginary time, and w 3 , w 4 on another, distant, slice along the cylinder of finite circumference L. This geometry is shown in Figure 1. In other words, w 12 and w 34 will be fixed, while w 13 and w 24 will be large and vary. In this limit, it will then be possible to compare the expansion (32) with the results of transfer matrix calculations, and identify, in particular, the set S.
Let us now be more precise. We set which means that the points w 1,2 and w 3,4 are a certain distance 3 2a apart on the vertical axis, l is the horizontal distance (imaginary time) between the two groups, and on top of this we have the centre of mass of w 3,4 shifted by x. A short calculation then gives We can then expand this from (28): where we have set The bracket [1 + O(ξ,ξ)] contains now contributions from the conformal blocks and contributions from the hyperbolic functions in the conformal map. The expansion (35) is the crucial tool that we will use systematically in our analysis below. In the following, we will sometimes use the short-hand notation We now discuss this in more detail.

Remark 4.
We also see that if we exchange w 1 and w 2 in (32), the leading contributions for a given ∆,∆ is multiplied by (−1) ∆−∆ . Hence primary fields with odd integer spin should contribute an opposite weight upon making this exchange.
For future reference, the definition of the channels is Henceforth, when denoting the probabilities P a1,a2,a3,a4 , the four labels specifying the partition P = {a 1 , a 2 , a 3 , a 4 } refer to the points z 1 , z 2 , z 3 , z 4 in that order. Clearly, then, P aaaa should have the same spectrum (and structure constants) in all channels, while, for instance, P aabb should have the same spectrum (and structure constants) in the t-and u-channels, while the spectrum should be different in the s-channel.

Exponents
Contrarily to what is implied in [1], the exponents of percolation-and more generally the Q-state Potts model in the FK cluster representation-are essentially known (with the exception of certain "exotic" exponents, see [33][34][35][36]). This knowledge relies on two stages. First, the transfer matrix sectors of the Potts model can be described in terms of standard modules of the affine TL algebra, as described in section 3.4. Second, the continuum limits of these objects are known in the form of spectrum generating functions within the corresponding CFT, as we now review. This is of course not yet the solution of the s-channel conundrum, but since we are able to formulate the four-point functions in terms of the FK transfer matrix (see section 4.3 and Appendix A.1), the results on the generating functions will narrow down the set of states than can possibly be part of the spectrum S. Extensive numerical analysis-corroborated by the solvability of a few special cases (see Appendix B)-will then lead to the results that we give in section 5.
The local FK connectivities in the geometrical Potts model and their evolution along the cylinder are described by a transfer matrix or, in the familiar extreme anisotropic limit, a Hamiltonian. Both transfer matrix and Hamiltonian exhibit the same conformal content-that is, eigenstates associated with local CFT operators, and the corresponding conformal weights h,h, together with their multiplicities. It is convenient to encode the latter into spectrum generating functions. Using for instance the Hamiltonian language and setting with λ adjusted so that the sound velocity is unity as usual, we define the generating function of levels (eigenenergies of H) and lattice momentum P as as traces of lattice operators, with the scaling limit [44] Tr Here ε 0 is the (non-universal) ground state energy per site in the limit N → ∞. The scaling limit is defined by taking N, β R , β I → ∞ while keeping the modular parameters 4 q(q) = exp − 2π N (β R ± iβ I ) (with β R,I real and β R > 0) finite. The parameters β R and β I define the size of the system in the two principal directions of the torus, while the trace ensures the periodic boundary conditions in the imaginary time direction. We recall that N = 2L is the number of sites in the system, often referred to as the length of the spin chain in the Hamiltonian limit. In other words, only even chains are relevant in our problem. On the right-hand side of (40), L 0 andL 0 are of course Virasoro generators, while c denotes the central charge.
The generating function (40) calculated in the modules W j,e 2iK is [44,45] 4 Here and elsewhere a notation of the type q(q) means that q is given by the first expression on the right-hand side (the one with a + sign), andq by the second expression (with a − sign). where is the (inverse of) the generating function for integer partitions, and η(q) is Dedekind's eta function. Instead of (4) we shall find it convenient to parameterise the number of states in the Potts model by so that q = e iπ m+1 . Note that to access the generic case (q not a root of unity) we do not restrict m to be integer, as would be the case for the minimal models. The corresponding central charge is then and we also use the Kac table parameterisation of conformal weights In "usual" CFT the labels (r, s) are positive integers, but as for the parameter m we shall here allow them to take more general values, as is already evident from (41). To make contact with standard references, it is also convenient in the following to introduce the Coulomb gas coupling constant g = m m+1 and the background electric charge e 0 = 1 m+1 . The operator associated with the order parameter has conformal weight h 1/2,0 [15,46], the primordial example of an "unusual" h rs , with r here being non-integer and s non-positive. This operator belongs to the generating function F 0,−1 .

Remark 5.
To compare with [1] one must set in their equation (1.1) β 2 = m m+1 (so that 1 2 ≤ β 2 ≤ 1), and q = Q. Moreover, the conventions used in their paper for the exponents are switched with respect to ours. In other words, they call ∆ sr what we call h rs (or ∆ rs ).
Restricting now to the cases of interest with momentum K = πp/M and M |j gives On top of this we also have to consider the generating function of levels in W 0,−1 , which reads Finally we need the generating function for the quotient module W 0,q 2 . The twist e iK = q corresponds in our notation to K π = e 0 , so we have first The subtraction necessary to obtain the module W 0,q 2 = W 0,q 2 /W 1,1 leads to the expression for the generating function of the corresponding levels: The subtraction can actually be implemented term by term. Introducing the characters of the so-called Kac modules-which are Verma modules where a single singular vector at level rs has been removed-, we haveF To summarise, corresponding to the modules (22) we expect (and confirm below) that the set of exponents contributing to the four-point connectivities is encoded into Moreover, we shall see below that only j even contributes, which cannot be foreseen at this stage. We note that for generic values of Q or q (i.e, with m irrational), there are no coincidences of exponents in the different sectors (generating functions)-this is also true on the lattice, where there are no coincidences of eigenvalues. Moreover, in a given sector, no two exponents differ by integers.

Numerical validation of the generating functions
It appears useful at this stage to test the internal coherence of the ingredients brought together this far. On one hand, in section 3.4 we have related the sectors of the Potts model transfer matrix T to certain standard modules, W j,z 2 and W 0,q 2 , of the affine TL algebra. On the other hand, we have just given their corresponding spectrum generating functions, F j,z andF 0,q 2 . This means that a numerical diagonalisation of T in the various sectors should produce-after a proper extrapolation to the continuum limit L → ∞the primaries and descendents (with multiplicities) of these generating functions. We are not aware of a previous careful study that this is indeed so.
To this end, we first outline in Appendix A.3 the extraction of the eigenvalue spectrum of T in the various sectors relevant for the Potts model. Fixing the values of the momentum and the conformal spin is a non-trivial operation that is expounded in Appendix A.4. Once this has been done, we can examine the spectrum of T ; this is done first for a generic value of Q in Appendix A.5.1, and then for a few non-generic values: Q = 4 in Appendix A.5.2, and Q = 2 in Appendix A.6. When combined, these three cases permit us to test examples of Verma modules with zero, one, or infinitely many singular vectors.
In all cases we find that the agreement with the expected spectrum generating functions is perfect. In the generic case, we are able to see descendents up to level 6 for the identity operator, and up to level 3 for other operators. Moreover, the set of primaries fully agree with the expectations from the affine TL algebra. In the Ising case we are able to follow the first 29 scaling levels and observe descendents up to level 9 for both operators (I and ) in the even sector, with an agreement better than 10 −4 for almost all scaling dimensions. Moreover, the degeneracy observed for each "completed" level is in perfect accord with the spectrum generating functions.
The techniques used in the numerical analysis may be of independent interest and can be consulted in the appendices (see also Appendix A.7 for a few practical remarks).

The numerical algorithm
The geometrical setup for four-point functions is shown in Figure 1. As stated earlier, our lattice discretisation consists in embedding a periodic square lattice G = (V, E) of width L Potts spins into the cylinder, with the edges E being either horizontal or vertical with respect to the figure (axial geometry). We possess two different strategies for obtaining numerical results for the correlation functions.
The first strategy gives access to the most general FK correlation functions, namely the 15 different P a1,a2,a3,a4 . It is practically feasible up to size L = 7, after a considerable numerical effort. It applies for both generic and non-generic values of Q, meaning that in the latter case it can determine the indecomposable structure of correlation functions.
The second strategy applies to a smaller set of correlation functions, namely the four order-parameter correlators G aaaa , G aabb , G abba and G abab . Its advantage is that it gives access to larger sizes, in practice up to L = 11, at a much smaller computational expenditure than the first method. However, it applies only to generic values of Q and, at least in its present form, cannot determine the Jordan block structure at non-generic Q-values.
We now briefly outline the two methods, while relegating all technical details to Appendices A.1-A.2.

First method
It is shown in the appendix that all P a1,a2,a3,a4 can be computed, for fixed values of the distances a, x and l, via a suitable modification of the FK representation of the transfer matrix in which certain clusters (viz., the ones touching one or more of the points w 1 , w 2 , w 3 , w 4 ) carry specific marks. The spectrum of this modified transfer matrix is contained within that of the original one, namely the one described in section 3.4 in terms of affine TL representations. The spectrum can be proven to be real, so we can order the distinct eigenvalues as Λ 0 > Λ 1 > · · · > Λ i > · · · . The correlation function then takes the following form, for generic values of Q, where the amplitudes A i = A i (a, x, L) are to be determined. The corresponding expression for non-generic values has the same form, but with the replacement whenever the eigenvalue Λ i is associated with a Jordan block of rank r i . In the latter case the generalised amplitudes A i (a, x, L) are again independent of l.

Remark 6.
There is of course an exact degenerescence of the scaling dimension of a CFT operator with non-zero spin ∆ −∆ and that of its conjugate (i.e., obtained by the exchange ξ →ξ). This is prefigured in the lattice discretisation by the exact degenerescence of the eigenvalues corresponding to eigenstates of T with opposite non-zero lattice momenta ±m (cf. Appendix A.4). Because of the regrouping of degenerate eigenvalues in (53) it should thus be remembered to divide the amplitude A i of such states by a factor of two when comparing to the CFT predictions (see Appendix B for many examples of this phenomenon).
For finite L the set of eigenvalues of T -and hence the number N ≡ i r i of (generalised) amplitudes to be determined-is finite. It follows that the form (53) is an exact expression, not merely an asymptotic expansion. Therefore, to determine the amplitudes A i -or the generalised amplitudes A (j) i for the cases with Jordan blocks-for given separations a, x and size L, it suffices in principle to numerically determine the spectrum {Λ i }, compute the correlation functions P a1,a2,a3,a4 for N different values l, and to invert the linear system (53).
In practice, of course, things are more complicated, and several remarks must be made (see Appendix A.7). The most important of those is that the magnitude of the terms in (53) decreases exponentially fast, in particular when N , and hence l, is large. Therefore both {Λ i } and P a1,a2,a3,a4 must be computed to an exceedingly high numerical precision. For instance, our most demanding computation (see section 5.2.2 for details) required a 4000-digit numerical precision.
Another remark is that when N 1 we often wish to determine only the first "few" amplitudes (corresponding to i ≤ some i max ). This can be done by using the expression (53), truncated to the first i max terms, as an asymptotic expression, i.e., by solving for A i the linear system provided by the numerically computed P a1,a2,a3,a4 with l = l min , l min + 1, . . . , l min + i max , where l min is taken sufficiently large. One then has to carefully check that the desired A i are stable, within the desired numerical precision, to small changes in l min .

Second method
Our other method applies to the computation of the order-parameter correlators G a1,a2,a3,a4 , defined in (7). This requires another variant of the FK transfer matrix with marked clusters, as described in details in Appendix A.2. The number of different marks allowed must be chosen as the number of different symbols among a 1 , a 2 , a 3 , a 4 , and the dimension of the transfer matrix grows with this number. In practice we have employed two different marks, in order to gain access to the four correlators G aaaa , G aabb , G abba and G abab .
The spin operator O a (σ k ), defined in (5), can be expressed within this basis and has essentially the effect of attributing the label a to the spin situated at vertex k ∈ V . Our geometrical setup is such that vertices w 1 and w 2 belong to the same time slice, while w 3 and w 4 belong to another time slice, with the relative positions within these two time slices being specified by Figure 1. We henceforth denote the spin operators simply by O a k , for k = 1, 2, 3, 4, and keep implicit their point of insertion in the relevant time slices.
Let v i | and |v i denote the left and right eigenvectors of T . In the case of simple eigenvalues we then have But even for the generic values of Q that we consider here, some of the eigenvalues are degenerate, due to symmetries of the lattice and of the order parameter symbols a k . In that case we denote the multiplicity of Λ i by d i , and we endow the corresponding eigenvectors with an extra label, |v i,j , where j = 1, 2, . . . , d i . The left and right eigenvectors can be obtained efficiently within an iterative diagonalisation scheme, such as the Arnoldi method (see again Appendix A.7 for details). Left and right vectors will obviously be orthogonal if they correspond to different eigenvalues, but the Arnoldi method does not guarantee orthogonality within the degenerate subspaces. It is however possible to perform an additional diagonalisation step that will ensure that the orthogonality holds with respect to both labels: Remark 7. Some readers are likely to be well acquainted with the representation theory of the TL algebra, in which "geometrical" scalar products (see e.g. [47,48]) are introduced between basis states (often called link patterns in the context of the loop representation). These scalar products "count" loops and clusters formed by the gluing of the states, leading to Q-dependent results. We must therefore stress that the scalar products appearing in this section are simply the standard Euclidean scalar products between ordinary vectors. Moreover, all eigenvectors turn out to be real, so there is no issue of complex conjugation.
We now claim that the amplitudes A i for the generic case without Jordan blocks can be obtained as It is crucial for the validity of this result that the orthogonalisation in degenerate subspaces has been performed (see Appendix A.2). We have performed extensive checks that the first and second methods give the same results, in situations where they are both applicable. The advantage of the second method is that it is numerically much more efficient, and hence gives access to larger sizes L. The sources for this gain in efficiency are explained in Appendix A.2.
Formula (57) has a nice geometrical interpretation that makes direct contact with Figure 1. Indeed the numerator of the formula (read from right to left) and the figure (read from left to right) are completely analogous: 1. The propagation from the free boundary condition at imaginary time t → −∞ to the time slice containing w 1 and w 2 corresponds to the production of the ground state |v 0 .
2. The first two operators are then inserted by O a1 O a2 .
3. The piece |v i,j v i,j | corresponds to the projection on a definite state appearing in the s-channel of the four-point function.
4. This is followed by the insertion of the two remaining operators, at w 3 and w 4 .
5. The projection on v 0 | matches the propagation to the free boundary condition at the other extreme of the cylinder (t → +∞).
We stress that the validity of (57) depends crucially on the orthogonality (56).

Continuum limit
To extract the continuum limit (L → ∞) of the amplitudes it is important to be able to associate each A i with a definite field in the continuum limit. 5 For instance, one question that one might want to answer is what would be the amplitude contributing to a given primary Φ ∆∆ , i.e., to identify the A i that will converge to the amplitude of Φ ∆∆ in the expansion (35). Several remarks are in order in this respect. Obviously, it is the ratio between two amplitudes-rather than each amplitude taken individually-that is universal and hence related to CFT. Therefore we assume tacitly in what follows that each amplitude of interest is measured via its ratio to the one that gives the leading contribution to the considered correlation function. Moreover, the conformal mapping to the cylinder implies that we should correct the raw lattice result by a conformal factor, namely the powers of sin 2πa L appearing in (35). Once this has been done, our main claim is that we are capable of analysing the numerical results so as to establish the convergence where the conformal amplitude has been defined in (37). More generally, we can obtain the corresponding results for subdominant contributions from the conformal blocks, corresponding to the amplitude multiplying a term of the type ξ NξN in the square bracket of (35). This is interpreted as the (total) amplitude of the descendents at level N ,N of the primary Φ ∆∆ . The challenge involved in making this identification is to make sure that we possess enough information about the lattice model to unambiguously associate a given field in the continuum limit with its "corresponding" eigenvalue of the transfer matrix in finite size L. In the lattice model, each A i is unambiguously associated with the eigenvalue Λ i . A careful study of the transfer matrix (see Appendix A.4) enables us to attribute to each eigenvalue three labels , k, m, formally restricting to a representation denoted V ,k,m . The first label gives the number of propagating FK clusters, so in the notation of the standard modules W j,z 2 we have j = for all = 1, while = 1 corresponds to j = 0. The second label k is directly related with the momentum of through-lines, via z 2 = e 2iπk/j . And finally the third label m is the lattice momentum that gives directly the conformal spin, h −h = m, at least if L is large enough to accommodate the desired spin.
While certainly very helpful, the three labels , k, m are not quite enough to determine which conformal block to associate with A i , nor at which level N ,N . Roughly speaking, the trouble is that the k'th smallest scaling dimension in the continuum limit does not necessarily correspond to the k'th largest eigenvalue of T in finite size L. While this is certainly true for L 1, there are numerous crossovers in finite size, and these have to be monitored carefully in order to make to correct identifications. How we overcame this delicate problem is described in Appendix A.5.
Finally, once the finite-size approximation A i to a given CFT amplitude has been determined, for several different sizes L, the numerical value of the latter is determined by finite-size extrapolation techniques. This is again discussed in Appendix A.5. A large number of concrete applications of the entire method can be found in section 5 and (with more details provided) in Appendix B.

Checks
Our approach, being based on properties of the lattice model, requires a careful control of the continuum limit. There are several aspects to this. The most obvious one is that, since we are studying four-point functions of a CFT, we should, ideally, have all ratios |z ij | 1 (where all distances are measured in units of the lattice spacing). On the cylinder, we have chosen to take points far apart along the cylinder axis, but placed, pairwise, on identical imaginary time slices. Since the cylinder widths L are limited for technical reasons, this means that |w 12 |, |w 34 | will be limited, in fact, to a few lattice spacings. 6 We first observe that the dependence A(x) of the amplitudes on the shift x in the space-like direction (see Figure 1) between the two groups of points is in fact trivial. Taking into account that all amplitudes have been normalised as ratios with respect to the leading one, as well as remark 6, it is seen from (35) that where s = ∆ −∆ denotes the conformal spin (which coincides with the lattice momentum, m = s).
An alternative means of deriving (59) goes through the inspection of (57). Imagine evaluating the first scalar product in the numerator in a geometry where the cylinder has been rotated by the amount −x. This rotation will re-align the second pair of operators O a3 O a4 with the first pair O a1 O a2 , like in the computation of A(0). The ground state v 0 | is obviously rotationally invariant, but an intermediate state |v i,j of lattice momentum m = 0 is not, and will therefore pick up a corresponding phase factor under the rotation. Summing this over the degenerate contributions ±m reproduces (59).
We have checked numerically on explicit examples that (59) holds true exactly in finite size. As a matter of fact, determining numerically the dependence A(x) is a convenient means of establishing the lattice momentum m of a given state, complementary to the techniques explained in Section A.4.
A maybe more subtle aspect is that the lattice observables are not in general pure scaling fields. This means that the conformal field whose four-point functions we want to study, is identified on the lattice as the Potts spin operator only up to additional corrections ("higher (or excited) spin operators"), whose contributions become negligible only when all distances are once again much larger than the lattice spacing: put otherwise, measured four-point functions on the lattice are a mixture of four-point functions of pure scaling fields.
For our approach to be useful, it is necessary to perform many tests in order to control these potential drawbacks. As discussed extensively in Appendix B, we have checked that, for the sizes we were able to access: • The mixture of excited spin operators can be neglected (see Appendix B.1); and • The values of the extrapolated amplitudes A Φ ∆∆ -see (37) Moreover, even for operators higher in the spectrum S where the lattice determination of amplitudes may not have fully converged, our approach, combined with the algebraic understanding of transfer matrix sectors, indicates unambiguously which coupling constants will remain non-zero in the scaling limit, even if error bars on their extrapolated values are not negligible. Figure 2: FK cluster configurations that contribute to the correlation functions P abab and P abba .

5.2
The s-channel spectrum of P abab or P abba : F j,e 2iπp/M (M |j, j≥2), j even We first study the cases where (at least) two distinct clusters are forced to propagate between the two distinguished time slices in Figure 1. Specifically, P abab is the probability that points w 1 , w 3 and w 2 , w 4 respectively belong to the same clusters. This quantity is also called P 2 in [1]. Similarly, P abba is the probability that w 1 , w 4 and w 2 , w 3 belong to the same clusters. These two correlation functions are depicted in Figure 2.

Results in finite size
It is evident from Figure 2 that the leading term in this sector should be the term corresponding to the propagation of two different clusters, that is, four cluster boundaries. The affine TL modules W j,z 2 (or their continuum counterparts F j,z 2 ) correspond to 2j through-lines, so the propagation of two clusters must have a contribution with j = 2. The corresponding generating function of levels has two sectors, depending on whether a pair of boundaries going around the system picks up a phase z 2 = 1 or z 2 = −1. No other choice is possible since for two pairs of boundaries (picking up a phase z 4 ), we do not want a phase. 7 Hence we expect the contribution of modules We have first checked that all the eigenvalues associated with these two modules contribute to the probabilities P abab and P abba for all finite sizes. This is however not all. Sectors with a higher number of clusters than the two imposed by the choice of indices might also be thought to contribute to these correlation functions. One could think of several mechanisms for such contributions. First, there might be more clusters, distinct from the two imposed by the boundary conditions, that "by chance" connect the two time slices. Second, the two clusters might have more complicated topologies, with for instance one of them (say, the one containing points w 1 and w 3 in the left part of Figure 2) starting out at one insertion point (here w 1 ), and wrapping all around the other cluster (containing points w 4 and w 2 ), before arriving at its terminal point (here w 3 ).
We have found it difficult to provide a convincing argument that certain subclasses of configurations will necessarily lead to further contributions to the correlation functions, in terms of the modules W j,z 2 ; we think there are underlying symmetry and branching rules considerations that may answer this riddle on general grounds, and that we do not yet control. Fortunately the numerical results are completely clear. We find that, for all finite sizes, all eigenvalues associated with the modules W j,e 2iπp/M , M |j, with j even, and only those, contribute to the probabilities P abab and P abba . For j = 4 for instance, this allows contributions from the momentum sectors z 2 = 1, exp(iπ/4), exp(iπ/2), exp(3iπ/4), and thus the following modules, in addition to those of (60), Note that for a given width L, the maximum value of j is bounded from above, j ≤ L. We have checked that, as L increases, higher values of j start contributing to the probabilities, provided that the separation 2a between the insertion points is sufficiently large. More precisely, we have observed that: • For L = 5 and separation 2a = 2, all the eigenvalues with j = 2, 4 and none of the eigenvalues with j = 0, 1, 3 contribute to the probabilities.
• Still for L = 7, but diminishing to separation 2a = 2, the contributions from j = 6 disappear. Diminishing further to 2a = 1, the contributions from j = 4 disappear as well.
The above result is corroborated by a closer study of the spectrum of the transfer matrix described in section 4.3.2, namely the one that produces the correlation functions of order parameter operator. Its eigenvalues are precisely those corresponding to the modules (22) with j even, while those corresponding to j odd are not observed at all. Motivated by the above list of observations, we conjecture that in fact a given sector j ∈ 2N contributes only when 2a ≥ j/2.

Non-zero coupling to the sector j = 6
The computation with L = 7 and 2a = 3 establishing that the sector j = 6 does contribute to the correlation functions P P is the most demanding among all of those made for this paper. For the benefit of readers interested in computational aspects (and those wanting to scrupulously assess the validity of our conclusions) we wish to describe it in some more detail-other readers may wish to skip this section and resume the reading below.
This computation was performed for the generic value Q = 3 2 . There is a total of 3 932 distinct eigenvalues in the sectors with j = 0, 2, 4, 6, corresponding to all possible momenta. Using the methods of Appendix A.4 we can classify them in sectors V ,k,m corresponding to propagating clusters with cluster momentum k and lattice momentum m. Ordering all the eigenvalues as Λ 1 > Λ 2 > · · · > Λ 3932 , with Λ 1 being the dominant eigenvalue in V 0 (i.e., the ground state), the dominant eigenvalues in the sectors V 1 , V 2 , V 4 and V 6 -which obviously have vanishing momenta, k = m = 0-are respectively Λ 2 , Λ 5 , Λ 205 and Λ 2390 .
Suppose first that we considered some correlation function coupling to all of these eigenvalues, and we wished to isolate the amplitude corresponding to Λ 2390 by using the first method of Appendix A.1. The ratio r = Λ 2390 /Λ 1 3.451 · 10 −4 , and since we need to determine 2 390 coefficients A i in (53) we will need the same number of equations, obtained by choosing the distance l = l min + 1, . . . , l min + 2390. We would need (at the very least) to take l min = 100 in order to be in the asymptotic regime. Then, since r 2490 2.6 · 10 −8821 , we see that the terms entering (53) would differ by almost nine thousand orders of magnitude, so allowing some margin for numerical instabilities we would have to compute the correlation function for (at least) 2 500 different values of l to a numerical precision of (at least) 10 000 digits. This task is hopelessly impossible, given that the transfer matrix of Appendix A.1 is of dimension ∼ 10 6 in this case.
To do better, we need to consider a particular well-chosen combination of correlation functions, designed so that it decouples from a sufficient number of low-lying states in the spectrum. The quantity is a such a good combination. On symmetry grounds, it decouples from the sectors with odd momenta k. The term P aaaa is rather easily checked to pick up contributions from the sectors V 1 , V 2 and V 4 (and maybe higher values of ), whereas P aabb couples in addition to V 0 . Meanwhile, P abab and P abba get contributions only from V 2 and V 4 (and maybe higher values of ); it is indeed clear that since these terms impose two long clusters (see Figure 2) they cannot couple to V 0 and V 1 . The surprising property of P * is now that, with the above choice of the two coefficients in its definition, all contributions from V 1 and V 2 disappear from the combination. In other words, P * couples to V 0 , V 40 , V 42 , and maybe V ,k with higher values of and even k. We have dim V 0 = 232, dim V 40 = 190 and dim V 42 = 182. Therefore we shall be able to settle whether there is a non-zero amplitude for the 6-cluster sector provided we can look beyond the first 604 eigenvalues.
To that end, we have computed P * for l = 100, 101, . . . , 900 to be on the safe side. Noting that r 900 1.3 · 10 −3116 we have performed the computations to a numerical precision of 4 000 digits. This required about 3 × 10 5 hours of single-processor CPU time. The conclusion is that we have unambiguously established that P * picks up non-zero contributions from the first few eigenvalues in each of the sectors V 60 , V 62 and V 64 , with amplitudes in the range ∼ 10 −10 . While these numbers may seem small, they follow the clear trend (observed in all cases) that the non-zero amplitudes decay exponentially with the index of the corresponding eigenvalue. Moreover, these amplitudes are numerically stable towards changing l min throughout the range l min ∈ [100, 200]. We have also checked that the absolute values of contributions which are genuinely supposed to be zero (such as those from sectors V 61 , V 63 and V 65 ) come out numerically as 10 −500 , which is fully compatible with the above estimates of the required numerical precision.

Exponents
Introducing our usual notation F j,z 2 , the spectrum in the sector with j = 2 propagating clusters is encoded in the generating functions This corresponds to the conformal weights given by (46): Note that for the first part of the spectrum, h −h is an even integer, while it is an odd integer for the first part. We will denote these two contributions by 2S and 2A respectively, where S stands for symmetric and A for antisymmetric. Going back to an earlier remark, this means that we should have, for primary fields We thus have, for the 2S part while for the 2A part Similarly, in the sector with j = 4 propagating clusters we have The corresponding conformal weights from (46) are but the fourth set is identical to the second one. As the number of clusters increases, so does the number of allowed sectors. Our finite-sizes observations are clearly in favour of the extension of this pattern, invariably with all even values of j. Clearly, we find exponents (h r,s , h r,−s ) in the s-channel with r ∈ Z and s ∈ 2Z, and r ∈ Z + 1/2, s ∈ 2Z. We will refer to these exponents as sets S Z,2Z and S Z+1/2,2Z , in analogy with [1]. These sub-spectra arise from the modules W j,1 and W j,−1 with j even. But we find that the s-channel, in finite size at least, contains many more exponents, arising from phases z 2 = 1, −1 and, in terms of exponents, corresponding to rational values of the first label with higher denominators, such as those with first label e + 1/4. The next key question is whether some sort of simplification might occur in the scaling limit-for instance, whether some sectors that contribute to the probabilities in finite size might do so with amplitudes that go to zero as L → ∞. We have seen absolutely no evidence of this. To make the point as clear as possible, we illustrate it on the case of the antisymmetric combination P abab − P abba .

Amplitudes and the antisymmetric combination P abab − P abba
The antisymmetry of the combination implies that only modules with z j = −1 contribute, which translates into primaries with h −h an odd number-what we have called earlier the j even, A sectors.
To make things concrete, we can take for instance Q = 1/2 (so m is irrational), in which case we find We observe that the field with (h 1/2,−2 + 2, h 1/2,2 ) has total dimension larger than (h 3/2,−2 , h 3/2,2 ). Therefore, at momentum 3, the field (h 3/2,−2 , h 3/2,2 ) will be the first contribution, and so will be (h 1/2,−2 , h 1/2,2 ) at momentum one. It is therefore very easy to identify the corresponding contributions to the four-point function: Since m is irrational, there is no mixing in the conformal mapping, and we have on the cylinder To restate the obvious, what we do then is measure the combination of probabilities on the left, identify the various terms on the right (via the exponential l-dependence of ξ,ξ), and account for the geometrical factors (the powers of 4 sin 2 2πa L ) to extract, for a given sizes L, an estimate of the amplitudes.

The (1/4, ∓4) amplitude
We give in figure (3)   While the amplitude is small in general, it is found to become large-nay divergent-for two special values: There are several ways to understand this. We will discuss a CFT analysis in the conclusion. From the lattice point of view, the divergence arises because the transfer matrix exhibits a Jordan cell in the lowest level of W 4,±i . This Jordan cell arises from representation theory of the Jones algebra for q = e iπ/8 , q = e 3iπ/8 . To illustrate this, take for instance the case q = e 3iπ/8 . The module W 2,−1 becomes reducible for this value of q, and admits a sequence of submodules as represented in Figure 4. The presence of submodules W 4,i , W 4,−i (in particular) suggests 8 that excited states in W 2,−1 (a module with two propagating clusters) mix with states in W 4,i , W 4,−i within the module involving four propagating clusters. This mixture leads to Jordan cells in the transfer matrix. As shown in (54)-and further discussed in Appendix B.4 in the case Q = 0-a Jordan cell in turns translates into a contribution to the correlation function that is linear in (imaginary) time on the cylinder. This, finally, corresponds formally to an infinite amplitude.

5.3
The s-channel spectrum of P aaaa : F j,e 2iπp/M (M |j, j≥2), j, jp/M even and F 0,−1 We now turn to P aaaa : this is the probability that all four points belong to the same cluster. It is called P 0 in [1]. For all finite sizes, we find that the modules W j,e 2iπp/M with M |j, j ≥ 2 and j even contribute when jp/M is even: this corresponds to the sectors with an even number of clusters propagating, and values of z Figure 4: Sub-module structure of W 2,−1 for q = e 3iπ/8 . Note the appearance of sub-modules isomorphic to W 4,±i , which lead to glueing of standard modules into bigger, indecomposable modules, and Jordan cells for the transfer matrix.
obeying z j = 1, what we have called earlier the j even, S sectors. Geometrically, these contributions arise from configurations where for instance the points 1, 3 and 2, 4 are joined by two clusters which are only connected outside of the interval between their two (imaginary) time slices. On top of this, we also have the contribution where the four points belong to a single cluster arising between their two (imaginary) time slices. As discussed earlier, having a cluster propagating along the cylinder does not imply that there are boundaries around the cluster. The corresponding module of the Jones algebra is thus not a module with j = 1: rather, it occurs as W 0,−1 , i.e., as a module with no through-lines, but for which non-contractible loops (which would cut the connection between 1, 2 and 3, 4) are forbidden. Like for P abab and P abba we find that all eigenvalues in these modules do contribute in finite size, and that none of the amplitudes seem to vanish as L → ∞. This suggests that the spectrum of critical exponents is given by F j,e 2iπp/M (M |j, j≥2), j, jp/M even and F 0,−1 .
In the two clusters (j = 2) sector, this leads to while in the four clusters (j = 4) sector we find These two contributions occur as well in S Z+1/2,2Z . New contributions appear for higher even values of j.
For instance we find also On top of this we have the 'one-cluster sector', which is described by F 0,−1 (i.e., non-contractible loops are killed). This corresponds to the set of conformal weights which is also in S Z+1/2,2Z .

5.4
The s-channel spectrum of P aabb : F j,e 2iπp/M (M |j, j≥2), j, jp/M even, F 0,−1 and F 0,q 2 The quantity P aabb is the probability for two "short clusters" (as opposed to the "long clusters" shown in Figure 2): points 1, 2 belonging to one cluster, points 3, 4 to the other. It is called P 1 in [1]. We find that all the eigenvalues occurring in P aaaa also contribute to P aabb . On top of these, we also find the eigenvalues from the module W 0,q 2 = W 0,q −2 . This module corresponds to a sector with no (forced) propagating cluster, which is obtained simply by giving non-contractible loops their bulk weight. As usual now, none of the corresponding amplitudes seem to vanish in the limit L → ∞. The operator content from W 0,q ±2 involves diagonal primaries, with weights where Of course this is the same set as the set after a shift of the electric charge. We will denote this set as S d Z,1 .

Summary
We can now summarise our spectra in the s-channel where we have allowed j to take positive or negative values, since the sets of exponents are invariant under Recall also that p, M are coprime integers, and that the value p = 0 in particular is allowed. The case p M = 1 2 appears already in [1]. Note that these are the generic results, i.e., those valid for m irrational. Some contributions vanish for special values of Q, such as Q = 0, Q = 2 and Q = 4 (see Appendix B) and in some cases Jordan blocks appear.
The spectra in the other channels follow from simple geometrical considerations: An important property of our spectra in the case of P aabb is that only states with positive conformal weights propagate along the cylinder: no "effective central charge" appears, despite the non-unitarity of the CFT. This is contrast with what would be observed, for instance, in the case of minimal models corresponding to m + 1 ≡ p p , p, p integer, where the effective ground state with c eff = 1 − 6 pp would appear. It is our understanding that a similar phenomenon takes place in the conjectured expressions of [1].
6 Comparison with results in [1] The comparison with the proposal in [1] requires some discussion, since the authors in this reference did not, in particular, provide conjectured results for P aaaa . The simplest quantity to consider is P abab − P abba ≡ P 2 − P 3 in the notations of that reference. Indeed, from eq. (3.2) in [1] we see that, in their notations, R 2 − R 3 = λµ(P 2 − P 3 ). The spectrum in the s-channel of P 2 − P 3 = P abab − P abba , according to our analysis, is made of the fields (h e+p/M,−j ; h e+p/M,j ) for e ∈ Z, with j even and pj/M an odd integer. Note that all these fields have h −h odd. In [1], meanwhile, the spectrum is S Z+1/2,2Z (after switching indices in [1] to make their conventions the same as ours), restricted like for us to odd spin h −h. So for instance the field with weights (h 1/4,∓4 ; h 1/4,±4 ) for which we have seen that the amplitude is generically non-zero, is absent in the solution proposed in [1]. This suggests that their solution is, generically, not the correct one, and that an infinity of fields is missing in their proposal. We illustrate this qualitatively in Figures 5,6. Meanwhile, it is fascinating to compare results for amplitudes that are predicted in [1] and which are also found to occur in our analysis. A good example of this is the first amplitudes for the sector with j = 2, namely A Φ h 3/2,−2 ,h 3/2,2 and A Φ h 1/2,−2 ,h 1/2,2 . The bootstrap in [1] produces amplitudes which are in fact simply related with those of Liouville field theory at c < 1, and thus admit analytical expressions [49,50]. In particular, their conjecture is where β 2 = m m+1 , . This conjecture reproduces results which are believed to be exact at Q = 0, 3, 4-the result for Q = 0 is discussed in our Appendix B.4; the result for Q = 4 follows from a work by A. Zamolodchikov (as discussed in [1]), and the result for Q = 3 is unpublished work of R. Santachiara.
Numerical results for this ratio are given in figure 7. They are intriguingly close -after reasonable extrapolation-to the formula (84). The agreement is worse near Q = 4, but as commented elsewhere in this paper, this discrepancy can possibly be attributed to the presence of a marginal operator affecting corrections to scaling. We do not know whether (84) might actually be exact, or whether it is just very close to the exact result. Numerics, at this stage, do not really allow us to settle this issue.
Meanwhile, the uncertainty of the numerical determination shown in Figure 7 can be estimated from the difference between the extrapolations through even and odd system sizes L. Given that this uncertainty is (for most values of Q) comparable to the distance to the conjectured result (84) is certainly a strong motivation for further improving the numerical algorithm and gain access to a few more sizes. This could  Figure 5 maybe be achieved if one could impose the sector and momentum constraints within our scalar product method (see Appendix A.2).

Conclusion
We believe that the numerical and algebraic evidence presented in this paper invalidates the results in [1]. This is a very intriguing conclusion, since, in particular, the authors of [1] presented Monte Carlo simulations of four-point functions in the plane that were in good agreement with their bootstrap prediction. It is possible that the conjecture in [1], while not the correct answer to the problem of describing geometrical correlations in the Potts model, is indeed a solution to the bootstrap, and moreover captures numerically the essential features of the four-point functions, failing only at an accuracy, or for values of the cross-ratio z, not accessible using the Monte-Carlo approach. If this is the case, this raises several questions, in particular about the number of possible solutions to the bootstrap, 9 and what, if anything, is truly described by the proposal in [1].
To shed more light on this issue, an obvious route is to build four-point functions following the methodology in [1] but based on our spectra. This is quite challenging technically, because of the large number of primary fields with dimensions of the same order of magnitude we would have to involve. Another, more fundamental aspect worth mentioning is that, in our spectrum, many of the conformal weights have degenerate values, with singular conformal blocks. It is not clear whether the regularisation procedure used in the bootstrap approach [1] is actually the relevant one for the Q-state Potts model. This, we believe, could be answered by numerical studies in the spirit of the present paper and [53].
A particularly intriguing fact is that we found numerically a ratio A Φ h 3/2,−2 ,h 3/2,2 /A Φ h 1/2,−2 ,h 1/2,2 which is not incompatible with the proposal in [1]; see Figure 7. It could be that the solution to the bootstrap relevant for the Potts model involves for this ratio a value close to the one in [1] and yet different, over the whole range Q ∈ [0, 4]. It could also be that the amplitudes in [1]-which, to the best of our understanding, are actually given by standard formulae for Liouville at c < 1, naively extended to the case of fields with h =h-are exact, but that something has to be added.
Adding "something" to the spectrum in [1] is definitely necessary if one wishes to avoid correlation functions with many singularities as Q is varied. To see why this is the case, we consider the contributions to the antisymmetric combination of probabilities: including now higher order terms in the conformal blocks Note that we used here conformal blocks where the dependency z −2h 1/2,0 (resp.z −2h 1/2,0 ) has been factored out. When Q → 4 cos 2 3π 8 ≡ Q * , we find that h 3/2,2 → h 1,2 , a degenerate value. The conformal block in the four-point function coming from (h 3/2,−2 , h 3/2,2 ) = (h 3/2,−2 , h 1,2 ) = − 1 32 has a null-state at level 2 for thez components, with weights The appearance of the null-state means that the conformal block F h 3/2,2 has a pole of the form we see that the amplitude of the singular term in the bracket in (85) is Meanwhile we observe that the weights for the singular term coincide with the weights from the 1/4 sector: Recall that we have found numerically that the amplitude of this field also has a simple pole when Q → Q * : so the amplitude of the second singular term in (85) is For technical reasons, we normalise all quantities by A Φ h 1/2,−2 ,h 1/2,2 (this amplitude is not expected to be singular), so we set We find numerically that on "resonance", there is a Jordan cell of rank two mixing the two terms, but no singularity. This means that we should have the condition To put things differently, the appearance of a null-state in the conformal block F h 3/2,2 leads to a divergence in the four-point function (assuming the formula for A Φ h 3/2,−2 ,h 3/2,2 given earlier is correct). To cancel this divergence, a contribution A Φ h 1/4,−4 ,h 1/4,4 is necessary. Moreover, this contribution must exhibit a simple pole, as we have observed numerically. It is possible that this picture generalises, with singularities in the proposal of [1] exactly cancelled out by the additional terms we find in our lattice analysis. This will be discussed elsewhere [54].
To conclude this paper, we re-iterate the remark that the eigenvalues contributing to the probabilities are (a subset of) those appearing in the Potts model partition function [7]. While this would be a well expected fact for a model defined locally such as the Ising model or any kind of height model, this is not so obvious in our case. Indeed, in a model where correlations are defined non-locally there is no clear connection between the partition function and at least some of the observables. To give a simple example, we know well that the probability that two points are connected with a cluster allowing two independent paths involves exponents not present in the partition function [33,35]. The fact that no such exponents are needed for the P a1a2a3a4 suggests that these are rather close to "ordinary" observables, and that we may be able to understand them in terms of fully local operators acting on the space of states. One of the main "elementary" mysteries in this description is why only sectors with an even number of clusters contribute. We believe that thinking more deeply about algebraic aspects of the problem on the lattice will shed some light on this question.

A Computing four-point functions on the cylinder
In this appendix we gather all the algebraic and numerical technology that enables us to compute four-point functions in the FK cluster model on a cylinder.
The geometrical setup is shown in Figure 1. For the sake of simplicity we shall define the model on an axially oriented square lattice of width L spins. The four points are inserted in two different time slices, separated by l lattice spacings, with points w 1 , w 2 in the first slice and points w 3 , w 4 in the second slice.
The algebraic tool is the transfer matrix T based on the join-detach algebra, which we have already briefly reviewed in section 3.1. To adapt it to the computation of correlation functions the crucial point is the representation of the algebra, i.e., the choice of the set of states on which T acts. While the partition function Z can be simply computed by taking these states to be set partitions of L points (see [55,56] for details) we shall need to endow these partitions with various kinds of marks.
We are interested in two different situations that present some subtle differences. First we describe how to compute directly the 15 different FK cluster correlation functions P a1,a2,a3,a4 , defined by (6), that lie at the heart of the numerical method outlined in section 4.3.1. Second, we go into the details of the alternative method of section 4.3.2 which is based on computing scalar products involving the order parameter operators (5).
A.1 Computation of P a 1 ,a 2 ,a 3 ,a 4 Let w k (with k = 1, 2, 3, 4) be four marked points on the cylinder, as shown in Figure 1. Let P = {a 1 , a 2 , a 3 , a 4 } be a set partition of the four points, defined by the corresponding integer labels a k . We wish to construct a transfer matrix that builds the statistical weight W P corresponding to a correlation function in which the marked points k and belong to identical (resp. different) FK clusters if a k = a (resp. a k = a ). For instance, the choice P = {1, 2, 1, 2} means that w 1 and w 3 belong to the same cluster, while w 2 and w 4 also belong to a common cluster which is different from the first one, as shown in the left panel of Figure 2.
On the cylinder there are B 4 = 15 different correlation functions, where the Bell number B N denotes the number of partitions of an N -element set. On the strip there would be only C 4 = 14 correlation functions, where C N denote the Catalan numbers.
We consider the square lattice of width L wrapped on a cylinder. The sites within each row are labelled by i = 0, 1, . . . , L − 1, with the labels considered modulo L by the periodic boundary conditions. To make contact with Figure 1, we let w 1 and w 2 correspond to sites 0 and 2a (with a ∈ N/2) in the row labelled by the transfer matrix "time" t 1 = 0, while w 3 and w 4 correspond to sites x and x + 2a in the row t 2 = l. A cluster containing (at least) one of the points w k is called a marked cluster.
To allow the marked clusters to wind around their insertion points, we consider a larger piece of the lattice, going from row t 0 = −M to row t 3 = l + M . We impose free boundary conditions on the two extremites of the finite cylinder defined by t ∈ [t 0 , t 3 ]. The desired correlation functions (or probabilities) are then given by the limit For practical purposes, to obtain P P (l) to a given numerical precision, it suffices to take M sufficiently large, so that the results for M and M + 1 coincide to the chosen precision. A more precise criterion for the choice of M can be given once the spectrum of the transfer matrix is known (see section A.7.2). The transfer matrix acts as usual on states {s 1 , s 2 , . . . , s L } which are certain set partitions of L points, but now endowed with suitable markings. We have s i = s j if and only if sites i and j are seen to be in the same FK cluster at a given time t-by this we mean that i and j are connected by a piece of FK cluster on the partly constructed cylinder t ∈ [t 0 , t]. The symbols s i can take the values 1, 2, . . . , L for unmarked clusters, and1,2,3,4 for marked clusters containing one of the points w k . If a marked cluster contains more than one marked point, the chosen symbol is the lowest one, in order to avoid any redundancy in the correspondence between cluster connectivities and states. Unmarked clusters are indistinguishable.
As shown in (11) the transfer matrix can be written as a product of elementary operators H i = I + vJ i and V i = vI + D i that add respectively a horizontal edge between sites i and i + 1 (mod L), and a vertical edge on top of site i. We have v = e K − 1, with K the Potts coupling, and the critical value on the square lattice is v c = √ Q. The join and detach operators, J i and D i satisfy the join-detach algebra with relations (12). We now describe a modified representation of this algebra that properly takes into account the possibility of having marked clusters.
A clusters is said to be left behind at site i if D i detaches a singleton in the set partition. When this happens, D i applies a weight Q, irrespective of whether the cluster being left behind is unmarked or marked.
Remark 8. In some applications it is natural to give weight 1 to marked clusters. In particular, this is mandatory in the limit Q → 0, since otherwise all correlation functions would be identically zero. However, for the time being we choose to give weight Q to any cluster, including the marked ones, since then 15 r=1 W Pr = Z(Q, v), the partition function on a cylinder of circumference L and length l +2M . The quantity Z(Q, v) can easily be obtained by independent means [57], and the sum rule then provides a valuable check of the correctness of the algorithm.
To apply the initial condition, we start from the all-singleton state {1, 2, . . . , L} at time t 0 . The final condition at time t 3 is to project out any state in which the desired connections between marked clusters have not been achieved. To be precise, if there exists two distinct labelsk =¯ in the connectivity state such that a k = a , then the state must be discarded. Any state that survives this projection is attributed a weight Q per cluster (again irrespectively of whether it is unmarked or marked). The weighted sum over retained states produces the sought-for W P .
Marked symbols are introduced in the time evolution by joining site i = 0 (resp. i = 2a) to the marked symbol1 (resp.2) at time t = t 1 , and by joining site i = x (resp. i = x + 2a) to the marked symbol3 (resp. 4) at time t = t 2 .
To impose the desired four-point connectivity, we modify the action of the join operators J i as follows: • When J i joins an unmarked cluster to a marked cluster with symbolk, the result is a marked cluster with the same labelk.
• J i can join two marked clusters with labelsk and¯ only if a k = a . In that case, the result is a marked cluster with the lowest label min(k,¯ ).
The detach operators D i need a more subtle modification. We describe the marked points w 1 and w 2 as unseen (i.e., not yet visited by the transfer process) when the time t < t 1 . Similarly the marked points w 3 and w 4 are unseen when t < t 2 . The required modification is: • The clusterk is allowed to be left behind only if the label a k is different from the labels of any other marked cluster in the connectivity state, and also different from the labels of all yet unseen marked points.

A.1.1 Example
Consider the case P = {a k } = {1, 2, 1, 2} of two propagating clusters. At times t 1 < t < t 2 , the points w 1 and w 2 have been seen, so the connectivity state can contain the symbols1 and2. In fact, both symbols must be present, because the J i operator cannot join1 and2 (since we have a 1 = a 2 ); moreover none of them can be left behind, because points P 3 and P 4 are unseen. Indeed, a 1 = a 3 then prevents1 from being left behind, and a 2 = a 4 prevents2 from being left behind. At later times t 2 < t < t 3 , all points have been seen. Symbols1 and3 can join, because a 1 = a 3 (and similarly for2 and4). But before this happens, none of the clusters1 and3 can be left behind (by the rule on D i ). When1 and3 join, the resulting cluster carries the symbol1 (by the second rule on J i ). After this happens,1 can be left behind (by the rule on D i , because3 is no longer used in the connectivity state).

A.1.2 Checks
In addition to extensive checks for small systems, where all configurations can be drawn by hand, we have checked that: Moreover, for the situation with shift x = 0 (see Figure 1) it is non-trivial from the point of view of the transfer algorithm that the following lattice symmetries hold true: 3. The four probabilities in which three points are in the same cluster are all equal.
4. The six probabilities in which two points are in the same cluster and the other two are in two distinct clusters, are equal two by two.
A more restricted set of lattice symmetries holds true also for x = 0; for instance the four probabilities in which three points are in the same cluster, are equal two by two. We have checked this as well.

A.1.3 Sample results
To help the readers desirous of writing their own implementation, we here give some sample results for P P (l) for a system with Q = 1 and size L = 5. We take distance l = 3 between the time slices and place the points in each slice as nearest neighbours (2a = 1), with no offset between the two slices (x = 0). The set of 15 probabilities (up to the symmetries described above) are then: It was necessary to take M = 200 to achieve 50 correct digits. Note that our code is written so that the number of digits of numerical precision can be adjusted to any desired value. As a final check we consider the limit l L. Then we expect where p is the finite probability that two nearest neighbours are in the same cluster. All other P Pr (l) will be negligible in that limit. By taking L = 5 and l = 384 (with still M = 200) we have verified that this is indeed the case, and we estimate p = 0.76315602507834269413 .
A.1.4 Case of Q → 0 The limit Q → 0 (with v = √ Q → 0 as well) is somewhat special from the point of view of normalisations, since then Z(Q, v) = 0. The partition function which is used to normalise the weights W P and turn them into probabilities, as in (95), is then taken as Z = W aaaa , the number of spanning trees containing all four marked points. Moreover, any marked cluster is assigned a weight 1 (instead of Q in the general case), whereas unmarked clusters still have the weight Q = 0, which implies that such clusters are disallowed. G a 1 ,a 2 ,a 3 ,a 4 The second numerical method presented in section 4.3.2 provides a means of computing one by one the amplitudes appearing in the order parameter correlation functions G a1,a2,a3,a4 defined in (7). This method is of a very different nature than the one just presented (i.e., for the computation of P a1,a2,a3,a4 ), since it does not compute the four-point correlator directly on a cylinder of finite length, but rather goes directly for the asymptotic quantities A i . Thus, the appearance of |v 0 and v 0 | in (57) amounts effectively to taking the limit M → ∞ of the distance to the boundary conditions, and the piece |v i,j v i,j | projects directly on an intermediate state in the s-channel, which is equivalent to taking the limit l → ∞.

A.2 Computation of
As a prerequisite to this approach-and in order recover in particular the same results as with the first method of section 4.3.1-it is however necessary to produce yet another representation of the join-detach algebra. In particular, we shall need the action of the order parameter operators O a to be well-defined for generic values of Q ∈ R. By (5), this hinges on giving a well-defined meaning-and in particular, a meaning that extends to non-integer Q-to the operator δ σ k ,a that fixes the value of the spin at point w k to be equal to the label a. Note that such extensions of order-parameter related quantities have played an important role in several recent pieces of work involving the present authors [11][12][13][19][20][21].
As stated around (9) we are here only interested in correlation functions G a1,a2,a3,a4 employing at most two distinct symbols, a and b. We shall therefore take T to act on basis states which are set partitions of L points with two marked values, denoted a and b. It is crucial that these values are now different by definition. This provides a difference with the computation of P a1,a2,a3,a4 , where two differently marked clusters could be joined under some circumstances. More precisely, the rules are now the following: • The operator δ σ k ,a (the non-trivial part of O a ) acts at the point w k by transforming an unmarked cluster touching that site into a marked cluster of label a. If the cluster is already marked with label b, O a acts as the identity times δ a,b .
• The join operator J i acts as usual on two unmarked clusters. When acting on an unmarked cluster and a marked cluster of label a, the result is a marked cluster of label a. Finally, when acting on two marked clusters of labels a and b, J i acts as the identity times δ a,b .
• The detach operator D i transforms site i into an unmarked singleton. If i was already a singleton beforehand, a weight Q is applied if the corresponding cluster is unmarked, and 1 if it is marked.
To compute the scalar products in (57) it is convenient to produce all the intervening vectors in the same space. In particular, the ground state eigenvector |v 0 is written within the space of set partitions with two marked values, although it is easily seen that its component along any state containing marked clusters is zero.

A.2.1 Orthogonalisation
It has already been stated in the main text that the scalar product method relies on the left and right eigenvectors being orthogonal, even within degenerate subspaces. Standard numerical methods for nonsymmetric matrices, such as the Arnoldi algorithm, do not immediately produce the eigenvectors in this form. Rather, within each degeneral subspace corresponding to the eigenvalue Λ i , the scalar product of eigenvectors come out as v i,j |v i ,j = α where the nombers α (i) j,j can be viewed as the elements of some matrix α (i) . However, if we replace the right eigenvectors by the linear combinations it is easy to see that we obtain the orthonormality v i,j | v i ,j = δ j,j provided we take β (i) = α (i) −1 . In practice, the size of degenerate subspaces is rather small (of dimension 1, 2, 4 or 8 in the problem at hand), and so any elementary method of producing the inverse matrix α (i) −1 will solve the problem conveniently.

A.2.2 Checks
We have made extensive checks that the A i obtained from the scalar product method are identical to those obtained from the more involved first method based on (53), provided one takes into account the linear relations (9), and forms the sum over orthogonalised degenerated subspaces in (57). We have verified that the scalar product method also gives the correct amplitudes of G aaaa in the simpler case where the states of T are constrained to have only one marked value a. It seems likely-although we have not actually tried this-that it will also extend to the most general G a1,a2,a3,a4 provided the number of marked values is (at least) equal to the number of different symbols a k in the correlation fucntion.

A.3 Spectrum of T
At many occasions throughout this work we need to obtain the eigenvalues Λ i of the transfer matrix. This is obviously a much easier problem than obtaining the correlation functions, and has been discussed in many places, so we can be rather brief.

A.3.1 First representation
The Λ i are related to the asymptotic decay of the two-point functions, so in analogy with section A.1 they can be obtained by using the representation of the join-detach algebra in the regime t 1 < t < t 2 , where the boundary conditions t 1 → −∞ and t 2 → ∞ have been pushed to the extremities of the cylinder. More precisely, we are interested in the sector with propagating FK clusters, so we can build on the same representation as in section A.1, but with states employing distinct marked symbols. The action of the elementary operators J i and D i must ensure that distinct clusters propagate along the cylinder, and the rules match those of section A.1 for t 1 < t < t 2 : • J i cannot join two marked clusters corresponding to different labels.
• D i cannot leave behind any marked cluster.
• The different labels are treated as indistinguishable.
The last rule means that if the first site in the cluster with markk is denoted i k , we quotient the set of states in order to impose i 1 < i 2 < . . . < i .
We have checked that the spectrum of this transfer matrix coincides with that of the loop model (3) in the affine TL representations (22), where as usual j = for j ≥ 2, while the first and the third terms in the direct sum correspond to = 0 and = 1 respectively.

A.3.2 Second representation
It is also of interest to express T in the representation of section A.2. There are now precisely two different marks, and the rules read: • J i cannot join two marked clusters corresponding to different labels.
• D i can leave behind a marked cluster (with weight 1).
• The different labels are treated as distinguishable.
We find that the spectrum of this transfer matrix reproduces precisely (22), but only for j = 0, 2, 4, . . ., including the structure W 0,q ±2 ⊕ W 0,−1 with two direct summands for j = 0 and the absence of (at least the first few) odd terms (j = 1, 3). At present it is not obvious to us what is the precise decomposition of this representation in terms of simple affine TL modules, but we take the observations just mentioned as a first sign that our main result on the s-channel of four-point functions might have a natural interpretation within this representation of the join-detach algebra. We therefore suggest that it might be worthwhile to study more precisely the algebra obtained from generators O a (σ i ), J i and D i , which contains the usual join-detach algebra as a sub-algebra.

A.4 Momentum sectors of T
As discussed in section 4.3.3, it is important to be able to associate the eigenvalues Λ i of T with the conformal properties that emerge in the continuum limit. To this end we shall need to attach to each Λ i more information than just the number of propagating FK clusters that it corresponds to.
It is most practical for our purposes to consider here the loop model, in its representation of standard modules W j,z 2 =e 2iK within the affine TL algebra T a N (n). We denote in this section by T N,2 the corresponding transfer matrix on link patterns with N = 2L strands (recall that L is the number of Potts spins in a row) and 2 through-lines. Each through-line picks up a factor z (resp. z −1 ) when it traverses the periodic boundary condition towards the right (resp. left).
Our goal is to make apparant two more quantum numbers: the momentum k of the through-lines, and the lattice momentum m. The former comes directly from the quantisation of z, and the latter is equal, in finite size, to the conformal spin: m = h −h. We shall obtain each momentum sector by transforming T N,2 into an appropriate matrix of smaller dimension.
To this end we proceed as follows. The dimension of W j,z is N N/2−j . Within this space we choose one representative state for each orbit of the cyclic group C L generated by u 2 , where u is the affine TL shift operator. Note that [T N,2 , u 2 ] = 0 for the Potts-model transfer matrix. The orbit length corresponding to a representative state |s is denoted g s , so the orbit can be written |s , u 2 |s , u 4 |s , . . . , u 2(gs−1) |s .
Obviously g s is a divisor of L, and the weighted sum over representative states (with each |s being weighted by its orbit length g s ) equals dim W j,z .
Example 1. For L = 3 and = 0 the dimension of W j,z is 6 3 = 20 (we do not need to take the quotient yet), and there are 8 distinct orbits (namely 6 with g s = 3, and 2 with g s = 1). For L = 8 and = 4 the dimension of W j,z is 16 4 = 1820, and there are 224 orbits with g s = 8, 6 orbits with g s = 4, and 2 orbits with g s = 2.
The motion within orbits gives rise to another momentum variable, distinct from the twist z of the affine TL algebra. Recall that the twist variable is z = e iπk/ for each through-line traversing the seam towards the right, and since the total phase factor for a turn of all 2 through-lines must be z 2 = 1 we have k = 0, 1, . . . , − 1. Similarly we now wish to construct the sector of lattice momentum ω m ≡ (ω) m , where ω ≡ e 2πi/L , by attributing a weight ω m to each translation of one Potts spin (hence two TL loop strands) within the orbits of C L . Since the total phase factor for a rotation through L spins must be (ω m ) L = 1 we have m = 0, 1, . . . , L − 1. We stress that k is related to the number of FK clusters , whereas m is related to the system size L.
We now wish to construct a restriction T N,2 ,m to given values of the momentum labels k ∈ Z and m ∈ Z L . This means that the spectrum of T N,2 (n, z) must be decomposed as the union of the spectra of T N,2 ,m , including multiplicities: spec T N,2 (n, z) = L−1 m=0 spec T N,2 ,m (n, z, ω m ) .
It is practical for us to denote the set of eigenvalues of T N,2 ,m (n, z, ω m ), with specified values of the momenta z = e iπk/ and ω m = e 2iπm/L , simply as V ,k,m . We denote the revelant dimensions as D = dim W j,z and d ,k,m = dim V ,k,m , where d ,k,m will be determined below. We should of course have for any = 0, 1, . . . , L and k ∈ Z , in accordance with (102). Our approach is to construct, for each m, a D × d ,k,m matrix S in , and a d ,k,m × D matrix S out (z, ω), such that T N,2 ,m (n, z, ω m ) = S out (z, ω)T N,2 (n, z)S in .
In the matrix S in , each state |s in the restricted space is mapped to the corresponding representative state in the full space, with a Boltzmann weight equal to g s . In other words, if the representative state |s is ordered as the j'th basis state in the restricted space and as the i'th basis state in the full space, then (S in ) ij = g s . All other matrix elements are zero.
In the matrix S out (z, ω), each state |t in the full space is identified as |t = u 2k |s , where |s is the representative state corresponding to |t (i.e., with |s living in the restricted space) and k the number of double shifts (by convention, towards the left) necessary to bring |t into the representative form. Notice that under these shifts, it is possible that a number of through-lines p will cross the periodic boundary condition (towards the left). We must keep in mind that the convention for T a N (n) is that each through-line that crosses the periodic boundary condition towards the right (resp. left) acquires the weight z (resp. z −1 ). Therefore, the Boltzmann associated with bringing |t into the representative form |s is (ω m ) k z −p /g s . Thus, if |t is ordered as the j'th basis state in the full space and its representative |s is the i'th basis state in the restricted space, then All other matrix elements are zero. A crucial point is that the orbit lengths must be compatible with the momentum that we impose. To be precise, the restricted states corresponding to given values of the labels k, m are those with orbit lengths g s satisfying (k − m)g s = 0 mod L .
The dimension d ,k,m is precisely determined by the number of restricted states satisfying this constraint.  (103) we notice that for any k ∈ Z , the eight values of m ∈ Z L corresponds to 8 cases where g s = 8 is allowed, 4 cases where g s = 4 is also allowed, and 2 cases where g s = 2 is also allowed, i.e., we recognise the values of the g s themselves in this count. Thus 8 × 224 + 6 × 4 + 2 × 2 = 1820, as it should. By the same reasoning, (103) is verified for arbitrary values of L and .
We have written an algorithm that produces the basis change matrices S in and S out (z, ω) very efficiently. The crux is obviously to deal with the ordering of the states and the identification of the relevant orbits. We stress that the construction with basis change matrices S in and S out (z, ω) is perfectly compatible with iterative diagonalisation schemes (such as the Arnoldi method), which are particularly efficient for dealing with transfer matrices that allow for a sparse matrix factorisation.
Finally we note that the spaces V ,k,m are invariant upon simultaneously changing the signs of both momenta (modulo and L, respectively). Thus We have checked this exhaustively for L = 5 and L = 6.

A.4.1 Checks
To check the algorithm, we have performed extensive tests on the case N = 2L = 10. We have first checked that [u 2 , T N,2 (n, z)] = 0 for = 0, 1, . . . , L and any values of n = √ Q and z. The sector decomposition was checked in details for a generic value Q = 3 2 . Our general result is that the union of V ,k,m where k = 0, 1, . . . , − 1 and m = 0, 1, . . . , L − 1 produces, for each , precisely the spectrum of the FK transfer matrix T N,2 with the correct multiplicities for degenerate eigenvalues. The advantage of this construction is thus twofold: it suffices to diagonalise much smaller matrices, and at the same time we can associate the quantum numbers ( , k, m) with each eigenvalue.
• For = 0 we formally set k = 0. Without taking theW j,z 2 quoitient, we have compared the 42 eigenvalues obtained from numerically diagonalising the FK transfer matrix with those obtained from the momentum sector decomposition of the full T a N (n) module. This allows us to assign the correct momentum label to each FK eigenvalue. We find that V 0,0,0 has 10 eigenvalues, while V 0,0,1 and V 0,0,2 each have 8. Taking into account the ±m degeneracies, this gives all the required 10 + 2 × (8 + 8) = 42 eigenvalues indeed.
• For = 1 we also formally set k = 0. We find that V 1,0,0 has 52 eigenvalues, but only 21 distinct eigenvalues; of these 5 are fourfold degenerate and the remaining 16 are twofold degenerate. Each of V 1,0,1 and V 1,0,2 contains 50 eigenvalues, but only 25 are distinct; each of these are twofold degenerate. This accounts for the required 52+2×(50+50) = 252 eigenvalues of the affine TL module. The overall twofold degeneracy is explained by the fact, that since there are no winding loops (recall d = 0) and the square lattice is selfdual, the T a N (n) module in fact decomposes in two isomorphic representations each of dimension 1 2 2L L = 126. However, this still leaves 5 degenerate eigenvalues in V 1,0,0 . These degeneracies are due to the choice of the square lattice, which is not only symmetric under rotations but also under reflections (i.e., the symmetry is the dihedral group D 5 , which is larger than the cyclic group C 5 ).
• The case = 5 is obviously very degenerate, as will generally be the case when = L.

A.5 Spectrum of T and the CFT limit
We wish to verify the CFT interpretation of the spectrum of the FK model transfer matrix T in the sectors V ,k,m . As above, denotes the number of propagating clusters, k is the twist parameter for the affine TL through-lines (with allowed values k = 0 for = 0, and k ∈ Z for ≥ 1), and m is the lattice momentum (with allowed values m ∈ Z L for a periodic strip of width L Potts spins). Moreover we have the symmetry (107).
As we shall see, the conformal interpretation of the labels k, m is as follows. The twist label k fixes the first Kac label in φ r,s to be r = k/ + e, where e is an integer. The momentum label m detemines the conformal spin, s = h −h = m mod L. Below we give detailed evidence that if the true value of the conformal spin is too large to be accommodated in a given size L (that is, |h −h| > L/2), there will in general be an appropriate eigenvalue with m = (h −h) mod L that nicely fits into the finite-size scaling formulae.
In the sequel we consider strips of widths L ≥ 5. For technical reasons, it is feasible to diagonalise T with respect to all quantum numbers ( , k, m) only for L ≤ 7, but such sizes are insufficient to numerically determine the scaling dimensions h +b for a significant number of low-lying excitations. However, we are able diagonalise T in each sector V ≡ k,m V ,k,m for higher values of L, by simply imposing the number of through-lines (and for = 1, also setting to zero the fugacity of non-contractible loops), without having to compute the basis change matrices S in and S out appearing in (104). Concretely, we have done so for = 0 and L ≤ 14, for = 1 and L ≤ 13, and for = 2 and L ≤ 12.
For a fixed value of , the eigenvalues are labelled Λ (i) (L), where we have supposed the ordering Λ (1) > Λ (2) > · · · . Note in particular that we disregard multiplicities (the eigenvalues with m = 0 are two-fold degenerate, as are those with m = 0 and k = 0). We shall refer to i as the "rank" of the eigenvalue within V . From the eigenvalues we can compute the effective scaling dimension x = h +h via The practical problem is that the rank i = i L of a given scaling level in V ,k,m will depend on L. We would of course expect i to stabilise at some constant value i ∞ for L sufficiently large, for the level that becomes the i ∞ 'st lowest-lying excitation in the conformal limit. The trouble is that for all but the few lowest excitations, i ∞ turns out to be comparable to or larger than the values of L that can be accessed numerically. As explained above, we however know the values of i L that are compatible with given (k, m) for L ≤ 7. Fortunately, by plotting x (i) (L) against 1/L, and fitting the available sizes to a polynomial in 1/L, we can gradually reconstruct the sequence i L by adding one size at a time, by carefully monitoring the quality of the fit and checking whether the extrapolation produces a reasonable value. In this way we obtain excellent fits that extend in general from L min = max (5, 2m) to the largest accessible size L max . The outcome of this procedure can be verified by checking that, for any given L, each value of i is attributed to one and only one scaling level. Moreover, the extrapolated scaling dimension x (i∞) should make sense in the CFT, being in particular compatible with the values of (k, m) that we have thus inferred. For m > 3 we cannot proceed in this way, because i L is only defined for L ≥ L min = 2m and we only know the sector labels for L ≤ 7. In this case, we can usually work the other way around, either guessing the possible values of i L from those which are not used by other levels, or by carefully monitoring the fits starting at the highest values of L (where i L can be supposed to be constant, or almost constant) and gradually proceeding to include also lower sizes and in the same time determining the corresponding i L .

A.5.1 A generic case
In Table 1 we report the outcome of these computations within V 0 , where we have taken an arbitrary and generic value Q = 1 2 . We show the sequence (i 5 , i 6 , . . .) for each level, the corresponding representation V ,k,m , and the numerically extrapolated scaling dimension x ≡ x (i∞) . Comparing this to the possible analytical values 10 allows us to identify the corresponding scaling field in the CFT limit. In the latter we use the notation L −n L −n to denote any highest-weight representation at level (n,n) in the holomorphic (resp. antiholomorphic) Verma module. For instance, L −2 denotes here what is usually called either L −2 or (L −1 ) 2 , or a linear combination of those terms. In a generic module we should then expect a multiplicity corresponding to p(n) × p(n), where p(n) denotes the number of integer partition of n. This count can however be smaller due to the presence of null vectors. Since the method just outlined for identifying the exponents may be of more general interest, we now give a detailed example of its application. The final fit, using all ten sizes L = L min , . . . , L max and a polynomial extrapolation of order 9 is shown in Figure 8. It leads to a very precise estimate x ≈ 4.50369, which is in excellent agreement with the exact scaling dimension 4.50378 · · · of the field φ 3,1 × φ 3,1 . Moreover, we can check that for the three smallest sizes, the eigenvalues with (i 5 , i 6 , i 7 ) = (11,14,16) all have lattice momentum m = 0, so the CFT field must indeed have conformal spin h −h = 0.
The contents of Table 1 can be quickly summarised by saying that we have numerically observed within V 0 three different primaries, along with their corresponding Kac modules (where the generic Verma module of φ r,s has been quotiented by one singular vector at level rs). More precisely: • For the identity operator I we have observed all descendents up to and including the total level N ≡ n +n = 5, and some of the descendents at level 6. This operator is seen to be degenerate at level 1.
• For the energy operator ε we have observed all descendents up to and including level N = 3, and some of the descendents at level 4. This operator is degenerate at level 2, as witnessed by the fact that we only observe a single state (and not p(2) = 2 states) at level (n,n) = (2, 0).
• For the second energy operator ε we have observed only the primary at level N = 0, due to the high scaling dimension x = 4.503782 · · · . Table 2 gives the results of a similar investigation for the sector V 1 . We here observe numerically two different primaries, along with their corresponding Verma modules. More precisely: • For the magnetisation operator σ we have observed all descendents up to and including the total level N = 3, and some of the descendents at level 4.
• For the second magnetisation operator σ we have observed all descendents up to and including level N = 2.
Remark 9. We stress that no degenerate states appear in these magnetisation operators, and accordingly their Verma modules are generic. This is in a nutshell why the determination of four-point functions in the bulk case is so difficult: the absence of degenerate states implies that we cannot write down differential equations satisfied by the correlation functions.
Finally, Table 3 shows our results for the sector V 2 . We see here the beginning of the conformal towers of the primaries φ e,2 × φ e,−2 for e = 0, 1 2 , 1, 3 2 , 2. As for V 1 , there are no singular vectors in this case neither.
A.5.2 The special case Q = 4 In Appendix B we show that the case Q = 4 can be compared with an exact solution. It is therefore a particularly important benchmark for the numerical method. Moreover, it is well known that the Q = 4 case is hampered by slow convergence, due to the logarithmic corrections produced by a marginally irrelevant operator. To make sure that our generic analysis applies in this case as well, we shall need the same kind of tables as above with this value of Q. They are obtained using the same methods as before, and are shown in Tables 4-6.
The precision of the extrapolated exponents x suffers somewhat from the logarithmic corrections to scaling. But assuming that the operator content of the generic case carries over, the assignment of sector labels is nevertheless certain for the levels shown in the tables. To this end, it is particularly helpful that the observed exponent difference between a descendent operator and its corresponding primary is determined with considerably better precision than the exponents themselves.
This phenomenon is vividly illustrated in Table 4. For example, the primary ε has the exact scaling dimension 1 2 , but the numerically measured exponent of 0.62 is very imprecise. However, the corresponding    3  3  3  3  3  3  3  3  3   Similarly, the exponent of the primary ε is measured as 2.35 instead of the exact value of 2. This is per se a catastrophic lack of precision-but it actually makes it easier to correctly identify the descendent L −1 ε , whose measured exponent comes out as 3.36 with an almost perfect integer gap! Without this phenomenon one could easily have mixed up the numerical value 3.36 with the other candidate exact value of 7 2 . In a similar fashion we have the Tables 5-6 for the sectors V 1 and V 2 , respectively. Note that the rank of eigenvalues whole momentum labels are fixed by analytic continuation-using the PT symmetry (107)-for small values of L are shown as tiny numbers in the tables. In the previous Tables 1-3 (for Q = 1 2 ) we have left blank such entries, although they can be determined in those cases as well.

A.6 Ising model
Finally we discuss the case of a unitary minimal model, namely the Ising model (Q = 2). A special case of the first of the identities (9) then reads [14] G aaaa = σ 1 σ 2 σ 3 σ 4 = P aaaa + P aabb + P abba + P abab , where we have used the short-hand notation σ i = Qδ Si,+ − 1 for the usual FK spin operator O a (σ i ) (5) and S i = ± are the four Ising spins. Inserting the definition of σ i , we can express G aaaa in terms of the probabilities P (S 1 , S 2 , S 3 , S 4 ) of having fixed values of the S i : It is straightforward to write a transfer matrix that computes the probabilities P (S 1 , S 2 , S 3 , S 4 ) in the Ising spin representation, by simply projecting on the required values of S i at the position of each operator. Doing this we have verified the above identity to 4000 decimal places for various systems. Since this relates very non-trivially the probabilities in the FK representation-whose computation is intricate, as we have seen in section A.1-to those in the Ising representation, this provides a strong test of the correctness of our transfer matrix setup.
We can now compute G aaaa at larger sizes and analyse it in terms of the eigenvalues of the Ising spin transfer matrix. Since two spin operators are inserted simultaneously, the propagating states should only be those of the Z 2 -even sector, and we have verified that this is indeed the case. j Added: Note that the Z 2 -even sector is the simple module X 0,i = W 0,q 2 =i − W 3,1 + · · · which is the 'top' corresponding to the V ,k,m (i 5 , i 6 , . . . , i 12 ) x Identification 5 6 7 8 9 10 11 12 Numerics Exact of scaling field V 100      Identification  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  Numerics  Exact  of scaling field  V000  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  φ1,1 × φ1,1 ≡ I  V000  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2 1 + 10 −12 1 φ2,1 × φ2,1 ≡ ε V001, V002  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  2.0000002 Table 7: Conformal spectrum in the Z 2 -even sector for the Ising model (Q = 2). The rank i k of the eigenvalues refer to this sector only, and not to the full spectrum of the corresponding FK model. leftmost diagram in Figure 9 (see [27] for details). In the conformal limit we expect the propagating states to be the identity, φ 1,1 × φ 1,1 ≡ I, and the energy operator, φ 2,1 × φ 2,1 ≡ ε. We postpone the further discussion of these results to Appendix B. However, to refine the analysis and parallel the discussion given above for generic Q and the case Q = 4, we shall again need to establish the precise correspondence between the finite-size and the conformal spectra. This is done in Table 7, where the rank of eigenvalues now refer only to the Z 2 -even part of the spin representation (and not to a sector of the full FK transfer matrix spectrum).
Note that the spin transfer matrix does not contain sufficient information to attribute a momentum label to the eigenvalues. However, since the Ising spectrum is included in the FK spectrum we can still rely on the sizes L = 5, 6, 7 to assign momentum labels to each level.
Doing this we encounter an interesting phenomenon: In some cases a given eigenvalue corresponds to two different momentum labels. For instance, the 3rd eigenvalue belongs simultaneously to V 001 and V 002 for L = 5, 6, 7; the 4th eigenvalue belongs to V 002 and V 003 for L = 6, 7; and the 7th eigenvalue belongs to V 001 and V 003 for L = 7. We recall that we count here only distinct eigenvalues, so statements of this type mean that the eigenvalue is degenerate (in addition to the usual degeneracy coming from the sign of the momentum), with different momentum identifications for each of the degenerate states. The momentum labels m assigned to the lowest-rank eigenvalues for L = 5, 6, 7 are shown in Table 8.

A.7.1 Diagonalisation of T
Our diagonalisation of the transfer matrix T is based on its decomposition (10) as a product of sparse matrices. Indeed, the elementary operators J i and D i have at most one non-zero entry per column. It is therefore feasible to compute w = T v-i.e., the action of T on a vector of weights v-without ever storing T and working only on its non-zero entries. It follows that it is highly efficient to diagonalise T by iterative methods that require only the operator w = T v and not T itself. Among such methods, we have found that the Arnoldi method is well suited for our situation, where T is a non-symmetric real matrix. We have used the C++ interface for the Fortran library Arpack [58] for producing both eigenvalues and eigenvectors, in cases where standard numerical precision (16 digits) is sufficient. This applies in particular to the scalar product method of section A.2. For the more direct method of section A.1, a considerably higher numerical precision is needed. To this end we have used the CLN package [59] (again written in C++) that performs floating point operators to any desired numerical precision. We are greatly indebted to Christian R. Scullard who has provided a pure C++ version of Arpack with templates that are compatible with CLN; this is unpublished work, but it has been described in recent publications on a different subject [60][61][62].

A.7.2 Boundary conditions
In the first method (see section 4.3.1) we compute the probabilities P P corresponding to FK cluster correlation functions on the cylinder (see Figure 1). As described in section A.1 this is done in practice by expressing the probability as the the ratio of partition functions W P on finitely long cylinders, where free boundary conditions have been imposed at a distance M from the insertion points of the cluster operators.
To ensure that the result does not depend on M , the latter much be chosen large enough. Let d denote the number of decimal digits in the desired numerical precision. Let Λ I and Λ ε be the largest and nextlargest eigenvalues of T within the sector V 0 ; it is easy to compute their values beforehand. We now claim that we must chose M so that (Λ ε /Λ I ) M ≤ 10 −d .
As the notation shows, these two eigenvalues are associated, throughout the critical regime 0 ≤ Q ≤ 4, with the scaling levels that determine respectively the identity I and the energy operator ε. It should be noticed that Λ ε thus defined is most definitely not the next-largest eigenvalue of T , if we take into account all sectors. This latter eigenvalue would rather be Λ σ , the largest eigenvalue within the sector V 1 , associated with the scaling level of the leading order parameter operator σ, which has the smallest non-zero scaling dimension in the CFT. In other words, we have Λ I > Λ σ > Λ ε . The point of the claim made above is then that outside the region t ∈ [t 1 , t 2 ] of the cylinder (see Section A.1 for notations) that contains the operator insertions, there are no constraints enforcing the propagation of clusters. This implies that the dominant decaying mode is given by the ratio Λ ε /Λ I , as can easily be checked by direct inspection of the numerically computed probabilities.

A.7.3 Attainable sizes
The largest sizes L attainable in the numerical work are essentially limited by the exponentially growing dimensions of the corresponding transfer matrices. For the first method (see section 4.3.1) there is however an additional complication, since we do not obtain the amplitudes A i directly but acquire them by solving the linear system (53).
Suppose first that there are K eigenvalues coupling to P P , and call Λ max and Λ min the largest and smallest of those. If all computations are done to d decimal digits, and we wish to obtain the final values of the amplitudes to d 0 digits, then we must chose d large enough that (Λ min /Λ max ) K > 10 d0−d . Indeed, we will need the separation between operators to take the values l = 1, 2, . . . , K in order to have enough equations to solve the linear system for A i . And assuming all A i to be of the order unity, the smallest term on the right-hand side of (53) should not be below the level of numerical precision.
Similar considerations can be made for the case where K is so large that we are unable (or do not want) to determine all the amplitudes. In that case, we truncate the sum in (53) at some i max , and provide i max different values of l. Since the formula is no longer exact, we must make sure to use it in the proper asymptotic regime. We therefore take l = i min , i min + 1, . . . , i max for some suitable large i min . In practice we have chosen i min = 100 and repeated the computation for a slightly higher value (e.g., i min = 120), in order to check that the A i were unchanged to the desired numerical precision.
In our largest computations, for L = 7, we determined more than 500 amplitudes in this way. This computation needed as much as d = 4000 digits of numerical precision and required the work of ∼ 100 computers simultaneously for a period of ∼ 3 months. See section 5.2.2 for further details on this computation.
For the second scalar product method (see section 4.3.2) it was enough to work with standard double precision (d = 16 digits). This method is altogether more efficient and could be done up to L = 11 in a much more reasonable time (at most a few days on a single computer).

A.7.4 Jordan blocks
With the first method it is possible to determine the Jordan block structure of correlation functions, by using a fit of the form (54). To do this in practice requires some reasonable prior knowledge about where to look for the Jordan blocks, and about their expected order r i . The analytical understanding of indecomposability in the Temperley-Lieb algebra has grown steadily over recent years, including in the affine TL case [26], so the search for Jordan blocks is certainly not devoid of any guidance. However, one can also take a more naive numerical approach to the problem.
Imagine that we are interested in some specific value Q = Q 0 , where indecomposability is known to occur. Consider a correlation function P P for which we have previously established, for generic Q, which transfer matrix sectors V ,k,m occur with a non-zero amplitude. At the non-generic Q 0 , some of the eigenvalues from different sectors will collide, as can be seen by explicit diagonalisation using the methods of section A.4 (or predicted analytically). The number of eigenvalues that collide at a given value Λ i is a natural guess for the largest possible rank r i of the Jordan block corresponding to the generalised amplitude (54). Trying the corresponding Ansatz for different values of i min (see above), we can examine whether the amplitudes are consistently determined (and non-zero) and arrive at a quantitative understanding of the Jordan block structure.

A.7.5 Extrapolation of amplitudes
Whichever method is used for obtaining the eigenvalue amplitudes A i on the cylinder, the results need to be extrapolated to infinite size (L → ∞) in order to be compared with predictions of conformal field theory. We shall see in Appendix B that in many cases the finite-size results are quite far away from their CFT limits, more often than not by a factor in the range 2-5. Obtaining reliable extrapolations is important, in particular because we only possess rather small sizes (L ≤ 11). Moreover, most amplitudes exhibit mod 2 parity effects, so that even and odd L must be treated separately.
Some guidance can be taken from our results for the Ising model where larger sizes can be obtained (L ≤ 16). We find the following method to provide quite accurate results. First we produce finite-size results with 2a L equal to or close to a constant, typically 1 2 . Since it turns out that some amplitudes vanish exactly, by symmetry reasons, for 2a L = 1 2 precisely, we often take 2a = (L − 1)/2 for odd L, and 2a = (L − 2)/2 for even L. These results are then corrected by a conformal factor, namely the powers of sin 2πa L appearing in (35). Finally, we extrapolate the corrected results, separately for odd and even L, by fitting all available data against a polynomial in 1/L (or in some cases leaving out the first or first few data points). The reliability of the extrapolation can then be estimated by comparing the two parities, which are expected to give identical results.

A.8 The split-attach algebra
The scalar product method described in section 4.3.2 and appendix A.2 requires us to find the left eigenvectors of T . One way to attain this is to take the transpose of the right eigenvectors of the transposed matrix T t .
The scalar product method described in section 4.3.2 and appendix A.2 requires us to find the left eigenvectors of T . One way to attain this is to take the transpose of the right eigenvectors of the transposed matrix T t .
In the join-detach algebra for the FK cluster model, T acts on basis states which are set partitions of L points. The basic building block of the sparse matrix factorisation scheme (10) is the action of the elementary join and detach operators, J i and D i , on these states. If this was presented in explicit matrix form, it would of course be trivial to take the transpose of the corresponding matrices. The point is, however, that when we use an iterative diagonalisation scheme (such as the Arnoldi method; see section A.7.1), it is inefficient to represent J i and D i as explicit matrices. Rather, we just need to provide their action on any state: for any in-state v in , the definition of J i (or D i ) in terms of set partitions provides the corresponding out-state v out = J i v in (or D i v in ). To treat the transpose in a similar setting means that we should instead answer the question: given v out , what are the possible v in that could lead to it (and with which transition weights)?
We therefore formally define the split operator S i = J t i and the attach operator A i = D t i as the transpose of the join and detach operator, respectively, and provide now their transition elements in terms of the basis states.

A.8.1 Without marked clusters
We first consider the analogue of the usual join-detach algebra, that is, without any marked clusters.
The operator S i performs a split between sites i and i + 1. In particular, S i is zero on states in which sites i and i+1 are not connected (i.e., in the same block of the set partition). If the two sites are connected, S i acts by breaking up the block containing i and i + 1 in all possible ways. To be precise, assume the block consist of k ≥ 2 points, ordered cyclically: i 1 < i 2 < · · · < i k < i 1 (the inequalities being considered modulo L), with i 1 = i and i 2 = i + 1. Place one "separator" between i 1 and i 2 . Place a second separator at the position of any of the above inequalities, i.e., between i and i +1 for any = 1, 2 . . . , k, with the subscripts considered mod k. Note that this includes the case where the two separators are at identical positions. Then break up the block by cutting it at the positions of the two separators. This produces the k possible out-states of S i , which each occur with transition weight 1. Note that the block is being broken into precisely two pieces in the k − 1 cases where the two separators are different, 11 whereas nothing happens (i.e., S i acts as the identity) in the remaining case where the separators coincide.
The operator A i attaches a singleton at site i to certain other blocks of the set partition, as we now describe. In particular, A i is zero on states in which site i is not a singleton. If i is an singleton, A i leaves the state unchanged with weight Q, and connects i to any other "visible" block in the set partition with weight 1. To determine whether a block is visible by i, draw the set partition as a hypergraph with vertex set 1, 2, . . . , L embedded in the half-infinite cylinder (concretely, in Figure 1 this would be the part of the cylinder situated to the left of the current time slice, with the vertices on the slice). If a block is adjacent to the face where i resides, the block is said to be visible. Another, less formal way to state the same, is that we connect i to any block that is not nested within another block, seen from the position of i.
It is straightforward (albeit somewhat laborious) to verify that the split-detach algebra thus defined satisfies the algebraic relations which are the transpose of those defining the join-detach algebra:

A.8.2 With marked clusters
We now provide the definition of the split-attach algebra in the representation with marked clusters having a fixed spin label. This is analogous to section A.2.
The split operator S i acts as zero, unless sites i and i + 1 belong to the same block (which can be unmarked or marked). Otherwise: • If i and i + 1 belong to the same unmarked block, S i leaves the block unchanged or splits it into two unmarked blocks, as before, with a total of k possibilities for a block of size k.
• If i and i + 1 belong to the same marked block, S i leaves the block unchanged or splits it into two blocks, of which one is unmarked and the other keeps the same mark as the original block. Both choices for which of the two blocks should carry the mark are realised.
The attach operator A i acts as zero, unless site i is an unmarked singleton. Otherwise: • With weight Q, the operator A i leaves the singleton block at i unchanged and unmarked.
• With weight 1, it marks site i by each possible colour of the mark which is used in none of the other blocks in the partition.
• With weight 1, it attaches site i to each of its visible blocks (which can be marked or unmarked).
The two first rules cover the cases where the corresponding detach operation is trivial, in the sence that D i would turn a (marked or unmarked) singleton into an unmarked singleton. The third rule covers the non-trivial case, where the corresponding detach operation acts on a (marked or unmarked) block of size ≥ 2.
With these modifications one can again check that the defining relations (115) are satisfied, but within the larger representation that allows clusters to be marked.

B Basic checks
In this appendix we check our general method for computing the amplitudes A i against a series of exact results. First we discuss the easy case of two-point functions and quantify the conformal content of the lattice spin operator. Next we compare our results for the s-channel spectrum of four-point functions against exact solutions for Q = 2, Q = 4, and the limit Q → 0.

B.1 Lattice observables and scaling fields
As discussed in the main text, one of the basic difficulties in our problem is that lattice observables correspond, in the scaling limit, to combinations of scaling fields weighed with appropriate powers of the cut-off. The question we want to investigate briefly here is how this might affect the measurement of amplitudes.
According to (8), the two-point function of the lattice spin operator becomes, in the geometrical formulation, proportional to the probability that two points belong to the same cluster. The leading dimension controlling the large distance behavior of this probability is h 1/2,0 and thus we expect, in the scaling limit, to have on the cylinder We have here used the same notation (33) as for the four-point function, so w 1 = ia, w 3 = i(a + x) + l, and w 13 = −ix − l. As before, we set ξ = e −2π(l+ix)/L , and the above dependence on ξ,ξ follows from a reasoning similar to the one leading to (34). We find numerically that only the sector F 0,−1 contributes, with As is usually the case in lattice models, the order operator on the lattice is not purely represented by a conformal field of weight (h 1/2,0 , h 1/2,0 ). In general, one expects instead that this field is only the first in a sum of the type where z,z are the complex coordinates corresponding to σ i , we have introduced the notation Φ r,s for conformal chiral fields with weight h rs , and stands for the lattice cutoff ( was set equal to unity so far).
In the limit where L, l → ∞ (that is, L/ , l/ → ∞ if L, l are measured in units of length) the contributions of the excited (higher) spin fields will scale away. For finite values of the parameters, they are unavoidably there. Note that This is larger for small values of m, i.e., smaller values of Q: therefore, the closer we will get to Q = 0, the faster these contributions will decrease with ( /L on the cylinder). We do not know how the Q-dependent amplitudes A e in (118) behave a priori, but we have studied numerically the two-point function in order to understand the amount of "mixing" of the order operator with the leading correction at weight h 3/2,0 . We have checked that this mixing is small and decreases significantly with increasing L.
As an example, we give here results for Q = 3/2-a case for which m is irrational, and all operators in the spectrum can be uniquely identified. The dimension of the order parameter is h 1/2,0 = 0.0583892, and by (116) the two-point function in the conformal limit expands as P scaling aa ∝ ξ h 1/2,0ξ h 1/2,0 1 + 0.116778(ξ +ξ) + 0.0136372ξξ + 0.0652078(ξ 2 +ξ 2 ) + . . . , where all numerical constants have been given to six-digit precision. The leading terms at momentum h −h = 1 and h −h = 2 are easily identified as the ground states of the corresponding momentum sectors. With no knowledge of the possible mixing with the term ξ h 3/2,0ξ h 3/2,0 one would look for the ξξ term in (120) in the lattice data as the contribution of the first excited state in the sector at vanishing momentum, but a careful analysis of the amplitudes as well as the scaling dimensions shows this is not what happens. On the lattice, the two point function contains, in addition to the terms in (120), a term in the bracket proportional to (ξξ) h 3/2,0 −h 1/2,0 . This follows from (118), which leads to and, after mapping on the cylinder, to The amplitude ratio A 1 /A 0 is non-universal, but, for a given lattice, it is a fixed quantity. The term L −4(h 3/2,0 −h 1/2,0 ) goes to zero as L goes to infinity, guaranteeing the disappearance of the unwanted terms in this limit. The question is how much this might affect the results for L finite. Measurement of the amplitudes-or rather the ratio with respect to the leading term-gives the following results: The extrapolation to L → ∞ was made via a polynomial of order 7 in 1/L, using all data points, since there are no discernable parity effects in this case. We observe that on one hand the expected ratios converge to their conformal values to a very good precision (4 or 5 digits). On the other, the one that is not expected-corresponding to the coupling to (h 3/2,0 , h 3/2,0 )-decreases fast, and converges to a value close to zero. 12 This example shows that one can obtain fine extrapolations, even though for the sizes considered in this paper (and the value Q = 3/2 taken here) the mixing of the field (h 3/2,0 , h 3/2,0 ) with the order parameter has a rather large relative coefficient of the order of 10 −2 .

B.2 The case Q = 2
We now consider the case Q = 2 (and later Q = 4) to check the consistency of our approach, and assess the quality of convergence of the amplitudes in the four-point functions. While in general, the geometrical probabilities cannot be expressed in terms of simple four-point functions in rational CFTs, the situation is better for these two values of Q. This has to do with the relationship [14] between correlation functions of spins in the Potts model and the geometrical objects, see eq. (9). We stress that for Q arbitrary, the left-hand sides of these equations are only formally defined: it is in fact the right-hand sides that give them a meaning. In the Q = 2 case the first relation (9a) reads simply G aaaa = P aaaa + P aabb + P abba + P abab , for Q = 2 . (123) On the other hand, recall from (7) that As discussed in A.6, for Q = 2, the order parameter Qδ σi,a − 1 coincides with the Ising spins S i = ±1, so this four-point function is nothing but the four-point function of the spin operator in the Ising model [63] σ(z 1 )σ(z 2 )σ(z 3 )σ(z 4 ) = 1 2 z 13 z 24 z 12 z 23 z 34 z 41 We also recall the structure constant C σσ = 1 2 , where denotes the energy operator; the latter appears as the non-trivial part of the scaling limit of S i S i+1 . The four-point function (125) involves two conformal blocks, corresponding to the fusion channels σσ ∼ 1 and σσ ∼ . Expanding the function G as in (28) gives then The normalised correlation function on the lattice should then be where s = sin 2πa L , and the extra terms arise from the conformal mapping as in (34).

B.2.1 Reduction from the generic case
We first checked-by computing the four-spin correlation function numerically for the ordinary (spin representation) Ising model on the cylinder-that the identity (123) holds exactly in finite size; details are given in section A.6. This computation is certainly a rather stringent test of the program that determines the P a1,a2,a3,a4 numerically. Note that the (large) set of eigenvalues that contribute to the geometrical correlations for generic values of Q reduces drastically-as expected-when we consider the combination (123). In algebraic terms, generically irreducible representations of the Jones algebra become reducible, and a complex structure of submodules of the relevant W j,z 2 develops. It turns out that only two simple quotients contribute to the Ising model: X 0,i and X 0,−1 , which are obtained as the the 'tops' of chains of modules according to the diagrams in Figure 9 (for a detailed discussion of the emergence of simple modules of the Jones algebra describing minimal models when q is a root of unity, see [27]). Figure 9: The simple modules X 0,i and X 0,−1 are the tops of these two diagrams representing the structure of the standard modules for q = e iπ/4 .
The continuum limit of the simple modules of the Jones algebra is fully expressed in terms of irreducible modules X r,s of the Virasoro algebra where the factor 2 indicates that, in fact, the representation splits into two isomorphic direct summands for the subalgebra generated by e i and u 2 , and X r,s denotes the irreducible Virasoro module with conformal weight h rs (and, e.g., character given by the well-known Rocha-Caridi formula [64]). We checked that the combination (123) receives contributions only from the tops in the diagrams in Figure 9 , which corresponds, in the continuum limit, to the required channels of the minimal Ising model. Note however that individual probabilities in (123), unlike their sum, do involve contributions outside the simple representations.

B.2.2 Numerical checks of the amplitudes
We now turn to the question of the convergence of the amplitudes measured on the lattice. Surprisingly, it seems this question has not been investigated much in the past (see [65] for early work). Since the Ising model is much easier to handle numerically than the general FK cluster model, we can study in this case much larger sizes (up to L max = 16), and explore in particular the nature of the convergence of the coefficients in the expansion (127).
We present in the following Figures 10-11 two different ways of handling the data. First, we consider the case where a is changed as L increases. For definiteness we consider 2a = L 2 when L is even and 2a = L−1 2 when L is odd. In both cases, sin 2πa L → 1 as L → ∞, but we should of course expect even/odd parity effects in L because of the different choices of 2a.
A priori the aim would be to compute the amplitude ratios, with respect to the ground state amplitude, corresponding to the five last terms in (127), namely the contributions to G aaaa of , (L −1 +L −1 ) , T +T , L −1 L −1 , and TT , respectively. However it turns out that (L −1 +L −1 ) and T +T correspond to eigenvalues that are exactly degenerate in finite size, so the corresponding contributions cannot be disentangeled. We thus list the amplitude ratios in (127) that we can access numerically along with their corresponding analytical predictions: L −1 L −1 : s 2 − 2s 4 + s 6 (129c) TT : The four panels of Figure 10 show the corresponding amplitudes ratios, where in each case we have divided the ratio by the expected analytical result just given. It is seen that in each case the result tends to 1 after extrapolation, thus confirming the whole analysis.
A number of remarks can be made about Figure 10. First, the extrapolations have to be carried out using some amount of common sense. In some cases the obvious non-monotonicity of the finite-size results (panel d) makes it clear that one should leave out the first few points from the fits. In other cases, (panel c) including all points in a high-order polynomial fit leads to the best results. Second, the comparison between fits through even and odd sizes appears to be a good measure for the "error bar" on the extrapolated value. Third, there is a tendency for the finite-size effects to grow as one goes to higher descendents. Thus, if we were limited to smaller sizes-such as L max = 11, as would be the case for the FK cluster model at generic values of Q-it should be expected that the last data point might very well be off the true extrapolated value by a factor of 2 or more.
It is also interesting to study what happens when a is fixed while L is increased. In this case, one does not expect the amplitude ratios to converge to the four-point CFT values reported in (129), since the latter require a, L 1 (in units of the lattice spacing) with fixed value of the ratio a L (and hence s fixed). To illustrate this, we consider for instance the case of 2a = 2 fixed, with L increasing-note that we do not expect parity effects using this protocol. To make the analysis comparable to the one above, we again divide the ratio by the CFT predictions (129). The results are shown in Figure 11, to be compared with the previous Figure 10. Clearly the convergence is no longer towards one. It is however remarkable that the numerical data still appear to converge to a finite value. For the energy operator and its first two descendents (panels a, b, c) we get fine extrapolations to the values 0.9292, 0.9289 and 0.9292 respectively, which are all compatible among themselves. For TT the extrapolations are more problematic, but still appear to converge to a finite value 0.53, which definitely appears to be different from the previous one.
We are not sure whether these values can be interpreted within CFT-neither whether the observed convergence is real or only apparent. Repeating the computations for 2a = 3, a similar extrapolation gives the consistent values 0.9686 and 0.9680 for the energy field and its first descendent, whereas TT gives another value 0.71.

B.3 The case Q = 4
As pointed out in [1], the Potts model at Q = 4 is particularly interesting, since in this case all geometrical correlation functions can be expressed in terms of spin correlation functions. This provides analytical results for all the P abcd , and thus provides a simple benchmark against which to test our approach. There is however a drawback to this idea: numerics for the Potts model at Q = 4 are known to be affected by stronger than usual corrections to scaling due to the presence of a marginal operator. As a result, convergence to the (c) (d) Figure 11: Same as Figure 10, but for a fixed distance 2a = 2 between each pair of operators. The panels are as in Figure 10. The extrapolations, shown as blue curves, now uses sizes L of both parities.
expected amplitudes appears to be slower than for other values of Q, especially small ones.
In any case, we start by observing that the Potts model partition function, can be expressed in the continuum limit as the sum [7,66], where the generating functions are obtained by taking the limit m → ∞ of our general result (46) In contrast with the case of generic Q [7], we see that only a very small subset of the F j,e 2iπp/M generating functions contribute. This corresponds to the fact-which we checked explicitly in finite size, and which has a simple representation theoretic interpretation; see below-that of all the usual Jones algebra modules, only those corresponding to z 2 = ±1 play a role at Q = 4.

B.3.1 Exact results via the Ashkin-Teller model
After these preliminaries, we now consider four-point functions per se. In [1], some of the probabilities at Q = 4 are given based on results in [14], combined with the careful reading of a paper by Al.B. Zamolodchikov [67]. Calling τ 1 , τ 2 the two spins of the Ashkin-Teller model, we introduce, following [67], the four quantities note that G here must not be confused with G in (23). It is then expected that 13 The quantities on the right-hand side are then calculated in [67]. We introduce first where z = z12z34 z13z24 as usual. We shall also need the Jacobi theta functions These two functions are used to define q(z) implicitly via We then set Finally we choose The claim is then that the z,z dependent part of the four-point function (the function G of section 4.1) is given by corresponding combinations of the functions G and R i . So for instance with What we must do then is to expand the combinations (138) in powers of z. This gives expressions of the form The point is that the set of weights (∆,∆) is the set of conformal weights appearing in the s-channel expansion of the four point function. (It is not necessarily the set of primary fields, as the exact result does not distinguish between primaries and descendents. Hence our use of the notation ∆ is slightly abusive when compared with the general case.) This must be the same set as the set of conformal weights contributing to the lattice observable. Hence, we must compare the sets occurring in (138) with the sets determined directly. These sets can be determined numerically, as for Q generic, and are found to be: Here the quotients of modules have a direct meaning in terms of subtracting eigenvalues common to various sectors in finite size. Like for other non-generic values of q, the modules contributing to the probabilities are no longer irreducible, and the Q = 4 results correspond to taking the simple irreducible tops. Indeed we find in general that Spec W 2k+1,1 ⊂ Spec W 2k,1 for k = 1, 2, . . .. We have moreover established numerically that Spec (P abab −P abba ) does not contain Spec W 4,−1 , which would otherwise have been a viable candidate in view of symmetries and the results for generic Q. For these reasons we conjecture that the complete result generalising (143) should in fact read: where we recall that ξ = e −2π(l+ix)/L ,ξ = e −2π(l−ix)/L .
The first and third terms in this expression are associated with the largest eigenvalue in F 0,−1 at zero momentum (this is denoted V 100 in section A.4). The second term also belongs to F 0,−1 , but at non-zero momentum (namely V 101 ). Finally, the fourth and fifth terms are the leading ones in F 2,1 (in respective momentum sectors V 200 and V 202 ). To associate these terms to definite eigenvalues, for each size L, requires going through the detailed analysis of section A.5. We moreover note that one has to be careful in the comparison, since on the lattice we do not individually observe the amplitudes of conjugate terms, such as the two terms in (ξξ) 1/16 (ξ +ξ), since they correspond to the same eigenvalue. Therefore, to make the comparison we must divide the numerically observed amplitude by two in such cases (cf. remark 6)-this has been done tacitly below and in subsequent similar cases. This leads to the following results for the second to fifth terms in (147): We see that the first two columns come out satisfactorily, and the ratios can be rather convincingly extrapolated to a number close to one, despite of the small number of sizes. For the higher-order terms (the last two columns of the table) the convergence is definitely slower, but still goes in the right direction. We have already seen (witness Figure 10 in the Ising case) that ratios of the order of 2 at small sizes (L ≤ 9 here) are not uncommon, and it is quite plausible that also these columns could be extrapolated to one, provided one could obtain a few more sizes. In conclusion, the test of P aaaa is fully compatible with the exact results.
The first term is the leading eigenvalue in F 0,−1 , corresponding in finite size to the largest eigenvalue in V 100 . The second and third terms are the two leading eigenvalues inF 0,1 , corresponding to the two largest eigenvalues in V 000 . The comparison with numerics-presented in the by now familiar format-here comes out as: It again appears convincing that the ratios will converge to one, i.e., that the numerical results confirm the analytical ones.
B.3.4 The combination P abab − P abba As a last example we discuss antisymmetric combination P abab − P abba of the probabilities of having two "long" propagating clusters, as shown in Figure 2.
The powers of q,q on the right-hand side are expected to be exactly the exponents appearing in the sum (142). Using the foregoing discussion and the expressions of R 1 , R 3 , we can on the other hand identify these exponents by performing a direct expansion whose first few terms are: This plays the role of the function G in our general discussion (34). It is easy to check that the two sets agree. Moreover, we have also checked that there is no gap in the spectrum, that is, all exponents predicted from the generating function (143) are indeed present with non-zero coupling constant.
To check the numerical values of the amplitudes, we need as usual to map the four-point function on the cylinder. We find the first few terms P abab − P abba ∝ (ξξ) 9/16 (ξ +ξ) + 9 4 − 9 sin 2 2πa L 4 (ξξ) 25 The first term corresponds to the lowest eigenvalue in the F 2,−1 sector with conformal spin h −h = 1, originating from V 211 in finite size. The second and third term correspond to the leading eigenvalue with h−h = 0 (resp. h−h = 2), originating from V 210 (resp. V 212 )-note that since the scaling dimension h+h = 25 8 is the same for these two terms their ordering in the finite-size data cannot be guessed straightaway, and one has to make use of the lattice momentum to reveal the correct conformal spin. Finally, the fourth term corresponds to the lowest eigenvalue in the sector with h −h = 3-we notice that this state is not present for L = 5, since this size is yet too small to accommodate the high value of the momentum. The numerical results for the last three terms being discussed (taking ratios with respect to the first term) now run as follows:

L
(ξξ) 25 . The exponents appearing in the various Ashkin-Teller correlators are thus in agreement with the spectra conjectured in [1] G : S Z,2Z, after the switch of r, s labels mentioned earlier in remark 5.

B.4 The case Q = 0
The case Q = 0 (or rather the limit Q → 0) is connected to the combinatorics of spanning trees and forests (see, e.g., [68][69][70]). It is interesting in the present context for two reasons. On the one hand, it gives rise to Jordan cells in the transfer matrix, and thus provides an opportunity to study their effect on the determination of amplitudes. On the other hand, it also turns out that the combination P abab − P abba can be explicitly calculated in this case, providing another non-trivial check of our method, this time in the region of small values of Q. While for Q generic, the two-point function of the Potts spin operator is proportional to the probability that the two points belong to the same cluster [see (8)], the limit Q → 0 needs to be handled with care, since the partition function itself vanishes. At leading order, the only terms left in the partition function are spanning trees: a naive definition-requiring that the two points still belong to the same cluster-would then lead to a trivially constant two-point function. To obtain non-trivial results it is better to change the normalisation by one factor of Q, that is, to redefine the partition function as the number of spanning trees. A natural redefinition of the two-point function will be given below; it follows by using the equivalence of the Q → 0 limit with symplectic fermions [71,72] and the theory of spanning trees.
The simplest, in fact, is to start by discussing the combination P abab − P abba . Indeed, consider four points labelled 1, 2, 3, 4 in the plane. As usual, we consider the square lattice in concrete calculations. We denote by N 13,24 (resp. N 14,23 ) the numbers of configurations of spanning trees where one tree connects points 1, 3 and a different tree connects points 2, 4 (resp. one tree connects points 1, 4 and a different ones points 2, 3). It is possible to show, by generalising Kirchoff's original discussion [73] along the lines of [74], that N 13,24 − N 14,23 = Det ∆ (12) (34) , where the right-hand side is the determinant of the lattice Laplacian after having removed the lines corresponding to points 1, 2 and columns corresponding to points 3, 4. Here, the Laplacian is defined as follows.
It is a matrix denoted ∆ whose linear size is the number of vertices of the graph, here a planar lattice, that we suppose loopless (i.e., no vertex is connected to itself) for simplicity. The diagonal elements ∆ ii are the number of edges incident on vertex i. The off-diagonal elements ∆ ij are equal to minus the number of edges connecting the vertices i and j. We also denote by ∆ (kl) the minor of ∆ obtained by erasing row k and column l. We recall that det ∆ = 0, since by definition the sum of all rows (or the sum of all columns) is zero. Moreover, by the Kirchhoff matrix-tree theorem [73], det ∆ (kl) is equal to the number of spanning trees on the graph. If we now go back to the Q → 0 limit of the Potts model and the definition (6) of the probabilities P a1,a2,a3,a4 in terms of the cluster expansion, we see that (155) is the leading contribution to P abab − P abba as Q → 0: in this limit indeed, any non-connected additional cluster gets cancelled by a power of Q, and so do cycles within the clusters (so each cluster is a tree indeed). Meanwhile, (155) can be calculated in the continuum limit, which is simply described by a pair of symplectic fermions (θ + , θ − ) with Euclidean action subject to summation over repeated indices, and with d +− = 1, d −+ = −1. The quantity in (155) is then nothing but θ + (1)θ + (2)θ − (3)θ − (4) and gives 15 [75,76] θ + (1)θ + (2)θ − (3)θ − (4) = ln z 14 z 23 z 13 z 24 2 . (157)

B.4.1 Two-point function
Before studying P abab − P abba , we can extract from this a useful definition of the two-point function as well. Indeed imagine sending point 1 to 3, and point 2 to 4. In this case, the configurations in N 14,23 become negligible, while those in N 13,24 describe a situation where 1 = 3 belongs to a tree and 2 = 4 to a different tree. We use this to define a two-point function: where N 1,2 is the number of configurations of spanning forests with two trees only, with points 1 and 2 belonging to different trees, while N 1 denotes the number of spanning trees. In the continuum limit we therefore have from (157) g 2 (z,z) = 2 ln |z/ | 2 , where is a short-distance cutoff. Note that cannot be eliminated by multiplicative renormalisation, as is usually done in CFT. We could, alternatively, set = 1, and say that the two-point function takes the form (159) up to a non-universal additive constant.
On the cylinder, we expect the following behavior in the conformal limit [77]: g 2 (w,w) = 2 ln L π sinh πw L + 2 ln L π sinh πw L .
We shall measure w, L in units of the lattice spacing. This means that will be a (non-universal) numerical constant of order unity. We now expand g 2 to see how it connects with the results from the transfer matrix: It is seen that g 2 consists of a linear term and a sum of exponentials. All the amplitudes can be compared with lattice calculations: the sum of exponentials arises from a sum over eigenvalues of the transfer matrix as usual, while the linear term arises because of the presence of a Jordan cell of rank two. Write (161) as where Λ 0 denotes the ground-state eigenvalue. We get from (161) that A = 4π L while B 0 = 4 ln L 2π , B 1 = −4, and B 2 = −2 (where the last term comes from n = 2 in the sum). However, if we change the normalisation of (162) so that A = 1 L , we get instead: As usual we can confront this with the numerical computations. We find that A = 1 L exactly in finite size, which motivates the above choice of normalisation. Moreover, g(1, 3) is found to couple only to a very small set of eigenvalues: apart from the rank-two Jordan cell parameterised by (B 0 , A), there are only two (resp. three) simple eigenvalues for L = 5 and 2a = 2 (resp. L = 7 and 2a = 3). Presented in the usual table form we find: Note that we have here set = 1 in the CFT result for B 0 in order to get an order of magnitude estimate. Since this term is not expected to be universal, the agreement for B 0 is here seen to be reasonable. The agreement for the other terms (B 1 , B 2 , B 3 ) is seen to be very good, even at these small sizes.
We now comment more generally on the structure of the two-point function in terms of the transfer matrix. The linear term in (161) indicates the presence of a Jordan cell of rank two (since we get a term linear in l), and that the two corresponding (pseudo) eigenvalues must be in the ground-state sector since there is no exponential decay. The sum of exponentials involves only integer conformal weights, and corresponds simply to coupling to descendants of the identity. Remarkably, note that we couple only to purely chiral or purely antichiral fields. This explains the scarcity of non-zero amplitudes observed numerically. It is also useful here to recall the structure of some of the modules for Q = 0. Since q is a root of unity, the W j,z 2 =e 2iK are reducible, and have a structure of submodules that depends on the values of j, z. For instance we get the structure shown in Figure 12.
Meanwhile, we can read the dimension content from the general formulae: Let us consider in more detail F 1,1 . This contains the exponents (h,h), here ordered by the total value of the scaling dimension h +h: We now go back to the combination P abab − P abba , which we write in terms of the variable z ≡ z12z34 z13z24 : and expanding this we get P abab − P abba ∝ ∞ n=1 z n +z n n .
Note that, if we exchange points 1 and 2, we have z = z12z34 z13z24 → z21z34 z23z14 = z z−1 . This corresponds to sending 1 − z → 1 1−z , and thus to changing the sign of P abab − P abba , as required from the geometrical interpretation. This does not mean that only odd spins appear in the expression (28), because in this case there are large degeneracies (in fact, only a few of the terms correspond to primary fields).
To compare with the numerical work, we rewrite (168) as P abab − P abba ∝ ∞ n=1 C n ξ n + h.c.
and choose the normalisation C 1 = 1, so that the amplitude ratios can be read off directly from the following coefficients. The first few then read explicitly As above, the numerical work was done for L = 5 and 2a = 2 (resp. L = 7 and 2a = 3). The contributing eigenvalues are then found to be exactly the same as for the two-point function g 2 (1, 2), corresponding to the terms n = 1, 2 (resp. n = 1, 2, 3) in (169). There is no longer any Jordan cell (the term n = 0). Note also that in both cases, we expect at all orders to have only chiral or antichiral contributions. The numerical results compared with the conformal predictions run as follows: Once again the agreement is very reasonable, given the rather small sizes. A peculiarity of this four-point function-as well as of the two-point function g 2 considered above-is that only fields with h = 0 orh = 0 appear. Since the spin h −h is limited to values 0, 1, . . . , L − 1 for a finite size L, this means only a very small number of eigenvalues contribute to the correlation function at finite L. This is quite different from the usual case, where having h −h fixed does not preclude eigenvalues with larger h +h from contributing. Of course, it could be that some eigenvalues would contribute in finite L, although their amplitude would tend to zero when L increases, because they would not be present in the continuum limit. But this does not seem to be the case here-a fact that can probably be proven exactly in finite size using the Laplacian on the cylinder.

B.4.3 The combination P aabb − P abba
It is also interesting to observe that, under crossing z 2 ↔ z 3 , ie z → 1 z , the combination P abab − P abba that we have just discussed becomes P aabb − P abba . We have studied the latter combination independently. From our analytical result (166) we now find P aabb − P abba ∝ ln z − ln(1 − z) + lnz − ln(1 −z) .
It is easy to ascertain that the correct sectors are observed in the numerical study. Indeed we observe now again a Jordan cell on the ground state, plus the same 2 (resp. 3) simple eigenvalues as before for the case L = 5 and 2a = 2 (resp. L = 7 and 2a = 3). The indecomposability parameter, i.e., the coefficient multiplying l in the correlation function, is found to be A = 1 L exactly in finite size, like it was for the twopoint function. Let us write the remaining terms, concerning the simple eigenvalues, as ∞ n=1 D n ξ n + h.c. The comparison with numerics can then be summarised as follows: It will not have escaped the reader's attention that the ratios are exactly the same as for the previously considered case P abab − P abba . This presumably extends to arbitrary L and implies that the symmetry z → 1 z is respected exactly on the lattice-something that we again suspect can be proven by a careful study of the Laplacian on the cylinder. Per se, this close relationship between the foregoing two four-point functions is pretty remarkable, since after all, in the lattice study of P abab we have two long clusters, while in P aabb we have two short ones.