Anisotropy in finite continuum percolation: Threshold estimation by Minkowski functionals

TL;DR

This paper investigates how anisotropy affects percolation thresholds in continuum models, using Minkowski functionals to estimate thresholds and analyzing finite versus infinite systems.

Contribution

It introduces anisotropic extensions to the Boolean model and develops geometric approximation formulas for percolation thresholds based on Minkowski functionals.

Findings

Percolation thresholds depend on anisotropy and can be accurately predicted with empirical parameters.

Explicit geometric approximations provide reliable lower bounds for thresholds.

Differences in thresholds for different directions vanish in infinite systems.

Abstract

We examine the interplay between anisotropy and percolation, i.e., the spontaneous formation of a system spanning cluster in an anisotropic model. We simulate an extension of a benchmark model of continuum percolation, the Boolean model, which is formed by overlapping grains. Here we introduce an orientation bias of the grains that controls the degree of anisotropy of the generated patterns. We analyze in the Euclidean plane the percolation thresholds above which percolating clusters in $x$- and in $y$-direction emerge. Only in finite systems, distinct differences between effective percolation thresholds for different directions appear. If extrapolated to infinite system sizes, these differences vanish independent of the details of the model. In the infinite system, the uniqueness of the percolating cluster guarantees a unique percolation threshold. While percolation is isotropic even for anisotropic processes, the value of the percolation threshold depends on the model parameters, which we explore by simulating a score of models with varying degree of anisotropy. To which precision can we predict the percolation threshold without simulations? We discuss analytic formulas for approximations (based on the excluded area or the Euler characteristic) and compare them to our simulation results. Empirical parameters from similar systems allow for accurate predictions of the percolation thresholds (with deviations of $<5\%$ in our examples), but even without any empirical parameters, the explicit approximations from integral geometry provide, at least for the systems studied here, lower bounds that capture well the qualitative dependence of the percolation threshold on the system parameters (with deviations of $5\%$--$30\%$). As an outlook, we suggest further candidates for explicit and geometric approximations based on second moments of the so-called Minkowski functionals.

Equations

$$\mathcal{P}(\theta ):=\frac{\Gamma (\frac{\alpha }{2}+1)}{2\sqrt{\pi }\Gamma (\frac{\alpha +1}{2})}\cdot {\cos }^{\alpha }(\theta ),\mathrm{for}\theta \in [-\frac{\pi }{2},\frac{\pi }{2}),$$

$$\varphi =1-\exp (-\gamma A)$$

$$C({\varphi }_c^{\mathrm{eff}})=\frac{1}{2}.$$

$$C(\varphi )$$

$${\varphi }_c^{\mathrm{eff}}-{\varphi }_c\propto L^{-1/\nu }$$

$${\varphi }_c^{\mathrm{eff},x}$$

$${\varphi }_c^{\mathrm{eff},y}$$

$${\varphi }_c^x\equiv {\varphi }_c^y.$$

$${\varphi }_c^{\mathrm{eff},x}$$

$${\varphi }_c^{\mathrm{eff},y}$$

$$A_{\mathrm{ex}}(K,Q)=V_{1,1}^0(K,Q)+A(K)+A(Q)$$

$$⟨A_{\mathrm{ex}}⟩={\iint }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}\mathrm{d}{\theta }_i\mathrm{d}{\theta }_j\mathcal{P}({\theta }_i)\cdot \mathcal{P}({\theta }_j)\cdot A_{\mathrm{ex}}(\vartheta ({\theta }_i)K,\vartheta ({\theta }_j)Q),$$

$$⟨A_{\mathrm{ex}}⟩=\frac{1}{2\pi }P(K)\cdot P(Q)+A(K)+A(Q).$$

$$A_{\mathrm{ex}}(\vartheta ({\theta }_i)R,\vartheta ({\theta }_j)R)=(a^2+b^2)\sin |{\theta }_i-{\theta }_j|+2ab\left(\cos |{\theta }_i-{\theta }_j|+1\right).$$

$$B_c={\gamma }_c\cdot ⟨A_{\mathrm{ex}}⟩,$$

$${\gamma }_c=\frac{B_c}{⟨A_{\mathrm{ex}}⟩}.$$

$$B_c=\frac{⟨A_{\mathrm{ex}}⟩}{A}\ln \frac{1}{1-{\varphi }_c},$$

$$\overline{\chi }(\gamma )=\gamma \left(1-\frac{\gamma }{2}⟨V_{1,1}^0⟩\right)\exp (-\gamma A),$$

$$⟨V_{1,1}^0⟩:={\iint }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}\mathrm{d}{\theta }_i\mathrm{d}{\theta }_j\mathcal{P}({\theta }_i)\cdot \mathcal{P}({\theta }_j)\cdot V_{1,1}^0(\vartheta ({\theta }_i)R,\vartheta ({\theta }_j)R),$$

$$\chi ({\gamma }_0)$$

$$\iff$$

$${\varphi }_0:=1-{\mathrm{e}}^{-{\gamma }_{0}A},$$

$$\mathcal{L}=\prod _{L_i}\mathrm{Prob}\left[{\varphi }_{c,i}^{\mathrm{eff},x}|{\varphi }_c,m^x\right]\cdot \mathrm{Prob}\left[{\varphi }_{c,i}^{\mathrm{eff},y}|{\varphi }_c,m^y\right]$$

$$\mathrm{Prob}\left[{\varphi }_{c,i}^{\mathrm{eff},z}|{\varphi }_c,m^z\right]=\frac{1}{\sqrt{2\pi }{\sigma }_{z,i}}\exp \left[-\frac{{\left({\varphi }_{c,i}^{\mathrm{eff},z}-{\varphi }_c-m^z\cdot L_i^{-1/\nu }\right)}^2}{2{\sigma }_{z,i}^2}\right],$$

$$\log \mathcal{L}=const.-\frac{{\chi }_x^2}{2}-\frac{{\chi }_y^2}{2}$$

$$\displaystyle\chi_{x}^{2}+\chi_{y}^{2}=\left(\left(\begin{array}[]{c}\phi_{c}^{\mathrm{eff},x}\\ \phi_{c}^{\mathrm{eff},y}\end{array}\right)-F(L^{-\frac{1}{\nu}})\right)^{t}\mathbf{Cov}^{-1}\left(\left(\begin{array}[]{c}\phi_{c}^{\mathrm{eff},x}\\ \phi_{c}^{\mathrm{eff},y}\end{array}\right)-F(L^{-\frac{1}{\nu}})\right)$$

