A critical exponent for shortest-path scaling in continuum percolation

Generating realizations of ${{X}^{\gamma }}$ amounts to finding the shortest path (if it exists) between o and ${{v}_{\gamma }}$ in ${{G}^{\gamma }}$. Our approach to simulating ${{X}^{\gamma }}$ is based on the method introduced in [5]. This is an extension of the classical Leath model, introduced in [12], where clusters are grown in shells. A similar approach to growing continuum percolation clusters was used in a three dimensional setting in [27].

Because we are interested in the asymptotic behavior of ${{X}^{\gamma }}$, we increase the speed of our algorithm by only considering values of ${{X}^{\gamma }}$ greater than 2. That is, we assume there are no points in ${{B}_{\gamma }}(o)\cap {{B}_{\gamma }}({{v}_{\gamma }})$, where ${{B}_{\gamma }}(v)=\{x\in {{\mathbb{R}}^{d}}:\left\|x-v\right\|\lt \gamma \}$.

We 'grow' clusters from the 'seed' points o and ${{v}_{\gamma }}$ in a sequential fashion by generating Poisson processes in unexplored regions within distance γ of existing points. We are able to proceed in this sequential manner due to the independence property of the Poisson process. We begin by generating points in ${{B}_{\gamma }}(o)\backslash {{B}_{\gamma }}({{v}_{\gamma }})$, which form the initial points of the first cluster, and points in ${{B}_{\gamma }}({{v}_{\gamma }})\backslash {{B}_{\gamma }}(o)$, which form the initial points of the second cluster. We mark the points o and ${{v}_{\gamma }}$ as explored and mark the newly generated points as unexplored, as they have not yet been explicitly considered by our algorithm. We then proceed to sequentially examine the unexplored points, which are ordered by their chemical distance from the seed point of their cluster. If two points have the same chemical distance, we choose one at random. We examine an unexplored point, v, by generating a Poisson process in ${{B}_{\gamma }}(v)$, the open ball centered at v with radius γ. We delete all new points in this ball whose Euclidean distance to an explored point is less than γ. We add the surviving new points to the appropriate clusters. If the two clusters connect, we stop the algorithm and return ${{X}^{\gamma }}$. If it is no longer possible to add a point (i.e., we have explored all regions within distance γ of an existing point and there are no surviving new points), we stop the algorithm and return the empty set, ∅. By choosing the new points to explore based on their chemical distance to the seed of their cluster, we are essentially using a variation of Dijkstraʼs algorithm, i.e., breadth-first search; see [28]. This means that, if the two clusters connect, we can immediately determine the shortest path between o and ${{v}_{\gamma }}$.

We use two sequences of sets, ${{({{U}_{m}})}_{m\geqslant 0}}$ and ${{({{E}_{m}})}_{m\geqslant 0}}$, to describe the unexplored and explored points generated by our algorithm. The set U_m describes the unexplored points remaining at the mth step and the set E_m describes the points that have already been explored at the mth step. These sets consist of triples of the form $(v,c,t)\in {{\mathbb{R}}^{d}}\times \{1,2\}\times \mathbb{N}$. The Euclidean coordinates of a point are given by v. If a point is in the cluster containing the origin then c = 1 and if a point is in the cluster containing ${{v}_{\gamma }}$ then c = 2. The chemical distance from a point to the originating point of its cluster is given by t.

Because the running time of the algorithm is potentially unbounded, we introduce a stopping condition as in [5]. If the sizes of the clusters grow beyond a certain threshold without meeting, we terminate the algorithm and return ∅. We measure the sizes of the clusters using the sequence ${{({{L}_{m}})}_{m\geqslant 0}}$, where L_m is the minimal chemical distance from a point in U_m to the originating point of the cluster containing it. We denote the threshold by ${{L}_{{\sf max} }}.$ Note that this stopping condition means we are simulating ${{X}^{\gamma }}$ conditional on ${{X}^{\gamma }}\leqslant 2{{L}_{{\sf max} }}+1.$

The algorithm is described precisely in the following.

