Four-point boundary connectivities in critical two-dimensional percolation from conformal invariance

We conjecture an exact form for an universal ratio of four-point cluster connectivities in the critical two-dimensional Q-color Potts model. We also provide analogous results for the limit Q → 1 that corresponds to percolation where the observable has a logarithmic singularity. Our conjectures are tested against Monte Carlo simulations showing excellent agreement for Q = 1, 2, 3.


Introduction
The study of the geometry of random two-dimensional fractals has revealed the emergence of a profound mathematical connection between probability theory and stochastic processes [1][2][3][4] on one hand and quantum field theory together with conformal symmetry on the other [5][6][7][8][9][10][11]. Historically, a number of exact results were derived for the fractal dimensions of two-dimensional critical clusters in basic models of statistical mechanics such as percolation or Ising, built on the seminal contribution [12] and the so-called Coulomb Gas approach [13]. A deeper insight on how conformal invariance could be relevant to describe geometrical observables followed after [14] when J. Cardy derived, with methods borrowed from (boundary) conformal field theory [15,16], an exact formula for the probability that at least a cluster should span the two horizontal sides of a rectangle in critical percolation [17]. The success of this approach suggested that a geometrical problem, such as critical percolation, could be solvable in two dimensions due to the infinite-dimensional nature of the conformal group [12].
However, at the same time, it was noticed how the conformal algebra associated to geometrical phase transitions could be more subtle [18]. In particular as a conformal field theory, critical percolation should have vanishing central charge (denoted by c) since its partition function does not depend on finite size effects [19,20]. However rightly at c = 0 the stress-energy tensor is a null field and the field theory if not trivial, cannot be unitary. Absence of unitarity has serious consequences on the Operator Product Expansion (OPE) -1 -

JHEP12(2018)131
and ultimately produces logarithmic singularities in the four-point functions [21,22]. Later, it was conjectured in [23] that the OPE of two chiral fields with scaling dimension h = 0 at c = 0 should have the following expansion (z ∈ C) lim z→0 φ(z)φ(0) = 1 z 2h 1 + 2h b z 2 (t(0) + log(z)T (0)) + . . . , (1.1) where T (z) is the null stress energy tensor and t(z) was called its logarithmic partner. The parameter b in eq. (1.1), termed the indecomposability parameter, is a universal number characterizing the c = 0 Conformal Field Theory (CFT) [23,24]. Its name, in particular, stems from the fact that the fields T and t span a Jordan cell of dimension two which makes the conformal dilation operator non-diagonalizable. CFTs that are built upon indecomposable representations of the Virasoro algebra are called logarithmic [25]. They are supposed to be ubiquitous in the study of random clusters and disordered two-dimensional systems [26,27]. For detailed studies in higher dimensions, see also [28]. During the last decade a lot of effort has been put in the classification of logarithmic CFTs with special success on finite domains; see [29] and references therein. However not many exact correlation functions have been explicitly calculated and tested in statistical mechanics. Important exceptions are G. Watts result [30] and other generalizations of Cardy crossing formula on polygonal domains, such as hexagons or octagons [31][32][33][34][35]. In particular, logarithmic singularities in crossing probabilities are hidden into higher-point correlation functions [35,36]; for instance the six-point functions of the field φ 1,2 in the notations of eq. (3.3). Such a field has vanishing scaling dimensions at c = 0 and its fourpoint function cannot be logarithmic [14], cf. eq. (1.1). We should also mention that the study of logarithmic conformal field theories in the bulk is considerably harder than on a finite geometry, due to the constraints of crossing symmetry. Recent developments for three [37][38][39] and four-point functions [40], are based on a conformal bootstrap approach to Liouville theory for c < 1.
In this paper we complement those existing results, introducing and exactly determining a geometrical observable that explicitly shows logarithmic behavior at criticality. We focus on the Q-color Potts model [41] on a bounded domain and construct a ratio between four-point cluster connectivities, see figure 1. We then follow closely [14] and symmetry arguments to obtain a fully analytic expression for such a quantity in terms of Virasoro conformal blocks [12]. The main technical assumption of our approach is that the functions in eq. (2.3) solve a third order differential equation that is associated to a null vector for a field with non-zero scaling dimension. Null vector decoupling for a non-unitary CFT is not granted a priori and requires care. However for the specific case considered here, it was proven self-consistent [23] at c = 0 and led to the conjecture b = −5/8 in eq. (1.1) for boundary percolation. The same value for b was later re-derived in [42] by algebraic means and it is now accepted as a universal number characterizing this universality class [27]. Our results in section 4 for boundary percolation could be then considered both a nontrivial application and a further test of the ideas in [23]; in particular the OPE (1.1). As noticed in [27,43], our geometrical observable is also logarithmic at Q = 2, i.e. in the (extended [44]) Ising model. The latter was already analyzed in [45] by the same authors Figure 1. Four points x 1 , x 2 , x 3 and x 4 are marked on the boundary of a simply connected domain D embedded into a regular two-dimensional lattice. The function P (12) (34) is the probability that x 1 and x 2 are connected into one FK cluster while x 3 and x 4 are connected into a different FK cluster. Analogously P (14) (23) is the probability that x 1 and x 4 are connected into one FK cluster while x 2 and x 3 are connected into a different FK cluster. Finally P (1234) is the probability that all the points belong to the same FK cluster. The non-normalized probability measure of any configurations (i.e. of any random graphs G) is given in eq. (2.2).
but we recast it here in a more general context. We also provide numerical checks of all our results through high-precision Monte Carlo simulations and further extend the analysis in [45] to the three-color and four-color Potts model. These cases, although conceptually analogous and technically simpler than Q = 1 and Q = 2 have been not considered before. The outline of the rest of the paper is as follows. In section 2 we introduce the Qcolor Potts model and the geometrical observable R. In section 3 we show how this ratio of cluster connectivities can be obtained from conformal invariance and derive explicit analytic expressions in section 4. The comparison of the CFT predictions against Monte Carlo simulations is addressed in section 5. After the conclusions, three technical appendices complete the paper.