$$\displaystyle F(L^{-1/\nu})=\left(\begin{array}[]{c c c}{1}&L^{-1/\nu}&0\\ {1}&0&L^{-1/\nu}\end{array}\right)\cdot\left(\begin{array}[]{c}\phi_{c}\\ m^{x}\\ m^{y}\end{array}\right)\;,$$

$$\mathbf{Cov}=\mathrm{diag}({\sigma }_{x,1}^2,\dots ,{\sigma }_{x,i}^2,\dots ,{\sigma }_{y,1}^2,\dots {\sigma }_{y,i}^2,\dots ).$$

Full text

††: \JSTAT

Keywords: anisotropy, continuum percolation, percolation threshold, uniqueness of percolating cluster, excluded area, Euler characteristic

Percolation is a geometric phase transition in which formerly local clusters spontaneous form a system spanning cluster. Dating back to the work by Flory [1] and Stockmayer [2] and later by Broadbent and Hammersley [3], percolation theory describes in general the critical behavior of connected clusters, either on lattices or in the continuum [4, 5]. It is a general concept that can be applied to a multitude of very different physical systems and complex networks [6].

In this paper, we are interested in the interplay of topology and geometry, more precisely, of percolation and anisotropy. We investigate continuum percolation, that is, percolation on a continuous space, in statistically homogeneous (also called stationary) two-phase random media. More specifically, we study a benchmark model of continuum percolation [7], the so-called Boolean model. It consists of overlapping grains that are distributed randomly in space. The intensity $\gamma$, which is also called number density, gives the expectation for the number of grain centers per unit area. It is one of the parameters of the system.

The overlapping grains form clusters of connected components. Starting from a dilute arrangement of grains, the mean size of these clusters grows with increasing area fraction of the grain phase, called the occupied area fraction. Below the critical area fraction, the so-called percolation threshold, these clusters are almost surely (i. e., with probability one) bounded. At occupied volume fractions above the percolation threshold, almost surely exactly one unbounded cluster emerges that spans the whole system. It is called the percolating cluster. Percolation is a geometric phase transition with universal scaling of geometrical properties, for example, the size of the percolating cluster.

Anisotropy appears in many man-made or natural materials, from geological formations [8] over fiber composites and laminates [9] to porous media [10]. Moreover, anisotropy manifests itself in different forms. A heterogeneous medium might be isotropic w.r.t. the distribution of the area but anisotropic w.r.t. of the orientation of the interface, or vice versa [11]. Anisotropic two-phase random media can be modeled by the Boolean model with an orientation bias of the grains [11, 12].

In these systems, percolation exhibits an interesting qualitatively different behavior in finite samples or in the idealized infinite limit. Although there is no phase transition in a finite system, an effective percolation threshold can be defined for the ensemble. As expected, it is anisotropic. In other words, the system is more likely to percolate in the preferred direction than in the perpendicular direction.

However, this difference vanishes in the thermodynamic limit. Even the most anisotropic system simultaneously percolates in all directions. The percolation threshold is isotropic, see also [13, 14], which we here link to the uniqueness of the percolating cluster. There have already been numerous studies of anisotropic continuum percolation, some using analytic calculations [e. g., 15, 16, 17] others are based on Monte Carlo simulations [e. g., 13, 18] or experiments [e. g., 19]. There are also detailed studies of anisotropic lattice percolation [20, 21, 22].

While the isotropy of the percolation threshold is independent of the details of the model, the value of the percolation threshold is a nonuniversal constant, that is, it depends on the details of the system and specifically on its anisotropy. For Boolean models an increasing degree of anisotropy typically also increases the percolation threshold of the grain phase [14, 23, 18]. Because of the complexity of the geometric problem, analytic results for the percolation thresholds of Boolean models are hard to derive. Precise numerical estimates require complex and time-consuming calculations. Explicit formulas for estimates of percolation thresholds help to adjust models to experimental systems and predict their physical properties. We determine, from simulations, the percolation thresholds for a score of models with varying degree of anisotropy. How well can we predict the dependence of the critical intensity on the model parameters?

First, we compare the simulation results to the well-known excluded-area approximation, which uses an empirical parameter [24, 14]. Then, we examine a purely geometric lower bound without any empirical parameter, which is based on a topology index from integral geometry, the Euler characteristic [25]. The simulation results suggest that the lower bound qualitatively captures well how the critical intensity depends on the model parameters. We also suggest new candidates of purely geometric approximations based on further integral geometric shape descriptors.

In Section 1, we simulate anisotropic Boolean models (see Section 1.1) and examine the behavior of continuum percolation both in finite samples and extrapolated to infinite systems. For finite systems, we define the effective percolation threshold in Section 1.2 and extrapolate it to infinite system sizes in Section 1.3. There we also provide several heuristic explanations of the empirical finding of an isotropic percolation threshold. In Section 1.4, we examine how the value of the percolation threshold changes with the anisotropy of the system.

In Section 2, we study analytic approximations of the percolation thresholds based on the excluded area (Section 2.1), the Euler characteristic (Section 2.2), or the second moments of Minkowski functionals (Section 2.3). In Section 3, besides concluding remarks, we given an outlook to open problems and generalizations as well as to how our results suggest techniques for future research [26].

1 Anisotropic “percolation” behavior in finite samples and isotropic thresholds in infinite systems

A percolating cluster in an infinite system is defined as an unbounded cluster; that is, a cluster which is not the subset of a compact observation window. The percolation threshold in the intensity $\gamma_{c}$ is the intensity at which a percolating cluster emerges, more precisely, it is the infimum of all intensities $\gamma$ at which a percolating cluster appears (almost surely). Alternatively, the percolation threshold $\phi_{c}$ in the occupied area fraction can be considered (if the mean area of a typical grain is strictly positive). A phase transition occurs if the percolation threshold takes a non-trivial value $\gamma_{c}\in(0,\infty)$ (e. g., for Boolean models with a bounded grain size). Above the percolation threshold, the probability that the cluster spans an observation window converges in the limit of an infinitely large window to unity.

