Percolation of sites not removed by a random walker in $d$ dimensions

TL;DR

This paper investigates how a random walk's site removal process affects percolation in high-dimensional lattices, revealing a sharp percolation threshold in dimensions above two and continuous behavior at the lower critical dimension.

Contribution

It introduces a detailed study of correlated site percolation caused by random walks, identifying the percolation threshold and critical exponents across dimensions 2 to 6.

Findings

Percolation probability approaches a step function as system size increases.

Percolation threshold u_c is approximately 3 for dimensions 3 to 6.

Correlation length diverges with exponents consistent with ν=2/(d-2).

Abstract

How does removal of sites by a random walk lead to blockage of percolation? To study this problem of correlated site percolation, we consider a random walk (RW) of $N=uL^d$ steps on a $d$-dimensional hypercubic lattice of size $L^d$ (with periodic boundaries). We systematically explore dependence of the probability $\Pi_d(L,u)$ of percolation (existence of a spanning cluster) of sites not removed by the RW on $L$ and $u$. The concentration of unvisited sites decays exponentially with increasing $u$, while the visited sites are highly correlated -- their correlations decaying with the distance $r$ as $1/r^{d-2}$ (in $d>2$). Upon increasing $L$, the percolation probability $\Pi_d(L,u)$ approaches a step function, jumping from 1 to 0 when $u$ crosses a percolation threshold $u_c$ that is close to 3 for all $3\le d\le6$. Within numerical accuracy, the correlation length associated with percolation diverges with exponents consistent with $\nu=2/(d-2)$. There is no percolation threshold at the lower critical dimension of $d=2$, with the percolation probability approaching a smooth function $\Pi_2(\infty,u)>0$.

Equations

$$N_{\mathrm{dist}}=A_{d}N,$$

$$p=\prod _i{p}_i=\exp \left(-\sum _i\frac{n_i}{L^d}\right)$$

$$=\exp \left(-\frac{N_{\mathrm{dist}}}{L^d}\right)=\exp (-A_{d}u).$$

$$⟨v(\vec{x})v(\vec{y}){⟩}_c\sim 1/|\vec{x}-\vec{y}{|}^{d-2}.$$

$$\nu =2/(d-2),\mathrm{for}3\le d\le 6.$$

$$|u_{\ast }-u_c|\sim L^{-1/\nu }.$$

$$\mathrm{slope}\approx (L/b{)}^{1/\nu }.$$

$$N_{\mathrm{dist}}=A_{2}N,\mathrm{with}{A}_2=\pi /\ln N.$$

$$p=\exp \left(-\frac{\pi u}{2\ln L}\right).$$

$${\Pi }_2\approx \mathrm{Erf}(1/\sqrt{u}).$$

$$P(x_0|x,t)=\frac{2}{L}\sin \left(\frac{\pi x_0}{L}\right)\sin \left(\frac{\pi x}{L}\right)\exp (-D{\pi }^{2}N/L^2),$$

Full text

I Introduction

In the simplest (Bernoulli) site or bond percolation problem Stauffer and Aharony (1991); Gremmett (1999) sites or bonds of a regular $d$-dimensional lattice are independently occupied with probability $p$. For an infinite system, there is a sharp percolation transition point $p_{c}$, such that for $p>p_{c}$ there exists an infinite cluster spanning the system. Close to $p_{c}$, many geometrical and physical properties become singular, as expressed by universal power-law dependencies on $|p-p_{c}|$. For example, the typical linear extension of finite clusters of connected sites, indicated by the correlation length $\xi$, diverges as $\xi\sim|p-p_{c}|^{-\nu}$. The universal critical exponent $\nu$ depends only on space dimension $d$, and is well known for Bernoulli percolation in all $d$ (see, e.g., Ref. Koza and Poa (2016)). A much studied problem in mathematical literature, percolation has also been used to model a variety of physical systems for the onset of connectivity and flow, e.g., for current passing through a random resistor network.

An early application of percolation is to gels formed by random crosslinking of polymers Adam et al. (1990). Gelation has acquired new interest in the context of reversible accumulations of nonspecific biological molecules into liquid like droplets with important functions as in transcription regulation Sabari et al. (2018). The reverse process of gel degradation is now also of relevance. In principle, the removal of connections, rather than their addition, does not qualitatively change the percolation picture. For example, in the process of hydrogel degradation Metters et al. (2000a, b, 2001); Jahanmir et al. (2018) connections are severed mostly uniformly in space, and the measured elastic and rheological properties Diederich et al. (2017); Schultz et al. (2012) resemble the gel formation process in reverse Larsen and Furst (2008).

While in Bernoulli percolation, the elements are added or removed randomly and independently, new behavior emerges if subsequent events are correlated. An extreme case is that of explosive percolation when the events are specifically chosen so as to delay the percolation transition DeSouza and Nagler (2015). A different extreme consists of removing entire straight lines from the system in each step Kantor (1986), e.g., by drilling through a solid sample Hilário et al. (2015), or by having very elongated obstacles (fibers) influencing molecular diffusion between cells Gomez et al. (2019), and also corresponds to a completely different universality class Hilário et al. (2015); Schrenk et al. (2016); Grassberger (2017); Grassberger et al. (2017). A variant of the latter is removal of sites or bonds performed by a meandering random walk (RW). This models a simplified version of a degradation process in which a single enzyme, or possibly a few enzymes, travel through a gel, breaking the crosslinks they encounter Berry et al. (2000); Fadda et al. (2003). In an early numerical study of this problem, which they named “random walk decay,” Banavar et al. Banavar et al. (1985) considered properties of the clusters of vacant sites, unvisited by the RW, on square and cubic lattices. In a later study, Abete et al. considered percolation of the vacant bonds on a cubic lattice Abete et al. (2004), finding numerically the threshold and several critical indices of the problem, which they called “pacman percolation.”