Algorithm 1 (Simulating Xγ).
1. Set ${{U}_{0}}=\varnothing$ and ${{E}_{0}}=\{(o,1,0),({{v}_{\gamma }},2,0)\}$. Set m = 0.
2. Generate $Y\sim {\sf Pois} ({{\kappa }_{d}}{{\gamma }^{d}})$ points, where ${{\kappa }_{d}}$ denotes the volume of ${{B}_{1}}(o)$ in ${{\mathbb{R}}^{d}}$, in ${{B}_{\gamma }}(o)$ and remove all points that also lie in ${{B}_{\gamma }}({{v}_{\gamma }})$. Add the remaining points to U0 with c = 1 and t = 1.
3. Generate $Y\sim {\sf Pois} ({{\kappa }_{d}}{{\gamma }^{d}})$ points in ${{B}_{\gamma }}({{v}_{\gamma }})$ and remove all points that also lie in ${{B}_{\gamma }}(o)$. Add the remaining points to U0 with c = 2 and t = 1.
4.Choose a point $p=(v,c,t)$ from the current set of unexplored points, Um, with minimum chemical distance t. Set ${{U}_{m+1}}={{U}_{m}}\backslash \{p\}$, ${{E}_{m+1}}={{E}_{m}}\cup \{p\}$, and ${{L}_{m}}=t$.
5. Let $p_{1}^{o}=(v_{1}^{o},c_{1}^{o},t_{1}^{o}),\ldots ,p_{J}^{o}=(v_{J}^{o},c_{J}^{o},t_{J}^{o})$ be the points in Em that are not in the same cluster as p. Let ${{\mathcal{J}}_{S}}=\{j:v\in {{B}_{\gamma }}(v_{j}^{o})\}$ be the set of indices associated with those points from $\{v_{1}^{o},\ldots ,v_{J}^{o}\}$ that can be connected to v by an edge with length less than γ. If ${{\mathcal{J}}_{S}}\ne \varnothing$ terminate the algorithm and return ${{X}^{\gamma }}=t+{{t}_{{\sf min} }}+1$, where ${{t}_{{\sf min} }}={\sf min} \{t_{j}^{o}:j\in {{\mathcal{J}}_{S}}\}$.
6. Generate $Y\sim {\sf Pois} ({{\kappa }_{d}}{{\gamma }^{d}})$ points, ${{\tilde{v}}_{1}},\ldots ,{{\tilde{v}}_{Y}}$, distributed uniformly in ${{B}_{\gamma }}(v)$. Set y = 1.
7. Let $v_{1}^{e},\ldots ,v_{K}^{e}$ be the first components of the points in Em. If ${{\tilde{v}}_{y}}\notin \bigcup _{k=1}^{K}{{B}_{\gamma }}(v_{k}^{e})$, set ${{U}_{m+1}}={{U}_{m+1}}\cup \{({{\tilde{v}}_{y}},c,t+1)\}$.
8. If $y\lt Y$, set $y=y+1$ and repeat from step 7.
9. If ${{L}_{m}}\gt {{L}_{{\sf max} }},$ terminate the algorithm and return ∅.
10. If ${{U}_{m+1}}\ne \varnothing$, set $m=m+1$ and repeat from step 4. Otherwise, return ∅.

In practice, the data structures used are chosen carefully in order to improve the performance of the above algorithm. In particular, the performance of steps 3 and 5 is improved considerably by using a data structure for the family of sets ${{({{E}_{m}})}_{m\geqslant 0}}$ that orders the points according to their Euclidean coordinates. For instance, one can discretize the Euclidean space into boxes and use an array structure to access all points in a given box; see, for example, the approaches taken in [9, 10, 27, 29].

