Exact site-percolation probability on the square lattice

In the present case, a natural set of subproblems is the computation of F_n,m (z) for a given occupation pattern in the mth row. For n = 2, for example, either the left site, the right site or both sites are occupied:

An all empty row cannot occur in spanning configurations. The mth row can only represent a spanning configuration if the (m − 1)th row is either or , hence

Complementing the equations for and , we can write this as matrix-vector product

For general values of n we can write

$$\begin{equation}{F}_{n,m}^{\sigma }(z)=\sum\limits _{{\sigma }^{\prime }}{T}_{\sigma ,{\sigma }^{\prime }}(z){F}_{n,m-1}^{{\sigma }^{\prime }}(z),\end{equation} \tag{ 8 }$$

where σ and σ' denote possible configurations that can occur in a row of width n. Each σ in row m represents all compatible configurations in the n × m lattice 'above' it. Appropriately we refer to σ as signature of all these configurations.

The idea of the algorithm to compute ${F}_{n,1}^{\sigma }(z)$ for all signatures σ and then use (8) to iteratively compute ${F}_{n,2}^{\sigma }(z)$, ${F}_{n,3}^{\sigma }(z)$, etc, thereby reducing the two-dimensional problem with complexity O(2^nm ) to a sequence of m one-dimensional problems of complexity O(λⁿ ) each.

This particular form of a dynamic programming algorithm is known as transfer matrix method. When transfer matrices were introduced into statistical physics by Kramers and Wannier in 1941 [29], they were meant as a tool for analytical rather than numerical computations. Their most famous application was Onsager's solution of the two dimensional Ising model in 1944 [30]. With the advent of computers, transfer matrices were quickly adopted for numerical computations, and already in 1980, they were considered 'an old tool in statistical mechanics' [31]. Over the decades, the method has been applied to a wide range of models, see [32, 33] and references therein for a review. Most of these applications focused on the enumeration of finite objects that exists in an infinite lattice, like polygons of a given circumference or self avoiding walks of a given length. The idea to use transfer matrices to compute spanning probabilities in finite lattices emerged only recently, and not in statistical physics but in theoretical computer science. For their analysis of a two-player board game, Yang et al [23] developed a dynamic programming algorithm to compute numerically exact values of the spanning probability for a given value of the occupation probability p. We build on their ideas for the algorithm presented here.

Signatures are the essential ingredients of the algorithm. How do we code them? And what is their number? Both questions can be answered by relating signatures to balanced strings of parentheses.