Independently, the mathematical community has also studied aspects of the above problem. While physicists focused on geometry and critical properties of clusters, primarily in $d=2$ and $d=3$ dimensions, mathematicians debated the very existence of a percolation threshold. For a random walk of length $N$ to cover a finite fraction of the $M=L^{d}$ sites of a hypercubic lattice, the number of steps $N$ must be proportional to $M$. The question of whether percolation of vacant sites stops for $N/M=u$ larger than a critical $u_{c}$ was addressed in several mathematical works, first for $d\geq 7$ Sznitman (2008) and later for $d\geq 3$ Sidoravicius and Sznitman (2009); Sznitman (2010). In the mathematical literature, this problem is referred to as “percolation of vacant sites of interlacement.” Although, the bounds on $u_{c}$ established in this literature were extremely broad, and in $d=3$ spanned many order of magnitude Balázs (2015), we shall see that $u_{c}$ is approximately 3 for $3\leq d\leq 6$.

More precisely, in this work we perform a detailed numerical study of site percolation on $d$-dimensional hypercubic lattices of linear dimensions $L$ ($0\leq~{}\{x_{1},x_{2},\dots,x_{d}\}\leq~{}L-1$) and volume $V=a^{d}M=L^{d}$ (with lattice constant $a=1$). A RW starts at an arbitrary initial position on the lattice and performs $N=uM=uL^{d}$ steps, obeying periodic boundary conditions in all directions, i.e., $x_{i}=L$ coincides with $x_{i}=0$ for $i=1,\dots,d$. This scaling of $N$ maintains a fixed fraction of vacant sites in $d\geq 3$. (The case of $d=2$ is discussed separately in Sec. VI.) We then define a configuration as “spanning” (“percolating”) if a continuous path of vacant sites exists between boundaries at $x_{d}=0$ and $x_{d}=L-1$. Note that while checking for this percolation condition, boundaries in directions $1,\dots,d-1$ are assumed to be periodic; this is sometimes referred to as “helical boundary conditions” Acharyya and Stauffer (1998). By considering a large number of realizations for each $L$ and $u$, we determine the spanning probability $\Pi_{d}(L,u)$. The fraction of vacant sites $p$ is a monotonically decreasing function of $u$. In $d\geq 3$ there is a sharp percolation threshold $u_{c}$, such that in the $L\to\infty$ limit the spanning probability becomes a step function with $\Pi_{d}=1$ for $u<u_{c}$ and $\Pi_{d}=0$ for $u>u_{c}$. (And in terms of $p$, this corresponds to a threshold $p_{c}$.) The relation between $u$ and $p$ for $d\geq 3$ is discussed in detail in Sec. II.

Unlike Bernoulli percolation, the sites removed by a RW are highly correlated. In Sec. III we demonstrate that the correlations between vacant sites decay as a power-law, proportional to $1/r^{d-2}$ with their separation $r$. Such correlated percolation has been argued to belong to a universality class characterized by an exponent $\nu=2/(d-2)$ for divergence of the correlation length Weinrib and Halperin (1983); Weinrib (1984). We set up to test this behavior in dimensions $3\leq d\leq 6$. Section IV details our numerical study of percolation in $d=3$, while Sec. V extends the exploration to $d=4$, 5, and 6. There is no percolation threshold in $d=2$, resulting in a smooth limiting function for spanning probability as demonstrated in Sec. VI. The limits of this function for very short and very long walks are discussed in Sec. VII. We conclude in Sec. VIII with some topics for future investigation.

II Statistics of sites not visited by a random walk

As a random walker performs $N=uL^{d}$ steps on a finite hypercubic lattice of volume $M=L^{d}$ (“box”) with periodic boundary conditions, the fraction of unvisited or vacant sites $p$ decreases with increasing $u$. If instead of performing an $N$-step RW, we randomly and independently select $N$ of sites on a lattice, the occupied sites would be random, uncorrelated, with some of them multiply occupied. More precisely, the fraction of vacant sites would be $p=\exp(-u)$, as indicated by the bottom dash-dotted line in Fig. 1, with their positions uncorrelated as in usual (Bernoulli) percolation. However, the obvious correlations in the positions of the random walker generate a different dependence $p(u)$ as $L\to\infty$, which also remains nonzero for any finite $u$.

The fractal dimension of a RW is 2, leading to qualitatively different behaviors in $d=2$ that are separately discussed in Sec. VI. The results described in this section, as well as in Sec. III, thus pertain only to dimensions $d\geq 3$. While the RW on an infinite hypercubic lattice does intersect itself, it can be shown rigorously that the number of distinct visited sites of a long walk of $N$ steps grows as Vineyard (1963); Montroll and Weiss (1965); Rubin and Weiss (1982)

[TABLE]

