N-cluster correlations in four- and five-dimensional percolation

Abstract

We study N-cluster correlation functions in four- and five-dimensional (4D and 5D) bond percolation by extensive Monte Carlo simulation. We reformulate the transfer Monte Carlo algorithm for percolation [Phys. Rev. E 72, 016126 (2005)] using the disjoint-set data structure, and simulate a cylindrical geometry L^{d-1}\times \infty, with the linear size up to L=512 for 4D and 128 for 5D. We determine with a high precision all possible N-cluster exponents, for N\text{}=\text{}2 and 3, and the universal amplitude for a logarithmic correlation function. From the symmetric correlator with N\text{}=\text{}2, we obtain the correlation-length critical exponent as 1/\nu \text{}=\text{}1.4610(12) for 4D and 1/\nu \text{}=\text{}1.737(2) for 5D, significantly improving over the existing results. Estimates for the other exponents and the universal logarithmic amplitude have not been reported before to our knowledge. Our work demonstrates the validity of logarithmic conformal field theory and adds to the growing knowledge for high-dimensional percolation.

Graphical abstract

Keywords

critical exponents / percolation / logarithmic conformal field theory / Monte Carlo algorithm

1 Introduction

Percolation [1] is a cornerstone of the theory of critical phenomena [2], and a central topic in probability theory [3, 4]. The bond percolation corresponds to the Q\text{}\to \text{}1 limiting case in the context of the Fortuin–Kasteleyn cluster representation of the Q-state Potts model [5, 6], and provides a simple yet vivid illustration of many important concepts for the latter. For a seminal review, see Ref. [7]. In two dimensions (2D), the algebraic use of symmetries — lattice duality [8], Yang–Baxter integrability [9, 10] and local conformal invariance [11, 12] — lead to a host of exact results for 2D systems including percolation. The bulk critical exponents β = 5/36 (for the order parameter) and ν = 4/3 (for the correlation length) are predicted by Coulomb-gas arguments [13], conformal field theory [14] and stochastic Loewner evolution theory [15], and are rigorously proved in the specific case of triangular-lattice site percolation [16]. Above the upper critical dimensionality d_u = 6, the mean-field values β = 1 and ν = 1/2 are believed to hold [17-19]. For dimensions , exact values of the critical exponents are still unavailable, and their estimates rely on numerical methods or perturbative methods [20]. Monte Carlo (MC) simulation remains a primary numerical method [21-23].

Apart from β and ν, there are many other critical exponents for percolation. At percolation threshold p_c, a variety of fractal dimensions are used to characterize the power-law scaling of the sizes of percolation clusters, the hulls, the external perimeters, the backbones and the shortest paths, etc. [2, 24], and a set of exponents is defined to account for the algebraic decay of connectivity probabilities (correlation functions) that two far-away regions are connected by a number of mono- or poly-chromatic paths [16, 25, 26]. Even in 2D, exact values of some of these exponents are still unknown. Recently, a family of N-cluster correlation functions has been studied in the framework of the logarithmic conformal field theory (LCFT) [27-30]. For integer N\ge 2, one considers connectivity probabilities that N distinct clusters propagate from a small neighborhood {\mathcal{V}}_i to another one {\mathcal{V}}_j far away, and for each N, constructs a family of correlation functions from the representation theory of the symmetric group. At p_c, these correlation functions decay algebraically as functions of the distance r, governed by a set of N-cluster exponents, and the amplitudes can exhibit rich behavior under rotations. In addition, it is predicted that a certain combination of correlation functions with N\text{}=\text{}2 depends logarithmically on distance, instead of as the usual power law, and the universal amplitude is closely related to the logarithmic coupling or indecomposability parameter in LCFT [31-34]. In 2D, the exact values of N-cluster exponents have been identified for N=\mathrm{2,3,4} with the help of MC simulation, and in 3D, high-precision numerical estimates are available for N=\mathrm{2,3} [30]. Similar results have also been obtained for the universal logarithmic amplitude. In 2D and 3D, the universality of the logarithmic amplitude was checked via simulations on different lattices.

In this work, we extend Ref. [30] to four and five dimensions. In MC study of percolation, one usually measures correlation functions in the torus geometry, i.e., a d-dimensional hypercube with periodic boundary conditions in each of the d directions [35]. One major limitation of this choice is that the required computer memory grows rapidly with system linear size as ~L^d, making it hard to simulate large system for large d. For example, simulation of a 4D (5D) hypercube of linear size L=192 (68) needs more than five gigabyte (GB) memory. On the other hand, a high-precision determination of percolation thresholds and of critical exponents would request data for large system sizes. We alleviate this problem by adopting the transfer Monte Carlo algorithm (transfer MC algorithm) in Ref. [36] to effectively simulate an infinitely long cylinder of size L^{d-1}\times \infty in d dimensions, in which each layer corresponds to a (d-1)-dimensional hypercube with periodic boundary condition.

The main strategy of the transfer MC algorithm, which can be regarded as a variant of the celebrated Hoshen-Kopelman (HK) algorithm [37], is to iteratively add a layer of L^{d-1} lattice sites during each MC step. Only the information about connected components, i.e., to which cluster each site belongs, of the current and newly added layers is stored in computer memory. The reduction of computer memory from \mathcal{O}(L^d) to \mathcal{O}(L^{d-1}) enables one to simulate much larger systems, L(\text{4}D)=512 with two GB memory and L(\text{5}D)=128 with four GB memory. Moreover, the disjoint-set data structure, a simple and well-known data structure in computer science [38], can be well implemented to efficiently update the connectivity information when occupied bonds are sequentially added to the cylindrical system [39, 40]. Previously, the disjoint-set data structure has been widely used in simulation of percolation [41-43]. Thus, we reformulate the transfer MC algorithm in Ref. [36] using the disjoint-set data structure.