Four-point boundary connectivities in the Q-color Potts model
The Q-color Potts model [41] is defined by the Hamiltonian (J > 0) where the spin variable s(x) takes only positive integer values up to Q, and the sum extends over next-neighboring sites on a certain bounded domain D embedded into a twodimensional regular lattice. The boundary conditions for the spin are free. The Potts partition function Z(Q) = {s(x)} e −H Q admits a well known graph expansion, the socalled Fortuin and Kasteleyn [46,47] representation. Let p = 1 − e −J the probability of drawing a bond between two next neighboring lattice sites in D, then it turns out, up to a multiplicative constant,

JHEP12(2018)131
In eq. (2.2) above, n b (resp.n b ) is the number of occupied (resp. empty) bonds in the domain D. Connected components, including isolated points, inside a graph G are called clusters (FK clusters). Each graph contains N c clusters in which the Potts spins are forced to have the same color, hence the factor Q Nc . When Q = 2, eq. (2.2) is the high-temperature expansion of the Ising model. Although the partition function Z(Q) can be defined for any complex Q, in this paper we will consider only positive integer values including however Q = 1 that corresponds to the percolation problem. In particular, we will be interested in determining boundary connectivities in the Q-color Potts model. Connectivities in the Q-color Potts model are probabilities that a certain set of points marked by x 1 , . . . , x n are partitioned into FK clusters. A non-normalized probability measure for the allowed graph configurations is given by m(G) = p n b (1 − p)n b Q Nc according to eq. (2.2). The normalized probability measure for the graphs would be of course Z −1 m(G) as Z(Q) = 1 only at Q = 1; the normalization factor is however not essential here since only ratios of connectivities will be considered. Moreover, if we focus on configurations in which n points are on the boundary, it can be shown [48] that the number of linearly independent connectivities is also the number of non-crossing non-singleton partitions of a set of n elements (also known as Riordan numbers [49]). In geometrical terms, linear independent boundary connectivities correspond indeed [48] to configurations where none of the boundary points belongs to an isolated FK cluster (i.e. it is disconnected from all the other points). In particular if n = 4 there are only three linearly independent fourpoint boundary connectivities that are schematically represented in figure 1. With obvious notations such four-point functions will be denoted by P (12) (34) , P (14) (23) and P (1234) , see again figure 1. Notice that, due to the planarity of the domain D, we have P (13)(24) = 0. This important simplification does not occur when the four-points are in the bulk; in such a case the number of linear independent four-point connectivities is four.
For 1 ≤ Q ≤ 4, the Potts model undergoes a second order ferromagnetic phase transition for a critical reduced inverse temperature J = J c . The ferromagnetic phase transition is instead of the first order for Q > 4. In geometrical terms, at J > J c (p > p c ) there is a finite probability that any bulk point will be connected to the boundary of D. Such a probability vanishes as a power law as J → J + c with a critical exponent that coincides with the one of the one-point function of the order parameter; for instance [50].
At J = J c , and in the scaling limit when the mesh of the lattice is sent to zero, connectivities, although strictly speaking vanishing, are conjectured to be conformally covariant [17]. It is then useful to define the dimensionless conformal invariant ratio R = P (14)(23) P (14)(23) + P (12)(34) + P (1234) .
Let us aso clarify the statement that R is a conformal invariant quantity. If x j ≡ (s j , t j ) belongs to D we introduce complex coordinates w j = s j + it j ∈ C. Then the bounded domain D (for instance the unit disk) can be mapped conformally through the mapping z = z(w) into the upper half plane. When the points w j are on the boundary of D, they are mapped on the real axis and we can always choose z 1 < z 2 < z 3 < z 4 . Conformal invariance -4 -