with subleading corrections of $O(N^{1/2})$ in $d=3$, and $O(\ln N)$ in $d=4$. The corrections do not increase with $N$ for $d\geq 5$, as self-intersections of remote parts of RWs become negligible. The coefficient $A_{d}$ in Eq. (1) depends on lattice type and space dimension. It is the inverse of the mean number $B_{d}$ of visits of a random walker to its initial position Aldous (1983); Nemirovsky et al. (1990); Brummelhuis and Hilhorst (1991) (see also chapter 3 in Ref. Hughes (1995)). [Also, the mean number of RW steps required to visit all sites of a finite box of size $M$ is $B_{d}M\ln M$ Brummelhuis and Hilhorst (1991) (see also Ref. Weiss et al. (1985)).] A detailed procedure for calculation of both $A_{d}$ and subleading corrections for various lattices was outlined by Montroll and Weiss Montroll and Weiss (1965), while an efficient numerical method can be found in Ref. Griffin (1990). The values of these coefficients for hypercubic lattices are $A_{3}=0.659$ Montroll and Weiss (1965), $A_{4}=0.807$, $A_{5}=0.865$ and $A_{6}=0.895$ Griffin (1990). (Here and thereafter, the accuracy of numbers without error bars is a single unit of the last digit or better.)

Sites visited by a RW are strongly correlated. In an infinite space, an $n$-step RW explores a volume of radius $r\sim n^{1/2}$ in which, for $d\geq 3$, the density of visited sites will be $n/r^{d}\sim~{}1/r^{d-2}$. Thus, if a particular lattice site on a hypercubic lattice is visited by a RW, then the probability of finding another visited site at a large separation $r$ will be $C_{d}/r^{d-2}$, where the lattice- and dimension-specific prefactor $C_{d}$ is of order of unity. Consider a random variable $w(\vec{x})$ which is 1, if a site at position $\vec{x}$ has been visited by a RW, and is 0 otherwise. Probabilities, of RW visiting positions $\vec{x}$ and $\vec{y}$ are correlated, such that $\langle w(\vec{x})w(\vec{y})\rangle/\langle w(\vec{x})\rangle\approx C_{d}/|\vec{x}-\vec{y}|^{d-2}$, where $\langle\rangle$ denotes an average over realizations of RWs.

To examine the problem on a finite hypercube with periodic boundary conditions, we first create a RW on an infinite lattice, then tile the space with boxes (hypercubes) of size $M=L^{d}$. The boxes are then cut out and superposed; a procedure that we refer to as “folding.” For $d\geq 3$, and for finite $u=O(1)$ and large $L$, the walk on an infinite lattice will have typical size (end-to-end distance) $\sim L^{d/2}\gg L$. In an infinite space most of the boxes will be empty, and the RW will visit of order $L^{d-2}$ distinct boxes. In the $i$th visited cube there will be $n_{i}$ distinct visited sites, and $n_{i}$ will be of order $L^{2}$. Thus, the fraction of unvisited (vacant) sites in a cube is $p_{i}=1-n_{i}/L^{d}\simeq\exp(-n_{i}/L^{d})$, as $n_{i}/L^{d}\sim 1/L^{d-2}$ is very small. The “folding” process from infinite space into a single periodic box effectively removes correlations between far away boxes, and we can treat the configurations arising from superposition of distinct boxes as uncorrelated. Therefore, the probability that a particular point in the periodic box has not been visited is

[TABLE]

This asymptotic (very large $L$) result has been rigorously proven in Ref. Brummelhuis and Hilhorst (1991).

In numerical simulations with moderate values of $L$ ($\sim$10–100), the finite-$N$ corrections to Eq. (1) are noticeable, especially in $d=3$ where there is a non-negligible probability for a RW that exited a box (on the infinite lattice) to return to it. Yet correlations of occupied sites even in adjacent boxes do not exceed $C_{d}/L^{d-2}$, and are much smaller for nonadjacent boxes (which are the majority). Thus, in the “folding” process we superimpose practically independent configurations. Also, the fraction of visited sites $n_{i}/L^{d}\sim 1/L^{d-2}$ is still significantly smaller than 1. The deviations from Eq. (II) for moderate $L$ can be corrected for by replacing $A_{d}$ with an effective $A_{d}(L)$. We examined numerically the functions $p(u)$ in $d=3$, 4, 5, and 6, and in all cases the results for various $L$ could be fitted extremely well by a pure exponential $p=\exp[-A_{d}(L)u]$. For $d=3$ we checked such dependence for $L=4$, 8,…,128. The results presented in Fig. 1 fit pure exponentials at the accuracy level of $\chi^{2}\sim 10^{-7}$. The slopes on the semilogrithmic plots depend on $L$ but converge very fast to the known value of $A_{3}$ indicated by the top line on Fig. 1. Of course, both this and the following statements, about “purely exponential behavior,” should not be taken in a strict mathematical sense: We know, that for a finite periodic box of size $M$, a walk of $B_{d}M\ln M$ steps will “completely occupy” the box. So, in $d=3$ and $L=128$, this will occur for $u>B_{3}\ln 128^{3}\approx 22$, far beyond the limits of applicability of the above discussion.

In $d\geq 4$ the distant parts of RW rarely intersect and therefore $A_{d}(L)$ converges very fast to its asymptotic value $A_{d}$. As an example, in Fig. 2 we show the $d=4$ case, where for lattices of sizes $L=4,8,16$, and 32 we generate $10^{4}$ configurations per data point to measure the relation between $u$ controlling the chain length and the fraction of unvisited sites $p$. As before, all curves are straight on semi-log scale, i.e., $p$ is an exponentially decaying function of $u$. The rate of that decay quickly converges to a constant, leading to $p=\exp(-A_{4}u)$ for large $L$.

We repeated these calculations also for $d=5$ and $d=6$. There was practically no $L$ dependence of $A_{d}(L)$ and the coefficients reached their asymptotic values already for small $L$.

III Correlated percolation

