Exact densities of loops in O(1) dense loop model and of clusters in critical percolation on a cylinder

The subject of this letter, O(1) dense loop model (DLM), is a particular case of O(n) loop models, a class of lattice models of statistical physics formulated in terms of ensembles of paths on the lattices. Having connections with many other models they sometimes provide an alternative convenient language for the analysis. An idea of representing the partition function of the Ising model as a sum over sets of weighted contours comes back to Peierls [1]. The summation over contours with weight n assigned to loops appeared from the polygonal representation of the partition function of the random cluster model [2, 3], which, in turn, is related to the Q-state Potts model with $\sqrt{Q}=n$. Also, a connection of the O(n) loop model with the O(n) vector model, from which the former inherited the O(n) name, suggests that the former can be used to predict the critical behavior of the latter [4]. The language of O(n) loop models turned out especially efficient within the framework of Coulomb gas and conformal field theory (CFT) [5, 6], while their scaling limit fits naturally into the Schramm–Loewner evolution picture [7].

Here we consider the O(1) DLM formulated as a measure on paths on the two-dimensional square lattice. A path passes through every bond exactly once, and two paths meet at every site without crossing each other, see figure 1. All path configurations have equal weights.

Zoom In Zoom Out Reset image size

To construct configurations by local operations we place a vertex at every lattice site, in which two pairs of paths at four incident bonds are connected pairwise in one of two possible ways shown in figure 2. Both vertices are assigned the unit weight.

Zoom In Zoom Out Reset image size

The lattice we consider here is bounded in one spacial direction and unbounded in the other. Specifically, it is a strip of the square lattice, infinite in the direction that coincides with one of the lattice directions, referred to as vertical, and finite in the other (horizontal) direction with even number L = 2N of sites. Periodic boundary conditions are implied in the horizontal direction, i.e. the strip is rolled into a cylinder. Under the uniform measure on paths only finite closed loops present on such a cylinder with probability one, each loop having the weight n = 1.

This model has been intensively studied on its own and especially in view of its connection with the critical bond percolation problem, which is a particular instance of the random cluster model related to a formal Q → 1 limit of the Potts model. To go from loops to percolation we construct a new square lattice of 45 degrees rotated orientation, for which the original lattice is the so called medial graph. To this end, we put sites of the new lattice to the center of every second face of the original lattice in a staggered way, connecting them by bonds passing through the nearest sites of the original lattice as shown in figure 1. The periodic boundary conditions for the original lattice suggest that the rotated lattice is also rolled into the cylinder. Then, we consider the bond percolation on the lattice constructed. Specifically a bond of the rotated lattice is said to be open, if it is between the loop arcs and closed if it crosses them, see figure 3. All the four vertices have equal weights, i.e. open and closed bonds have equal probabilities p = 1/2, which is the critical point of the bond percolation on the infinite square lattice.

Zoom In Zoom Out Reset image size

The studies of percolation having been continuously conducted since the late fifties of the last century culminated in plenty brilliant results, see [8, 9] and references therein. In particular, the connection of percolation with O(1) DLM, in turn related to the exactly solvable six-vertex model, was found especially useful for calculating some observables and critical exponents [10]. In the context of this letter we mention the calculation of the density of critical percolation clusters on the infinite plane lattice, i.e. in the L → ∞ limit of our cylinder, performed in a seminal paper [11]. The result was obtained in the form of an integral, of which the approximate numerical value was provided. The exact value of the integral was later presented in [12] together with numerical evidences of universality of finite size corrections to this quantity occurring in confined geometries. An explicit form of the finite size corrections to the density of critical percolation clusters on the strip and on the cylinder were conjectured in [13] using arguments based on the Coulomb gas technique and CFT. Specifically, the corrections come from the conformal anomaly, the value of which was found from the mapping of the O(1) DLM to the Coulomb gas [5].

In the O(1) DLM language the quantity related to the density of percolation clusters is the density of loops. In fact, the latter can be considered as a combination of two quantities, which can be studied separately. Indeed, there are two types of loops on the cylinder, contractible and non-contractible, which do not and do wind around the cylinder respectively. Hence, we will be interested in the average numbers of loops of both types per site of the lattice. We will use the notations ν_c(L) and ν_nc(L) for the densities of contractible and non-contractible loops respectively.

To explain the relation of loop densities to the densities of percolation clusters, we note that every contractible loop is either circumscribed on a percolation cluster that does not wrap around the cylinder or is inscribed into a circuit inside a percolation cluster. The latter loop can also be thought as circumscribed on the dual percolation cluster on the dual rotated lattice. The critical point is self-dual. This means that the average numbers of percolation clusters and of dual percolation clusters are equal, and so are the average numbers of the circumscribed and the inscribed loops. Thus, the average number of non-contractible loops per unit length of the cylinder is twice the average number of the critical percolation clusters on the rotated lattice. Since the rotated lattice contains twice less sites per unit length of the cylinder than the original lattice, the density of percolation clusters not wrapping around the cylinder coincides with ν_c(L).

Also every percolation cluster that wraps around the cylinder is bounded by a pair of non-contractible loops and every non-contractible loop runs along the boundary of such a cluster. Thus, similarly to the above, we argue that ν_nc(L) coincides with the average number of critical percolation clusters wrapping around the cylinder per site of the rotated lattice.

How the mean density of non-contractible loops in O(n) model on a cylinder depends on loop weights was studied in [14]. There, the conformal anomaly being a function of the weights of contractible and non-contractible loops was obtained from the finite size correction to the energy of the XXZ chain with twisted boundary conditions [15–17] divided by a coefficient termed the sound velocity [18]. This allowed a determination of the finite size corrections to the loop densities.