In a finite sample, no infinite cluster can appear. The percolating cluster is therefore defined as a cluster that connects opposite sides of the sample. In other words, a cluster that percolates in $x$- or $y$-direction spans the system horizontally or vertically, respectively. Figure 1 shows a finite sample in which two clusters span the system in $x$- but none in $y$-direction. Phase transitions occur only in the thermodynamic limit. However, there is an apparent percolation threshold close to which the probability for the appearance of a percolating cluster steeply increases.

If the elongated particles have an orientation bias in $x$-direction, a finite anisotropic sample percolates more likely in $x$- than in $y$-direction. In other words, the occupied volume fraction at which the probability that the system percolates in $x$-direction with a given probability is lower than the corresponding volume fraction in $y$-direction. However, we show in this section that the difference vanishes in the limit of an infinite system size, no matter how anisotropic the system may be. We relate this observation to the uniqueness of the percolating cluster.

1.1 Anisotropic Boolean models

The (isotropic) Boolean model is not only a benchmark model of continuum percolation but also a popular model of heterogeneous materials [27, 28]. Examples include fractured materials or hydrating cement-based materials [29], sedimentary rock [30, 31] and wood composites [32].

Figure 2 visualizes the definition of the Boolean model. First, points are randomly placed in space distributed uniformly. They form a completely independent point process, the so-called Poisson point process. The number of points in a finite observation window follows a Poisson distribution. In the second step, each point is decorated with a grain. The shape and orientation follows the so-called grain distribution. We here only consider stationary processes, that is, statistically homogeneous models for which both intensity and orientation distribution of the grains are spatially constant.

We use grains with a constant shape, namely, rectangles with a fixed size and aspect ratio. To explore a range of different degrees of anisotropy in our simulations, we use a parametric orientation distribution from [12]. The anisotropy can be tuned by adjusting the orientation distribution of the grains. The orientation of a grain is defined by the angle $\theta$ between its main axis and the $x$-axis. Its probability density function is

[TABLE]

see figure 3. The parameter $\alpha$ controls the degree anisotropy of the system: $\alpha=0$ results in a uniform distribution and hence to an isotropic system; the larger $\alpha$ the stronger the anisotropy; we define that $\alpha=\infty$ corresponds to a full alignment with the $x$-axis.

The expected area fraction of the grain phase is called the occupied area fraction $\phi$. It is independent of the orientation distribution and given by

[TABLE]

where $A$ is the area of a single grain [33]. A dependency on the intensity always corresponds to an dependency on the occupied area fraction.

1.2 Connectivity and effective percolation thresholds

Percolation in an infinitely large system is a geometric phase transition. Above the critical occupied area fraction $\phi_{c}$ there is a percolating cluster with probability one, and below this threshold there is almost surely no percolation of the grains. The probability that there is a percolating cluster as a function of the occupied area fraction $\phi$ is called the connectivity $C$. It is a step function that changes its functional value at $\phi_{c}$ from zero to unity, see figure 4.

We are interested in the behavior of finite systems. In this case, there is no phase transition. Depending on statistical fluctuations, a system spanning cluster may appear or it may not appear. For all occupied area fractions $0<\phi<1$ there is a finite probability that there is a percolating cluster. The connectivity is expected to be a smooth and to take on values $0<C<1$.

The connectivity can be estimated from simulated Boolean models simply by the fraction of samples that contain a percolating cluster. We use polygons for an exact representation of the Boolean model to avoid strong pixelation errors that would occur in pixelated Boolean models even for fine pixelations (because unconnected clusters can become connected or vice versa). We simulate anisotropic Boolean models at different expected area fractions and for various sizes of the simulation box, and determine in each case the connectivity. Figure 4 shows (for a Boolean model of aligned rectangles with aspect ratio $1/4$), as expected, a smooth increase of the connectivity for an increasing occupied area fraction $\phi$.

In figure 4, we distinguish the connectivity in $x$- and $y$-direction in a Boolean model with fully horizontally aligned rectangles ($\alpha=\infty$). The connectivity $C^{x}$ is the probability that there is a cluster that spans the system in $x$-direction, and accordingly the connectivity $C^{y}$ considers clusters that span the system in $y$-direction. More particles are needed for vertical than for horizontal percolation, and a fluctuation that spans the system horizontally is more likely than vertically. The system is more likely to percolate horizontally than vertically (at least close to the percolation transition). Even for relatively large system sizes, there are occupied volume fractions at which the system most likely percolates in $x$- but not in $y$-direction.

Since there is no phase transition in a finite system, the percolation threshold $\phi_{c}$ is not well-defined. Instead an effective percolation threshold $\phi_{c}^{\mathrm{eff}}$ must be defined. A simple and intuitive definition is via

[TABLE]

A trivial estimate inverting the estimated connectivity as a function of the occupied area fraction would not be robust against statistical fluctuations. A more robust estimate has been defined in [34]. The derivative of the connectivity ${\mathrm{d}C}/{\mathrm{d}\phi}$ as a function of the occupied area fraction $\phi$ is the probability that the system starts to percolate at $\phi$. Another possibility to define an effective percolation threshold is $\phi_{c}^{\mathrm{eff}}:=\int_{0}^{1}\phi({\mathrm{d}C}/{\mathrm{d}\phi})\rmd\phi$ [4]. Because of the phase transition, the behavior of the system is close to the percolation threshold and for large systems dominated by the universal scaling. If the derivative of the connectivity ${\mathrm{d}C}/{\mathrm{d}\phi}$ is approximated by a Gaussian distribution, the connectivity can be approximated by an error function with a constant offset and prefactor:

[TABLE]

Fitting this rescaled error-function to the connectivity as a function of $\phi$, where $\Delta$ and $\phi_{c}^{\mathrm{eff}}$ are the fit parameters, see figure 4, provides a numerically robust estimate of the effective percolation threshold $\phi_{c}^{\mathrm{eff}}$. In the limit of infinite system size, the error function converges to a step function and $\phi_{c}^{\mathrm{eff}}\rightarrow\phi_{c}$.