To quantify the correlations between sites visited by the RW following the “folding” from infinite space, let us examine the two-point correlation function $\langle w(\vec{x})w(\vec{y})\rangle$. We argued previously that the RW segments from distinct boxes are almost uncorrelated, after being “folded” into the single box with periodic boundary conditions. A randomly selected site $\vec{x}$ in the periodic system has a probability $\langle w(\vec{x})\rangle=1-p=1-\exp(-A_{d}u)$ of being occupied. This point may have been visited a few times either by the same RW segment, or by RW segments belonging to distinct boxes before “folding.” Either way, the probabilities of finding a site at position $\vec{y}$ that has been visited by the segment of RW in the same box (before “folding”) is given by $C_{d}/|\vec{x}-\vec{y}|^{d-2}$, and it is modified by the number of distinct boxes to which the segments belonged before “folding.” Thus, the correlated part of $\langle w(\vec{x})w(\vec{y})\rangle$ will decay as $1/|\vec{x}-\vec{y}|^{d-2}$, as long as $|\vec{x}-\vec{y}|$ is significantly smaller than the box size $L$. Of course, the chances that position $\vec{y}$ also has been visited is dominated by segments of RW that belong to other boxes before the “folding.” This constant background part must be subtracted when measuring the cumulant, leading to $\langle w(\vec{x})w(\vec{y})\rangle_{c}=F(u)/|\vec{x}-\vec{y}|^{d-2}$, where $F(u)$ is a smooth function of order unity for $u=O(1)$. (In $d=3$, an equivalent argument can be found in Ref. Fadda et al. (2003).)

Since we are interested in the properties of vacant sites, we define a related variable $v(\vec{x})=1-w(\vec{x})$, with $\langle v(\vec{x})\rangle\equiv p$. Products of these variables at distinct positions have the property $\langle v(\vec{x})[v(\vec{y})+w(\vec{y})]\rangle=\langle v(\vec{x})\rangle=p$ and $\langle w(\vec{x})[v(\vec{y})+w(\vec{y})]\rangle=\langle w(\vec{x})\rangle=1-p$. Since our system is translationally invariant and symmetric under inversion, we have $\langle v(\vec{x})w(\vec{y})\rangle=\langle w(\vec{x})v(\vec{y})\rangle$. Thus, by subtracting the two previous equalities from each other, we find $\langle v(\vec{x})v(\vec{y})\rangle=\langle w(\vec{x})w(\vec{y})\rangle+2p-1$. Consequently, the vacant sites exhibit the same decay of cumulants as visited sites, i.e.,

[TABLE]

The two-point correlation function captures only one aspect of the interesting information about the system. The RW in a periodic box has the distinct property of being a single cluster (noting periodic boundary conditions). While vacant sites can form many clusters, they tend to aggregate into one large cluster. The unusual properties of both small and “infinite” clusters have been studied by several authors Banavar et al. (1985); Teixeira (2009, 2011); Benjamini and Sznitman (2008); Abete et al. (2004).

In “usual” (Bernoulli) percolation there are no correlations between occupied sites or bonds. This problem has a lower critical dimension of $d=1$, where $p_{c}=1$, and an upper critical dimension of $d_{c}=6$ Toulouse (1974), above which mean-field behavior is expected; e.g., with the correlation length diverging with exponent $\nu_{\rm B}=1/2$. Weinrib analyzed stability of the Bernoulli percolation universality class to correlations Weinrib (1984) (following a similar treatment for critical phase transitions Weinrib and Halperin (1983)). By appealing to a generalized Harris criterion Harris (1974), he demonstrated that short-range correlations, as well as power-law correlations decaying as $1/r^{a}$ with $a>d$, do not modify the universality class of Bernoulli percolation. However, for $a<d$, the relevance of the correlations is determined by the extended Harris criterion Weinrib and Halperin (1983): If $a\nu_{\rm B}-2<0$, then the correlations are relevant. Equation (3) shows that vacant sites have a power law correlation with $a=d-2$. The quantity $(d-2)\nu_{\rm B}-2$ is $-1.12$, $-0.62$, $-0.29$ and 0, for $d=3,4,5$, and 6, respectively Koza and Poa (2016). This expression becomes positive for $d\geq 7$.

Percolation with long-range correlations is a well-researched subject Coniglio and Fierro (2009). (For a more recent study see Ref. Gori et al. (2017).) In most physical contexts, correlations are generated by thermodynamic systems such a critical Ising model. It has been shown Weinrib and Halperin (1983); Weinrib (1984) that, for power-law correlations, the exponent characterizing divergence of the correlation length equals $\nu=2/a$, which in our case is

[TABLE]

Abete et al. Abete et al. (2004) studied percolation of vacant bonds on lattices of up to $60^{3}$ and found $\nu=1.8\pm 0.1$, consistent with $\nu=2$ expected from Eq. (4). In this work we consider site percolation on large lattices, and confirm Eq. (4) in $d=3,4$, and 5. In $d=d_{c}=6$ power-law correlations $\sim 1/r^{4}$ represent a marginal perturbation, and both Eq. (4) and standard Bernoulli percolation lead to $\nu=1/2$.