Unfortunately the CFT related arguments applied to models in confined geometries were suitable only for obtaining at most sub-leading terms of the asymptotics of the mean cluster size, while the exact formulas still remained off the scope of this approach. The possibility of the next step opened in the beginning of 2000s. Then, a burst of interest to the O(1) DLM was ignited by an observation by Razumov and Stroganov of a nice combinatorial structure of the ground state of the XXZ chain and the six-vertex model at a specific combinatorial point [19]. A connection of the results of [19] to the O(1) DLM was pointed at in [20]. A number of relations between O(1) DLM, the six vertex model, the XXZ model, the fully packed loop model and alternating sign matrices came from the studies of this subject [21–23]. In particular, several sum rules for the components of the ground state eigenvector of the O(1) DLM transfer matrix and its generalizations were obtained [24–30]. Also statistics of several observables describing connectivity of boundary points [31–33] and loop embeddings [34] on finite lattices were studied in lattices with different boundary conditions like cylinder, strip, etc. For many of them exact formulas were either conjectured or proved. However, to our knowledge, the simplest quantities like the mentioned densities of loops aka the densities of critical percolation clusters are not yet in this list.

In this letter we fill this gap. We obtain the exact formulas for ν_c(L) and ν_nc(L) for any even L. To this end we exploit the connection between the free energies of the O(n) DLM and the six-vertex model. The latter can be found as the largest eigenvalue of the corresponding transfer-matrix, of which the derivatives with respect to the fugacities of contractible and non-contractible loops yield the mean values of interest. The eigenvalue satisfies the Baxter's T–Q equation as well as the conjugated T–P equation, which are the functional relations between the eigenvalue and two polynomials Q(x) and P(x) having zeroes on the roots of two systems of Bethe equations. Both equations were solved for Q(x) and P(x) in the so called stochastic point, corresponding to the O(1) DLM, by Fridkin, Stroganov, Zagier (FSZ) in [35]. Furthermore, simultaneous differentiation of the transfer-matrix eigenvalue expressed in terms of Q(x) and P(x) and of the Wronskian relation between them allows one to exclude unknown terms and to express the derivatives of the eigenvalue at the stochastic point via the derivatives of the known Q(x) and P(x) in their arguments [36]. The rest of the work and the main technical challenge of this letter is a reduction of the complicated expressions obtained in terms of the hypergeometric functions evaluated at special values of the parameters and of the argument to a manageable rational form. It is performed with the use of Kummer's theorem for hypergeometric function and its contiguous generalizations. This is the program we complete in the next sections. Similar calculations were done also in [37] with the eigenvalue of the XXZ chain and in [38] in context of the raise and peel model [39].

In this section we present the formulas obtained and compare them with the results mentioned above. Stating the results, to avoid an alternate use of the parameters N and L, we always use the parameter N implying that L = 2N.

2.1. Contractible loops

The density of contractible loops is obtained in the form.

$$\begin{align}\hfill {\nu }_{\mathrm{c}}\left(2N\right)& =\frac{3{\Gamma}\left(\frac{N}{2}\right){\Gamma}\left(\frac{3N}{2}+\frac{1}{2}\right)}{4{\Gamma}\left(\frac{3N}{2}\right){\Gamma}\left(\frac{N+1}{2}\right)}+\frac{{\pi }^{2}{2}^{-2N}{3}^{2-3N}{\Gamma}\left(3N\right)}{{\Gamma}{\left(\frac{N}{2}+\frac{1}{6}\right)}^{2}{\Gamma}{\left(\frac{N}{2}+\frac{5}{6}\right)}^{2}{\Gamma}\left(N\right)}-\frac{5}{2}\hfill \\ \hfill & =\frac{1}{8},\frac{17}{160},\frac{913}{8960},\frac{3953}{39\enspace 424},\frac{14569}{146\enspace 432},\frac{3945\enspace 737}{39\enspace 829\enspace 504},\dots \hfill \end{align} \tag{ 1 }$$

In the second line we show the numerical values of the quantity obtained for N = 1, ..., 6. One can see that the densities are rational numbers. Indeed, using the reflection formula for gamma functions the lhs of (1) can be recast in the form of an explicit rational function of N

$$\begin{equation*}{\nu }_{\mathrm{c}}\left(2N\right)=\frac{{2}^{-2\left(N+1\right)}{3}^{2-3N}\left(2-{\left(-1\right)}^{N}\right)\left(3N-1\right)!}{\left(N-1\right)!{\left(\frac{5}{6}-\frac{N}{2}\right)}_{N}^{2}}+\frac{3{\left(\frac{N+1}{2}\right)}_{N}}{4{\left(\frac{N}{2}\right)}_{N}}-\frac{5}{2},\end{equation*}$$

written in terms of the factorials and the Pochhammer symbols (a)_n = a(a + 1)...(a + n − 1). Having no oscillatory factors the formula (1) is more suitable for the asymptotical analysis. Using the Stirling formula we obtain an asymptotic expansion, which in the first three orders is

$$\begin{equation}{\nu }_{\mathrm{c}}\left(2N\right)=\frac{3\sqrt{3}-5}{2}+\frac{1}{4\sqrt{3}}{\left(2N\right)}^{-2}-\frac{23}{48\sqrt{3}}{\left(2N\right)}^{-4}+O\left({N}^{-6}\right).\end{equation} \tag{ 2 }$$

2.2. Non-contractible loops

The density of non-contractible loops we obtained is