We adopt this procedure, but we use the intensity $\gamma$ as an input instead of the occupied area fraction $\phi$. In a finite system, there are statistical fluctuations in $\phi$. To avoid systematic errors, we replace the volume fractions in (4) by the corresponding intensities.

The definition of the effective percolation threshold via (4) has the advantage that it can be easily interpreted and implemented. However, there are two major disadvantages. First, unnecessary samples are simulated especially at high thresholds; they are computationally expensive to simulate because of the large number of grains. Second, the precision is limited by unknown systematic errors. Although the error function is a good approximation of the connectivity as a function of the occupied area fraction, high precision estimates of the connectivity reveal in finite systems statistically significant deviations from an error function. The statistical errors in our simulation are of the order $\mathcal{O}(10^{-4})$, whereas systematic errors are of $\mathcal{O}(10^{-3})$. There are more efficient algorithms using inhomogeneous Boolean models for more precise estimates of the percolation threshold of homogeneous models [e. g., 35, 36, 37]. In future research, these algorithms could be generalized to anisotropic systems. For our analysis of anisotropic continuum percolation, the procedure based on (4) is sufficient.

An error function is fitted to each curve in figure 4. The resulting estimates of the effective percolation thresholds are depicted by vertical dashed lines. The effective percolation thresholds in $x$-direction are distinctly smaller than in $y$-direction for all simulated systems sizes $L$. However, the differences in the effective percolation thresholds depend on the system sizes. The differences for horizontal and vertical directions become smaller with increasing system size.

1.3 Isotropic percolation threshold

Because percolation is a phase transition, the scaling of the effective percolation thresholds with the system size is universal for large system sizes. The scaling exponent does not depend on the grain distribution, it is even the same for continuum percolation and lattice percolation; see [e. g., 38, 39]:

[TABLE]

with the universal critical exponent $\nu=\frac{4}{3}$ [4].

The percolation threshold $\phi_{c}$ can be estimated from the effective thresholds $\phi_{c}^{\mathrm{eff}}$ by an extrapolation from the finite system sizes $L$ to $L\rightarrow\infty$. Figure 5 shows for an isotropic Boolean model with rectangles with aspect ratio $b/a=\frac{1}{4}$ the numerical estimates of $\phi_{c}^{\mathrm{eff}}$, which are derived from fits of error functions to the connectivity (in $x$-direction). It also plots the extrapolation to $\phi_{c}$ and the corresponding bands of one or two standard deviations.

For anisotropic models, we have to consider that the proportionality constant is nonuniversal. It varies for the effective percolation thresholds $\phi_{c}^{\mathrm{eff},x}$ in $x$- and $\phi_{c}^{\mathrm{eff},y}$ in $y$-direction. Therefore, we a priori distinguish also in the limit of infinite system size the percolation thresholds in $x$- or in $y$-direction. In other words, we estimate the thresholds in $x$- or in $y$-direction separately:

[TABLE]

However, we will show that even for the most anisotropic systems $\phi_{c}^{x}\equiv\phi_{c}^{y}$.

We simulate Boolean models with rectangles for different values of $\alpha$, that is, for different degrees of anisotropy. Their aspect ratios ${b}/{a}$ are fixed within one model, but we simulate several models with different ratios ${b}/{a}\in\{1,\frac{3}{4},\frac{1}{2},\frac{1}{4}\}$ and anisotropy parameters $\alpha\in\{0,1,3,15,\infty\}$. The models thus vary from isotropic orientation distributions to full alignment. For each system, we determine the connectivity $C$ as a function of the intensity for differently large simulation boxes with linear system sizes $L\in\{40a,50a,60a,70a,75a,80a,85a,95a\}$. As discussed above, the more anisotropic a Boolean model the larger the simulation boxes must be so that the connectivity is accurately described by an error function (so that the finite size scaling can be applied). For each system, we fit error functions separately to the both the connectivity in $x$- and in $y$-direction. The fit is applied to the five largest system sizes which are simulated for the according Boolean model. We thus estimate the effective percolation thresholds $\phi_{c}^{\mathrm{eff},x}$ and $\phi_{c}^{\mathrm{eff},y}$. Next we extrapolate the effective percolation thresholds to the limit of infinite system size using (7) with fit parameters $\phi_{c}^{x}$ and $m^{x}$ or $\phi_{c}^{y}$ and $m^{y}$.

Figure 6 shows for the six exemplary Boolean models the effective percolation thresholds $\phi_{c}^{\mathrm{eff},x}$ and $\phi_{c}^{\mathrm{eff},y}$ as a function of the rescaled system size. It also plots the extrapolation and the corresponding bands of one or two standard deviations. Within statistical significance, the percolation threshold is isotropic. Even for the most anisotropic systems, we find $\phi_{c}^{x}\equiv\phi_{c}^{y}$. Similarly, Balberg et al. [13, 14] found isotropic percolation thresholds in two- and three-dimensional anisotropic systems of sticks.

This isotropy of the percolation threshold, which means that the Boolean model either percolates simultaneously in all directions or not at all, can be expected for Boolean models for at least three reasons: First, consider the extreme case of aligned rectangles. A dilation can transform a Boolean model with aligned squares into a Boolean model with aligned rectangles of arbitrary aspect ratios. However, the affine transformation neither changes the area fraction nor the topology of the model. Therefore, the percolation threshold for aligned rectangles must be independent of the aspect ratio. Moreover, because the model can either be dilated in $x$- or $y$-direction, the percolation threshold must be isotropic. Balberg and Binenbaum [13] provide a similar argument for the anisotropic systems of sticks. They also discuss a renormalization-group argument why in an anisotropic square lattice an isotropic percolation threshold is expected.

Another plausibility argument for the isotropic percolation threshold is that, as indicated in Section 1.2, percolation in $x$-direction is almost surely not possible along a straight line but only along a winded path. This requires that also grains with different $y$-coordinates are connected.

The most fundamental reason for the isotropy of the percolation threshold is that it is closely related to the uniqueness of the percolating cluster in continuum percolation. Meester and Roy [41] have proven for Boolean models of spheres with varying size that there can be at most one unbounded component, but the results can be generalized to grains with general shape, see also [7]. Either with probability one or zero there is either exactly one or no percolating cluster, in other words, the percolating cluster is (almost surely) unique. Similar proofs could be derived for other percolation models [e. g., 42, 43].