In a system of infinite size the spanning probability is a step function, jumping between 1 and 0 as the percolation threshold is crossed. However, for finite $L$ it exhibits a smooth crossover between these values. It has been shown for Bernoulli percolation Langlands et al. (1992); Cardy (1992); Langlands et al. (1994); Watts (1996) that for large $L$ the value of spanning probability $\Pi_{d,{\rm B}}(L,p)$ at $p=p_{c}$ reaches a universal value independent of microscopic lattice details or the consideration of bond or site percolation. Thus, the critical spanning probability $\Pi^{c}_{d,{\rm B}}\equiv\lim_{L\to\infty}\Pi_{d,{\rm B}}(L,p_{c})$ is a universal number characterizing a percolation universality class. For example, $\Pi^{c}_{2,{\rm B}}=1/2$ in two-dimensional Bernoulli site and bond percolation with free boundary conditions. However, the value of the spanning probability for finite $L$, and, consequently, its universal limit of $\Pi^{c}_{d,{\rm B}}$, does depend on the macroscopic definition of spanning, such as requiring spanning in several directions simultaneously or using periodic versus free boundary conditions. Such dependence on macroscopic definitions has been observed by several authors Stauffer et al. (1994); Hovi and Aharony (1996); Lin et al. (1998); Acharyya and Stauffer (1998). The value of $\Pi^{c}$ thus may serve as an extra indicator of differences between universality classes. In the problem of percolation of vacant sites of a RW, we can similarly define $\Pi_{d}^{c}\equiv\lim_{L\to\infty}\Pi_{d}(L,u=u_{c})$ characterizing this type of percolation.

IV Percolation in $d=3$

In $d=3$ we considered site percolation along the $x_{3}$ direction for cubic lattice sizes $L=4,8,16,\dots,512$. For every $L$ and $u$ we generated many realizations of RWs of $N=uL^{3}$ steps. Spanning probability $\Pi_{3}(L,u)$ was calculated by averaging $10^{4}$ realizations for each $u$ for $L\leq 64$, and 4000 configurations for each $u$ for $L\geq 128$. The limiting factor in the computations was the large lattice size, and correspondingly long walks reaching $N=5\times 10^{5}$, necessitating long times required to process each configuration. (The previous study by Abete et al. Abete et al. (2004) reached lattice sizes $L=60$ and considered bond percolation.) Figure 3 depicts the spanning probability as a function of $u$. As expected, the transition becomes sharper with increasing $L$. The curves also exhibit a very strong drift toward larger values of $u$ with increasing $L$.

There are several methods to determine the critical value $u_{c}$ for percolation. The low accuracy of the results prevents us from using points $u=u_{*}$ of maximal slope of $\Pi_{d}(L,u)$ as estimates of $u_{c}$. However, we can examine the $L$-dependence of values $u_{*}$ for which $\Pi_{d}(L,u_{*})=c$. Independently of the choice of $c$, we expect $\lim_{L\to\infty}u_{*}=u_{c}$. Moreover, for large $L$, we expect

[TABLE]

For $d=3$, Fig. 4 depicts the dependence of successive estimates of $u_{c}$ on $L^{-1/2}$, as suggested by Eq. (5) with anticipated $\nu\approx 2$. (Note that for small $L$, the estimates of $u_{c}$ do not follow the asymptotic form of Eq. (5), and are sometimes even nonmonotonic functions of $L$.) All four lines extrapolate to $u_{c}=3.15\pm 0.01$.

An additional set of estimates (open circles) is obtained by looking at points of intersection $u_{*}$ of two sequential curves of $\Pi_{d}$: E.g., for $L_{1}$ and $L_{2}$ we may look for $\Pi_{d}(L_{1},u_{*})=\Pi_{d}(L_{2},u_{*})$ and study the resulting $u_{*}$ as a function of $(L_{1}L_{2})^{-1/2\nu}$. For $d=3$, when $L_{1}=L$, then $L_{2}=2L$. This sequence of estimates leads to the same $u_{c}$. This $u_{c}$ corresponds to $p_{c}=0.125\pm 0.001$, which is rather close to the threshold of $p_{c}=0.139\pm 0.001$ found for bond percolation in Ref. Abete et al. (2004). The fact that the site percolation threshold found in our work is smaller than the bond percolation threshold of Ref. Abete et al. (2004) is somewhat surprising.

By following the value of $\Pi_{3}(L,u_{*})$, at points $u_{*}$ where curves $\Pi_{3}(L,u)$ for successive $L$ intersect, we estimate the critical spanning probability $\Pi^{c}_{3}=0.04\pm 0.01$. This value is slightly larger than $\Pi^{c}_{3}\approx 0.032$ found by Abete et al. Abete et al. (2004) in bond percolation; both are significantly smaller than the analogous number 0.513 Acharyya and Stauffer (1998) in Bernoulli percolation.

Close to the percolation threshold, the correlation length $\xi$ diverges as $|u-u_{c}|^{-\nu}$. For finite system size $L$, as long as $L>\xi$ the percolation probability resembles that of an infinite system, i.e., $\Pi_{d}(L,u)\approx 1$ or $\approx 0$, for $u<u_{c}$ or $u>u_{c}$, respectively. For $\xi\gtrsim L$ the value of $\Pi_{d}(L,u)$ decreases from close to 1 to near 0 as $u$ increases. Therefore, the transition region approximately appears when $b|u-u_{c}|^{-\nu}>L$, where $b$ is a numerical prefactor. Thus, the width of the transition region scales as $\Delta u\approx(L/b)^{-1/\nu}$. Since $\Pi_{d}$ changes between 1 and 0 in the transition region, we expect the absolute value of the slope in that region to be $\approx 1/\Delta u$ or

[TABLE]

(Here and thereafter we disregard the negative sign of the slope.)

In $d=3$ we measured the absolute value of the maximal slope of the curves in Fig. 3 for each $L$, and the results are presented in Fig. 5. According to Eq. (6), the slope of this relation on a logarithmic plot should be $1/\nu$. We observe a slow and noisy increase of the slope with increasing $L$. The inverse of the extrapolated slope produces an estimate of $\nu=2.04\pm 0.08$ that is larger than $\nu=1.8\pm 0.1$ that was found on significantly smaller systems in Ref. Abete et al. (2004) and agrees well with the value $\nu=2$ predicted by Eq. (4). Instead of using the largest slopes for each $L$, we could have concentrated on significantly smaller slopes of the curves at $u=u_{c}$. Such an approach may provide a useful estimate of possible systematic errors. However, our data are too noisy and inaccurate in this area to produce reliable results.