$$\begin{align}\hfill {\nu }_{\mathrm{n}\mathrm{c}}\left(2N\right)& =\frac{{2}^{2\left(N-2\right)}{\Gamma}\left(N\right)}{N{\pi }^{2}{\Gamma}\left(3N\right)}\left({3}^{3N}{\Gamma}{\left(\frac{N}{2}+\frac{1}{6}\right)}^{2}{\Gamma}{\left(\frac{N}{2}+\frac{5}{6}\right)}^{2}-\frac{12{\pi }^{2}{\Gamma}{\left(\frac{3N}{2}\right)}^{2}}{{\Gamma}{\left(\frac{N}{2}\right)}^{2}}\right)\hfill \\ \hfill & =\frac{1}{8},\frac{11}{320},\frac{421}{26\enspace 880},\frac{1403}{157\enspace 696},\frac{4189}{732\enspace 160},\frac{952\enspace 067}{238\enspace 977\enspace 024},\dots \hfill \end{align} \tag{ 3 }$$

where we again show the first six rational values. The explicitly rational expression of ν_nc(2N) is

$$\begin{equation*}{\nu }_{\mathrm{n}\mathrm{c}}\left(2N\right)=\frac{3{\times}{2}^{2\left(N-1\right)}\left(N-1\right)!}{N\left(3N-1\right)!}\left({3}^{3N-2}\left(2+{\left(-1\right)}^{N}\right){\left(\frac{5}{6}-\frac{N}{2}\right)}_{N}^{2}-{\left(\frac{N}{2}\right)}_{N}^{2}\right).\end{equation*}$$

The Stirling formula applied to (3) to three leading orders yields

$$\begin{equation}{\nu }_{\mathrm{n}\mathrm{c}}\left(2N\right)=\frac{1}{\sqrt{3}}{\left(2N\right)}^{-2}-\frac{17}{18\sqrt{3}}{\left(2N\right)}^{-4}+\frac{1021}{216\sqrt{3}}{\left(2N\right)}^{-6}+O\left({N}^{-8}\right).\end{equation} \tag{ 4 }$$

2.3. Comparison to earlier results

As one could expect, the leading order term in the asymptotic expansion (2) coincides with the asymptotic value of the critical percolation cluster density on the infinite plane lattice obtained by Temperley and Lieb [11] and promoted to the exact numerical value in [12]. The finite size correction to this value obtained in [13] referred to all (both contractible and non-contractible) loops. It is to be compared to the sum of the sub-leading term in (2) and the leading term in (3). Our result differs by the factor of 2 from that in [13]. This difference stems from the distinction of the length scales between the original and the rotated lattice. Specifically, the result in [13] followed from the CFT prediction [40, 41] for the form of universal finite size correction to the specific free energy on the cylinder of circumference L

$$\begin{equation*}{f}_{L}={f}_{\infty }-\frac{\pi c}{6}{L}^{-2},\end{equation*}$$

where, when applied to the lattice, the length is supposed to be measured in the lattice spacings. As in our case the lattice spacing of the rotated lattice is $\sqrt{2}$, the length L should be replaced by $\sqrt{2}N$, which explains the discrepancy.

The conformal anomaly term used in [14] to evaluate the density of non-contractible loops in O(n) DLM referred to the original lattice. Indeed, the leading term in (3) exactly coincides with the value that follows from the formula in [14] with parameters corresponding to unit fugacities of both contractible and non-contractible loops. The sub-leading term of (2) can also be obtained from [14] by differentiation of the anomaly term with respect to the fugacity of contractible loops.

The solution is based on the relation between the O(1) DLM and the asymmetric six-vertex model via the directed loop model introduced as follows [3]. We give an orientation to every loop, which now can be either clockwise of anti-clockwise. To this end, we put an arrow in one of two possible directions to each arc within the vertices in figure 2. Thus we obtain eight vertices of the directed loop model. To obtain six vertices of the six-vertex model out of them we put arrows on the bonds incident to every site in directions consistent with the directions of the arcs, as shown in figure 4, ignoring the arc connectivities. Then, two vertices of the six-vertex model will be the sums of pairs of vertices of the directed loop model and the other four will be one-to-one. To be short, a summation over the loop orientations within the directed loop model leads us to the undirected DLM, while the summation over the arc connectivities, with information about the arrow directions kept, yields the six-vertex model.

Zoom In Zoom Out Reset image size

We assign the following vertex weights using the prescription of [42].

$$\begin{align*}\hfill {a}_{1}& =z\enspace {\mathrm{e}}^{\mathrm{i}\frac{\varphi }{L}},\hspace{25.0pt}{a}_{2}=z\enspace {\mathrm{e}}^{-\mathrm{i}\frac{\varphi }{L}},\hfill \\ \hfill {b}_{1}& ={\mathrm{e}}^{-\mathrm{i}\frac{\varphi }{L}},\hspace{25.0pt}{b}_{2}={\mathrm{e}}^{\mathrm{i}\frac{\varphi }{L}},\hfill \\ \hfill {c}_{1}& =z{q}^{1/2}+{q}^{-1/2},\hspace{25.0pt}{c}_{2}={q}^{1/2}+z{q}^{-1/2}.\hfill \end{align*}$$

Here z is an auxiliary spectral parameter that will be set equal to one in the end. At z = 1, the weight choice ensures that the contractible and non-contractible loops come with weights

$$\begin{equation*}w={2}_{q}=q+{q}^{-1}\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad v=2\enspace \mathrm{c}\mathrm{o}\mathrm{s}\enspace \varphi ,\end{equation*}$$

respectively. In particular we obtain the unit weights, w = v = 1, in the so called stochastic point

$$\begin{equation}q={\mathrm{e}}^{\frac{\mathrm{i}\pi }{3}},\hspace{25.0pt}\varphi =\pi /3,\end{equation} \tag{ 5 }$$