Occupied sites in a signature can belong to the same cluster. Either because they are neighbours within the signature or because they are connected through a path in the lattice beyond the signature. In any case the set of signature sites that belong to the same cluster has a unique left to right order. This allows us to encode cluster sites using left and right parentheses (figure 1). We represent the leftmost site of a cluster by two left parentheses $\left(\left(\right.\right.\;\;\;$, the rightmost site by two right parentheses $\left.\left.\right)\right)$ and any intermediate cluster site by $\left.\right)\left(\right.\;$. An occupied site that is not connected to any other site within the signature is represented by ().

Zoom In Zoom Out Reset image size

The resulting string is a balanced parentheses expression (with some white spaces from unoccupied sites), i.e., the total number of left parentheses equals the total number right parentheses, and counting from left to right, the number of right parentheses never exceeds the number of left parentheses. The latter is due to the topological constraints in two dimensions (two clusters cannot cross). Identifying all sites of a cluster in a signature means matching parentheses in a balanced string of parentheses, a task that can obviously be accomplished in time O(n).

Since we take into account only configurations that contain at least one spanning cluster, we need to identify clusters that are connected to the top row. We could follow [23] and use brackets [[, ]], ][ and [] to represent those clusters in a signature, thereby maintaining the concept of balanced strings of parentheses. It is more efficient, however, to use only a single symbol like || for sites that are connected to the top. We can think of all || sites being connected to each other by an all occupied zeroth row. Identifying the left- and rightmost || site in a signature can again be done in time O(n).

The total number of signatures for spanning configurations depend on the width n and the depth m. Some signatures can occur only in lattices of a minimum depth. The most extreme case in this respect are nested arcs, represented by signatures like

which can only occur if m ⩾ ⌈n/2⌉. If m is larger than this, however, the total number of signatures will only depend on n. In the appendix A we will prove that this number is given by

$$\begin{equation}{S}_{n}=\sum\limits _{r=0}^{\lfloor \frac{n-1}{2}\rfloor }\left(\genfrac{}{}{0pt}{}{n+1}{2r+2}\right)\left(\genfrac{}{}{0pt}{}{2r+2}{r}\right)\frac{3}{r+3}\equiv {M}_{n+1,2},\end{equation} \tag{ 9 }$$

where M_n+1,2 is the second column in the Motzkin triangle (A.8). The first terms of the series S₁, S₂, ... are

$$\begin{equation*}1,3,9,25,69,189,518,1422,3915,10\,813,29\,964,83\,304,232\,323,649\,845,\dots \end{equation*}$$

We also show in the appendix A, that S_n grows asymptotically as

$$\begin{equation}{S}_{n}\sim \frac{\sqrt{3}{3}^{n+3}}{\sqrt{4\pi {n}^{3}}}\left(1-\frac{159}{16n}+\frac{36\,505}{512{n}^{2}}+O\left(\frac{1}{{n}^{3}}\right)\right).\end{equation} \tag{ 10 }$$

The central idea of the transfer matrix method (8) is to transform a two dimensional problem into a sequence of one dimensional problems, as a result of which the complexity reduces from O(2^nm ) to O(m2ⁿ S_n ). This idea can be applied again. By building the new row site by site (figure 2), we can subdivide the one dimensional problem into a sequence of n zero dimensional problems. This reduces the complexity from O(m2ⁿ S_n ) to O(nmS_n ).

Zoom In Zoom Out Reset image size

In this approach, the new row is partially filled up to column c. The corresponding signatures contain a kink at column c. Each signature is mapped to two signatures, depending on whether the added site at columns c + 1 is empty or occupied (figure 2).

The number of signatures S_n,c depends on n and c. We can think of a partially filled row as a full row of width n + 1, where an extra site has been inserted between positions c and c + 1, this extra site being empty in all signatures. This gets us

$$\begin{equation}{S}_{n}\leqslant {S}_{n,c}\leqslant {S}_{n+1},\end{equation} \tag{ 11 }$$

where the first inequality reflects the fact that the site at position c has only two neighbours (left and up), and therefore less restrictions to take on 'legal' values. The actual value of S_n,c depends on c, see table 2.

The number of signatures is paramount for the complexity of the algorithm. Let ${S}_{n,{c}^{\star }}\equiv {\mathrm{max}}_{c}\,{S}_{n,c}$. Figure 3 shows that ${S}_{n,{c}^{\star }}\simeq 0.350\times {2.85}^{n}$. In particular, ${S}_{n,{c}^{\star }}$ exceeds the 'giga' threshold (10⁹) for n > 20. This is significant because we need to store pairs (σ, F^σ ) for all signatures. For ${S}_{n,{c}^{\star }}$ above 10⁹, this means hundreds of GBytes of RAM. In this regime, space (memory) rather than time will become the bottleneck of our algorithm.

Zoom In Zoom Out Reset image size

It turns out, however, that we can reduce the number of signatures. Consider a signature in row m that contains the pattern . All signatures in row m + 1 that are derived from it, do not depend on whether the single () is actually occupied. Hence we can identify all signatures with signatures and store only one pair (σ', F^σ + F^σ ') instead of two pairs (σ, F^σ ) and (σ', F^σ '). For the same reason, we can also identify with and with and store only one of each pair.

In addition, the symmetry of the problem allows us to identify each signature of a full row (c = n) with its mirror image.

Let ${\tilde{S}}_{n,c}$ denote the number of signatures that haven been condensed by these rules, and let ${\tilde{S}}_{n,{c}^{\star }}={\mathrm{max}}_{c}\,{\tilde{S}}_{n,c}$. Figure 3 shows that ${\tilde{S}}_{n,{c}^{\star }}\simeq 0.196\times {2.61}^{n}$. The number of condensed signatures exceeds the 'giga' threshold only for n > 23.

The elementary step of the transfer matrix algorithm is the addition of a new lattice site in a partially filled row (figure 4). This is implemented in a subroutine extend(σ, c) that takes a signature σ and a column c as input and returns a pair σ₀, σ₁, where σ₀ (σ₁) follows from σ by adding an empty (occupied) site in column c. The rules for these computations are depicted in table 3. The signatures σ₀ and σ₁ are of course condensed before they are returned. The subroutine extend is then called from within a loop over rows and columns of the lattice (figure 5).

Zoom In Zoom Out Reset image size

Table 3. Extension rules of the transfer matrix algorithm. Each cell shows the resulting configuration if the added site is empty (top) or occupied (bottom). Some update rules are non-local: ^a || may only be deleted, if there is at least one other || elsewhere in the signature, ^b deleting $\left(\left(\right.\right.$ or $\left.\left.\right)\right)$ means shortening of the corresponding cluster, ^c the whole cluster has to be connected to the top, and ^d two clusters merge. Configurations with must never occur.

Zoom In Zoom Out Reset image size

The algorithm maintains two lists L_old and L_new of configurations, where a configuration consists of a signature σ and the corresponding generating function F^σ . L_old is initialised with the 2ⁿ − 1 first row configurations, in which each occupied site is connected to the top row by definition. L_new is initially empty.

A configuration (σ, F^σ ) is removed from L_old and σ₀, σ₁ = extend(σ, c) is computed. Then the configurations (σ₀, F^σ ) and (σ₁, zF^σ ) are inserted into L_new. This is repeated until L_old is empty. At this point, L_old and L_new swap their names, and the process can be iterated to add the next site to the lattice.

If the lattice has reached its final size nm, the total generating function is the sum over all F^σ for (σ, F^σ ) ∈ L_old.

With its six symbols , ||, $\left(\left(\right.\right.\;\;\;$, $\left.\right)\left(\right.\;\;\;$, $\left.\left.\right)\right)$ and (), a signature can be interpreted as a senary number with n digits. For n ⩽ 24, this number always fits into a 64-bit integer. Since

$$\begin{equation}{2}^{64}={3520\,522\,010\,102\,100\,444\,244\,424}_{6},\end{equation} \tag{ 12 }$$

we can also represent signatures of size n = 25 by 64-bit integers, provided we choose the senary representation wisely: $\left.\right)\left(\right.={4}_{6}$ and $\left.\left.\right)\right)={5}_{6}$.