JHEP12(2018)131
implies that R calculated on the upper half plane is only a function of the anharmonic ratio η = z 21 z 43 z 31 z 42 (2.4) where z ij = z i − z j . We can then determine R on the original domain D simply replacing z j = z(w j ) (j = 1, . . . , 4) into the expression for R(η). In the next sections, we will obtain exactly R. The conjecture for eq. (2.3) is valid for all the integer values Q = 1, 2, 3, 4 and will be eventually tested against Monte Carlo simulations. Anticipating the content of the remaining sections, the reader will find a plot of our theoretical predictions for the universal ratio in eq. (2.3) in figure 4.
It is important to observe that, although their leading short-distance singularities are the same, four-point boundary connectivities cannot be obtained from the knowledge of the boundary correlation functions of the Potts order parameter. A paradigmatic example [48] is Q = 2, where the unique boundary four-point function of the spin is a linear combination of the three connectivities in figure 1. Three-point boundary connectivities for the Q-color Potts model (and in particular at Q = 3) could be determined instead from the known solutions [51,52] of the boundary bootstrap equations for the Minimal Moldels.

Duality and conformal symmetry
Boundary-condition-changing operators and duality. The quantum field theory that describes the critical large-distance fluctuations of the two-dimensional Q-color Potts model is a conformal field theory (CFT) [53] with central charge In principle p ∈ R, however since Q is integer and 1 ≤ Q ≤ 4, we only consider the cases p = 2 (percolation), p = 3 (Ising), p = 5 (three-color Potts model) and p → ∞ (four-color Potts model). Notice that at Q = 1 the central charge in eq. (3.1) is zero, a well-known circumstance of critical percolation that in particular implies that in such a case the CFT is non-unitary (the only unitary CFT with zero central charge is indeed trivial). This lack of unitary reflects itself into the presence, as we will discuss, of logarithmic singularities [18,[21][22][23] in the four-point connectivities. Let us now briefly review some basics of (boundary) CFTs in two dimensions [15,16]. When a scaling field is inserted at the boundary of a planar domain D, its scaling dimensions are eigenvalues of the operator L 0 [15]. The scaling dimensions of a primary field (for a definition see [54]) φ r,s sitting at the boundary are then given by Figure 2. Four boundary partition functions (from left to right) Z * αβαβ , Z * αβγδ , Z * αβαγ and Z * αβγα . Different colors correspond to different fixed boundary conditions for the (dual) Potts spin. The following relation holds Z * αβαβ + Z * αβγδ = Z * αβαγ + Z * αβγα , showing that only three of them are linearly independent. They are related [48,55] through a duality transformation to the three linearly independent connectivities in figure 1, calculated with free boundary conditions. At the critical point and in the scaling limit, they are also proportional, to the four-point functions of the bcc φ αβ . and it turns out that for our purposes r, s are positive integers. In such a case the fields φ r,s are also dubbed degenerate [12] and their correlation functions satisfy partial differential equations of degree rs. The operator content of the CFT depends on the boundary conditions. For the Q-color Potts model natural boundary conditions for the spin variable on D are either free or fixed to a definite color α = 1, . . . , Q. Remarkably, conformal symmetry is however also compatible with inhomogeneous boundary conditions that are associated to the insertion of scaling fields at the boundary called boundary-condition-changing (bcc) operators [16]. For instance if the value of the spin on the boundary switches from α to β = α nearby x, such a discontinuity in the boundary conditions is realized, in the scaling limit, by the insertion of a scaling field φ αβ (x). In particular, the bcc operator φ αβ can be identified with a field φ r,s , whose scaling dimensions are given in eq. (3.3); we will discuss which field in the next subsection. Before we remind, as pointed out first in [14], how correlation functions of bcc operators can be related to connectivities in the Q-color Potts model. The argument exploits the duality transformation of the Potts partition function on a planar domain D; in particular the mapping also requires a transformation of the lattice. A duality transformation [48] relates partition functions on the dual lattice with fixed boundary conditions and in the ordered phase (J * ≥ J c ) to connectivities calculated with free boundary conditions on the original lattice in the disordered phase (J ≤ J c ).For a definition of the dual of an FK cluster and the lattice-dependent mapping J * (J), we refer to sectionII of [41]. There are four different dual partition functions, involving four boundary points, that we could consider for J * ≥ J c . They are denoted by Z * αβαβ , Z * αβγδ , Z * αβαγ and Z * αβγβ and represented in figure 2. Regions coloured differently at the boundary of D correspond to regions with different fixed boundary conditions for the dual spins. It is also important to remark that the following relation [48,55] holds: showing that actually only three of them are linearly independent. Among the duality relations [48,55], involving the four-partition functions above and the three linearly independent connectivities with free boundary conditions in -6 -JHEP12(2018)131 figure 1, we will only need (see the discussion below eq. (3.11)) the following In eq. (3.4), A is a lattice-dependent normalization constant that however does not depend on x 1 , . . . , x 4 . eq. (3.4) was obtained in [48] directly in the scaling limit, but it can also be understood as follows. Perform an FK graph expansion of the partition function Z * αβαβ , then the dual FK clusters cannot connect regions where the spins are fixed to have different colors at the boundary. We can distinguish three cases: • Dual graph configurations contain at least a dual cluster connecting the two regions with boundary conditions β (there is an horizontal dual crossing). Applying a duality transformation these configurations are in one-to-one correspondence with the ones that contribute to P (12) (34) .
• Dual graph configurations contain at least a dual cluster connecting the two regions with boundary conditions α (there is a vertical dual crossing). Applying a duality transformation these configurations are in one-to-one correspondence with the ones that contribute to P (14)(23) , see figure 3a.
• Dual graph configurations do not contain any cluster that connects regions on the boundary with the same color (there are no dual crossings). Applying a duality transformation these configurations are in one-to-one correspondence with P (1234) .
Notice that cannot be simultaneous horizontal and vertical crossings; this possibility was instead investigated in [30]. Summing over the three possibilities we obtain eq. (3.4). The partition function Z * αβαβ (x 1 , x 2 , x 3 , x 4 ) is in turn proportional [14,16] in the scaling limit to the four-point function with A constant and P t ≡ P (12)(34) + P (14)(23) + P (1234) . Eq. (3.5) will be used to derive the expansion in conformal blocks in eq. (3.13) for the universal ratio R in eq. (2.3).
Conformal blocks and the universal ratio R. Following [14], we identify the bcc operator φ αβ with the field φ 1,3 , whose scaling dimensions are given (cf. eq. (3.3)) by The identification holds for any integer 1 ≤ Q ≤ 4. Let us now consider the boundary fourpoint function of the field φ 1,3 ; as discussed in section 3 we work on the upper half plane, denoted hereafter by H. The four points are then ordered on the real axis and chosen such that z 1 < z 2 < z 3 < z 4 . Exploiting global conformal symmetry on the upper half plane (i.e. SL(2, R) Moebius transformations), the four-point function of φ 1,3 can be written as The function G(η) solves [12] an Ordinary Differential Equation (ODE) of degree 3 that can be obtained by the condition of decoupling of the null-vector at level three in the Verma module of φ 1,3 . Although null vector decoupling for a non-unitary CFT requires care, for the specific case of operators on the first column of the Kac table, this appears consistent with all available predictions based on chiral logarithmic CFTs [23,24]. It has been also employed in [56] for an analogous geometrical problem in percolation. The derivation of such a differential equation is standard, the reader can consult for instance [57]. It turns out and h depends on p as in eq. (3.6) whereas p is related to the central charge and Q by eq. (3.1). The equation above is of Fuchsian type with regular singular points in η = 0, 1 and ∞, therefore we can write its Frobenius series near η = 0 as G ρ (η) = η ρ ∞ k=0 a k η k where conventionally we set a 0 = 1. The series has radius of convergence |η| < 1 however this is enough for our purposes since by a global conformal transformation we can set z 1 = 0, z 2 = η, z 3 = 1 and z 4 = ∞, and therefore recalling that the four points are ordered on the boundary η ∈ [0, 1]. Notice also that eq. (3.8) is symmetric under the transformation η → (1 − η). The exponent ρ solves the indicial equation and the three roots coincide with the scaling dimensions h 1,1 , h 1,3 and h 1,5 given in eq. (3.3). This is of course expected since [12] the roots of eq. (3.9) are the scaling dimensions of the leading singularities produced in the OPE lim z 1 →z 2 φ 1,3 (z 1 )φ 1,3 (z 2 ) as it can be seen from eq. (3.7); in particular h 1,5 = 3h + 1. When the differences between the roots of the indicial equation are not integer and always for the largest root, the Frobenius series are directly the Virasoro conformal blocks of the CFT. Denoting by F c ρ the four-point conformal block where we fixed the external legs to be the field φ 1,3 (of dimension h(c)) and the internal field has scaling dimension ρ, we have (3.10) To lighten the notation we did not show the dependence of G ρ from c (alias h). The coefficients in the series expansion in eq. (3.10) can be calculated from Zamolodchikov recursive formula [58], see appendix B, and provide a further verification of the solution of eq. (3.8). Notice that, compared to [58], we have omitted the pre-factor η −2h in the definition of the conformal block in eq. (3.10).
When the indicial equation have roots ρ 1 and ρ 2 , such that ρ 1 − ρ 2 = ν is a positive integer but no solution is logarithmic the Frobenius series may not necessarily coincide with the conformal blocks. The two power series can indeed mix at order ν.
Since the conformal blocks are directly the contribution in eq. (3.7) of the OPE channels, we will propose an identification of the universal ratio R in eq. (2.3) in terms of them. The identification is based on the following observations.
Obs. 1: The three linearly independent connectivities P (12) (34) , P (14) (23) and P (1234) are in the scaling limit proportional to the four-point function of φ 1,3 . Therefore, once the common prefactor in eq. (3.7) has been factored out, they are linear combinations of the three conformal blocks.
Obs. 2: The field φ 1,2k+1 when inserted at the boundary point x anchors k FK-clusters [59,60]. In particular, if we consider the limit x 1 → x 2 in P (14)(23) we generate configurations where two distinct FK clusters meet at x 2 and are separated by a dual FK cluster (see the previous subsection and figure 3b). Such a cluster configuration is associated to the insertion at x 2 of the field φ 1,5 and therefore the leading singularity for η → 0 of P (14) (23) has to be the same as the one of the conformal block F c 3h+1 . However 3h + 1 > h > 0 and we conclude then that no other conformal blocks can enter P (14)(23) except the one of φ 1,5 . This observation identifies (apart from an overall constant) the numerator in eq. (2.3) as (1 − η) −2h G 3h+1 through eq. (3.10).
Obs. 3: Consider the OPE of two bcc operators φ αβ as they appear in eq. (3.5), it has the structure φ αβ · φ βα = 1 + X + . . . , (3.11) where 1 is the identity field and X denotes a field with scaling dimension larger than zero that is compatible with the boundary conditions. Certainly such a field cannot be φ αβ , since this would imply a discontinuity of the boundary conditions that is not allowed by the OPE in eq. (3.11), see also [48] for analogous arguments for kink fields in the bulk. Therefore we are led to the conclusion that the conformal blocks that enter the denominator in eq. (2.3) can only be F c 0 (i.e. the identity conformal block) and F c 3h+1 (i.e. the φ 1,5 conformal block).
then it easy to verify, substituting into eq. (3.7) that the symmetry under the exchange x 1 → x 3 requires the function S to satisfy the functional equation S(η) = S(1 − η). Obviously S is a solution of eq. (3.8) and from the fusion matrix given in [53] (see eq. (5.11) there) it can be moreover verified that it exists only a linear independent function S, constructed from the conformal blocks of the identity and φ 1,5 that satisfies such a property (the other can be chosen a linear combination of the conformal blocks of φ 1,3 and φ 1,5 ).
In summary the universal ratio R in eq. (2.3) will be where A Q imposes R(1) = 1 that should be clear from the geometrical interpretation in figure 1.

Analytic expressions for the connectivities
We are ready to present explicit results for the ratio R in eq. following, are shown in figure 4. Finally, for more technical details, we invite to read the three final appendices. Q = 1; percolation. We have (cf. eq. (3.1)) p = 2 and h 1,3 = h = 1/3, h 1,5 = 3h+1 = 2.
The solutions G 1/3 (η) and G 2 (η) are free of logarithms and can be easily determined using the method described in appendix A. It turns out and we actually generated O(10 5 ) terms in all the power series. The conformal blocks F 0 1/3 (η) and F 0 2 (η) are obtained through eq. (3.10) as it can be verified from appendix B. The conformal block of the identity field is singular at c = 0 [23] and needs to be regularized. In particular, eq. (3.12) does not hold verbatim at c = 0; although it might be still true in the limit Q → 1. Here, we proceed pragmatically finding the Frobenius solution with ρ = 0 of eq. (3.8), that is symmetric under the exchange η → (1 − η); such a solution is logarithmic. Using the method described in appendix A and in particular the normalization b 0 (σ) = σ we obtaiñ Notice also that the second linearly independent solution symmetric under the transformation η → (1−η) can be chosen G 1/3 (η)+G 1/3 (1−η). However such a function cannot enter into S(η) since it contains a subleading singularity η 1/3 . In conclusion, we conjecture that at Q = 1, S(η) =G 0 (η), given in eq. (4.3).
It is also important to mention that the coefficient of the logarithm in eq. (4.3) is related to the Gurarie-Ludwig [23] indecomposability parameter b. It was indeed argued that CFTs at c = 0 could be characterized by a universal number b appearing in the regularized OPE of a chiral field with itself, see eq. (1.1). In particular if h is the (chiral) scaling dimension of such a field, the leading small η behavior of the functionG 0 (notice the definition of the prefactor in eq. (3.7)) has to be [23] Comparing with eq. (4.3) for h = 1/3 we obtain b = −5/8. This is indeed the same value of the indecomposability parameter that was argued to describe critical boundary -11 -

(5.2)
In particular, the points z 1 = 0, z 3 = 1 are mapped through eq. (5.2) into the midpoints of an equilateral triangle with length-side two and vertices at w(−1/3) = −1, w(1/3) = 1 and w(∞) = −i √ 3. The image of the point z 2 moves therefore along the triangle between w(z 1 ) = 0 and w(z 2 ) = e −iπ/3 . Symmetries of the triangle are also taken into account in order to enhance the statistics. The algorithm employed is the Swendsen-Wang cluster algorithm [63] giving direct access to the FK clusters. The random number generator is given in [64] and the number of samples collected is up to 10 10 for the largest sizes considered. As the size is increased all crossing events become rarer and this happens in a more severe way for higher values of the parameter Q. This can obviously be traced back to the leading scaling dimension h 1,3 setting the dimensions of the numerator and -14 -

JHEP12(2018)131
denominator in eq. (2.3). Its value gets bigger as Q is increased. This has limited the maximal size of the triangular lattice for Q = 3, 4 to L = 129.
The data from the simulations are presented in figure 5. We first note that the simulation data fall on the predicted curves with good and increasing accuracy as the system size increases. In figure 6 we show the deviations from the CFT conjectures for the different sizes considered. The main plots show the absolute values of deviations (notice the logarithmic scale) while the insets refer to the signed relative deviations. As it can be observed, a non-trivial behavior emerges. Some features are easily explained: the deviations near the endpoints, especially apparent in the relative deviations plot, are obviously due to short-distance lattice effects, indeed their extent is confined to smaller regions as the size is increased. Extrapolation of the ratio R (eq. (2.3)) in thermodynamic limit would require a theory of finite size corrections which is largely unknown and beyond the purposes of this paper. In order to give a quantitative assessment of the convergence to our theoretical predictions we employ the L ∞ norm in the space of functions. The distance between the Monte Carlo data R MC Q (η) and the CFT prediction R CFT Q (η) is then defined by The data obtained will be extrapolated to the thermodynamic limit by fitting the L dependence through the following Ansätze: a power law of the form a PL + b PL L γ PL and only at Q = 4 (see later) a logarithmic scaling of the form a log + b log log(L) γ log . Equipped with these tools we discuss now the different values of Q considered. Since they will display qualitatively different behaviors we will separate the Q = 1, 2, 3 cases from Q = 4. Q = 1, 2, 3 cases. Percolation (Q = 1) exhibits a remarkable convergence to the conjectured formula (eq. (4.5)). The numerical data approach the thermodynamic limit in a smooth way (see figure 5). Due to this regularity we can also try to inspect the pointwise convergence. As a case study, we concentrate on the point η = 1/2. If we extrapolate to the thermodynamic limit the numerically obtained values by fitting them with a power law model R MC Q=1 (1/2) = a PL + b PL L γ PL we get a PL = 0.11766(4) fully consistent with the conjectured expression that yields R Q=1 (1/2) = 0.117680185 . . . (see figure 7 left panel). This is the same kind of agreement obtained in the most precise tests of CFT (see for example [65] where a result correct to six significant digits for the coefficient of universal area distribution of clusters in percolation is obtained). Let us now turn to the more relevant global convergence: as can be seen in the right panel of figure 7 the distance d ∞ strikingly converges, with a power law decay, to a value compatible with zero within small errors. Now we turn to the cases Q = 2 and Q = 3. Concerning pointwise convergence, we observe that the curves become more complex (especially for Q = 3), alternating regions that overshoot or undershoot the predicted curve. In addition these regions move as the size is increased making more difficult a phenomenological description, lacking the support of a theory of finite size corrections for the problem under scrutiny. When we examine instead L ∞ convergence (shown in figure 8) the value of the distance extrapolates to zero with high precision, supporting a convincing convegence to our formulas in eqs. (4.10) and (4.14). The result for the Q = 3 Potts model also suggests that small subleading corrections to the scaling should be taken into account. In conclusion, the Monte Carlo study leaves little room for doubts about the validity of the CFT predictions also at Q = 2, 3.
Q = 4 Potts model. The four-color Potts model presents different features. Indeed if we look at figure 6 we see that relative differences (lower part of the plot) are considerably larger than on the other cases. The fact that they appear drifting toward zero could be reassuring, however, looking at the upper part of the plot (absolute value of deviation in logarithmic scale), we observe that convergence, if present, to eq. (4.17) is actually very slow. In order to assess the issue quantitatively, we inspected the behavior of the distance d ∞ in eq. (5.3). Now the choice of the correct model for finite size corrections becomes crucial. Differently from the previous cases, for the Q = 4 Potts model [66][67][68][69] logarithmic corrections to the scaling are also expected; see also [38]. An additional source of complications might be moreover related to the marginality of the boundary operator φ 1,3 rightly at Q = 4. The data, shown in figure 9, are reasonably described (χ 2 per dof approximately 1.5) by the logarithmic fit a log + b log log(L) γ log with a log = (1.7 ± 2.6) · 10 −3 , compatible with zero and a slow convergence to the theoretical curve. We also tried to fit the  data with a power law decay (χ 2 per dof approximately 0.5); the two fitting curves in figure 9 are barely distinguishable. The power law fit yields an O(1) term (a PL = (8.9 ± 0.7) · 10 −3 ) that, although small, is about one order of magnitude larger than at Q = 1, 2, 3 and non-zero within the confidence interval.
We believe that both power law and logarithmic corrections should be present. Disentangling these two effects is notoriously a difficult task, requiring larger system sizes. Possible strategies we envisage for settling the question could be: the search of a geometrical model within the same universality class of FK clusters at Q = 4 but free of logarithmic corrections; a more efficient sampling of crossing events to access larger sizes, or to undertake the full analytical study of finite size corrections. All these paths are certainly worth further investigation. 3)) between the measured R MC Q and the CFT prediction as a function of L for Q = 4. The data can be fitted both with a logarithmic scaling a log + b log log(L) γ log and with a power law model a PL + b PL L γPL . We obtain a log = (1.7 ± 2.6) · 10 −3 (compatible with the theory) and a PL = (9.0 ± 0.7) · 10 −3 ; see main text for a discussion. The shaded area are the errors in the determination of a PL and a log .