while the quantities of interest are given by derivatives

$$\begin{equation}{\nu }_{\mathrm{c}}={\left.w\frac{\mathrm{d}}{\mathrm{d}w}\right\vert }_{w=1}{f}_{L}\left(w,1\right),\hspace{25.0pt}{\nu }_{\mathrm{n}\mathrm{c}}={\left.v\frac{\mathrm{d}}{\mathrm{d}v}\right\vert }_{v=1}{f}_{L}\left(1,v\right)\end{equation} \tag{ 6 }$$

of the specific free energy f_L (w, v). The latter is equal to the logarithm of the largest eigenvalue Λ(z) of corresponding row-to-row transfer-matrix normalized to the number of sites L in a horizontal row

$$\begin{equation*}{f}_{L}\left(w,v\right)={L}^{-1}\enspace \mathrm{ln}\enspace {{\Lambda}}_{\mathrm{max}}\left(1\right).\end{equation*}$$

For the solution of the six vertex model we refer the reader to [10, 43, 44]. To summarize, the transfer matrix describing the transition between two subsequent horizontal rows of vertical bonds can be written in a basis of configurations with fixed positions of up-arrows. Since the action of the transfer matrix preserves the number of up- and down-arrows, the eigenspace is a direct sum of invariant subspaces indexed by the number M of up arrows in a horizontal row, which can take any integer value within the range 0 ⩽ M ⩽ L. For M fixed, we use the Bethe ansatz to diagonalize the transfer matrix and obtain the eigenvalues in the form

$$\begin{equation}{\Lambda}\left(z\right)={z}^{L}\enspace {\mathrm{e}}^{\mathrm{i}\varphi }\prod\limits _{i=1}^{M}\frac{1+\left(z+{2}_{q}\right){z}_{i}\enspace {\mathrm{e}}^{-\frac{2\mathrm{i}\varphi }{L}}}{z-{z}_{i}\enspace {\mathrm{e}}^{-\frac{2\mathrm{i}\varphi }{L}}}+{\mathrm{e}}^{-\mathrm{i}\varphi }\prod\limits _{i=1}^{M}\frac{{2}_{q}z+1+{z}_{i}z\enspace {\mathrm{e}}^{-\frac{2\mathrm{i}\varphi }{L}}}{{z}_{i}\enspace {\mathrm{e}}^{-\frac{2\mathrm{i}\varphi }{L}}-z},\end{equation} \tag{ 7 }$$

evaluated at numbers z₁, ..., z_M being the roots of the Bethe ansatz equations (BAE)

$$\begin{equation}{z}_{j}^{L}={\left(-1\right)}^{M-1}\prod\limits _{l=1}^{M}\frac{1+{\mathrm{e}}^{-\frac{2\mathrm{i}\varphi }{L}}{2}_{q}{z}_{j}+{\mathrm{e}}^{-\frac{4\mathrm{i}\varphi }{L}}{z}_{l}{z}_{j}}{1+{\mathrm{e}}^{-\frac{2\mathrm{i}\varphi }{L}}{2}_{q}{z}_{l}+{\mathrm{e}}^{-\frac{4i\varphi }{L}}{z}_{l}{z}_{j}},\quad i=1,\dots ,M.\end{equation} \tag{ 8 }$$

For further convenience we make a variable change

$$\begin{equation*}z=\frac{u-q}{1-qu},\hspace{25.0pt}{z}_{i}={\mathrm{e}}^{\frac{2\mathrm{i}\varphi }{L}}\frac{{u}_{i}-q}{1-q{u}_{i}},\end{equation*}$$

to arrive at the following form of the eigenvalue

$$\begin{equation}{\Lambda}\left(\frac{u-q}{1-uq}\right)={\left(\frac{u-q}{1-uq}\right)}^{L}\frac{{\mathrm{e}}^{\mathrm{i}\varphi }}{{\left(-q\right)}^{M}}\prod\limits _{i=1}^{M}\frac{u{q}^{2}-{u}_{i}}{u-{u}_{i}}+\frac{{\left(-q\right)}^{M}}{{\mathrm{e}}^{\mathrm{i}\varphi }}\prod\limits _{i=1}^{M}\frac{u{q}^{-2}-{u}_{i}}{u-{u}_{i}}\end{equation} \tag{ 9 }$$

and the system of BAE for the parameters u₁, ..., u_M to be substituted to (9)

$$\begin{equation}{\mathrm{e}}^{2\mathrm{i}\varphi }{\left(\frac{{u}_{i}-q}{1-q{u}_{i}}\right)}^{L}={\left(-1\right)}^{M-1}\prod\limits _{j=1}^{M}\frac{{q}^{2}{u}_{j}-{u}_{i}}{{q}^{2}{u}_{i}-{u}_{j}},\quad i=1,\dots ,M.\end{equation} \tag{ 10 }$$

Every solution of this system yields a particular eigenvalue. A straightforward approach would be to find the solution of (8) corresponding to the largest eigenvalue, to substitute it into (7), and to differentiate. Of course, finding a solution of nonlinear algebraic system explicitly is not a manageable problem. One, however, can get around this problem by writing an equation directly for the eigenvalue. The necessary technique is based on the T–Q equation and its FSZ solution in the stochastic point.

Alternatively, the eigenvalue problem can be rewritten as a single functional relation for the polynomial

$$\begin{equation*}Q\left(u\right)=\prod\limits _{i=1}^{M}\left(u-{u}_{i}\right),\end{equation*}$$

with the roots being a particular solution of (7). To this end we note that the eigenvalue (7) must be polynomial in z of degree at most L. Hence the quantity