We remark that, while the transfer-matrix technique [44] — a powerful research tool in statistical mechanics — also studies a cylindrical geometry, the transfer MC algorithm is a MC sampling method by definition. In the former, all possible configurations are summed up when one goes from one layer to the next, so the results are exact for the given size L once the probabilities have exponentially converged. The price is that the required computer memory grows exponentially fast as L increases, and thus the use of the transfer-matrix method is normally restricted to two-dimensional systems.

We extensively simulate critical 4D and 5D bond percolation, and by finite-size scaling analysis, determine with a high precision all possible N-cluster exponents for N=2 and 3 and the universal amplitude for the logarithmic correlation function. For N=2, the exponent for the symmetric correlation function reduces to the two-arm exponent, which for the case of percolation, is also related to the red-bond exponent y_red and the correlation-length exponent ν. In the renormalization group treatment, the latter further relates to the thermal renormalization exponent as y_{\text{t}}=1/\nu. We obtain the correlation-length exponent as 1/\nu \text{}=\text{}1.4610(12) for 4D and 1/\nu \text{}=\text{}1.737(2) for 5D, which are consistent with and significantly improve over the most recent results y_{\text{t}}(\text{4}\text{D})=1/\nu =1.459(6) and y_{\text{t}}(\text{5}\text{D})=1.747(5) in Ref. [45]. The reliability of our estimates and quoted errors are carefully examined, and it is suggested that the thermal exponent y_t in 5D is unlikely to be the central value 1.747 reported in Ref. [45]. In addition, our estimates agree well with the preliminary results y_{\text{t}}(\text{4}\text{D})=1.453(37) and y_{\text{t}}(\text{5}\text{D})=1.741(9) [46], which are obtained from the {\varphi }^3 computations up to five-loop order.

The remainder of this work is organized as follows. Section 2 describes in detail the reformulation of the transfer MC algorithm, and Section 3 defines N-cluster correlation functions and the logarithmic correlation functions. The simulation details and the fitting method are explained in Section 4, and the results are presented in Section 5. Finally, a discussion is given in Section 6.

2 Algorithm

In this section we shall describe in detail the transfer MC algorithm in Ref. [36] in the language of bond percolation on a d-dimensional hypercubic lattice of side length L. Analogous procedures can be readily obtained for site percolation. The main strategy of the algorithm is sketched in Fig. 1. A layer of L^{d-1} lattice sites is added during each MC step, and repeating this operation leads to an infinitely long cylinder L^{d-1}\times \infty. We refer to the direction along the cylinder as the transfer direction t and the perpendicular ones as the spatial directions x, and thus each lattice site is specified by its coordinate (x, t).

Each bond is occupied with probability 0\le p\le 1, and two sites connected through a chain of occupied bonds are said to be in the same connected component, which is also called a percolation cluster. Occupied bonds are randomly placed with probability p within the new layer t+1 and between the tth and (t+1)th layers. The connectivity of the (t+1)th layer is solely determined by the connectivity of the tth layer, and the newly placed bonds. Therefore, we only need to store the connectivity of the two most recent layers. This reduces the memory cost from \mathcal{O}(L^d) to \mathcal{O}(L^{d-1}), enabling one to simulate much larger systems.

Note that the cylinder is actually half-infinite, since the tth-layer connectivity is affected only by historical layers . In other words, the current layer t is a free surface. To study bulk behavior, we simultaneously simulate two such cylinders, and then “glue" the two free surfaces into a bulk system by adding a further layer of randomly occupied bonds in-between, as illustrated in Fig. 1. A caution is that before the “gluing", the connectivities of the free surfaces should be backed up to keep growing the cylinders.

The disjoint-set data structure, a simple tree-like data structure, supporting two simple operations — Union and Find, can be used to maintain the connectivity information. For convenience, in the actual coding we specify a lattice site by an integer x\in [\mathrm{1,2}V] with x\le V for the tth layer and x>V for the (t+1)th layer, where V=L^{d-1} is the volume of each layer. As shown in Fig. 2, each site x has a “parent" lattice site x^{\prime }, indexed by x^{\prime }=parent(x). All sites in the same percolation cluster form a tree graph, which is uniquely identified by the tree root. A lattice site is a root iff it points to itself, x=parent(x). Thus, given any site x, the label of the percolation cluster can be easily found by following index parent(x) until the root is reached, as illustrated by operation Find(x) in Fig. 2.

For the transfer MC algorithm, occupied bonds are added along the transfer direction t sequentially. When a bond is added between site x and y, we apply the operation Union(x, y) to change the connectivity information: making one of the two roots be the “parent" site of the other. Meanwhile, a trick is adopted: the new root is set to be the “parent" sites of both x and y, as illustrated in Fig. 2. The core ingredient for using the disjoint-set data structure is then the Union–Find algorithm [38]. We mention that the tree-like data structure ignores the cycle information about multi-connectivities in percolation clusters, and thus cannot be directly applied if there are bond-deletion operations [35].

The computational complexity of the Union–Find method mainly depends on operation Find(x) for finding the label of a cluster, as determined by the depths of trees. There exist many techniques that can be utilized to avoid deep trees, such as path compressing and union by tree depth or tree size [39, 40]. Certain combinations of these techniques are guaranteed to provide near-constant-time complexity on average to merge two clusters [39, 40]. In our implementation of Union(x, y), we simply make the root of y be the “parent" site of the root of x, which is convenient for propagating the connectivity information from the tth to the (t+1)th layer. Further, at the end of constructing connectivity for an entire layer, operation Flat(x) is taken to make all trees fully flat, i.e., with minimal depth. In the simulation, we find that it suffices to take the simple implementation in Fig. 2.

With these operations, we reformulate the transfer MC algorithm into the following steps.

Step 1: Construct clusters in the new layer. For each of the two cylinders, we take a new layer of V=L^{d-1} lattice sites, sequentially visit each pair of neighboring sites and place an occupied bond with probability p, and construct percolation clusters within the layer. The pseudocode is shown in Step 1, where “rand( )" draws a uniform random number in [\mathrm{0,1}).