The lists L maintained by the algorithm contain up to ${\tilde{S}}_{n,{c}^{\star }}$ configurations. An efficient data structure is mandatory. We take advantage of the fact that the signatures can be ordered according to their senary value. This allows us to use an ordered associative container like set or map from the standard C++ library. Those guarantee logarithmic complexity $O(\mathrm{log}{\tilde{S}}_{n,{c}^{\star }})$ for search, insert and delete operations [34]

Equipped with a compact representation for the signatures and an efficient data structure for the list of configurations, we still face the challenge to store the large coefficients A_n,m (k) of the generating functions ${F}_{n,m}^{\sigma }$. These numbers quickly get too large to fit within the standard fixed width integers with 32, 64 or 128 bits (figure 6).

Zoom In Zoom Out Reset image size

One could use multiple precision libraries to deal with this problem, but it is much more efficient to stick with fixed width integers and use modular arithmetic: with b-bit integers, we compute ${U}_{n,m}^{(i)}(k)\equiv {A}_{n,m}(k)\,\mathrm{mod}\,\,{p}_{i}$ for a set of prime moduli p_i < 2^b such that ∏_i p_i > max_k   A_n,m (k). We can then use the Chinese remainder theorem [34, 35] to recover A_n,m (k) from the ${U}_{n,m}^{(i)}(k)$. This way we can trade space (bit size of numbers) for time (one run for each prime modulus).

Figure 7 shows the actual memory consumption to compute F_n,n with b-bit integers. Note that data structures like red-black-tress require extra data per stored element. Unfortunately the product of all primes below 2⁸ is not large enough to compute the coefficients in F_n,n via modular arithmetic and the Chinese remainder theorem if n ⩾ 19. Hence eight-bit integers are essentially useless.

Zoom In Zoom Out Reset image size

Extrapolation of the data in figure 7 shows that the computation of F_n,n for n ⩽ 21 can be done with 16-bit integers on machines with 256 GB of RAM. Larger lattices require 384 GB (n = 22) or 1 TB (n = 23).

There is an alternative approach to determine F_n,m with less memory. Instead of enumerating the coefficients A_n,m (k), one can compute F_n,m (z) for nm integer arguments z (again modulo a sufficient number of primes) and then use Lagrange interpolation to recover F_n,m . This approach requires only a single b-bit integer to be stored with each signature. The prize we pay is that now we need nm computations, each consisting of several runs for a sufficient number of prime moduli. This was the approach we used to compute F_22,22 on 30...60 machines with 256 GB each.