V Percolation in $d=4$, 5, and 6

In higher dimensions $d$, we followed a similar strategy to the one described in the previous section. For $d=4$ and 5, we generated RWs of $N=uL^{d}$ steps starting from $L=4$ and then doubling $L$ until $L=64$ for $d=4$ or $L=32$ for $d=5$. (With increased $d$, we had to limit the maximal size $L$.) For each $L$, from repeated tests of spanning along $x_{d}$, we determined $\Pi_{d}(L,u)$. Results of these measurements are depicted in Fig. 6.

The “drift” of the curves to the right in $d=4$, and especially in $d=5$, is significantly smaller than in $d=3$, allowing a rather accurate determination of the percolation threshold. In $d=4$ we find $u_{c}=2.99\pm 0.01$, which corresponds to occupation fraction $p_{c}=0.0898\pm 0.0007$ of unvisited sites. This value is roughly half of the percolation threshold of 0.197 for regular (Bernoulli) site percolation on a hypercubic lattice in $d=4$ Mertens and Moore (2018). In $d=5$ we estimate $u_{c}=3.025\pm 0.008$, which corresponds to the critical fraction of vacant sites $p_{c}=0.0730\pm 0.0006$. This result is again about half of the percolation threshold of 0.141 for Bernoulli site percolation on a hypercubic lattice in $d=5$ Koza and Poa (2016). It is difficult to determine the critical spanning probability $\Pi_{d}^{c}$ with any accuracy. However, by following the values of the intersections of $\Pi_{d}(L,u)$ curves with sequential values of $L$, we estimate that $\Pi_{4}^{c}$ is between 0.1 and 0.2, while $\Pi_{5}^{c}$ is about 0.4.

We next estimate values for the critical exponent $\nu$ from the dependence of slope on $L$ near the transition point. In $d=4$, we were able to measure accurately both the maximal slope and the slope at the estimated $u_{c}$. Differences between these two methods provide us with an estimate of the possible systematic error. We find that $\nu=1.0\pm 0.1$, which agrees with $\nu=1$ expected from Eq. (4), and is very different from the Bernoulli value of $\nu_{\rm B}=0.685$ Koza and Poa (2016). In $d=5$ the largest slope practically coincides with the slope at $u_{c}$ (except for $L=4$). We detect a slight dependence of the effective exponent $\nu_{\rm eff}$ on $L$, and the extrapolated value is $\nu=0.65\pm 0.03$, where the error bars reflect the uncertainty in the extrapolation procedure. This agrees excellently with $\nu=2/3$ expected from Eq. (4), and differs from the regular (Bernoulli) percolation exponent of $\nu_{\rm B}=0.57$ Koza and Poa (2016).

In $d=6$, we generated random walks of length $uL^{6}$ for $L=4,6,8,11$, and 16. This required dealing with lattice sizes $M\approx 1.7\times 10^{7}$ and RWs reaching $6\times 10^{7}$ steps. Figure 7 depicts the spanning probability $\Pi_{6}(L,u)$ as a function of $u$ for several values of $L$. The maximal rate of decrease of all the curves is close to the point where $\Pi\sim 1/2$ or slightly higher. There is almost no drift in the curves with increasing $L$, and we estimate $u_{c}=3.10\pm 0.05$, which corresponds to critical fraction of vacant sites $p_{c}=0.062\pm 0.003$. The latter is significantly smaller than the threshold 0.109 for Bernoulli site percolation Koza and Poa (2016). The intersection points between $\Pi_{6}(L,u)$ for successive $L$ drift strongly upward, leading to a rather large estimate of $\Pi_{6}^{c}\sim 0.9$ for the critical spanning probability. Since $d=d_{c}=6$ is the common upper critical dimension Toulouse (1974), both the regular percolation theory and Eq. (4) posit $\nu=1/2$. Unfortunately, the maximal lattice size of $L=16$ computationally available to us, is too small to determine $\nu$ with any degree of accuracy. In addition to the obvious limitations of small $L$, we note that at $d_{c}$ there are logarithmic corrections in cluster size distributions Essam et al. (1978); Nakanishi and Stanley (1980) that further complicate detection of the trends. A straightforward power-law fit in the range $4\leq L\leq 16$ produces $\nu_{\rm eff}\approx 0.6$. Due to the above mentioned limitations, we believe that our numerical result does not contradict the expected value of $\nu=1/2$.

VI Spanning probability in $d=2$

The behavior of a two-dimensional RW of $N=uL^{2}$ steps on the $L\times L$ square lattice is quite different from the higher dimensions discussed previously. For $u=O(1)$, the root mean-squared end-to-end distance of a walk is of the order of the linear size $L$, rendering the imagery of multiple “foldings” of a much longer RW in Sec. II inapplicable. As in higher dimensions, we will test for spanning along the $x_{2}$ direction, while assuming periodicity in the $x_{1}$ direction. Clearly, for $u\ll 1$ a typical walk is simply too short to block percolation of vacant sites along $x_{2}$, while for $u\gg 1$ a single “circumnavigation” of the walk in $x_{1}$ direction will almost certainly block the percolation in $x_{2}$ direction. In earlier work, Banavar et al. Banavar et al. (1985) studied properties of clusters of vacant sites in a system of this type. They were primarily interested in the fractal and fracton dimensions in the regime where the RW covers a finite fraction of a lattice, i.e., for $u$ of order of unity.