Step 2: Add the new layer to the cylinder. The new layer of lattice sites, where the percolation clusters are already constructed, is added to the cylinder by randomly placing bonds between the tth and (t\text{}+\text{}1)th layers. The pseudocode is given in Step 2 and is obviously similar to Step 1. A key ingredient is to incorporate the connectivity information of the tth layer into that of the (t+1)th layer, which is naturally realized by our implementation of operation Union(z\mathrm{,}x). This leads to an important property that every site in the (t+1)th layer points to some site also from the (t+1)th layer. Therefore, after Step 2, all the sites from the (t+1)th layer alone form a self-contained disjoint-set data structure, and the tth-layer connectivity becomes obsolete and can be simply discarded.

Step 3: Print the connectivity onto the tth layer. The self-contained connectivity of the newly added (t+1)th layer is printed onto the tth layer, as in Step 3. By repeating Steps 1, 2, 3, one can grow a half-infinite cylinder along the transfer direction.

Step 4: Glue two free surfaces into a bulk system. The layer in the front of each of the two cylinders is a free surface, and the study of bulk behavior can be achieved by gluing the two free surfaces into a bulk system. After Step 3, the “gluing" operation can be readily applied to the two (t+1)th layers, without affecting the tth-layer connectivity. By a slight modification of Step 2, one obtains two bulk layers (the former surface layers) on which the N-cluster correlation functions can be sampled.

The above-formulated transfer MC algorithm can be further optimized from several aspects. For instance, instead of adding layer by layer, one can grow the cylinder by adding site by site, analogous to the sparse-matrix factorization in the traditional transfer-matrix technique. This can further save the computer memory, particularly if one is only interested in the surface properties. For high-dimensional percolation, the occupation probability p is small near criticality, and one can apply the trick of cumulative probability to use fewer random numbers [35].

3 Sampled quantities

We use the transfer MC algorithm of Section 2 to simulate bond percolation on 4D and 5D hypercubic lattices, with periodic boundary conditions in each of the perpendicular directions. Measurements take place in the bulk layers, and all sites considered below are within the same bulk layer. The observables are the same as those in Ref. [30] and shall be explained for completeness.

Let {\mathcal{V}}_i\equiv (i_1\mathrm{,}{i}_2\mathrm{,}\dots \mathrm{,}{i}_N) denote N lattice sites in a small neighborhood. We usually take their positions to be aligned, r_{i_{m+1}}=r_{i_m}+\delta, with m\text{}=\text{}\mathrm{1,2,}\dots \mathrm{,}N-1. For \left|\delta \right|=1, i_m and i_{m+1} are nearest neighbors. Let another site set {\mathcal{V}}_j\equiv (j_1\mathrm{,}{j}_2\mathrm{,}\dots \mathrm{,}{j}_N) be distant from {\mathcal{V}}_i by r=r_j-r_i, with r=\left|r\right|\gg 1. We consider events in which N distinct percolation clusters propagate from {\mathcal{V}}_i to {\mathcal{V}}_j, i.e., each cluster connects a site in {\mathcal{V}}_i to another site in {\mathcal{V}}_j. There are N! such events, symbolically represented as and for N=2, , , , , and for N=3, etc. For instance, is the probability for the event that i₁ and j₁ are connected by one cluster, and i₂ and j₂ are connected by another cluster. Event differs from event in the pairing between sites in {\mathcal{V}}_i and {\mathcal{V}}_j, i.e., i₁ is connected to j₂ and i₂ is connected to j₁ for . Within each bulk layer, the probability for each event is sampled, and \delta is perpendicular to r with \left|\delta \right|=1.

For N=2 and 3, according to the LCFT theory, appropriate linear combinations of the probabilities [,, etc.] give access to the operator content of the underlying field theory [27, 29]. More precisely, each combination corresponds, in the continuum limit, to the two-point function of an operator. This correspondence relies on the local {\mathcal{S}}_N symmetry between the N spins of {\mathcal{V}}_i (or {\mathcal{V}}_j), and the {\mathcal{S}}_Q symmetry of the Q-state Potts model. Note that {\mathcal{S}}_Q is subtly non-trivial, since percolation is not Q=1 but rather Q\to 1. The definitions of observables acting on N=2 and N=3 spins are recalled below. Each of them corresponds, technically, to a pair of Young diagrams for {\mathcal{S}}_N and {\mathcal{S}}_Q [29].

Consider first observables describing the propagation of N=2 clusters. There are two different combinations, corresponding to the symmetric and antisymmetric Young diagrams of {\mathcal{S}}_2,

In the continuum limit these correlation functions correspond to the two-point functions of two operators {\mathcal{O}}_{2\text{s}} and {\mathcal{O}}_{2\text{a}} respectively. We use different fonts to distinguish between the probabilities of events which are directly measured in the numerical work, and certain combinations thereof which are found to have particular scaling forms in the continuum limit. Below, we also use the term observable to describe a two-point function. The scaling dimensions of these operators in 2D CFT are available [29, 30].

For N=2, we could also define a logarithmic correlation function F as [27]:

\begin{array}{}\text{(2)}& F(r)=\frac{{\mathbb{P}}_0(r)+{\mathbb{P}}_1(r)-{({\mathbb{P}}_{\ne })}^2}{P_{2\text{s}}(r)}\sim \delta \ln (r)\mathrm{,}\end{array}

where is the probability that each of the four specified points belongs to a different percolation cluster; {\mathbb{P}}_1 is the probability that the points belong to three different clusters, one of which propagates from one site in {\mathcal{V}}_i to another site in {\mathcal{V}}_j, i.e., . Note that increases with r and converges to {({\mathbb{P}}_{\ne })}^2 for r\to \infty, where {\mathbb{P}}_{\ne } is the probability that the two points in {\mathcal{V}}_i belong to different percolation clusters. The composite observable F is expected to behave logarithmically as in Eq. (2), with δ a universal factor according to LCFT [27].

For N=3 clusters, the relevant combinations are

where {\mathbb{P}}_{N\circ } (with subscript \circ =\text{s}\mathrm{,}\text{m}\mathrm{,}\text{a}) refers to the symmetric, mixed and antisymmetric Young diagrams of symmetry {\mathcal{S}}_3.