Conclusions
In this paper we constructed an universal ratio R, see eq. (2.3), that involves four-point boundary connectivities of FK clusters in the two-dimensional Q-color Potts model. Exploiting lattice duality and conformal symmetry we conjectured an exact expression for R at criticality for any integer values 1 ≤ Q ≤ 4. In particular we considerably expanded the study in [45], to the three-color and four-color Potts model. Remarkably we also provided a conjecture for R in the percolation problem that corresponds to the limit Q → 1. Our theoretical results are plotted in figure 4.
The percolation case is particularly interesting since critical properties are described by a non-unitary CFT with vanishing central charge. Non-unitary extensions of minimal conformal models [12] are notoriously hard to address theoretically and few exact correlation functions have been obtained during the years. In particular earlier studies focused on generalizations of Cardy formula [14].
We calculated exactly four-point functions at c = 0 of an operator with non-vanishing scaling dimension and interpreted them as cluster connectivities in critical percolation. In particular, the sum of connectivities in figure 1 furnishes a fully explicit example of a logarithmic singularity at c = 0. Consistently with previous analysis [42] and the original proposal [23] we also re-obtained the value b = −5/8 for the indecomposability parameter of boundary percolation. This is another direct, although not easily accessible numerically, verification of the indecomposability parameter b in boundary percolation. We checked extensively our conjectures with high-precision Monte Carlo simulations on a triangular lattice, confirming both universality of the ratio in eq. (2.3) and its remarkable agreement with the predictions of conformal invariance for integers 1 ≤ Q ≤ 3. For the case Q = 4 the agreement is confirmed provided the finite size corrections are assumed to scale logarithmically with the system size [66]; a detailed study of finite size corrections at Q = 4 is a problem that we left for future work.