The isotropy of the percolation threshold follows by a heuristic argument from the uniqueness of the percolating cluster. If we assume that the unbounded cluster spans the system only in $x$- but not in $y$-direction, two percolating clusters that are not connected to each other can appear with finite probability, see figure 7. Assume that the percolating (dark blue) cluster is bounded in $y$-direction, that is, the cluster lies completely below some bound (maximum of red line). Because of the statistical homogeneity of the system, the probability that there is another percolating cluster in the domain above the bound (colored cyan in figure 7) should be nonzero. However, this would lead to a contradiction of the uniqueness of the unbounded cluster. Therefore, we conjecture that in a statistically homogeneous system with almost surely at most one unique unbounded component this percolating cluster must be unbounded both in $x$- and in $y$-direction, and the percolation thresholds must coincide:

[TABLE]

1.4 Threshold estimates for models with varying degrees of anisotropy

The percolation threshold is by plausibility argument and empirical findings confirmed to be isotropic. So the finite size scaling from (7) can be replaced by

[TABLE]

where the percolation threshold in the infinite system $\phi_{c}$ is isotropic. It is because of the nonuniversal prefactors $m^{x}$ and $m^{y}$ that the effective percolation thresholds differ in $x$- and $y$-direction.

We simultaneously extrapolate the effective percolation thresholds $\phi_{c}^{\mathrm{eff},x}$ and $\phi_{c}^{\mathrm{eff},y}$ to the same $\phi_{c}$ using a maximum likelihood fit described in A. The final results for the isotropic percolation thresholds $\phi_{c}$ are plotted in figure 8.

If the grains are aligned, the percolation threshold does not depend on the aspect ratio of the grains, because in this special case changes in the aspect ratio of the single grains correspond to dilations of the Boolean model, as mentioned above. The latter do not change the topology and therefore not the percolation threshold. For all other orientation distributions of the grains, the percolation threshold decreases with a decreasing ratio of smaller to larger side length. Put differently, the more elongated a rectangle with rotational degrees of freedom the smaller the percolation threshold. In the limit of stick percolation, that is, for vanishing aspect ratio $b/a=0$, the percolation threshold vanishes as well $\phi_{c}=0$. Comparing different orientation distribution, we find in agreement with previous findings [14], the more anisotropic the orientation distribution, the larger is the percolation threshold.

2 Approximations and bounds on percolation thresholds

Numerical estimates of the percolation threshold can get computationally very expensive. However, applications might need the threshold for some specific model parameters for which no simulation results are yet available. A possibility to predict the percolation threshold without the need for expensive simulations would be of great importance; that is, an explicit formula approximating the percolation threshold. Such an explicit formula should be easy computable and as accurate as possible [44]. An example is the scaling relation for varying dimensions from [45], which is based on a lower bound [46] and a conjecture that no convex shape provides higher thresholds than the hyperspheres.

A common approximation is the excluded area approximation (or excluded volume in three dimensions) [24, 14, 47, 48, 49, 50, 51, 23]. It tries to predict percolation thresholds based on simulation results of similar models by approximating the dependence of the critical intensity on the model parameters by the change in the excluded area. The excluded area for convex objects can be expressed explicitly by a so-called mixed intrinsic volume from integral geometry. In section 2.1, we compare the approximation of the excluded area to the percolation thresholds for the models with different degrees of anisotropy discussed above.

In contrast to the excluded area approximation, which uses an empirical parameter, Mecke and Wagner [25] suggested a purely geometrical approximation. It uses an additive topological constant from integral geometry, the Euler characteristic. For a compact set, it is given by the number of connected components minus the number of holes. The critical intensity is approximated by the zero of its mean value (as a function of the intensity). The Euler characteristic is positive for intensities below this zero and negative above. Intuitively speaking, this indicates a change from a predominant role of single clusters to a network-like structure. Explicit formulas for the mean Euler characteristic are known for very general Boolean models [52], that potentially allow for estimates of the percolation threshold. At least for the parametric Boolean models studied here, the approximation is a lower bound, see Section 2.2. We expect it to be a lower bound for even more general Boolean models, because the Euler characteristic does not distinguish holes in connected clusters from the percolating void phase.

Moreover, we find that the trend, i. e., the qualitative behavior, is very well described by the zero of the Euler characteristic. For not too small aspect ratios of the rectangles, this finding allows us to very accurately estimate the percolation thresholds of Boolean models if the percolation threshold is known for squares that follow the same orientation distribution as the rectangles.

The Euler characteristic belongs to a class of so-called Minkowski functionals (or intrinsic volumes) from integral geometry [53, 33, 11]. In two-dimensions, the two other Minkowski functionals are area and perimeter. We investigate the potential of further approximations of the percolation threshold that do not use an empirical parameter, namely, the extremal points of the variances and covariances of the Minkowski functionals. In Section 2.3, we given an outlook on candidates for close upper bounds on the percolation thresholds.

2.1 Excluded area approximation

The excluded area is defined for two particles $K$ and $Q$. It is the area of the region defined by all those center positions of the second particle $Q$ which leads to intersection of the two particles [54]. In the following, the number of grains that intersect an object is called its number of bonds. In other words, the second particle $Q$ cannot enter this area without intersecting $K$. Figure 9 visualizes the excluded area of two rectangles.

The excluded area $A_{\mathrm{ex}}$ of two convex grains $K$ and $Q$ can be explicitly expressed using the so-called mixed functional $V^{0}_{1,1}$ from integral geometry [33],

[TABLE]

Note that our notation differs from [24], in which the same quantity is denoted by $\langle A\rangle$, but $\langle A_{\mathrm{ex}}\rangle$ indicates the total excluded area. Averaging $A_{\mathrm{ex}}(K,Q)$ with the orientation distribution $\mathcal{P}(\theta)$ of the grains $K$ and $Q$ yields the averaged excluded area for single grains

[TABLE]

where $\vartheta(\theta)$ denotes a rotation by the angle $\theta$.

For the isotropic orientation distribution $\mathcal{P}(\theta)\equiv\mathrm{const.}$, the average mixed intrinsic volume separates. In this special case, it is proportional to the product of the perimeters of the grains. For arbitrary dimensions, a general formula for the exclusion volume of identical randomly oriented convex particles is given is given in [45]. Therefore, the excluded area in two dimensions can, as it is well-known [55], be expressed by the areas $A$ and the perimeters $P$ of the single grains