At criticality, the N-cluster correlation functions {\mathbb{P}}_{N\circ } are expected to decay algebraically as r^{-2X_{N\text{o}}}, with (a priori) distinct and symmetry-dependent scaling dimensions, i.e., critical exponents X_{N\text{o}}. It is also predicted that under rotations of the relative angle between the two neighborhoods {\mathcal{V}}_i and {\mathcal{V}}_j, the amplitudes of the algebraically decaying functions exhibit nontrivial rotational dependences, in accordance with the corresponding conformal spins [29, 30].

The following is the list of sampled quantities,

• N=2: P_{2\text{s}}, P_{2\text{a}}, and F,

• N=3: P_{3\text{s}}, P_{3\text{m}}, and P_{3\text{a}}.

4 Simulation and fitting method

The simulation of the 4D and 5D bond percolation uses the transfer MC algorithm of Section 2 and is carried out at the percolation threshold, which is taken as p_c(\text{4}\text{D})=0.16013122 [47, 48] and p_c(\text{5}\text{D})=0.11817145 [47]. In 4D, we take system sizes L = 4, 6, 8, 10, 12, 14, 16, 18, 20, 24, 32, 40, 48, 64, 96, 128, 192, 256, 384, 512. For each half-infinite cylinder, more than 2\times {10}^9 surface layers are generated for each L\le 32, and at least 4\times {10}^7 layers for each L>64. In 5D, we take L = 4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 64, 80, 96, 128, and for each half-infinite cylinder, generate more than 5\times {10}^8 surface layers for each L\le 24 and at least 10⁷ layers for each L>24. Initial simulations of 10L layers for each half-infinite cylinder are discarded before measurements are taken. In total, about 2\times {10}^6 CPU hours \approx 228 CPU years are used.

Measurements are taken within each of the two finite bulk layers of L^{d-1} lattice sites. According to the finite-size scaling theory [49], we expect that at criticality, the N-cluster correlation functions of distance r behave as

\begin{array}{}\text{(3)}& P_{N\text{o}}(r\mathrm{,}L)\sim r^{-2X_{N\text{o}}}{\tilde{P}}_{N\text{o}}(r/L)\mathrm{,}\end{array}

where {\tilde{P}}_{N\text{o}} is a universal scaling function and X_{N\text{o}} is an N-cluster exponent (with subscript \circ \in {s, a} for N=2 and \circ \in {s, m, a} for N=3). For simplicity, we set r=L/2 so that P_{N\text{o}}(r\mathrm{,}L) only depends on the linear size L as

\begin{array}{}\text{(4)}& P_{N\text{o}}(L)\sim L^{-2X_{N\text{o}}}\mathrm{.}\end{array}

Note that one can in principle take r=aL with a constant. Nevertheless, while a smaller value of a would enhance the amplitude of P_{N\text{o}}, stronger finite-distance corrections would also occur. Meanwhile, the choice of r=L/2 benefits from the fact that the N-cluster correlators have equal contributions from both directions of connection in the periodic system.

We expect that the logarithmic function F(r=L/2) diverges logarithmically as F(L)\asymp \tilde{\delta }\ln (L). Nevertheless, since the correlators at distance r=L/2 in Eq. (2) are unavoidably affected by finite system sizes, as reflected by the finite-size scaling form in Eq. (3), it is not clear whether the amplitude \tilde{\delta } is universal and equals to the coefficient δ obtained by studying the r-dependence. Therefore, we carry out further extensive simulations to measure the universal amplitude δ in a new procedure, in which the system sizes are fixed to L=384 in 4D and L=96 in 5D. We take distances r=\mathrm{3,4,5,6,7,8,10,13,16} in 4D and r=\mathrm{3,4,5,6,7,8} in 5D. For each half-infinite cylinder, more than 3\times {10}^7 (6\times {10}^7) surface layers are generated in 4D (5D). The extra simulations alone took about 10⁶ CPU hours \approx 114 CPU years.

According to the least-squares criterion, we fit the MC data of the N-cluster correlation functions to

\begin{array}{}\text{(5)}& P_{N\text{o}}(L)=L^{-2X_{N\text{o}}}\left(a+b{L}^{y_1}\right)\mathrm{,}\end{array}

the logarithmic correlation function F(r=L/2) to

\begin{array}{}\text{(6)}& F(L)=a+\tilde{\delta }\ln (L)\mathrm{,}\end{array}

and logarithmic correlation function F(r) to

\begin{array}{}\text{(7)}& F(r)=a+\delta \ln (r)\mathrm{,}\end{array}

where a is a constant and b{L}^{y_1} accounts for the leading finite-size correction term with exponent .

As a precaution against other correction-to-scaling terms which we fail to include in the fitting ansatz, we impose a lower cutoff L\ge L_{\text{m}} (r\ge r_{\text{m}}) on the data points admitted in the fits, and systematically study the effect on the {\chi }^2 value when increasing L_{\text{m}} (r_{\text{m}}). Generally, we prefer fits corresponding to the smallest L_{\text{m}} (r_{\text{m}}) for which the goodness of fit is reasonable and subsequent increase in L_{\text{m}} (r_{\text{m}}) does not cause the {\chi }^2 value to drop by vastly more than one per degree of freedom. In practice, by “reasonable” we mean that {\chi }^2/\text{D}\text{F}\lesssim 1, where “DF” is the number of degrees of freedom.

The error of our estimates consists of two parts, the statistical error and the systematic error. The statistical error is the error of MC simulations due to the randomness of the sampling procedure. All the observables have a statistical error, which will enter into the fitting results. The errors listed in the tables for each individual fit are all statistical errors. The systematic error is due to finite-size corrections. To account for it, we perform fits with different values of L_m and y₁. The confidence interval of our final estimates are set to be the union of confidence intervals in all individual fits with different fit conditions.

5 Results

5.1 N=2 correlation functions