A Frobenius series
For mathematical details we refer to the classic volume [70]. Given the ODE in eq. (3.8) we write its truncated Frobenius series as G ρ (η) = η ρ N k=0 a k η k . We denote by L the action of the differential operator in eq. (3.8) then it turns out The k = 0 coefficient in eq. (A.1) fixes ρ throughout Once we fixed the normalization a 0 = 1, the remaining coefficients a 1 , . . . , a N in the power series are obtained solving a linear lower triangular system of equations. The solution is required in symbolic form to avoid numerical truncation errors when summing the series and can be obtained efficiently up to N = O(10 5 ) with Mathematica. In practice, the N equations for the unknowns a 1 , . . . , a N are obtained from (A.1) as The constant vector of the linear system is fixed by the normalization condition a 0 = 1. However, when there are roots in eq. (A.2) that differ by integers, the linear system (A.3) might be inconsistent [70]. In particular let ρ 1 > ρ 0 two roots of eq. (A.2) such that ρ 1 − ρ 0 = ν is a positive integer and denote by M (ρ 0 ) the N × N lower triangular matrix associated to eq. (A.3) for ρ = ρ 0 . For simplicity we assume that no other pair of roots are separated by integers. In this case the coefficient a ν will appear as a free variable (i.e. the element [M (ρ 0 ) ] νν = 0). Furthermore let M (ρ 0 ) ν the matrix obtained replacing the ν-th column of M (ρ 0 ) with the constant vector of the linear system. If the linear system is consistent (i.e. rk(M (ρ 0 ) ν ) = rk(M (ρ 0 ) )) then the particular solution will give the coefficients of the power series of G ρ 0 , normalized by a 0 = 1. The kernel of M (ρ 0 ) will be spanned by the coefficients of the power series G ρ 1 and necessarily a 0 = · · · = a ν−1 = 0; conventionally we choose a ν = 1.
In this way we generated all the linearly independent power series at Q = 3 (where a 3 is a free variable, G 0 is the particular solution and G 3 spans the kernel of M (0) ) and Q = 4 (where actually a 1 and a 4 are both free variables, G 0 is the particular solution and the kernel of M (0) is spanned by G 1 and G 4 ). At Q = 1 and Q = 2 the linear -19 -JHEP12(2018)131 system associated to eq. (A.3) can be inconsistent. For percolation this happens when ρ 0 = 0, ρ 1 = 3h + 1 = 2 and for Ising when ρ 0 = h = 1/2 and ρ 1 = 3h + 1 = 5/2; in both cases ν = 2. When the linear system is inconsistent the particular solution does not exist and it will be replaced by a Frobenius power series with a logarithmic singularity; G ρ 1 instead continues to span the kernel of M (ρ 0 ) and is free of logarithms.
The coefficients of the logarithmic solution are determined as follows (again we refer to [70] for a comprehensive discussion that includes the case of repeated roots in eq. (A.2)). We introduce a formal power series and solve the linear system in (A.3) choosing b 0 (σ) = (σ−ρ 0 ). In this way all the coefficients g k (σ) are analytic in the limit σ → ρ 0 . However since [70] b k (ρ 0 ) = 0 for k < ν, it follows that G ρ 0 determined through (A.4) is b ν (ρ 0 )G ρ 1 (recall that we chose a ν = 1 for G ρ 1 ). To obtain the linear independent solution associated to ρ 0 we observe that from the analyticity of the b k 's the differential operator L commutes with the derivative with respect to σ and it turns out The linear independent solution associated to the root ρ 0 is thenG ρ 0 (η) ≡ ∂ σ | σ=ρ 0 G σ , i.e.
where we defined β ≡ b ν (ρ 0 ) and c k = db k dσ σ=ρ 0 . It should be noticed that multiplying b 0 (σ) by any analytic function F (σ) that is O(1) at σ = ρ 0 , the above procedure produces an equally valid solution of eq. (3.8) that corresponds to a linear combination of eq. (A.6) and G ρ 1 ; in particular c 2 can be arbitrarily redefined. To efficiently generate the power series (A.6) we can first determine G ρ 1 then fix c 0 = 1 (that in turns fixes β) and choose with the choice F (σ) = 1. Finally we set up a recursion for the c k with k > 2 requiring eq. (A.6) to be a solution of eq. (3.8).