For the planar case, d = 2, the basic idea of the algorithm is illustrated in the following. We begin with the two original vertices (which have different symbols in order to distinguish the two clusters); see figure 1, left. The area around them is still unexplored. We first examine the left point, o. The pattern marks the area where we know there are no points. We generate a homogeneous Poisson process in the region of radius γ surrounding o. In this case, one point is generated; see figure 1, center. We mark the new points as belonging to cluster 1 and add the explored point to the set of explored points (marked in black). The same is done for the other original vertex, ${{v}_{\gamma }}$. Again, only one point is generated; see figure 1, right. We next generate points in the region surrounding the point in the first cluster of chemical distance 1 to the originating point. The region that has already been explored (i.e., in which the Poisson process has already been generated) is marked in gray; see figure 2, left. We mark this point as explored and add the new points to the list of points to be explored. We do the same for the point with chemical distance 1 to the originating point of the second cluster; see figure 2, center. Next, we consider one of the points in the first cluster that is of chemical distance 2 to its originating point. This point is within distance γ of an explored point in the second cluster, forming a path between o and ${{v}_{\gamma }}$ of chemical length 4; see figure 2, right. The algorithm terminates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As stated above, for $\gamma ={{\gamma }_{{\rm c}}}$, we expect that the chemical distance displays power law behavior. For d = 2, we provide numerical evidence of this in section 5. Thus, estimating the critical exponent λ amounts to estimating the parameters of a (right-truncated) power law distribution. There are two commonly used techniques: linear regression and maximum likelihood. It has been demonstrated in a number of studies that maximum likelihood is a more robust approach (see [30, 31]). In particular, almost all the standard assumptions about statistical error made for least squares estimation are violated, resulting in inaccurate estimators of λ and incorrect estimators of standard errors. We therefore follow the approach described in [30, 32]. That is, we assume that there exists an ${{x}_{{\sf min} }}$ such that, for all $x\gt {{x}_{{\sf min} }}$

$$\begin{eqnarray}&&\mathbb{P}({{X}^{\gamma }}=x)=c{{x}^{-\lambda }},\end{eqnarray} \tag{ 1 }$$

for some $c\gt 0$. Given a sample $X_{1}^{\gamma },\ldots ,X_{N}^{\gamma }$, generated according to the algorithm described in section 3.1, we determine λ and ${{x}_{{\sf min} }}$ as follows.

The random variable ${{X}^{\gamma }}$ conditioned on being less than or equal to ${{x}_{{\sf max} }}=2{{L}_{{\sf max} }}+1$ has probability mass function $\ell :\{{{x}_{{\sf min} }},\ldots ,{{x}_{{\sf max} }}\}\to [0,1]$ with

$$\begin{eqnarray}&&\ell (x)=\mathbb{P}({{X}^{\gamma }}=x\;|\;{{x}_{{\sf min} }}\leqslant {{X}^{\gamma }}\leqslant {{x}_{{\sf max} }})=\frac{{{x}^{-\lambda ({{x}_{{\sf min} }})}}}{\sum _{y={{x}_{{\sf min} }}}^{{{x}_{{\sf max} }}}{{y}^{-\lambda ({{x}_{{\sf min} }})}}}.\end{eqnarray} \tag{ 2 }$$

This yields the maximum likelihood estimator $\hat{\lambda }({{x}_{{\sf min} }})$, which is the solution of

$$\begin{eqnarray*}&&\frac{\sum _{y={{x}_{{\sf min} }}}^{{{x}_{{\sf max} }}}{{y}^{-\lambda }}{\rm log} y}{\sum _{y={{x}_{{\sf min} }}}^{{{x}_{{\sf max} }}}{{y}^{-\lambda }}}=\frac{1}{N}\mathop{\sum }\limits_{i=1}^{N}{\rm log} X_{i}^{\gamma }.\end{eqnarray*}$$

Thus, $\hat{\lambda }$ can be easily calculated via standard numerical procedures.

The choice of ${{x}_{{\sf min} }}$ involves a trade-off. If a too small value of ${{x}_{{\sf min} }}$ is chosen, it may not be true that equation (1) holds, even approximately. However, if ${{x}_{{\sf min} }}$ is too large, only a small fraction of the total sample will be larger than ${{x}_{{\sf min} }}$ and the estimator of λ will be inaccurate. In practice, we estimate $\lambda ({{x}_{{\sf min} }})$ for ${{x}_{{\sf min} }}={{x}_{0}},{{x}_{0}}+1,{{x}_{0}}+2,\ldots ,{{x}_{1}}$, where $0\leqslant {{x}_{0}}\lt {{x}_{1}}\leqslant {{x}_{{\sf max} }}-1$. We then evaluate the goodness-of-fit of each $\hat{\lambda }({{x}_{{\sf min} }})$ by computing the Kolmogorov–Smirnov statistic