$$\begin{equation*}T\left(u\right)={\Lambda}\left(\frac{u-q}{1-uq}\right){\left(1-uq\right)}^{L}{\left(-q\right)}^{M-L}\end{equation*}$$

is polynomial in u of at most the same degree. In terms of T(u) the formula (9) takes the form

$$\begin{equation}T\left(u\right)Q\left(u\right)={\mathrm{e}}^{\mathrm{i}\varphi }\phi \left(u{q}^{-1}\right)Q\left(u{q}^{2}\right)+{\mathrm{e}}^{-\mathrm{i}\varphi }\phi \left(uq\right)Q\left(u{q}^{-2}\right){\left(-q\right)}^{2M-L},\end{equation} \tag{ 11 }$$

where we introduced

$$\begin{equation}\phi \left(u\right)={\left(1-u\right)}^{L}\end{equation} \tag{ 12 }$$

and multiplied both sides of (9) by the denominator of its rhs. A condition of polynomiality of T(x) suggesting that the rhs of (11) is divisible by Q(u), i.e. the Bethe roots are zeroes of rhs of (11), is equivalent to the system (10) of BAE. The idea, however, is to attack the problem by solving the functional relation (11) for two unknown polynomials T(u) and Q(u).

Before going to the solution we introduce a conjugated problem for T(u) and another polynomial P(u) of degree (L − M). This problem arises if we solve the eigenproblem for the transfer-matrix in a different basis keeping track for positions of the down-instead of the up-arrows. This is equivalent to solving the original six-vertex model with the weights a₁, b₁, c₁ and a₂, b₂, c₂ exchanged respectively, which in practice is realized by the change φ ↔ −φ. Since we did nothing but the change of the basis, the solution of the problem with (L − M) down-arrows should produce the same eigenvalues Λ(z). Thus, for T(u) we obtain

$$\begin{equation}T\left(u\right)P\left(u\right)={\left(-q\right)}^{2M-L}\enspace {\mathrm{e}}^{-\mathrm{i}\varphi }\phi \left(u{q}^{-1}\right)P\left(u{q}^{2}\right)+{\mathrm{e}}^{\mathrm{i}\varphi }\phi \left(uq\right)P\left(u{q}^{-2}\right).\end{equation} \tag{ 13 }$$

Multiplying equations (11) and (13) by P(u) and Q(u) respectively, subtracting one from the other and analyzing the structure of zeroes of the terms of equation obtained, we arrive at the quantum Wronskian relation between Q(u) and P(u).

$$\begin{equation}\phi \left(u\right)=\frac{{\mathrm{e}}^{\mathrm{i}\varphi }Q\left(qu\right)P\left({q}^{-1}u\right)-{\mathrm{e}}^{-\mathrm{i}\varphi }Q\left({q}^{-1}u\right)P\left(qu\right){\left(-q\right)}^{2M-L}}{{\mathrm{e}}^{\mathrm{i}\varphi }{\left(-q\right)}^{2M-L}-{\mathrm{e}}^{-\mathrm{i}\varphi }}.\end{equation} \tag{ 14 }$$

Substituting this into either of (11) or (13) we also obtain

$$\begin{equation}T\left(u\right)=\frac{{\mathrm{e}}^{2\text{i}\varphi }Q\left({q}^{2}u\right)P\left({q}^{-2}u\right)-{\mathrm{e}}^{-2\mathrm{i}\varphi }Q\left({q}^{-2}u\right)P\left({q}^{2}u\right){\left(-q\right)}^{2M-L}}{{\mathrm{e}}^{\mathrm{i}\varphi }{\left(-q\right)}^{2M-L}-{\mathrm{e}}^{-\mathrm{i}\varphi }}.\end{equation} \tag{ 15 }$$

These two equations are the key formulas for our solution.

Now, we are in a position to proceed with the solution for the largest eigenvalue. This solution corresponds to the choice

$$\begin{equation}M=N=L/2.\end{equation} \tag{ 16 }$$

Thus, the factor (−q)^2M−L disappears from the equations (14) and (15), which then take exactly the form studied in [38]. In this case, the solutions of T–Q and T–P equations for the XXZ chain with anisotropy parameter Δ = −1/2 and twisted boundary conditions corresponding to the stochastic point (5) were found in [36]. Here we keep to the notations of [38] and refer the reader to the formulas in that paper. The solution for T(x) has the form

$$\begin{equation}T\left(x\right)={\left(1+x\right)}^{2N}\end{equation} \tag{ 17 }$$

that in particular suggests Λ(1) = 2^2N , which is simply the number of possible arrangements of vertices from figure 2 within one horizontal row. The polynomials Q(x) and P(x) are looked for in the form

$$\begin{equation}Q\left(x\right)={f}_{Q}\left(x\right){\left(1+x\right)}^{-2N},\hspace{25.0pt}P\left(x\right)={f}_{P}\left(x\right){\left(1+x\right)}^{-2N},\end{equation} \tag{ 18 }$$

where in our case it is convenient to represent the functions f_Q (x) and f_P (x) in terms of Gauss hypergeometric function (see [38])

$$\begin{align}\hfill {f}_{Q}\left(x\right)& =\frac{{\Gamma}\left(2/3\right)}{{\Gamma}\left(2/3-N\right)}\left(\frac{{\Gamma}\left(2/3\right)}{{\Gamma}\left(2/3+N\right)}{}_{2}{F}_{1}\left(1/3-N,-N,1/3;-{x}^{3}\right)\right.\hfill \\ \hfill & \quad \left.+\frac{{x}^{2}N{\Gamma}\left(-2/3\right)}{{\Gamma}\left(1/3+N\right)}{}_{2}{F}_{1}\left(2/3-N,1-N,5/3;-{x}^{3}\right)\right),\hfill \end{align} \tag{ 19 }$$