Since the fractal dimension of a RW is 2, it is not surprising that its behavior in $d=2$ exhibits important differences from higher dimensions. For example, the number of distinct sites visited by an $N$-step walk in $d=2$ on an infinite lattice is modified from Eq. (1) to the leading order as Dvoretzky and Erdős (1951); Erdős and Taylor (1960); Vineyard (1963); Montroll and Weiss (1965)

[TABLE]

Correspondingly, the number of visits of a random walker to its initial position increases logarithmically with $N$, namely $B_{2}=(\ln N)/\pi$, as opposed to a constant $B_{d}$ for $d>2$. Similarly, the number of steps required to visit all sites on a lattice of $M$ sites is $\sim M\ln^{2}M$, i.e. with an extra logarithm compared to higher dimensions. Thus, when considering RWs of $N=uL^{2}$ steps on a square lattice of $M=L^{2}$ sites, with $u=O(1)$, the density of occupied sites actually vanishes in the limit $L\to\infty$. This again justifies arguments used in Sec. II to demonstrate that, even for moderate $L$, there is a pure exponential dependence of the fraction of vacant sites $p$ on the parameter $u$.

The results of our numerical study of the dependence of $p$ on $u$ are presented in Fig. 8. We performed simulations for $L=4,8,16,\dots,512$. All curves are well approximated by a pure exponential, and appear as straight lines on the semi-logarithmic plots in Fig. 8. Steps visible for $L=4$ (and, to lesser extent, on $L=8$) are a consequence of $N$ being an integer, requiring downward truncation to integer of $uL^{2}=16u$ (or $64u$). The apparent slope $A(L)$ in Fig. 8 keeps decreasing with increasing $L$. In fact, even the product $A(L)\ln L$ has some residual dependence on $L$, but approaches $\pi/2$ for large $L$. We verified this convergence to five-digit accuracy. [The extra factor of 2 is due to the fact that to leading order $\ln N=\ln(uL^{2})\approx 2\ln L$.] Thus for large $L$,

[TABLE]

In the asymptotic limit of large $L$, this result has been proven in Ref. Brummelhuis and Hilhorst (1992). The choice of $u$ over $p$ as the control parameter is inconsequential in higher dimensions where the two quantities are related by a fixed function. In $d=2$, the relation between $p$ and $u$ in Eq. 8 depends on $L$, and using $u$ creates a somewhat different perspective.

Figure 9 depicts numerically calculated accurate values of $\Pi_{2}(L,u)$ obtained using large statistics ($10^{5}$ samples per data point). We studied the spanning probability on a square lattice for $L=4,8,16,\dots,512$; the limiting factor for the largest $L$ was the need to evaluate $\Pi_{2}(L,u)$ with high accuracy. As before, the steps seen on the graph for $L=4$ are a result of truncating $uL^{2}=16u$ to integer $N$. We immediately note that there is no sign of $\Pi_{2}(L,u)$ becoming a step function with increasing $L$: The two-dimensional problem does not have a nontrivial percolation threshold. On the other hand, we clearly see that for large $L$, the spanning probability approaches a smooth function: $\lim_{L\to\infty}\Pi_{2}(L,u)=\Pi_{2}(\infty,u)$.

It can be argued that the existence of a finite $\Pi_{2}(\infty,u)$ is a consequence of the reduced role of lattice spacing in $d=2$. For $d\geq 3$ the discreteness of the lattice plays a crucial role: For percolation in the $x_{d}$ direction, a RW that is a one-dimensional path must completely block passage between boundaries at $x_{d}=0$ and $x_{d}=L-1$. As this can only be achieved by creating a $(d-1)$-dimensional continuous surface separating the space into disconnected parts, the RW path must either be on a lattice, or endowed with some thickness $a$ to block a finite volume. By contrast in $d=2$, percolation, say, in the vertical direction, can be blocked by a path crossing the system horizontally, which can be accomplished by a RW of zero thickness, as long as it acquires an extension $R\sim aN^{1/2}$, comparable to the system dimension $aL$. The corresponding ratio $(R/L)^{2}$ is proportional to $u$ and independent of $a$, would then determine the finite probability $\Pi_{2}(\infty,u)$

Since the relation between $p$ and $u$ depends on $L$, as in Eq. (8), reexpressing $\Pi_{2}(\infty,u)$ as a function of $p$ leads to a function that depends on $L$. As depicted in Fig. 10, with increasing $L$ the corresponding $\Pi_{2}(L,p)$ shift toward $p=1$. Thus $p=1$ serves as a trivial percolation threshold in $d=2$.

VII Very short and long walks in $d=2$

VII.1 Limit of small $u$

For $u\ll 1$, the typical linear size of a RW, such as its end-to-end distance, is $\sqrt{u}L$. Since this is shorter than the linear size, $L$, of the lattice, a typical RW cannot block percolation in the vertical ($x_{2}$) direction, and $\Pi_{2}\approx 1$. The only reason for deviation of $\Pi_{2}$ from 1 is due to rare configurations that, stretching far beyond the typical end-to-end distance, circumnavigate the periodic box in the horizontal ($x_{1}$) direction, with end segments intersecting the starting segments to completely block percolation in the vertical direction. For example, for $L=64$ at $u=0.1$ the probability of absence of percolation is about one part per $10^{5}$, and Fig. 11 depicts one such rare event. Since such an event requires a RW of root-mean-square length $\sqrt{u/2}L$ in the horizontal direction to be stretched by at least $L$ in that direction, a lower bound for the probability of getting a nonpercolating configuration can be obtained by integrating the Gaussian probability distribution of the end-to-end distance from $L$ to infinity. The resulting integral can be expressed in terms of the error function, as