The N=2 correlation functions scale as P_{\text{2}\text{s}}\sim L^{-2X_{\text{2}\text{s}}} and P_{\text{2}\text{a}}\sim L^{-2X_{\text{2}\text{a}}}. The MC data for these correlation functions are shown in Fig. 3. We fit the data for P_2s and P_2a to Eq. (5) and the fitting results are shown in Table 1.

The P_2s data are well described by Eq. (5) with y₁ fixed to -4 and -3 in 4D and 5D respectively. In 4D, the data for L_{\text{m}}\ge 24 can be fitted without the finite-size correction term, i.e., b=0. From these fits we take our final estimate to be X_{\text{2s}}(\text{4}\text{D})=2.5390(12) and X_{\text{2s}}(\text{5}\text{D})=3.263(2).

The exponent X_2s is also called the two-arm exponent [16, 26]. For the case of percolation, it is further related to the thermal renormalization exponent y_{\text{t}}=1/\nu or the red-bond exponent y_red as X_{\text{2s}}=d-y_{\text{t}}=d-y_{\text{red}} [50, 51]. Thus, from the results for X_{\text{2s}} we obtain y_{\text{t}}(\text{4D})=1.4610(12) and y_{\text{t}}(\text{5D})=1.737(2). In comparison with the most recent results y_{\text{t}}(\text{4D})=1.459(6) and y_{\text{t}}(\text{5D})=1.747(5) in Ref. [45], our estimates have a much higher precision.

To examine the reliability of our final results, we plot P_{\text{2s}}{L}^{2X_{\text{2s}}} vs. L^-4 (L^-3) for 4D (5D) in Fig. 4. Three different values of X_{\text{2s}} are used, corresponding to our final quoted value, as well as those with three standard deviations 3\sigma away. Using y_{\text{t}}(\text{4D})=1.4610 and y_{\text{t}}(\text{5D})=1.737, approximately straight lines are produced for both 4D and 5D, and in the L\to \infty limit, the P_{\text{2s}}{L}^{2X_{\text{2s}}} values quickly converge to some constants. It can be seen that finite-size corrections are rather minor, particularly in 4D. By contrast, when using the values away by 3\sigma, the P_{\text{2s}}{L}^{2X_{\text{2s}}} data significantly bend upward or downward, illustrating the robustness of our results. In 5D, the black star data points correspond to y_{\text{t}}=1.743. If y_{\text{t}}=1.747 (X_{\text{2s}}=d-y_{\text{t}}=3.253), 5\sigma away from our estimate, were used, the bending-down curvature would be even more severe. Thus, even though our estimate y_{\text{t}}=1.737(2) and result y_{\text{t}}=1.747(5) in Ref. [45] are basically consistent, the central value y_{\text{t}}=1.747 is nearly excluded for 5D.

Compared to P_{\text{2s}}, the antisymmetric correlation function P_{\text{2a}}(L) decays much faster as a function of distance r, because of the cancellation between the connectivity probabilities and ; the P_{\text{2a}}(L) data are also less precise. The fitting results for P_{\text{2a}} are also shown in Table 1. In 4D, leaving y₁ as a free parameter in the fit shows that y₁ is around 2 and 3 with large uncertainties. From the fitting results with y₁ fixed to 2 and 3, we obtain X_{\text{2a}}(\text{4D})=3.14(3). For 5D, our final estimate is X_{\text{2a}}(\text{5D})=4.04(8). Following the same trick as for P_{\text{2s}}, the reliability of our results is examined in Fig. 5, which plots P_{\text{2a}}{L}^{2X_{\text{2a}}} vs. L^-3 for 4D and 5D using three different values of X_{\text{2a}}.

5.2 Logarithmic correlation functions

In the definition of the composite logarithmic correlation function F by Eq. (2), there exist cancellations among various connectivity probabilities in the numerator, as well as cancellations of the algebraic decay between the numerator and the denominator, leading to the logarithmic divergence as a function of distance r. The exact structure in Eq. (2) is derived directly in the thermodynamic limit. It is not a priori clear that the cancellations in F(r\mathrm{,}L), calculated at half linear size r=L/2, are the same to ensure a logarithmic finite-size dependence F(r=L/2)\sim \ln (L).

The data for F(L) are shown in Fig. 6, where the largest system size is rather limited — L(\text{4D})=32 and L(\text{5D})=16. Due to the cancellations, it is computationally expensive to have accurate MC data for larger L. Nevertheless, the approximately straight lines in the semi-log plot (Fig. 6) clearly show the logarithmic diverging behavior F\sim \ln (L). According to the least-squares criterion, the data are fitted to Eq. (6) and the results are shown in Table 2. From the fitting results, we take our final estimate of the amplitude \tilde{\delta } to be \tilde{\delta }(\text{4D})=1.11(5) and \tilde{\delta }(\text{5D})=0.76(4).

The data for F(r) are shown in Fig. 7, which clearly confirm the logarithmic diverging behavior F\sim \delta \ln r. They are fitted to Eq. (7) and the results are shown in Table 2. We take our final estimates to be \delta (\text{4D})=1.16(1) and \delta (\text{5D})=0.74(6), consistent with the results for \tilde{\delta }. In 4D, δ has a significantly smaller error than \tilde{\delta }.

To further check whether the two amplitudes always agree, we have also performed extensive simulations to measure \tilde{\delta } in 2D. In previous work we found \delta =1.12(3) [30], in good agreement with the analytical result \delta (\text{2D})=2\sqrt{3}/\pi [27]. We now obtain \tilde{\delta }=1.03(1) in 2D, indicating that the two amplitudes do differ in general. By contrast, the difference in 4D is at best marginally discernible, and in 5D not at all, to within the numerical precision of the data shown in Table 2.

5.3 N=3 correlation functions

The N=3 correlation functions include P_3s, P_3m, and P_3a, and the data are shown in Fig. 8. As N increases, the probability that N distinct clusters propagate from a small neighborhood to another one far away drops quickly. The N=3 correlation functions decay very rapidly as L increases, particularly the antisymmetric correlator P_3a. Because of cancellations of various connectivities in Eq. (3), the scaling of P_3m is effectively the sub-leading behavior of 3-cluster propagating probabilities, and that of P_3a corresponds to the sub-sub-leading behavior. As a result, it is very challenging to obtain meaningful data with reliable statistical errors for large N and/or large L.