$$\begin{align}\hfill {f}_{P}\left(x\right)& =\frac{{\Gamma}\left(2/3+N\right)}{{\Gamma}\left(2/3\right)}\left(\frac{{\Gamma}\left(2/3-N\right)}{{\Gamma}\left(2/3\right)}{}_{2}{F}_{1}\left(2/3-N,-N,2/3;-{x}^{3}\right)\right.\hfill \\ \hfill & \quad \left.+\frac{xN{\Gamma}\left(1/3-N\right)}{{\Gamma}\left(4/3\right)}{}_{2}{F}_{1}\left(1/3-N,1-N,4/3;-{x}^{3}\right)\right).\hfill \end{align} \tag{ 20 }$$

The final step is to calculate the derivatives (6) of the free energy with respect to fugacities w and v at v = w = 1. This corresponds to differentiation of the eigenvalue with respect to q and φ at the stochastic point. Thus, in terms of T(u) we have

$$\begin{equation}{\nu }_{\mathrm{c}}=\frac{1}{2}+\frac{{\left(1-{q}^{-2}\right)}^{-1}}{2N}\frac{\mathrm{d}\mathrm{ln}\enspace T\left(1\right)}{\mathrm{d}q},\hspace{25.0pt}{\nu }_{\mathrm{n}\mathrm{c}}=-\frac{1}{2\sqrt{3}N}\frac{\mathrm{d}\mathrm{ln}\enspace T\left(1\right)}{\mathrm{d}\varphi },\end{equation} \tag{ 21 }$$

where the coefficients of the derivatives come from the change of variables v and w in (6) to q and φ respectively.

To differentiate T(1) in q we use its expression given by rhs of (15). The result

$$\begin{equation*}\frac{\mathrm{d}\mathrm{ln}\enspace T\left(1\right)}{\mathrm{d}q}={2}^{-2N}\left(2A+B\right),\end{equation*}$$

consists of two parts. One is from the explicit dependence of rhs of (15) on q. The other is from an implicit dependence of Q(x) and P(x) on q via the q-dependence of Bethe roots. The first one is encapsulated in the letter A. It contains the derivatives with respect to q applied to the arguments of Q(q^±2 u) and P(q^±2 u). The coefficient 2 of A comes from the exponent of q^±2. The quantity A can explicitly be calculated to

$$\begin{align}\hfill A& =\frac{q{Q}^{\prime }\left({q}^{-2}\right)P\left({q}^{2}\right)-{Q}^{\prime }\left({q}^{2}\right)P\left({q}^{-2}\right)-{q}^{-1}Q\left({q}^{2}\right){P}^{\prime }\left({q}^{-2}\right)-{q}^{-1}Q\left({q}^{-2}\right){P}^{\prime }\left({q}^{2}\right)}{q-{q}^{-1}} \hfill \\ \hfill & =\frac{1}{\mathrm{i}\sqrt{3}}\left\{q{f}_{Q}^{\prime }\left({q}^{-2}\right){f}_{P}\left({q}^{2}\right)-{f}_{Q}^{\prime }\left({q}^{2}\right){f}_{P}\left({q}^{-2}\right)-{q}^{-1}{f}_{Q}\left({q}^{2}\right){f}_{P}^{\prime }\left({q}^{-2}\right)-{q}^{-1}\right.\hfill \\ \hfill & \quad {\times}\left.{f}_{Q}\left({q}^{-2}\right){f}_{P}^{\prime }\left({q}^{2}\right)\right\}-2N\left\{\left({f}_{Q}\left({q}^{-2}\right){f}_{P}\left({q}^{2}\right)-\frac{1+{q}^{-1}}{\mathrm{i}\sqrt{3}}{f}_{Q}\left({q}^{2}\right){f}_{P}\left({q}^{-2}\right)\right)\right\}\hfill \end{align}. \tag{ 22 }$$

Here and below we assume that e^iφ = q = e^iπ/3 and, hence, q + q⁻¹ = 1, $q-{q}^{-1}=\mathrm{i}\sqrt{3}$ and q³ = q⁻³ = −1, q² = q⁻⁴ = −q⁻¹ = −q⁵, etc. In particular, after differentiating we replace e^±iφ by q^±1 respectively. The second part, expressed via the letter B, is unknown. However, B can be found from yet another derivative of equation (14)

$$\begin{equation*}{\left.\frac{\mathrm{d}\phi \left(u\right)}{\mathrm{d}q}\right\vert }_{u={q}^{3}}=-A+B=0,\end{equation*}$$

which apparently vanishes on one hand and includes A and B on the other. Note that we first differentiate and then substitute u = q³ = −1. As a result we obtain

$$\begin{equation}\frac{\mathrm{d}\mathrm{ln}\enspace T\left(1\right)}{\mathrm{d}q}=3{\times}{2}^{-2N}A.\end{equation} \tag{ 23 }$$

Similarly for the derivatives in φ we have

$$\begin{equation*}\frac{\mathrm{d}\mathrm{ln}\enspace T\left(1\right)}{\mathrm{d}\varphi }=\frac{1}{\sqrt{3}}\left(-1+\frac{2C+D}{{2}^{2N}}\right),\hspace{25.0pt}{\left.\frac{\mathrm{d}\phi \left(u\right)}{\mathrm{d}\varphi }\right\vert }_{u={q}^{3}}=\frac{1}{\sqrt{3}}\left(D-C-{2}^{2N}\right)=0,\end{equation*}$$

with