Yet another method suggests itself if one is mainly interested in evaluating the spanning probability R_n,m (p) for a given value of p. Here we need to store only a single real number with each signature, and there is no need for multiple runs with different prime moduli. Quadruple precision floating point variables with 128 bit (__float128 in gcc) allow us to evaluate R_n,m (p) with more than 30 decimals precision. We used this approach to compute R_n,n (p) for n = 23, 24, see table 4.

Figure 8 shows the time to evaluate R_n,n (p) for given value p and with 30 decimals accuracy, measured on a laptop. The time complexity corresponds to the expected complexity $O({n}^{4}{\tilde{S}}_{n,{c}^{\star }})$. The n² iterations to build the lattice, combined with the O(n) complexity of the subroutine extend and the $O(\mathrm{log}{\tilde{S}}_{n,{c}^{\star }})=O(n)$ complexity of the associative container set yield a scaling factor n⁴. The number of signatures ${\tilde{S}}_{n,{c}^{\star }}\simeq 0.196\times {2.61}^{n}$ then provide the exponential part of the time complexity.

Zoom In Zoom Out Reset image size

Despite the exponential complexity, the algorithm is surprisingly fast for small values of n, even on a laptop. Systems up to 12 × 12 can be computed in less than one second, 16 × 16 takes 90 seconds. The largest system that fits in the 16 GByte memory of the author's laptop is 21 × 21. A single evaluation of R_21,21(p) on the laptop takes about 9 hours.

We have computed the generating functions F_n,m (z) for n ⩽ 22 and all m ⩽ n. Here we will briefly discuss some properties of the coefficients A_n,m (k). For simplicity we will confine ourselves to the quadratic case m = n.

A_nn (k) has a maximum A_max for k = k_max (n). Our data reveals that

$$\begin{equation}{k}_{\mathrm{max}}(n)=\begin{cases}\frac{1}{2}n(n+1)\quad \hfill & (n\leqslant 4)\hfill \\ \frac{1}{2}n(n+1)-1\quad \hfill & (4< n\leqslant 22)\hfill \end{cases}\end{equation} \tag{ 13 }$$

It is tempting to conjecture that ${k}_{\mathrm{max}}(n)=\frac{1}{2}n(n+1)-O(\mathrm{log}\,n)$.

A_max itself scales like $\simeq \hspace{-2pt}{2}^{{n}^{2}}$ (figure 6). This number determines how many prime moduli of fixed width we need to compute A_n,n (k) via the Chinese remainder theorem.

The width of the maximum at k_max scales like n⁻¹, as can be seen in the 'data collapse' when plotting $\frac{{A}_{n,n}(k)}{{A}_{\mathrm{max}}}$ versus n(k − k_max) (figure 9).

Zoom In Zoom Out Reset image size

The properties of various finite-size threshold estimators for two-dimensional percolation have been discussed by Ziff and Newman [36]. Here we will focus on estimators p_med(n) and p_cell(n), defined as

$$\begin{equation}{R}_{n,n}({p}_{\mathrm{m}\mathrm{e}\mathrm{d}})=\frac{1}{2}\end{equation} \tag{ 14a }$$

and

$$\begin{equation}{R}_{n,n}({p}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}})={R}_{n-1,n-1}({p}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}})\end{equation} \tag{ 14b }$$

The values for these estimators up to n_max = 24 are shown in table 4.

Both estimators are supposed to converge with the same leading exponent,

$$\begin{equation}{p}_{\mathrm{m}\mathrm{e}\mathrm{d}}(n)-{p}_{\text{c}}\simeq {c}_{\text{med}}{n}^{-1-1/\nu }\hspace{25.0pt}{p}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}(n)-{p}_{\text{c}}\simeq {c}_{\text{cell}}{n}^{-1-1/\nu },\end{equation} \tag{ 15 }$$

where $\nu =\frac{4}{3}$ is the critical exponent of the correlation length in two-dimensional percolation [37]. This scaling can clearly be seen in the data (figure 10).

Zoom In Zoom Out Reset image size

Any convex combination of finite size estimators is another finite size estimator. We can choose λ in λp_med(n) + (1 − λ)p_cell(n) such that the leading order terms of p_med(n) and p_cell(n) cancel:

$$\begin{equation}\lambda ={\left(1-\frac{{c}_{\text{med}}}{{c}_{\text{cell}}}\right)}^{-1}.\end{equation} \tag{ 16 }$$

The precise values of c_med and c_cell are not known, but their ratio follows from scaling arguments [36]:

$$\begin{equation}\frac{{c}_{\text{med}}}{{c}_{\text{cell}}}=-\frac{1}{\nu }=-\frac{3}{4}.\end{equation} \tag{ 17 }$$

Hence we expect that the estimator

$$\begin{equation}{p}_{\star }(n)\equiv \frac{4}{7}{p}_{\mathrm{m}\mathrm{e}\mathrm{d}}(n)+\frac{3}{7}{p}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}(n)\end{equation} \tag{ 18 }$$

converges faster than both p_med(n) and p_cell(n). This is in fact the case, as can be seen in figure 10.

Figure 10 also depicts the estimator p_pol(n), which is the root of a graph polynomial defined for percolation on the n × n torus [38–40]. Empirically, it converges like ∼n⁻⁴ in leading order. Jacobsen [24] computed this estimator for n ⩽ 21 and extrapolated the results to obtain

$$\begin{equation}{p}_{\text{c}}=0.59274605079210(2).\end{equation} \tag{ 19 }$$

This is the most precise value up to today. We have used it as reference p_c in figure 10 and we will also use it as reference in figures 12 and 13.

Extrapolation methods are commonly based on the assumption of a power law scaling

$$\begin{equation}p(n)={p}_{\text{c}}+\sum\limits _{k=1}{A}_{k}{n}^{{{\Delta}}_{k}}\end{equation} \tag{ 20 }$$

with exponents 0 > Δ₁ > Δ₂ > .... These exponents can be determined by the following procedure.

We start by considering the truncated form $p(n)={p}_{\text{c}}+{A}_{1}{n}^{{{\Delta}}_{1}}$. From each triple of consecutive data points p(n), p(n − 1) and p(n − 2) we can compute three unknowns p_c, A₁ and Δ₁. This gets us sequences A₁(n) and Δ₁(n) that we can use to compute Δ₁ = lim_n→∞Δ₁(n) and A₁ ≡ lim_n→∞ A₁(n). Figure 11 shows Δ₁(n) for p(n) = p_med(n) plotted versus n⁻¹ along with a 4th-order polynomial fit. The result ${{\Delta}}_{1}=-\frac{7}{4}$ is of course expected. It is stable under variations of the order of the fit and the number of low order data points excluded from the fit.

Zoom In Zoom Out Reset image size

We then subtract the sequence ${A}_{1}{n}^{{{\Delta}}_{1}}$ from p(n) to eliminate the leading scaling term ${n}^{{{\Delta}}_{1}}$. Considering again a truncated form $p(n)-{A}_{1}{n}^{{{\Delta}}_{1}}={A}_{2}{n}^{{{\Delta}}_{2}}$, we can get sequences A₂(n) and Δ₂(n) from three consecutive terms of this sequence. And again we extrapolate these sequences by fitting a low order polynomial in n⁻¹, thereby obtaining Δ₂ and A₂.

Obviously, this can be iterated to determine the exponents Δ_k one by one, although the stability of the procedure deteriorates to some extent with each iteration. Figure 11 shows the extrapolations for Δ₂, Δ₃ and Δ₄ for p_med with results

$$\begin{equation}{{\Delta}}_{2}=-2.7516\approx -\frac{11}{4}\quad {{\Delta}}_{3}=-3.7553\approx -\frac{15}{4}\quad {{\Delta}}_{4}=-4.7590\approx -\frac{19}{4}\end{equation} \tag{ 21 }$$

The results for p_cell are similar. These results provide compelling evidence that

$$\begin{equation}{{\Delta}}_{k}=-\frac{4k+3}{4}\hspace{25.0pt}({p}_{\mathrm{m}\mathrm{e}\mathrm{d}},{p}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}).\end{equation} \tag{ 22 }$$

These are the exponents that we will use in our extrapolations.

Applying the same procedure, Jacobsen [24] concluded that the scaling exponents for p_pol are

$$\begin{equation}{{\Delta}}_{k}=-2k-2\hspace{25.0pt}({p}_{\mathrm{p}\mathrm{o}\mathrm{l}}).\end{equation} \tag{ 23 }$$