B Recursive formula for the Virasoro conformal blocks
The Virasoro conformal blocks F c ρ (η) can be obtained directly using Al. Zamolodchikov recursive formula [58]. The formula is conveniently written in terms of the elliptic nome q = e iτ , where the modulus τ is related to the anharmonic ratio by Explicit expressions for R r,s and h r,s (c) (cf. eq. (3.1) and eq. (3.6)) can be found in [58] and we will not repeat them here. Since one is interested in generating a series expansion in q (and ultimately in η) of eq. (B.3) up to order N , the level of recursion is fixed by rs = N . Notice that in principle when the internal field is degenerate, i.e. ρ = h r,s the conformal block might be singular. For all the non-logarithmic cases examined in this paper, this does not happen as the corresponding factor R r,s in the numerator of eq. (B.3) vanishes as well. The use of the recursive formula for calculating Virasoro conformal blocks for rational values of the central charge c < 1 has been discussed also in [71,72]. Since The ODE that is associated to the percolation problem is obtained replacing h = 1/3 in eq. (3.8) and reads 4(2η−1)G(η)+(−6+8η−8η 2 )G(η)+3(η−1)η(2(2η−1)G (η)+3(η−1)ηG (η)) = 0. (C.1) Explicit form for functions entering the ratio R Q=1 . Remarkably, one can formally find two linearly independent hypergeometric solutions [73] (one long and the other short) L respectively, and F S , which is real, constitute an independent basis of the solutions of (C.1) in the range 0 ≤ η ≤ 1/2. The computation of F L (η) relies upon the evaluation of a 3 F 2 hypergeometric function with argument in the range [1, ∞) and it is stable for numerical evaluation (in its Mathematica implementation) so it will be preferred to explicit series expression ofG 0 (η) and G 2 (η) (and G 1/3 (η) too) that will be also given. This is possible due to the inversion formulas JHEP12(2018)131 z → 1/z see [74] that fix the value assumed by 3 F 2 on the branch cut (1, ∞). All in all the convention results in the function acquiring a value that equals the one below the branch cut 3 F 2 (a 1 , a 2 , a 3 ; b 1 , b 2 , z) ≡ lim →0 + 3 F 2 (a 1 , a 2 , a 3 ; b 1 , b 2 , ze i(2π− ) ) for real z > 1.
In order to retrieve the functionsG 0 (η), and G 2 (η) in section 4 from the above F L (η), F S (η) within the range 0 ≤ η ≤ 1 we have to proceed in the following way: we use linear combinations of F . (C.6) The functionG 0 (η) is obtained by imposing the function to be continuous, having vanishing first derivatives and continuous second derivative in η = 1/2 and setting the function to be one for η = 0. This yields the following expressioñ G 0 (η) = α 0 √ 3 sin 2π 9 − cos 2π 9 F where ζ = min(η, 1 − η). In order to reproduce G 2 (η) we have to impose the function to be equal to F S (η) for 0 ≤ η ≤ 1/2 and impose continuity of the function and the first two derivatives at η = 1/2. The outcome is (C.8) The power series expansions of eqs. (C.7)-(C.8) coincide with eq. (4.3) and eq. (4.2) respectively. Moreover since we can also calculate explicitly the value of this function in η = 1 we can fix the normalization constant A 1 for the ratio R Q=1 given in (4.5)