$$\begin{align}\hfill C& ={q}^{2}Q\left({q}^{2}\right)P\left({q}^{-2}\right)+{q}^{-2}Q\left({q}^{-2}\right)P\left({q}^{2}\right)\hfill \\ \hfill & ={q}^{2}{f}_{Q}\left({q}^{2}\right){f}_{P}\left({q}^{-2}\right)+{q}^{-2}{f}_{Q}\left({q}^{-2}\right){f}_{P}\left({q}^{2}\right),\hfill \end{align} \tag{ 24 }$$

coming from the explicit dependence on φ of the numerators of T(u) and ϕ(u) and D being the unknown part. There are also constants not included neither in C nor in D, which came from differentiating the denominators. Then we obtain

$$\begin{equation}\frac{\mathrm{d}\mathrm{ln}\enspace T\left(1\right)}{\mathrm{d}\varphi }=\frac{\sqrt{3}C}{{2}^{2N}},\end{equation} \tag{ 25 }$$

which agrees with formula obtained previously in [36] in context of XXZ chain with imaginary magnetic field.

In fact, equations (21)–(25) together with formulas (18)–(20) already provide explicit answers. However, they are not yet of satisfactory form being sums of hardly computable terms. For example, consider the four quantities f_Q (q^±2), f_P (q^±2) to be substituted to (22) and (24). Obtained from (19) and (20) by setting x = q^±2 each of them consists of two terms containing the hypergeometric function ₂ F₁(a, b, c; t) with some parameters a, b, c, evaluated at t = −x³ = −1. Likewise, by direct differentiation of (19) and (20) with the use of formula

$$\begin{equation}\partial \;{}_{2}{F}_{1}\left(a,b,c;t\right)\partial t=ab{c}^{-1}{}_{2}{F}_{1}\left(a+1,b+1,c+1;t\right)\end{equation} \tag{ 26 }$$

the derivatives ${f}_{Q}^{\prime }\left({q}^{{\pm}2}\right)$ and ${f}_{P}^{\prime }\left({q}^{{\pm}2}\right)$ necessary for (22) can be calculated to similar three term expressions. The formulas obtained in this way are not very informative. In particular, they are not suitable for the asymptotic analysis. Our next goal is to transform them to a more tractable form.

To this end, we need to evaluate the hypergeometric functions ₂ F₁(a, b, c; −1) that appear in f_Q (q^±2), f_P (q^±2) and in ${f}_{Q}^{\prime }\left({q}^{{\pm}2}\right),{f}_{P}^{\prime }\left({q}^{{\pm}2}\right)$. The well known Kummer's theorem [45] gives such an evaluation resulting in a ratio of gamma functions, provided that the relation c = 1 + a − b between the parameters holds. One can see that the parameters of hypergeometric functions in (19) and (20) satisfy relations 1 + a − b − c = ±1 shifted by ±1. In addition, according to (26), the derivatives ${f}_{Q}^{\prime }\left({q}^{{\pm}2}\right),{f}_{P}^{\prime }\left({q}^{{\pm}2}\right)$ in (22) will contain the hypergeometric functions with the argument t = −1 and parameters satisfying the original Kummer's relation as well as the relations 1 + a − b − c = ±2 shifted by ±2. Luckily, contiguous generalizations of the Kummer's theorem proved in [46] are applicable to these cases. They express the hypergeometric functions at t = −1 with parameters satisfying shifted relations as sums of ratios of gamma functions as follows

$$\begin{align}\hfill {}_{2}{F}_{1}\left(a,b,1+a-b+n;-1\right)& =\frac{{\Gamma}\left(1+a-b+n\right){\Gamma}\left(1-b\right)}{2{\Gamma}\left(a\right){\Gamma}\left(1-b+n\right)}\hfill \\ \hfill & \quad {\times}\sum\limits _{k=0}^{n}{\left(-1\right)}^{k}\left(\begin{matrix}\hfill n\hfill \\ \hfill k\hfill \end{matrix}\right)\frac{{\Gamma}\left(\frac{a}{2}+\frac{k}{2}\right)}{{\Gamma}\left(\frac{a}{2}+\frac{k}{2}-b+1\right)},\hfill \end{align} \tag{ 27 }$$

$$\begin{align}\hfill {}_{2}{F}_{1}\left(a,b,1+a-b-n;-1\right)& =\frac{{\Gamma}\left(1+a-b-n\right){\Gamma}\left(1-b-n\right)}{2{\Gamma}\left(a\right){\Gamma}\left(1-b-n\right)}\hfill \\ \hfill & \quad {\times}\sum\limits _{k=0}^{n}\left(\begin{matrix}\hfill n\hfill \\ \hfill k\hfill \end{matrix}\right)\frac{{\Gamma}\left(\frac{a}{2}+\frac{k}{2}\right)}{{\Gamma}\left(\frac{a}{2}+\frac{k}{2}-b+1-n\right)},\hfill \end{align} \tag{ 28 }$$

where $n\in {\mathbb{N}}_{0}$. For example, applying this to (19) and (20) we obtain