For p_⋆(n), we could not reliably compute exponents Δ_k . It can already be sensed from the curvature in the logarithmic plot (figure 10), that p_⋆(n) is not determined by a set of well separated exponents. This is the reason why we will not use p_⋆(n) for extrapolation in this contribution, despite its fast convergence.

A study of extrapolation methods in statistical physics [41] demonstrated that Bulirsch–Stoer (BST) extrapolation [42] is a reliable, fast converging method. It is based on rational interpolation. For a given sequence of exponents 0 < w₁ < w₂ < ... and n data points p(1), ..., p(n), one can compute coefficients a₀, ..., a_N and b₁, ..., b_M with N + M + 1 = n such that the function

$$\begin{equation}{Q}_{N,M}(h)=\frac{{a}_{0}+{a}_{1}{h}^{{w}_{1}}+{a}_{2}{h}^{{w}_{2}}+\cdots +{a}_{N}{h}^{{w}_{N}}}{1+{b}_{1}{h}^{{w}_{1}}+{b}_{2}{h}^{{w}_{2}}+\cdots +{b}_{M}{h}^{{w}_{M}}}\end{equation} \tag{ 24 }$$

interpolates the data: p(k) = Q_N,M (1/k) for k = 1, ..., n. The value of a₀ yields the desired extrapolation p(k → ∞).

The best results are usually obtained with 'balanced' rationals N = ⌈(n − 1)/2⌉ and M = ⌊(n − 1)/2⌋. This is what we will use for our extrapolations.

The BST method discussed in [41] uses exponents w_k = kω, where the parameter ω is chosen to match the leading order of the finite size corrections. Here we choose ${w}_{k}=-{{\Delta}}_{k}=\frac{4k+3}{4}$.

Having fixed the exponents and the degrees of nominator and denominator, we still have a free parameter: the number of data points to be included in the extrapolation. Since we are interested in the asymptotics n → ∞ we will of course use the data of the large systems, but we may vary the size n₀ of the smallest system to be included. We expect that the extrapolation gets more accurate with increasing n₀, but only up to a certain point after which the quality deteriorates because we have to few data points left. This intuition is confirmed by the data (figure 12).

Zoom In Zoom Out Reset image size

The accuracy of an extrapolation based on system sizes n₀, ..., n_max is gauged by the difference between extrapolations of data points for n₀ + 1, ..., n_max and n₀, ..., n_max − 1. These are the errorbars shown in figure 12.

It seems reasonable to focus on the regime of n₀ where the results form a plateau with small errorbars, and to take the mean of these values as the final estimate of p_c. Using the values for n₀ = 14, 15 (p_med) and n₀ = 15, 16 (p_cell) (see table 5) we get

$$\begin{equation}{p}_{\text{c}}=0.592746050786(3)\hspace{25.0pt}({p}_{\mathrm{m}\mathrm{e}\mathrm{d}}),\end{equation} \tag{ 25a }$$

$$\begin{equation}{p}_{\text{c}}=0.59274605060(8)\hspace{25.0pt}({p}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}).\end{equation} \tag{ 25b }$$

The first value deviates from the reference value (19) by 2 errorbars, the second by 2.5 errorbars.

The reference value p_c 19 was obtained with different extrapolation approach. To see how consistent our method is, we applied it also to the estimator p_pol(n) from Jacobsen [24, table 2]. Here we used the exponents (23) in the BST extrapolation.

Figure 13 shows, that the errorbars are smaller than the errorbars on p_med and p_cell (note the scaling of the y-axis). The estimator p_pol is also more robust under variation of n₀. Taking the average of the estimates for n₀ = 1, ..., 10 gives

$$\begin{equation}{p}_{\text{c}}=0.5927460507920(4)\quad ({p}_{\mathrm{p}\mathrm{o}\mathrm{l}}).\end{equation} \tag{ 26 }$$

This value agrees with Jacobsen's value (19) within the errorbars.

Zoom In Zoom Out Reset image size

We believe that the superiority of the estimator p_pol is due to the better separation of scales: Δ_k+1 − Δ_k = −1 for p_med, p_cell versus Δ_k+1 − Δ_k = −2 for p_pol.