$$\begin{eqnarray}&&D({{x}_{{\sf min} }})=\mathop{{\sf max} }\limits_{{{x}_{{\sf min} }}\leqslant x\leqslant {{x}_{{\sf max} }}}|F(x)-\hat{F}(x)|,\end{eqnarray} \tag{ 3 }$$

where, given the estimator $\hat{\lambda }({{x}_{{\sf min} }})$, F(x) is the theoretical distribution function based on equation (2) and $\hat{F}(x)$ is the empirical distribution function. We choose the values of ${{x}_{{\sf min} }}$ and $\hat{\lambda }({{x}_{{\sf min} }})$ corresponding to the smallest value of D.

We performed our simulation not only for $\gamma ={{\gamma }_{{\rm c}}}$, but also for other choices of γ. The resulting simulations provide strong evidence that polynomial decay of the conditional probability mass function of ${{X}^{\gamma }}$, given in (2), occurs only for $\gamma ={{\gamma }_{{\rm c}}}$. The results are shown in figure 4.

Zoom In Zoom Out Reset image size

We carried out a numerical investigation of the critical exponent for d = 2. The critical percolation threshold for this model, ${{\gamma }_{{\rm c}}}$, is estimated in [10] to be ${{\gamma }_{{\rm c}}}\approx 1.198\;468$; see, [36] and [37] for two additional papers where this threshold is estimated.

Using the standard Monte Carlo procedure, we generated a sample of size $N=1.5\times {{10}^{9}}$ with ${{L}_{{\sf max} }}=4000$. Generating this sample took approximately 26.7 d computer time. We estimated $\lambda ({{x}_{{\sf min} }})$ for ${{x}_{{\sf min} }}=3,\ldots ,7990$. According to the Kolmogorov–Smirnov statistic, given in (3), the optimal value of ${{x}_{{\sf min} }}$ is ${{x}_{{\sf min} }}=458$. This gives a rate estimate of $\hat{\lambda }(458)=2.1025$.

The estimates of $\lambda ({{x}_{{\sf min} }})$ for ${{x}_{{\sf min} }}=380,\ldots ,6000$ are shown on the left of figure 5. The estimates of $\lambda ({{x}_{{\sf min} }})$ for ${{x}_{{\sf min} }}=380,\ldots ,2000$ are shown on the right of figure 5. The values of $\hat{\lambda }({{x}_{{\sf min} }})$ appear fairly stable for at least the first 2000 values of ${{x}_{{\sf min} }},$ with the subsequent loss of stability due to the smaller sample sizes used for large values of ${{x}_{{\sf min} }}.$ Closer investigation of values in this region shows that $\hat{\lambda }({{x}_{{\sf min} }})$ fluctuates between 2.102 and 2.106 with the estimate obtained using the Kolmogorov–Smirnov statistic at the lower end of the range. Note that this degree of fluctation shows the sensitivity of the estimator to the choice of ${{x}_{{\sf min} }}.$

Zoom In Zoom Out Reset image size

Note that, prior to the simulations with ${{L}_{{\sf max} }}=4000$, we carried out simulations using ${{L}_{{\sf max} }}=3000$ with a sample size of $N={{10}^{9}}$. These simulations gave an optimal value of ${{x}_{{\sf min} }}=854$ with $\hat{\lambda }(854)=2.1064$. This suggests that there is considerable variability in the estimator and also sensitivity to the cut-off point ${{L}_{{\sf max} }}.$