For the symmetric correlator P_3s, we obtain data with maximum L=64 for 4D and 32 for 5D, with magnitudes in \mathcal{O}({10}^{-14}) and relative errors about 20%. The fitting results are shown in Table 3. In 4D, P_3s data for L\ge 8 can be well described with correction exponent y_1=-3 and the correction amplitude vanishes, b=0, for L_{\text{m}}\ge 10. In 5D, P_3s data can be described with correction exponent y₁ around -3. The fits with y_1=-3 or -5/2 give consistent results L_{\text{m}}\ge 8. From these fits we take our final estimate to be X_{\text{3s}}(4D) = 4.39(2) and X_{\text{3s}}(5D) =5.15(20), of which the reliability is tested in Fig. 9.

For P_3m, we have data with maximum L=28 for 4D and 16 for 5D, which has a magnitude in \mathcal{O}({10}^{-12}) and relative errors about 20%. The fitting results by Eq. (5) are shown in Table 3. In both 4D and 5D, P_3m data can be described with correction exponent y_1=-4. From these fits we take our final estimate to be X_{\text{3m}}(4D) = 4.9(3) and X_{\text{3m}}(5D) =6.55(20). In Fig. 10, we plot P_{\text{3m}}{L}^{2X_{\text{3m}}} vs. L^-4 for 4D and 5D using three different values of X_{\text{3m}}.

It is very difficult to obtain a reliable estimate of the exponent X_{\text{3a}} for the antisymmetric correlator P_3a, which decays extremely rapidly as L increases. For instance, P_{\text{3}\text{a}} in 5D decays algebraically with exponent 2X_{\text{3a}}\approx 16.0, leading to a tiny value \le {10}^{-14} already for L=16. In addition, as a residual effect that survives during cancellations of various 3-cluster connectivities, P_3a has the same statistical variance as the symmetric correlator P_3s, making it hard to reduce the statistical error of P_3a. We have the P_3a data only up to L=14 for 4D and 10 for 5D, which has a magnitude in \mathcal{O}({10}^{-11}) and relative errors about 10%. The 5D L=12 data point has a higher relative error as shown in Fig. 11. The fitting results by Eq. (5) are shown in Table 3. From these fits we take our final estimate to be X_{\text{3a}}(4D) = 6.9(2) and X_{\text{3a}}(\text{5D})=8.4. In 5D, since the number of data points is too small, we are not able to give a robust error bar for exponent X_{\text{3a}}. In Fig. 11, we plot P_{\text{3a}}{L}^{2X_{\text{3a}}} vs. L^-3 for 4D using three different values of X_{\text{3a}} in the upper panel, and P_3a vs. L directly for 5D in the lower panel.

6 Discussion

We study the N-cluster correlation functions for critical percolation in 4D and 5D. We reformulate the transfer MC algorithm in Ref. [36] by utilizing a disjoint-set data structure, and carry out extensive Monte Carlo simulations for cylindrical systems of size L^{d-1}\times \infty. The linear system size is up to L=512 in 4D and 128 in 5D. From finite-size scaling, we report the estimates of all the N=2 and 3 exponents and the universal logarithmic amplitude δ. Table 4 summarizes the d-dependence of these exponents and the amplitude δ.

For N=1, the only correlation function describes the probability that two distant sites are in the same percolation cluster, and the exponent X₁ is the magnetic dimension X_{\text{h}}=d-y_{\text{h}}, which has been extensively studied in Refs. [20, 47, 52, 53]. The symmetric exponent X_{\text{2s}} is the so-called two-arm exponent. For percolation, i.e., the Q\to 1 limit of the Q-state Potts model, it is identical to the thermal scaling dimension X_{\text{t}}=d-y_{\text{t}}. Thus, our estimate of X_{\text{2s}} also gives a determination of y_{\text{t}} for 4D and 5D, which significantly improves over the existing results. To our knowledge, the other exponents for N=2 and 3 and the universal amplitude δ have not been reported for 4D and 5D.

For spatial dimensions d>2, exact results are unavailable for critical phenomena, mainly due to the finiteness of the conformal group. The approach from logarithmic conformal field theory (LCFT) focuses on nonlocal observables rather than local ones. Combining the symmetry group {\mathcal{S}}_N of N clusters and the {\mathcal{S}}_Q symmetry for the Q-state Potts model, the LCFT approach gives exact structural properties on nonlocal connectivities and defines an infinite family of N-cluster correlation functions that have symmetry-dependent critical exponents. Together with Ref. [30], our work confirms the validity of LCFT and produces one of the scarce pieces of knowledges for high-dimensional percolation with [54].

Above the upper dimensionality d\ge d_{\text{u}}=6, the critical behavior is governed by the Gaussian fixed point in the framework of renormalization group. This implies that, at criticality, the random path connecting a pair of distant sites along a percolation cluster is effectively a simple random walk in the continuum limit. As a consequence, an N-cluster connectivity event, that there are N distinct clusters propagating from a neighborhood {\mathcal{V}}_i to another one far away {\mathcal{V}}_j, can be regarded as the event that there are N simple random walks between {\mathcal{V}}_i and {\mathcal{V}}_j. Since the random walks are independent to each other, we speculate that for d\ge 6, the symmetric exponent simply takes value X_{N\text{s}}=2N, supported by the known results X_1=2 and X_{\text{2s}}=4. Further, we speculate that, in the continuum limit, the cancellation of connectivity probabilities in correlators P_{N\text{o}} would be so complete that the N-cluster exponents X_{N\text{o}}, except X_{N\text{s}}, are ill-defined, denoted as symbol “/" in Table 4. To confirm/falsify these speculations, theoretical insights are desired since numerical study is extremely challenging.