[TABLE]

This is just an estimate, since besides circumnavigating the periodic cell horizontally, the RW must also self-intersect. Since for $u\ll 1$ nonpercolating configurations are rare, we sampled $10^{6}$ configurations for each data point for small $u$ with $L=512$. Below $u=0.1$ we could not find nonpercolating configurations. However, in the range $0.1\leq u\leq 0.3$ we got rather accurate values of $\Pi_{2}(512,u)$ that fit quite well to ${\rm Erf}(1/\sqrt{u})$, as expected in Eq. (9).

VII.2 Limit of large $u$

It is interesting to view the results of Fig. 9 on a semilogarithmic scale, as depicted in Fig. 12. Since each data point was obtained by averaging $10^{5}$ random configurations, the statistical errors exceed 10% for $\Pi_{2}$ below $10^{-3}$, and exceed 30% for $\Pi_{2}$ below $10^{-4}$ (comprising all fluctuations in the bottom right corner of the figure). We note that for $u>3$, the spanning probability for $\Pi_{2}\ll 1$ seems to decay exponentially (while maintaining some residual dependence on $L$). To understand this behavior of the spanning probability, we take a closer look at the shapes of percolating configurations. Figure 13 depicts four such examples, for different values of $u$, with the black squares indicating sites visited by the RW, while the white area corresponds to vacant sites. By construction, the black sites always form a single cluster, while the vacant sites can be split into many clusters. We notice that most of the vacant sites also form a single large cluster. For $u=1$, vacant sites percolate most of the time, and Fig. 13(a) represents a “typical” configuration. For $u=4$, 6, and 8, the percolating configurations are exceptional (dropping to below $10^{-6}$ for $u=8$). As $u$ increases, the rare percolating clusters start to resemble narrow white bands connecting the top and bottom boundaries, and the probability of such a rare configuration is estimated next.

Consider a two-dimensional RW of length $N=uL^{2}$ on a square $L\times L$ lattice with periodic boundary conditions. The probability distribution of the end point can be estimated from the diffusion equation $\partial P/\partial N=D\nabla^{2}P$, with diffusion constant $D=a^{2}/4$, where $a=1$ is the lattice spacing. Consider an extreme configuration where the percolating cluster of unvisited sites is a single straight line in the vertical direction. The allowed configurations of RW in that case are those that can fit between two vertical lines separated by distance $L$. The total number of $N$-step RWs is $4^{N}$. The fraction of configurations that fit within a strip of width $L$ can be found by solving the diffusion equation with absorbing boundary conditions. The diffusion equation can be separated into two independent parts: one for the vertical direction which imposes no limits on the walker, and one for the horizontal component with absorbing boundary conditions. With the horizontal coordinate denoted by $x$, the normalized eigenstates are $\Psi_{n}(x,t)=\sqrt{2/L}\sin(\pi nx/L)\exp(-Dn^{2}\pi^{2}N/L^{2})$. For large $u$, and hence $N=uL^{2}$, the solution is dominated by the term with $n=1$, leading to the probability

[TABLE]

to find a random walker that started at $x_{0}$ and arrived at $x$. To obtain the total survival probability, we need to integrate over $x$ and to average over the starting point $x_{0}$. (The vertical direction is completely free and does not affect the probability.) The final result is the ‘survival probability’ of $G\propto\exp(-D\pi^{2}N/L^{2})\propto\exp(-\pi^{2}u/4)$. (We omitted the numerical prefactor of this expression.) Note that our answer does not depend on $L$, and $G\sim\exp(-\pi^{2}u/4)$ provides a lower bound for $\Pi_{2}$ at large $u$. This result clearly underestimates $\Pi_{2}(L,u)$, as the percolating channel does not have to be a straight line; it can be undulating or inclined. We anticipate that these factors will modify the prefactor for $G$, but leave unchanged the leading exponential part that is depicted by the dashed line in Fig. 12. It clearly underestimates $\Pi_{2}$ by several orders of magnitude, but its slope (on the semilogarithmic plot) is close to the behavior of the numerical results for large $u$.

VIII Discussion

In this work we studied percolation of sites unvisited by RWs on a periodic lattice. We extended previous results by looking at larger lattice sizes in $d=3$ and by considering all dimensions $2\leq d\leq 6$. Our primary goal was to find the dimensionality dependence of such characteristics as the percolation threshold $u_{c}$, the critical exponent $\nu$, and the critical spanning probability $\Pi_{d}^{c}$. While our results agree well with the general theory of correlated percolation, much remains unknown. We concentrated on a single critical exponent $\nu$, and did not attempt to calculate additional critical exponents, e.g., describing cluster sizes and fractal dimensions, which may provide an interesting direction for further study. The RW creates rather unique constraints on cluster structure that may lead to interesting results not only for $d\leq 6$, but also in higher dimensions when the critical behavior is mean-field like.

The original impetus for the study of this type of percolation came from degradation of a gel. Percolation in such a gel is accompanied by dramatic changes in physical properties such as elasticity. However, theoretical studies (with the exception of Ref. Banavar et al. (1985)) have so far concentrated on purely geometrical quantities. Thus it would be interesting to extend the theoretical investigation of such correlated percolation to conductivity and diffusion.

Acknowledgements.

Y.K. thanks D. Gomez for stimulating discussions on related subjects. M.K. was supported by the National Science Foundation through Grant No. DMR-1708280. Y.K. was supported by the Israel Science Foundation Grant No. 453/17.