We also used splitting to generate samples of ${{X}^{{{\gamma }_{{\rm c}}}}}$, carrying out the splitting at the level 400. This produced a sample of ${{X}^{{{\gamma }_{{\rm c}}}}}$ conditioned on being larger than 800. The splitting point was chosen based on the standard Monte Carlo sample with ${{L}_{{\sf max} }}=3000$ and $N={{10}^{9}}$ which returned an optimal choice of ${{x}_{{\sf min} }}=854$. We calculated a sample of $N=1.5\times {{10}^{7}}$ with a splitting factor of s = 100, giving us a sample in the tail of equivalent size to the sample gained using standard Monte Carlo with $N=1.5\times {{10}^{9}}$. This took approximately 20.2 d to generate, which is roughly 75% of the time required by the standard Monte Carlo approach. The optimal value of ${{x}_{{\sf min} }},$ according to the Kolmogorov–Smirnov statistic, given in (3), is ${{x}_{{\sf min} }}=1065$. This gives a rate estimate of $\hat{\lambda }(1065)=2.1051$.

The estimates of $\lambda ({{x}_{{\sf min} }})$ for ${{x}_{{\sf min} }}=900,\ldots ,6000$ are shown on the left of figure 6 and the estimates of $\lambda ({{x}_{{\sf min} }})$ for ${{x}_{{\sf min} }}=900,\ldots ,2000$ are shown on the right of figure 6. As in the standard Monte Carlo case, the values of $\hat{\lambda }({{x}_{{\sf min} }})$ appear fairly stable for at least the first 2000 values of ${{x}_{{\sf min} }},$ with the subsequent loss of stability due to the smaller sample sizes used for large values of ${{x}_{{\sf min} }}.$ Closer investigation of values in this region shows that, at least for ${{x}_{{\sf min} }}\gt 1000$, $\hat{\lambda }({{x}_{{\sf min} }})$ fluctuates between about 2.1035 and 2.1055 with the estimate obtained using the Kolmogorov–Smirnov statistic at the higher end of the range.

Zoom In Zoom Out Reset image size

We then carried out an estimate of λ using the combined sample from the standard Monte Carlo and splitting approaches. This sample of ${{X}^{{{\gamma }_{{\rm c}}}}}$ conditioned on being larger than 800 was equivalent to that obtained using standard Monte Carlo with a sample size $N=3\times {{10}^{9}}$. The optimal value of ${{x}_{{\sf min} }},$ according to the Kolmogorov–Smirnov statistic, is ${{x}_{{\sf min} }}=991$. This gives a rate estimate of $\hat{\lambda }(991)=2.1045$.

The estimates of $\lambda ({{x}_{{\sf min} }})$ for ${{x}_{{\sf min} }}=900,\ldots ,6000$ are shown on the left of figure 7. The estimates of $\lambda ({{x}_{{\sf min} }})$ for ${{x}_{{\sf min} }}=900,\ldots ,2000$ are shown on the right of figure 7. As in both the standard Monte Carlo and splitting cases, the values of $\hat{\lambda }({{x}_{{\sf min} }})$ appear fairly stable for at least the first 2000 values of ${{x}_{{\sf min} }}.$ Closer investigation of values in this region shows that, as in the splitting case, for ${{x}_{{\sf min} }}\gt 1000$, $\hat{\lambda }({{x}_{{\sf min} }})$ fluctuates between about 2.1035 and 2.1055. Here, the estimate obtained using the Kolmogorov–Smirnov statistic is in the middle of the range.

Zoom In Zoom Out Reset image size

The four estimates we have obtained are $\hat{\lambda }=2.1025$ for the standard Monte Carlo sample with ${{L}_{{\sf max} }}=4000$ and $N=1.5\times {{10}^{9}}$, $\hat{\lambda }=2.1064$ for the standard Monte Carlo sample with ${{L}_{{\sf max} }}=3000$ and $N={{10}^{9}}$, $\hat{\lambda }=2.1051$ for the splitting sample with ${{L}_{{\sf max} }}=4000$ and $N=1.5\times {{10}^{9}}$, and $\hat{\lambda }=2.1045$ for the joint sample with sample size equivalent to a standard Monte Carlo sample of $N=3\times {{10}^{9}}$. We note that standard errors are not available for these estimators, which are also biased (though not as biased as estimates obtained using least squares). Nevertheless, these estimates broadly agree with the results for lattice percolation given in [5] and [7], as well as the estimate based on the results in [26].