[TABLE]

Also for rectangles with arbitrary orientations, there is an explicit formula for the mixed intrinsic volume in (11) [12, p. 82]. The excluded area of two rectangles $R$ with orientations $\theta_{i}$ and $\theta_{j}$ and side lengths $a$ and $b$ is given by

[TABLE]

This equation allows for an explicit expression approximating the percolation thresholds of Boolean models with rectangles.

The product of the critical intensity $\gamma_{c}$ and the average excluded area $\langle A_{\mathrm{ex}}\rangle$ is called the total excluded area. It is equivalent to the average number of bonds per object at the critical intensity $B_{c}$, because the average number of bonds of a grain is equal to the average number of intersecting grains and thus equal to the average number of particles within the excluded area [24]

[TABLE]

which was confirmed numerically [e. g., 24, 14, 50].

Balberg et al. [24] proposed that $B_{c}$ is, at least approximately, an invariant quantity for similar systems, following a similar argument from lattice percolation [56]. Actually, the average number of bonds per object $B_{c}$ is not a universal constant but depends on the grain shape and orientation distribution. However, it varies to reasonable approximation only slightly for similar systems [e. g., 24, 40]. If the average number of bonds per object at the critical intensity $B_{c}$ is empirically known (e. g., for a Boolean model of squares) and assumed to be the same for similar systems (e. g., Boolean models with rectangles following the same orientation distribution) the critical intensity of these similar models can be approximated by the explicit expresssion

[TABLE]

In numerous simulation studies, the assumption of a system invariant total excluded area is shown to be a useful approximation [49, 51, 23] for which the critical intensity scales correctly with the particle size [50] and which approximates the functional dependence of the percolation threshold on the anisotropy of the model [24, 14].

The excluded area for overlapping rectangles with anisotropic orientation distribution is also calculated in [23], where the orientations are random within the interval $[-\alpha,\alpha]$ (following [24]). The approximation is compared to Monte Carlo simulations of isotropic models. For increasing anisotropy, an increase of the percolation threshold is predicted, which we have confirmed by our simulations above.

Here we compare the excluded area approximation to the Monte Carlo estimates of the anisotropic parametric Boolean model studying the variety of models with different degrees of anisotropy from Section 1.

At the critical intensity, the average number of bonds $B_{c}$ for Boolean models can be derived from the percolation threshold $\phi_{c}$ using (2) and (15):

[TABLE]

where $A$ is the area of a single grain. Using the numerical results for the percolation thresholds $\phi_{c}$ from Section 1, we estimate $B_{c}$ for overlapping squares for the different orientation distributions. Using (16), we finally approximate the critical intensities of the overlapping rectangles with varying aspect-ratio.

Figure 10 compares this excluded area approximation to the simulation results from figure 8. For aligned rectangles, the independence of the percolation thresholds from the aspect ratio is correctly reproduced, because the excluded area of aligned rectangles $A_{\mathrm{ex}}(R,R)=4ab$ is proportional to the area of a single grain (which is here chosen to be unity). For the other anisotropic or isotropic Boolean models, the excluded area approximation seems to be a close upper bound for the rectangles if the empirical parameters of a square following the same orientation distribution are used. However, there are significant deviations from the prediction, for example, of the percolation threshold of randomly orientated rectangles with aspect ratio $\frac{1}{4}$.

2.2 Euler characteristic approximation

In contrast to the excluded area approximation that requires the input of an empirical parameter ($B_{c}$), Mecke and Wagner [25] suggested a purely geometrical approximation, namely, the zero of the mean Euler characteristic as a function of the intensity. They were motivated by the idea that a vanishing Euler characteristic implies a balance between clusters and holes. Moreover, for self-dual matching lattices, like the square-lattice, the zero of the Euler characteristic and the critical point coincide exactly [57]. They compared the critical intensity $\gamma_{c}$ to the zero of the Euler characteristic $\gamma_{0}$, for example, for overlapping discs and spheres, and suggested that $\gamma_{0}$ could be a close bound for $\gamma_{c}$; see also [58].

Mecke and Seyfried [59] found in Monte Carlo simulations of a large class of discretized Boolean models a good qualitative agreement, in that the dependence on the shape of the constituents is captured well, but they found also quantitative deviations. Neher et al. [44] analyzed two-dimensional lattice graphs, for which the Euler characteristic provides a good approximation of the critical point. The criterion has, for example, been applied to study the percolation in level sets of Gaussian random fields [60, 61], biopolymeric networks [62], flow in heterogeneous soil [63], or surface roughness of medical implants [64]. Moreover, the connection between percolation and the Euler characteristic triggered predictions of flow in porous media that are based on the Minkowski functionals [65, 66].

More precisely, we here refer to the density of the Euler characteristic $\overline{\chi}$, that is, before taking the limit of an infinitely large observation window the functional value is rescaled by the size of the window [33]. Explicit formulas are known for the density of the Euler characteristic of Boolean models as a function of the intensity [52]. The criterion is thus easily applicable to a great variety of continuum percolation models.

Here we calculate this approximation for the Boolean models with varying anisotropy that are studied above, and we compare it to the simulation results for the percolation thresholds from Section 1. We investigate how well the dependence of the threshold on the grain shape and orientation distribution is captured by the zero of the Euler characteristic. The density $\overline{\chi}$ of the Euler characteristic of a Boolean model with rectangles as a function of the intensity $\gamma$ is given by

[TABLE]

see [52], where $\langle{V}^{0}_{1,1}\rangle$ is the average of the mixed intrinsic volume ${V}^{0}_{1,1}$ of two rectangles,

[TABLE]

see (11) and (12). The zero of the Euler characteristic then follows in a straightforward way:

[TABLE]

which appears similar to the excluded area approximation in (16). The occupied area fraction $\phi_{0}$ at which the mean Euler characteristic changes its sign is according to (2) then given by

[TABLE]

where $A$ is the area of a single grain.