In a signature, the parentheses are not independent. Occupied sites that are adjacent in a row belong to the same cluster. We deal with this local connectivity by considering runs of adjacent occupied sites as single entities that we may call h(orizontal)-clusters. A row of width n contains at most ⌈n/2⌉ h-clusters. Different h-clusters may be connected through previous rows, and again we have a unique left-right order and the topological constraints that allow us to identify h-clusters with strings of balanced parentheses expression, using the set of symbols ((, )), )(and (). Any legal signature with r h-clusters corresponds to a balanced string of 2r parentheses and vice versa.

We write the occupation pattern of a row as n-bit word s ∈ {0,1}ⁿ , where 0 represents an empty site and 1 represents an occupied site. Let r(s) denote the number of h-clusters, i.e. runs of consecutive 1s, in s. Then the number of signatures is

$$\begin{equation*}\sum\limits _{s\in {\left\{0,1\right\}}^{n}}{C}_{r(s)}.\end{equation*}$$

This sum can be simplified by computing the number of n-bit words with exactly r h-clusters. For that we add a leading zero to the word. This ensures that every h-cluster has a zero to its left. The position of this zero and the position of the rightmost site of the h-cluster specify the h-cluster uniquely. Hence the number of n-bit words with r h-clusters is given by $\left(\genfrac{}{}{0pt}{}{n+1}{2r}\right)$, providing us with

$$\begin{equation}\sum\limits _{r=0}^{\lfloor \frac{n+1}{2}\rfloor }\left(\genfrac{}{}{0pt}{}{n+1}{2r}\right){C}_{r}\equiv {M}_{n+1},\end{equation} \tag{ A.3 }$$

where M_n is known as the nth Motzkin number [48, 49].

Up to this point we have ignored that our signatures have an additional symbol || for occupied sites connected to the top row. In a signature with r h-clusters, any number w ⩽ r of h-clusters can be connected to the top row:

$$\begin{equation}{S}_{n}\leqslant \sum\limits _{r=1}^{\lceil \frac{n}{2}\rceil }\left(\genfrac{}{}{0pt}{}{n+1}{2r}\right)\sum\limits _{w=1}^{r}\left(\genfrac{}{}{0pt}{}{r}{w}\right){C}_{r-w}.\end{equation} \tag{ A.4 }$$

This is only an upper bound because it does not take into account that h-clusters not connected to the top row cannot be connected to each other across a top row h-cluster. In the example shown in figure A1, the number of signatures compatible with the top row h-cluster is not C₆ = 132 but C₃ C₁ C₂ = 10.

Zoom In Zoom Out Reset image size

We need to consider the compositions that w top row h-clusters can impose upon the r − w other h-clusters. The composition of an integer n into k parts is defined as an ordered arrangement of k non-negative integers which sum to n [50]. The compositions of 4 into 2 parts are

$$\begin{equation*}(0,4),(1,3),(2,2),(3,1),(4,0).\end{equation*}$$

A number of w top row h-clusters can divide the r − w other h-clusters into w + 1 parts at most. This gets us

$$\begin{equation*}{S}_{n}=\sum\limits _{r=1}^{\lceil \frac{n}{2}\rceil }\left(\genfrac{}{}{0pt}{}{n+1}{2r}\right)\sum\limits _{w=1}^{r}\sum\limits _{\begin{subarray}{c}{i}_{k}\geqslant 0\\ {i}_{1}+\cdots +{i}_{w+1}=r-w\end{subarray}}{C}_{{i}_{1}}\dots {C}_{{i}_{w+1}}.\end{equation*}$$

Using Catalan's convolution formula [51] for 1 ⩽ k ⩽ n,

$$\begin{equation}\sum\limits _{{i}_{1}+\cdots +{i}_{k}=n}{C}_{{i}_{1}-1}\dots {C}_{{i}_{k}-1}=\frac{k}{2n-k}\left(\genfrac{}{}{0pt}{}{2n-k}{n}\right)\hspace{25.0pt}(1\leqslant k\leqslant n),\end{equation} \tag{ A.5 }$$

we can write this as

$$\begin{equation}{S}_{n}=\sum\limits _{r=1}^{\lceil \frac{n}{2}\rceil }\left(\genfrac{}{}{0pt}{}{n+1}{2r}\right)\sum\limits _{w=1}^{r}\frac{w+1}{2r-w+1}\left(\genfrac{}{}{0pt}{}{2r-w+1}{r+1}\right).\end{equation} \tag{ A.6 }$$

The sum over w can be done using binomial coefficient identities [52, table 174]. The result reads

$$\begin{align}\hfill {S}_{n}& =\sum\limits _{r=1}^{\lceil \frac{n}{2}\rceil }\left(\genfrac{}{}{0pt}{}{n+1}{2r}\right)\frac{3}{r+1}\left(\genfrac{}{}{0pt}{}{2r}{r-1}\right)\hfill \\ \hfill & =\sum\limits _{r=1}^{\lceil \frac{n}{2}\rceil }\left(\genfrac{}{}{0pt}{}{n+1}{2r}\right){C}_{r+1,r-1},\hfill \end{align} \tag{ A.7 }$$

where (A.7) follows from (A.2) with a = r + 1 and b = r − 1. This result can be categorised by considering the generalisation of (A.3). Analog to Catalan's triangle, Motzkin's triangle [53] is defined as

$$\begin{equation}{M}_{n,k}=\sum\limits _{r=0}^{\lfloor \frac{n-k}{2}\rfloor }\left(\genfrac{}{}{0pt}{}{n}{2r+k}\right){C}_{r+k,r}.\end{equation} \tag{ A.8 }$$

Hence S_n = M_n+1,2, the second column in Motzkin's triangle. This sequence is number A005322 in the on-line encyclopaedia of integer sequences.

Although Motzkin's triangle is a well known object in combinatorics, we could not find an explicit formula for the asymptotic scaling of M_n,2. Hence we compute it here, following the procedure outlined in [54].

Our starting point is the generating function for S_n = M_n+1,2 [53]

$$\begin{equation}F(x)=\sum\limits _{n > 0}{S}_{n}{x}^{n}=\frac{{\left(\left(\right.1-x-\sqrt{(1+x)(1-3x)}\right)}^{3}}{8{x}^{5}}.\end{equation} \tag{ A.9 }$$

The pole that determines the asymptotics of the coefficients is x = 1/3. Expanding F in terms of $\sqrt{1-3x}$ gives

$$\begin{equation*}F(x)=9-27\sqrt{3}{(1-3x)}^{\frac{1}{2}}+\frac{279}{2}(1-3x)-\frac{1485\sqrt{3}}{8}{(1-3x)}^{\frac{3}{2}}+\cdots \,.\end{equation*}$$

The coefficients of the binomial series that appear in this expansion are

$$\begin{equation}\begin{aligned}\hfill \left[{x}^{n}\right]{(1-3x)}^{\alpha }& ={(-3)}^{n}\left(\genfrac{}{}{0pt}{}{\alpha }{n}\right)={3}^{n}\left(\genfrac{}{}{0pt}{}{n-\alpha -1}{n}\right)\hfill \\ \hfill & ={3}^{n}\frac{(-\alpha )(1-\alpha )\dots (n-1-\alpha )}{n!}.\hfill \end{aligned}\end{equation} \tag{ A.10 }$$

For α = 0, 1, 2, 3, ..., the binomial series is finite and does not contribute to the asymptotic scaling. For α ≠ 0, 1, 2, 3, ..., the binomial coefficient can be written as

$$\begin{equation}\left(\genfrac{}{}{0pt}{}{n-\alpha -1}{n}\right)=\frac{{\Gamma}(n-\alpha )}{{\Gamma}(-\alpha ){\Gamma}(n+1)},\end{equation} \tag{ A.11 }$$

and we can use the asymptotic series for the Γ function

$$\begin{equation}\begin{aligned}\hfill \frac{{\Gamma}(n-\alpha )}{{\Gamma}(n+1)}& \sim {n}^{-\alpha }\left(\frac{1}{n}+\frac{\alpha (\alpha +1)}{2{n}^{2}}+\frac{\alpha (\alpha +1)(\alpha +2)(3\alpha +1)}{24{n}^{3}}\right.\hfill \\ \hfill & \left.\quad +\frac{{\alpha }^{2}{(\alpha +1)}^{2}(\alpha +2)(\alpha +3)}{48{n}^{4}}+O\left(\frac{1}{{n}^{5}}\right)\right).\hfill \end{aligned}\end{equation} \tag{ A.12 }$$

to obtain

$$\begin{equation}[{x}^{n}]F(x)\sim \frac{\sqrt{3}{3}^{n+3}}{\sqrt{4\pi {n}^{3}}}\left(1-\frac{159}{16n}+\frac{36\;505}{512{n}^{2}}+O\left(\frac{1}{{n}^{3}}\right)\right).\end{equation} \tag{ A.13 }$$