Given a spatial dimensionality , we conjecture that the symmetric exponents are superadditive, i.e., X_{N\text{s}}>X_{N^{\prime }\text{s}}+X_{N^{″}\text{s}} for any set of positive integers with N=N^{\prime }+N^{″}. This concavity property immediately yields another inequality X_{N\text{s}}>N{X}_1. The conjecture holds true for N=\mathrm{1,2,3} and d=\mathrm{2,3,4,5}, as shown in Table 4. Moreover, the ratio X_{\text{2s}}/X_1 and X_{\text{3s}}/X_1 increase as d decreases. The underlying argument is as follows. As d is lowered from d_{\text{u}}=6, critical percolation clusters become more and more compact [21, 35, 51, 53]. Consider the process of simultaneously growing N distinct clusters from a neighborhood of N lattice sites, the clusters are very likely to merge together if they do not die out, due to their geometric compactness/fatness. Thus, in the unlikely event that all the N clusters survive up to distance r, they would separate from each other as r increases. This leads to our expectation that {\mathbb{P}}_N\le {\mathbb{P}}_{N^{\prime }}\cdot {\mathbb{P}}_{N^{″}} for any positive integers N=N^{\prime }+N^{″}.

Finally, we remark on the d-dependent behavior of the universal logarithmic amplitude δ. Unlike the N-cluster exponents X_{N\text{o}} that are monotonically increasing functions of d, the δ value starts from 2\sqrt{3}/\pi \approx 1.103 in 2D, reaches a maximum 1.52(3) in 3D, and then decreases to 1.16(1) for 4D and 0.74(6) for 5D. The universal amplitude δ characterizes both the LCFT at Q=1 and the limit of conformal field theories when Q\to 1, and it is predicted that [27]

\begin{array}{}\text{(8)}& \delta =2\times \underset{Q\to 1}{{\displaystyle \lim }}\frac{X_2-X_{\text{t}}}{Q-1}\mathrm{,}\end{array}

where the two-arm exponent X₂ collides with the thermal exponent X_{\text{t}} in the Q\to 1 limit. On the complete graph, it has been rigorously proved [55] that as long as , the Fortuin–Kasteleyn cluster representation of the Q-state Potts model belongs to the same universality class as the mean-field percolation. We speculate that this would hold for dimension d>6. Thus, one has X_2=X_{\text{t}} and \delta =0 for d\ge 6, in line with the decreasing tendency of δ as d increases. The value \delta =0 might indicate that for d\ge 6, the logarithmic correlation function F(r) no longer has logarithmic r-dependence.

References

[1] Broadbent S. R. and Hammersley J. M., Percolation processes, Proc. Camb. Philos. Soc. 53(3), 629 (1957)

[2] Stauffer D. and Aharony A., Introduction to Percolation Theory, 2nd Ed., London: Taylor Francis, 1994

[3] Grimmett G. R., Percolation, 2nd Ed., Berlin: Springer, 1999

[4] Bollobás B. and Riordan O., Percolation, Cambridge: Cambridge University Press, 2006

[5] Kasteleyn P. W. and Fortuin C. M., Phase transitions in lattice systems with random local properties, J. Phys. Soc. Jpn. 26(Suppl.), 11 (1969)

[6] Potts R. B., Some generalized order-disorder transformations, Proc. Camb. Philos. Soc. 48(1), 106 (1952)

[7] Wu F. Y., The Potts model, Rev. Mod. Phys. 54(1), 235 (1982)

[8] Kramers H. A. and Wannier G. H., Statistics of the two-dimensional ferromagnet (Part I), Phys. Rev. 60(3), 252 (1941)

[9] Lieb E. H., Exact solution of the problem of the entropy of two-dimensional ice, Phys. Rev. Lett. 18(17), 692 (1967)

[10] Baxter R. J., Partition function of the eight-vertex lattice model, Ann. Phys. 70(1), 193 (1972)

[11] Belavin A. A., Polyakov A. M., and Zamolodchikov A. B., Infinite conformal symmetry in two-dimensional quantum field theory, Nucl. Phys. B 241(2), 333 (1984)

[12] Friedan D., Qiu Z., and Shenker S., Conformal invariance, unitarity, and critical exponents in two dimensions, Phys. Rev. Lett. 52(18), 1575 (1984)

[13] Nienhuis B., in: Phase Transition and Critical Phenomena, edited by Domb C., Green M., and Lebowitz J. L., Academic Press, London, 1987, Vol. 11

[14] Cardy J. L., in: Phase Transition and Critical Phenomena, edited by Domb C., Green M., and Lebowitz J. L., Academic Press, London, 1987, Vol. 11

[15] Lawler G. F., Schramm O., and Werner W., The dimension of the planar Brownian frontier is 4/3, Math. Res. Lett. 8(4), 401 (2001)

[16] Smirnov S. and Werner W., Critical exponents for two-dimensional percolation, Math. Res. Lett. 8(6), 729 (2001)

[17] Aharony A., Gefen Y., and Kapitulnik A., Scaling at the percolation threshold above six dimension, J. Phys. A 17, L197 (1984)

[18] Hara T. and Slade G., Mean-field critical behaviour for percolation in high dimensions, Commun. Math. Phys. 128(2), 333 (1990)

[19] Fitzner R. and van der Hofstad R., Mean-field behavior for nearest-neighbor percolation in d > 10, Electron. J. Probab. 22(0), 22 (2017)

[20] Gracey J. A., Four loop renormalization of ϕ³ theory in six dimensions, Phys. Rev. D 92(2), 025012 (2015)

[21] Wang J. F., Zhou Z. Z., Zhang W., Garoni T. M., and Deng Y. J., Bond and site percolation in three dimensions, Phys. Rev. E 87(5), 052107 (2013)

[22] Xu X., Wang J. F., Lv J. P., and Deng Y. J., Simultaneous analysis of three-dimensional percolation models, Front. Phys. 9(1), 113 (2014)

[23] Paul G., Ziff R. M., and Stanley H. E., Percolation threshold, Fisher exponent, and shortest path exponent for four and five dimensions, Phys. Rev. E 64(2), 026115 (2001)