Figure 11 compares $\phi_{0}$ from (23) to the estimates of the isotropic percolation thresholds $\phi_{c}$ from figure 8. At least for these Boolean models, the zero of the Euler characteristic is a lower bound on $\phi_{c}$. There is a significant quantitative deviation, but similar to the findings for discretized Boolean models [59], we find that the qualitative behavior of $\phi_{c}$, i. e., the dependence on anisotropy and grain shape, is well described by $\phi_{0}$.

To further investigate this accurate qualitative description, we plot in figure 12 the difference of the critical intensity of rectangles and the critical intensity of squares with the same orientation distribution. We compare these differences to the differences of the corresponding zeros of the Euler characteristic. For the Boolean models studied here, we find a remarkable agreement. Only for an intermediate regime of strongly preferred orientation, but no full alignment ($\alpha=25$), we find distinct (but still small) deviations. However, there is an excellent agreement of the difference in the critical intensities to the difference of the zeros of the Euler characteristic not only for the special case of full alignment but also for randomly oriented rectangles. If the critical intensity of squares is used as an empirical parameter, similar to the excluded area approximation, the zero of the Euler characteristic allows at least for these Boolean models a quite precise prediction of the percolation threshold, compare figure 10.

Like the critical average number of bonds $B_{c}$, the offset between the critical intensity $\gamma_{c}$ and the zero $\gamma_{0}$ of the Euler characteristic depends on the shape of the grains; it is different for randomly oriented squares, 0.19688(1) [67], and discs, 0.128085(2) [68]. Note, however, that assuming a constant offset to the number of grains per area of a single grain would lead to a nonphysical behavior in the limit of vanishing aspect ratios. For grains with unit area, both the critical intensity and the zero of the Euler characteristic trivially vanish in the limit of extreme anisotropy, i. e., for sticks, because of the choice of units. For example, for randomly oriented sticks of length unity, the zero of the Euler characteristic ($\gamma_{0}=\pi$) is as expected a lower bound on the critical intensity ($\gamma_{c}=5.6372858(6)$ [69]).

2.3 Second-moment approximations

The Euler characteristic is one of the so-called Minkowski functionals (or intrinsic volumes) from integral geometry [53, 33, 11]. In two-dimensions, the Minkowski functionals are the area $W_{0}$, the perimeter $W_{1}$, and a functional proportional to the Euler characteristic $W_{2}=2\pi\chi$. These three Minkowski functionals characterize, according to Hadwiger’s theorem, the complete “additive shape information” of a compact and convex set $K$ in the following sense: they form a basis of all functionals $F$ that are motion invariant, continuous, and additive [70, 33]. The latter means that $F(K\cup L)=F(K)+F(L)$ for two disjoint sets $K$ and $L$. The second moments of Minkowski functionals and generalizations thereof are known analytically [71, 72, 73, 74]. Remarkably, a behavior similar to thermodynamical quantities in statistical physics has been found for the second moments of Minkowski functionals [53, 75].

As an outlook, we here examine the prospects of further approximations of the percolation threshold based on the second moments of the Minkowski functionals without using empirical parameters. We compare the extremal points of the variances and covariances of the Minkowski functionals to the empirical estimates of the critical intensity.

An interesting example, which can serve as a motivation of this approach, is cell percolation in planar Poisson Voronoi tessellations. Voronoi cells of a Poisson point process are colored black with some probability $p$, and a percolating cluster is an unbounded and connected collection of black cells. The critical probability of this model of continuum percolation is $p_{c}=\frac{1}{2}$ [76]. Last and Ochsenreither [77] have proven that this coincides with the global maximum of the asymptotic variance of the Euler characteristic. Moreover, Hug et al. [73] found for overlapping discs that both the local minimum of the variance of the perimeter and the local minimum of the covariance of area and Euler characteristic are good approximations of the percolation threshold. For example, the local minimum of the asymptotic covariance $\sigma_{0,2}$ is found at $\gamma\approx 1.15$ compared to the critical intensity $\gamma_{c}\approx 1.13$ [36], where the unit of area is again defined by the grain.

In our outlook here, we compare for four different Boolean models with rectangles the extremal points of the variance-covariance structure of the Minkowski functionals to the percolation threshold. We investigate whether the second moments could allow for accurate predictions of the thresholds without the need for empirical parameters, especially testing whether the local minima of $\sigma_{0,2}$ or $\sigma_{1,1}$ are close to the percolation threshold for anisotropic models.

Figures 13 and 14 show the different asymptotic variances or covariances for a Boolean models with either squares or rectangles (with aspect ratio $\frac{1}{2}$). These are either fully aligned or follow an isotropic orientation distribution. For these four Boolean models very accurate estimates of the percolation thresholds can be found in literature [40, 78, 45, 67]. We use these precise estimates with at least four significant number of digits because we indeed find some remarkable agreement between the extremal points and the percolation threshold.

Like for the overlapping discs, a local minimum point of the variance of the perimeter is an excellent approximation of the percolation threshold for aligned squares or rectangles with aspect ratio $\frac{1}{2}$. The critical intensity of aligned squares is $\gamma_{c}=1.0982(3)$ [45], the local minimum point of the asymptotic variance of the perimeter is $\gamma_{m}[{b}/{a}=1]\approx 1.0994$, which is derived by numerically minimizing $\sigma_{W_{1}W_{1}}$ using [79]. While this is a remarkably good approximation, there is a statistically significant difference. Moreover, while the critical intensity must not depend on the aspect ratio of the aligned rectangles, see Section 1.3, the local minimum point of the asymptotic variance of the perimeter slightly depends on the aspect ratio. For an aspect ratio of $\frac{1}{2}$, the local minimum point is at $\gamma_{m}[{b}/{a}=\frac{1}{2}]\approx 1.0952$, and for an aspect ratio of $\frac{1}{10}$, it is $\gamma_{m}[{b}/{a}=\frac{1}{10}]\approx 1.0727$. The local minimum point of the variance is also neither an upper nor a lower bound for the percolation threshold. For aligned rectangles, the minimum points can be below the critical intensity. However, the minimum points are above the critical intensities for aligned squares or for the randomly oriented squares and rectangles, see figure 13.

Another candidate for an approximation of the percolation threshold is the local minimum point of the covariance of area and Euler characteristic. Figure 14 shows that it is a very good approximation for aligned squares and rectangles with aspect ratio $\frac{1}{2}$ and a close upper bound if the grains are randomly oriented.