$$\begin{align}\hfill {f}_{Q}\left({q}^{{\pm}2}\right)& =\frac{\pi }{\sqrt{3}}\frac{{\Gamma}\left(\frac{2}{3}\right)}{{\Gamma}\left(\frac{2}{3}-N\right)}\left[\frac{1}{{\Gamma}\left(\frac{2}{3}+N\right){\Gamma}\left(\frac{1}{3}-N\right)}\left(\frac{{\Gamma}\left(\frac{1}{6}-\frac{N}{2}\right)}{{\Gamma}\left(\frac{1}{6}+\frac{N}{2}\right)}+\frac{{\Gamma}\left(\frac{2}{3}-\frac{N}{2}\right)}{{\Gamma}\left(\frac{2}{3}+\frac{N}{2}\right)}\right)\right.\hfill \\ \hfill & \quad \left.-\frac{{q}^{\mp 2}}{{\Gamma}\left(\frac{1}{3}+N\right){\Gamma}\left(\frac{2}{3}-N\right)}\left(\frac{{\Gamma}\left(\frac{1}{3}-\frac{N}{2}\right)}{{\Gamma}\left(\frac{1}{3}+\frac{N}{2}\right)}+\frac{{\Gamma}\left(\frac{5}{6}-\frac{N}{2}\right)}{{\Gamma}\left(\frac{5}{6}+\frac{N}{2}\right)}\right)\right]\hfill \\ \hfill {f}_{P}\left({q}^{{\pm}2}\right)& =\frac{{\Gamma}\left(\frac{2}{3}+N\right)}{2{\Gamma}\left(\frac{2}{3}\right)}\left[\left(\frac{{\Gamma}\left(\frac{1}{3}-\frac{N}{2}\right)}{{\Gamma}\left(\frac{1}{3}+\frac{N}{2}\right)}+\frac{{\Gamma}\left(\frac{5}{6}-\frac{N}{2}\right)}{{\Gamma}\left(\frac{5}{6}+\frac{N}{2}\right)}\right)\right.\hfill \\ \hfill & \quad \left.+{q}^{{\pm}2}\left(\frac{{\Gamma}\left(\frac{1}{6}-\frac{N}{2}\right)}{{\Gamma}\left(\frac{1}{6}+\frac{N}{2}\right)}+\frac{{\Gamma}\left(\frac{2}{3}-\frac{N}{2}\right)}{{\Gamma}\left(\frac{2}{3}+\frac{N}{2}\right)}\right)\right],\hfill \end{align}$$

where we also used the reflection identity

$$\begin{equation}{\Gamma}\left(z\right){\Gamma}\left(1-z\right)=\frac{\pi }{\mathrm{sin}\enspace \pi z}\end{equation} \tag{ 29 }$$

to reduce the number of gamma functions. Similar though longer formulas can also be obtained for ${f}_{Q}^{\prime }\left({q}^{{\pm}2}\right)$ and ${f}_{P}^{\prime }\left({q}^{{\pm}2}\right)$. Substituting the eight obtained expressions into (22) and (24) we arrive at large sums of rational combinations of gamma functions. The remaining part of the work though tedious is straightforward. It consists in step by step simplification of the obtained expressions with the use of the reflection identity (29), the Legendre duplication formula

$$\begin{equation}{\Gamma}\left(z\right){\Gamma}\left(z+1/2\right)=\sqrt{\pi }{2}^{1-2z}{\Gamma}\left(2z\right)\end{equation} \tag{ 30 }$$

and the trigonometric identities. At every step we try to reduce the number of gamma functions and trigonometric functions appearing from the reflection identity. Finally, we arrive at the shortest form we are able to obtain, which being substituted to (23) and (25) and then to (21) gives us formulas (1) and (3).

To summarize we have obtained the densities of contractible and non-contractible loops on the square lattice and of critical percolation clusters on the 45 degree rotated square lattice, both rolled into a cylinder. To this end, we used the mapping of the DLM to the six-vertex model, and the method of solution of the T–Q and T–P equations proposed earlier by FSZ. The results must also be related to some correlation functions over the stationary state of the associated Markov chain like it was e.g. in context of its continuous-time analogue raise and peel model [38]. Which one, however, is not obvious and requires careful investigation in the more complex discrete framework of the O(1) DLM.

It has been shown that the leading orders of the asymptotic expansion reproduce previous results. In particular, the sub-leading order of the density of contractible loops and the leading order of the density of non-contractible loops are known to be universal having a meaning within the CFT. The exact formulas obtained here allow in principle deriving the asymptotic expansion to any finite order. It would be interesting to understand whether the higher order finite size corrections can be combined into quantities with any degree of universality [47].

The used technique must be applicable to a wider class of models. In particular the T–Q equation for the XXZ model with free boundary conditions related to the DLM on the strip was solved in [36] and a quantity similar to our (23) was found in terms of Q- and P-polynomials. No final explicit formulas like those presented here were however obtained. Another direct generalization is to consider the lattice with odd L. In this case one infinite path on the cylinder exists and no non-contractible loops present. The technique used here was extended to odd L in [37], where it was applied to the eigenvalue of the Hamiltonian of the XXZ model. No results for the discrete lattice model were yet obtained.

More ambitious problems would be an extension of the technique to obtain exact higher cumulants of loop densities, which were also asymptotically predicted from the CFT arguments in [13]. It would also be of interest to find extensions to loop models related with higher spin [48] or rank [49] integrable systems. Another challenge is the generalization of the technique used here to integrable models at the elliptic level like e.g. eight-vertex model and XYZ spin chain. In principle the parameters analogous to our q and φ of the six vertex model can also be introduced in that case [50] as well as the technique of T–Q equation is applicable. This makes potentially possible to evaluate the conjugate ground state observables as the derivatives of the largest eigenvalue with respect to these parameters. Also, there is an analogue of the combinatorial point [51] that shows similarly nice properties and admits obtaining exact analytic results about the ground state eigenvalue [52] and eigenvector [53, 54]. Still it is not known whether an analogue of the FSZ method exists in that case. Also there is no interpretation of the eight-vertex model in terms of loops. Is there an appropriate generalization of the loop picture for that case is an interesting question. The study of these issues in both analytic and combinatorial perspectives is the matter for further investigation.

The work is supported by Russian Foundation of Basic Research under Grant 20-51-12005.

All data that support the findings of this study are included within the article (and any supplementary files).