[24] Stanley H. E., in: Percolation Theory and Ergodic Theory of Infinite Particle Systems, edited by Kesten H., IMA Volumes in Mathematics and Its Applications Vol. 8, New York: Springer-Verlag, 1987

[25] Beffara V. and Nolin P., On monochromatic arm exponents for 2D critical percolation, Ann. Probab. 39(4), 1286 (2011)

[26] Aizenman M., Duplantier B., and Aharony A., Path-crossing exponents and the external perimeter in 2D percolation, Phys. Rev. Lett. 83(7), 1359 (1999)

[27] Vasseur R., Jacobsen J. L., and Saleur H., Logarithmic observables in critical percolation, J. Stat. Mech.: Theory Exp. 07, L07001 (2012)

[28] Vasseur R. and Jacobsen J. L., Operator content of the critical Potts model in d dimensions and logarithmic correlations, Nucl. Phys. B 880, 435 (2014)

[29] Couvreur R., Lykke Jacobsen J., and Vasseur R., Non-scalar operators for the Potts model in arbitrary dimension, J. Phys. A Math. Theor. 50(47), 474001 (2017)

[30] Tan X. J., Couvreur R., Deng Y. J., and Jacobsen J. L., Observation of nonscalar and logarithmic correlations in two- and three-dimensional percolation, Phys. Rev. E 99(5), 050103 (2019)

[31] Gurarie V. and Ludwig A. W. W., Conformal field theory at central charge c = 0 and two-dimensional critical systems with quenched disorder, arXiv: hep-th/0409105 (2004)

[32] Mathieu P. and Ridout D., From percolation to logarithmic conformal field theory, Phys. Lett. B 657(1–3), 120 (2007)

[33] Vasseur R., Jacobsen J. L., and Saleur H., Indecomposability parameters in chiral logarithmic conformal field theory, Nucl. Phys. B 851(2), 314 (2011)

[34] Gurarie V. and Ludwig A. W. W., Conformal algebras of two-dimensional disordered systems, J. Phys. Math. Gen. 35(27), L377 (2002)

[35] Huang W., Hou P. C., Wang J. F., Ziff R. M., and Deng Y. J., Critical percolation clusters in seven dimensions and on a complete graph, Phys. Rev. E 97(2), 022107 (2018)

[36] Deng Y. J. and Blöte H. W. J., Monte Carlo study of the site-percolation model in two and three dimensions, Phys. Rev. E 72(1), 016126 (2005)

[37] Hoshen J. and Kopelman R., Percolation and cluster distribution (I): Cluster multiple labeling technique and critical concentration algorithm, Phys. Rev. B 14(8), 3438 (1976)

[38] Galler B. A. and Fisher M. J., An improved equivalence algorithm, Commun. ACM 7(5), 301 (1964)

[39] Tarjan R. E. and Van Leeuwen J., Worst-case analysis of set union algorithms, J. Assoc. Comput. Mach. 31(2), 245 (1984)

[40] Tarjan R. E., A class of algorithms which require nonlinear time to maintain disjoint sets, J. Comput. Syst. Sci. 18(2), 110 (1979)

[41] Newman M. E. J. and Ziff R. M., Efficient Monte Carlo algorithm and high-precision results for percolation, Phys. Rev. Lett. 85(19), 4104 (2000)

[42] Newman M. E. J. and Ziff R. M., Fast Monte Carlo algorithm for site or bond percolation, Phys. Rev. E 64(1), 016706 (2001)

[43] Danziger M. M., Gross B., and Buldyrev S. V., Faster calculation of the percolation correlation length on spatial networks, Phys. Rev. E 101(1), 013306 (2020)

[44] Blöte H. W. J. and Nightingale M. P., Critical behaviour of the two-dimensional Potts model with a continuous number of states; A finite size scaling analysis, Physica A 112(3), 405 (1982)

[45] Koza Z. and Poła J., From discrete to continuous percolation in dimensions 3 to 7, J. Stat. Mech.: Theory Exp 2016(10), 103206 (2016)

[46] Borinsky M., Gracey J. A., Kompaniets M., and Schnetz O. (in preparation) (2020)

[47] Mertens S. and Moore C., Percolation thresholds and Fisher exponents in hypercubic lattices, Phys. Rev. E 98(2), 022120 (2018)

[48] Xun Z. P. and Ziff R. M., Precise bond percolation thresholds on several four-dimensional lattices, Phys. Rev. Research 2(1), 013067 (2020)

[49] Cardy J., Finite-size Scaling, Vol. 2, Elsevier, 2012

[50] Stanley H. E., Cluster shapes at the percolation threshold: and effective cluster dimensionality and its connection with critical-point exponents, J. Phys. Math. Gen. 10(11), L211 (1977)

[51] Xu X., Wang J. F., Zhou Z. Z., Garoni T. M., and Deng Y. J., Geometric structure of percolation clusters, Phys. Rev. E 89(1), 012120 (2014)

[52] Ballesteros H. G., Fernández L. A., Martín-Mayor V., Muñoz Sudupe A., Parisi G., and Ruiz-Lorenzo J. J., Measures of critical exponents in the four-dimensional site percolation, Phys. Lett. B 400(3–4), 346 (1997)

[53] Zhang Z. J., Hou P. C., Fang S., Hu H., and Deng Y. J., Critical exponents and universal excess cluster number of percolation in four and five dimensions, arXiv: 2004.11289 (2020)

[54] Araújo N. A. M., Grassberger P., Kahng B., Schrenk K. J., and Ziff R. M., Recent advances and open challenges in percolation, Eur. Phys. J. Spec. Top. 223, 2307 (2014)

[55] Luczak M. and Luczak T., The phase transition in the cluster-scaled model of a random graph, Random Structures Algorithms 28(2), 215 (2006)

[56] Guttmann A. J. and Jacobsen J. L., Lattice models and integrability: A special issue in honour of F. Y. Wu, J. Phys. A Math. Theor. 45(49), 490301 (2012)

Rights & permissions

Higher Education Press 2020