Only for one of the Boolean models, a vanishing correlation can serve as a good approximation of the percolation threshold: for randomly oriented squares the first zero of the covariance of perimeter and Euler characteristic ($\propto\sigma_{W_{1}W_{2}}$) seems to be in good agreement with the critical intensity, see figure 14. However, for randomly oriented rectangles there is a distinct difference. Moreover, for aligned squares or rectangles, as well as for overlapping discs, this zero crossing from negative to positive correlation does not exist.

3 Conclusion and Outlook

We have investigated percolation in anisotropic systems, both for finite simulation boxes and in the limit of infinite system size. We have especially elucidated how the anisotropy of the Boolean model, which is a geometric property, influences percolation, which is a topological property. There is both universal and non-universal behavior. The isotropy of the percolation threshold is of the first kind. In infinite systems, there is no difference in effective percolation thresholds for different directions. Independent of the details of the model, the difference between the effective percolation thresholds observed in finite systems vanishes in the thermodynamic limit. This has here been related to the uniqueness of the percolating cluster, see figure 7.

In a finite sample of an anisotropic Boolean model, a cluster is more likely to span the system in the preferred than in the perpendicular direction, see figure 1. This anisotropic percolation in finite samples is demonstrated by determining the effective percolation thresholds, see (3) and (4), by a fit of an error function to the connectivity, see figure 4. Although a distinct difference remains even for very large system sizes, the extrapolation to infinite system sizes shows that even the most anisotropic model percolates simultaneously in all directions, see figure 6. This finding we here refer to as the isotropy of the percolation threshold. Our results are in agreement with previous findings [13, 14, 80, 18].

For Boolean models with different degrees of anisotropy, we have estimated the isotropic percolation thresholds by a simultaneous fit of the finite size scaling in $x$- and $y$-direction, see (10)–(37). While the isotropy of the percolation threshold is universal, its value is non-universal, that is, it depends on the grain distribution (or the degree of orientation bias). In agreement with previous findings [14], the more anisotropic the orientation distribution, the larger is the percolation threshold, see figure 8.

There are efficient algorithms that use inhomogeneous Boolean models for high-precision estimates of percolation thresholds in homogeneous and isotropic models, e. g., see [35, 36, 37]. The isotropy of the percolation threshold implies that these algorithms can be directly generalized to anisotropic systems. Different directions of the gradient in the intensity should yield (within statistical errors) the same estimate of the percolation threshold. Such methods could strongly increase the precision of estimates of the percolation thresholds in anisotropic models.

In two dimensions, the percolation of grain or void phase are complementary phenomena, that is, there is (almost surely) a percolating cluster of the void phase at occupied area fractions below the threshold, but it vanishes (almost surely) above the threshold. We can therefore directly apply the numerical results and the approximations also to void percolation and the critical porosity.

The Minkowski functionals help to study the relation between geometry, e. g., quantified by area or perimeter, and topology, e. g., quantified by the Euler characteristic. They thus allow for insights into percolation as a topological phase transition. Their generalization, the mixed volumes, provide an explicit formula for the well-known excluded area approximation of percolation thresholds, see (11)–(16). If the percolation thresholds of squares with different orientation distributions are used as empirical parameters, we find that the excluded area approximation is in relatively good agreement with the percolation thresholds of the rectangles, see Fig. 10. For the systems studied here, it provides a close upper bound, which implies that the critical average number of bonds per grain slightly decreases for the rectangles compared to squares with the same orientation distribution.

In contrast to the excluded area approximation, the zero of the Euler characteristic serves as a purely geometric approximation of the percolation threshold. We show that it provides a lower bound that captures well the qualitative dependence on the system parameters, see figure 11. If an empirical parameter is introduced (similar to the excluded area approximation), it allows for quite precise predictions of $\gamma_{c}$ within a narrow range of parameters studied here, see figure 12.

Using the second moments of the Minkowski functionals, we suggest and compare new candidates for accurate approximations of the percolation thresholds, see Figs. 13 and 14. Some are found to be in surprisingly good agreement for the four systems that are tested here, for example, the minimum point of the covariance of area and Euler characteristic. More examples and further research are needed to test this suggestion. An important prospect for applications is that no empirical parameters would be needed because analytic formulas are known for the second moments of the Minkowski functionals [73, 74]. Similar techniques could also easily be applied to percolation on lattices for which the covariances of the Minkowski functionals can be calculated straightforwardly [81, 26]. Moreover, our results can also lead to more fundamental questions. Why is there such a good agreement between the percolation threshold and some of the extremum points of the second moments of the Minkowski functionals as a function of the intensity? Is it a coincidence or is there a connection between local fluctuations and the global topology?

We thank the German Research Foundation (DFG) for the Grants No. HU1874/3-2, No. LA965/6-2, and No. ME1361/11 awarded as part of the DFG-Forschergruppe FOR 1548 “Geometry and Physics of Spatial Random Systems.”

Appendix A Simultaneous extrapolation of vertical and horizontal effective percolation thresholds

Based on the finite size scaling in (10), we determine the isotropic percolation threshold $\phi_{c}$ via a simultaneous extrapolation of the effective percolation thresholds in $x$- and in $y$-direction. More precisely, we estimate the percolation threshold via a maximum likelihood fit of the parameters $\phi_{c}$, $m^{x}$, and $m^{y}$.

The effective percolation thresholds in $x$- and in $y$-directions could be correlated, but the correlation cannot be considered in this analysis. We neglect the correlation of $\phi_{c,i}^{\mathrm{eff},x}$ and $\phi_{c,i}^{\mathrm{eff},y}$ in the simulation and use the likelihood function

[TABLE]

with

[TABLE]

where for each system size $L_{i}$ the error of the effective percolation thresholds $\phi_{c,i}^{\mathrm{eff},x}$ and $\phi_{c,i}^{\mathrm{eff},y}$ are given by $\sigma_{x,i}$ and $\sigma_{y,i}$, respectively. The log-likelihood function can be written as

[TABLE]

with

[TABLE]

and

[TABLE]

where the covariance matrix $\mathbf{Cov}$ is assumed to be diagonal

[TABLE]

Therefore, the maximum likelihood estimation equals a minimum least square fit with block matrices.