Continuum percolation expressed in terms of density distributions

TL;DR

This paper introduces a new integral equation approach to analyze the connectivity in pairwise interacting systems, linking pair connectedness with density distributions and pair correlation functions, applicable across various dimensions and interaction types.

Contribution

It develops a novel formalism relating pair connectedness to density distributions, extending analytical solutions to higher dimensions and complex interactions.

Findings

Analytical solutions for pair connectedness in one-dimensional systems with specific potentials.

Extension of the formalism to higher-dimensional systems like the 3D ideal gas.

Framework accommodating external fields and long-range interactions.

Abstract

We present a new approach to derive the connectivity properties of pairwise interacting n-body systems in thermal equilibrium. We formulate an integral equation that relates the pair connectedness to the distribution of nearest neighbors. For one-dimensional systems with nearest-neighbor interactions, the nearest-neighbor distribution is, in turn, related to the pair correlation function g through a simple integral equation. As a consequence, for those systems, we arrive at an integral equation relating g to the pair connectedness, which is readily solved even analytically if g is specified analytically. We demonstrate the procedure for a variety of pair-potentials including fully penetrable spheres as well as impenetrable spheres, the only two systems for which analytical results for the pair connectedness exist. However, the approach is not limited to nearest-neighbor interactions in one dimension. Hence, we also outline the treatment of external fields and long-ranged interactions, and we illustrate how the formalism can applied to higher-dimensional systems using the three-dimensional ideal gas as an example.

Figures

Equations

$${\rho }^{2}P({\mathbf{r}}_i,{\mathbf{r}}_j)\mathrm{d}{\mathbf{r}}_i\mathrm{d}{\mathbf{r}}_j$$

$$g({\mathbf{r}}_i,{\mathbf{r}}_j)p({\mathbf{r}}_i,{\mathbf{r}}_j):=P({\mathbf{r}}_i,{\mathbf{r}}_j).$$

$$f({\mathbf{r}}_i,{\mathbf{r}}_j)=\exp (-\beta V({\mathbf{r}}_i,{\mathbf{r}}_j))-1,$$

$$g(i,j)=1+f(i,j)+(f(i,j)+1)\sum _{n=1}^{\infty }{\rho }^n{\beta }_n,$$

$$e(i,j)$$

$$=:$$

$$f^{†}(i,j)=e^{†}(i,j)f^{\ast }(i,j)=e^{\ast }(i,j)-1.$$

$$P(i,j)$$

$$=$$

$$P(i,j)$$

$$P(i,j)=C^{†}(i,j)+\rho \int \mathrm{d}k{C}^{†}(i,k)P(k,j).$$

$$P(r)=I(r)+\lambda {\int }_a^r\mathrm{d}xK(r,x)P(x),$$

$$P(r)=I(r)-\lambda {\int }_0^r\mathrm{d}xH(r,x)I(x)$$

$$H(r,x)=-\sum _{n=0}^{\infty }{\lambda }^n{K}_{n+1}(r,x)$$

$$K_1(x,y)$$

$$K_{n+1}(x,y)$$

$$H_{\rho }(r,x)=-K(r,x)+\lambda {\int }_x^r\mathrm{d}zK(r,z)H_{\rho }(z,x),$$

$$P(r)=p(r)g(r),$$

$$\omega (\tau )=\rho g(\tau ).$$

$${\omega }^{\prime }(\tau )=\omega (\tau )\mathcal{P}((0,\tau )\mathrm{empty}|\mathrm{particle}\mathrm{at}\tau ),$$

$$g^{(n+1)}(r_1,...,r_{n+1})=g^{(n)}(r_1,...,r_n)g^{(2)}(r_n,r_{n+1})$$

$$I_1(r):={\int }_0^r\mathrm{d}{r}_1{g}^{(3)}(0,r_1,r){\rho }^{(1)}(r_1)$$

$$=\rho {\int }_0^r\mathrm{d}{r}_1{g}^{(2)}(0,r_1)g^{(2)}(r_1,r).$$

$$I_2(r):={\rho }^2{\int }_0^r\mathrm{d}{r}_1{\int }_{r_1}^r\mathrm{d}{r}_2{g}^{(4)}(0,r_1,r_2,r)$$

$$={\rho }^2{\int }_0^r{\int }_{r_1}^r\mathrm{d}{r}_1\mathrm{d}{r}_2{g}^{(2)}(0,r_1)g^{(2)}(r_1,r_2)g^{(2)}(r_2,r).$$

$${\omega }^{\prime }(0,r)=g(0,r)-\rho {\int }_0^r\mathrm{d}xg(0,x){\omega }^{\prime }(x,r),$$

$$g(0,r)={\omega }^{\prime }(0,r)+\rho {\int }_0^r\mathrm{d}x{\omega }^{\prime }(0,x)g(x,r)$$

$$P(0,r)=$$

$$+$$

$$P(r)$$

$$=$$

$$I(r)=\rho {\int }_{r-d}^d\mathrm{d}x{\omega }^{\prime }(r-x)g(x)$$

$$K(r)$$

$$I(r)$$

Full text

I Introduction

Clustering of particles into connected aggregates is a process that occurs frequently in nature as well as in materials processing. The conditions under which a cluster becomes system spanning are of particular technological interest, as such a cluster might support mechanical stress (e.g. in the case of gels) or transport charges (e.g. in the case of conductive particles immersed in an insulating host matrix). The parameters associated to the emergence of a system-spanning cluster define the percolation threshold Bollobás and Riordan (2006); Stauffer and Aharony (2018), the calculation of which for different systems has been subject of a vast number of studies. Accordingly, a rich methodology has evolved ranging from simulation Seaton and Glandt (1987); Lee and Torquato (1988); Rintoul and Torquato (1997); Miller and Frenkel (2003); Consiglio et al. (2003) over liquid state theory Xu and Stell (1988); Coniglio et al. (1977); DeSimone et al. (1986); Kyrylyuk and van der Schoot (2008); Chiew and Stell (1989); Chatterjee (2000) to renormalization group techniques Cardy (1996), stochastic tools Meester and Roy (1996); Grimmett (1999), and conformal field theory Smirnov (2001); Cardy (1992). Nevertheless, exact (and non-trivial) results are only known for a couple of discrete systems, e.g. random graphs Bollobás (2001) or specific two-dimensional lattice systems Sykes and Essam (1964); Smirnov and Werner (2001).

In disordered systems such as complex liquids, the interplay between liquid structure and connectivity is non-trivial. As a consequence, theories that are general enough to be applicable to a variety of different systems but still allow for immediate computation and sensible estimates of the key quantities, are rare. Thus, the major part of recent work on percolation is dedicated to tailor-made approaches for specific systems, for instance rod-like systems Kyrylyuk and van der Schoot (2008); Drwenski et al. (2017); Meyer et al. (2015); Nigro et al. (2013); Kale et al. (2015); Mutiso et al. (2012); Jadrich and Schweizer (2011); Schilling et al. (2015), which have attracted particular attention due to their use as fillers in composite materials.

Although percolation itself is trivial in one-dimension if the connectivity range of a each individual particle remains finite Schulman (1983), resolving the distance dependence of the probability of two particles being part of the same cluster is not trivial. As this quantity exists in any dimension, an exact solution of a one-dimensional problem provides the perfect benchmark for a more general formalism. However, to our knowledge, even in one dimension the connectivity problem has been solved exactly only for two systems: the ideal gas of non-interacting particles Domb (1947) and the system of impenetrable hard rods Vericat et al. (1987); Drory (1997). These cases were cracked in completely different ways, each tailored to the specific system. However, both solutions can be obtained straightforwardly from the same integral equation as we show below.

In this work, we derive an exact integral equation for the pair connectedness for given arbitrary pair-interactions. In one-dimensional systems with nearest-neighbor interactions, this integral equation requires only the pair-density as input and no approximations. For higher dimensions and long-ranged interactions, a closure relation is required, however, for an intuitively accessible quantity. We demonstrate the virtue of this perspective for the three-dimensional ideal gas.

The fundamental aim of our considerations is to provide a framework that links connectivity properties to thermal distribution functions. As approximate pair distributions are known for many interaction potentials either in analytical form or from experiments and computer simulations, an immediate link to connectivity functions is of considerable practical value. We demonstrate how thermal distribution functions can be used as input to a computational scheme that yields the corresponding connectivity properties.

We revisit established integral equations and summarize the necessary basics in Section II. In Section III we work out our framework for nearest-neighbor interactions and discuss the two analytically known test cases of non-interacting ideal and impenetrable hard-core particles as well as a numerical example. Generalizations of our approach to external fields, long-ranged interactions and higher dimensions, including an analysis of three-dimensional fully penetrable spheres, are discussed in Section IV.

II Established Integral Equations

The quantity we focus on within this paper is the pair connectedness $P({\boldsymbol{r}_{i}},{\boldsymbol{r}_{j}})$ as, for instance, defined by Coniglio Coniglio et al. (1977) by demanding that

[TABLE]

describes the absolute probability to find particles $i$ and $j$ within the corresponding volume elements $\;\mathrm{d}{\boldsymbol{r}_{i}}$ and $\;\mathrm{d}{\boldsymbol{r}_{j}}$ belonging to the same connected cluster; $\rho$ is the number density of the system. Notice that we consider particles $i$ and $j$ connected if the distance between their assigned coordinates $|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|$ is smaller than a constant threshold $d$. In that sense, the coordinates $r_{i}$ and $r_{j}$ can be regarded as centers of spherical connectivity shells of diameter $d$ so that a connection corresponds to overlapping connectivity shells. This notion of connectivity is commonly referred to as Boolean model. Definition (1) implies $P({\boldsymbol{r}_{i}},{\boldsymbol{r}_{j}})\leq g({\boldsymbol{r}_{i}},{\boldsymbol{r}_{j}})$, where $g$ denotes the pair-distribution function. Furthermore, we define the connection probability $p$ that particles centered at ${\boldsymbol{r}_{i}}$ and ${\boldsymbol{r}_{j}}$ are part of the same cluster given their existence as

[TABLE]

It seems natural to assume that a complete description of the connectivity properties requires complete information on the thermodynamic equilibrium, i.e., density distributions to arbitrary order. Acquiring this critical knowledge is commonly subsumed as solving the thermal problem. In contrast to that, extracting the connectivity properties for given density distributions is referred to as the percolation problem.

Following reference Hansen and McDonald (1990), the pair-distribution function can be written as a diagrammatic density expansion using the Mayer $f$-bonds

[TABLE]

for an arbitrary pair potential $V$:

[TABLE]

where $\beta_{n}$ denotes the sum over the so-called irreducible cluster integrals of the second kind of order $n$; numbers $i$ stand for particle positions ${\boldsymbol{r}_{i}}$. Each diagram in $\beta_{n}$ contains two labeled white circles, $n$ black circles, and $f$-bonds between them such that there is no direct bond between the white circles, but if you drew an imaginary line between them, each diagram would be free of connecting circles. Removal of a connecting circle splits the diagram into two or more separate components. Schematic representations of $\beta_{1}$ and $\beta_{2}$ are depicted in Figure 1.

The pair connectedness can be treated with exactly the same expansion by modifying the $f$-bonds such that they distinguish between connected and disconnected pairs of particles. For Boolean models the Boltzmann factor $e(i,j)=f(i,j)+1$ can be split into a connected ($\dagger$) and a disconnected ($\ast$) part, respectively, using the Heaviside step function $\Theta$ Hill (1955):

[TABLE]

This translates into the corresponding Mayer bonds via

[TABLE]

Applying eqn. (6) to eqn. (4) yields an expansion of $g$ in terms of the connectivity bonds $f^{\dagger}$ and the blocking functions $f^{\ast}$. Finally, if the sum is restricted to diagrams that feature a connection between the white circles established purely via $f^{\dagger}$-bonds, we arrive at the density expansion of the pair connectedness function $P$:

[TABLE]

where $\beta_{n}^{\dagger}$ contains all diagrams of $\beta_{n}$ which after replacing each $f$ by either $f^{\dagger}$ or $f^{\ast}$ feature a path of $f^{\dagger}$-bonds connecting the two white circles. Naturally, $P(i,j)=g(i,j)$ follows for $|{\boldsymbol{r}_{i}}-{\boldsymbol{r}_{j}}|<d$ from eqn. (7). Eqn. (8) is the expansion that has to be reproduced by any alternative approach in order to be exact.

The approach to continuum percolation by Coniglio et al. Coniglio et al. (1977) starts off by dividing the diagrams in eqn. (8) into nodal and non-nodal parts. The latter constitute the direct connectivity $(C^{\dagger})$. Recognizing that the nodal-diagrams can be constructed as products of non-nodal diagrams results in the connectivity Ornstein-Zernike (OZ) relation

[TABLE]

This approach of casting the non-nodal diagrams bears the indisputable advantage, that due to its structural equivalence to the standard OZ relation, the entire toolbox from liquid state theory can be applied. However, akin to liquid state theory Hansen and McDonald (1990), the direct connectivity requires approximations or assumptions as it is still an infinite sum of arbitrarily complicated diagrams. Hence, eqn. (9) grants insight into the structure of $P$, but does not actually solve the problem unless the direct connectivity $C^{\dagger}$ function is known. Moreover, the percolation threshold is related to $C^{\dagger}(r\to\infty)$. Unfortunately, the closure relations employed in liquid state theory which allow for an analytic treatment tend to make assumptions specifically for this regime. The Percus-Yevick closure for instance assumes $C^{\dagger}(r>d)=0$ and can thus not be expected to be accurate in predicting the percolation threshold. We therefore use a slightly different approach.

On the diagrammatic level, the thermal and the percolation problem share the same complexity. Unfortunately, with the exception of a few simple systems there are no exact solutions to the thermal problem available to work with. However, there are many decent approximations as well as experimental data and simulation results, which can be used as a starting point. We are hence interested in a computational scheme to construct $P$ from accessible observables like the pair-distribution function or the nearest-neighbor distribution. Both functions inherently contain the diagrams $\beta_{n}$ in their corresponding diagrammatic expansions which we endeavor to exploit.

In our deliberations, Volterra equations of the second kind Tricomi (1985), i.e., equations of the form

[TABLE]

will play a key role. Here, $P$ is to be determined for given functions $I$ and $K$ (referred to as inhomogeneity and kernel, respectively) and a real number $\lambda$. Notice, that the lower boundary of the integral can always be chosen as zero by supplementing the kernel with a respective Heaviside function $\Theta(x-a)$. This way, eqn. (10) can always be cast in a form to which Laplace transform techniques can easily be applied. When $I$ and $K$ are specified, $P$ is obtained through simple numerics, or even analytically, if $I$ and $K$ are known analytically. If $I$ and $K$ are $L^{2}$-functions, the unique solution of eqn. (10) (except for functions that vanish almost everywhere) can be formally written as

[TABLE]

where $H$ denotes the resolvent kernel defined by

[TABLE]

with the iterated kernels $K_{n}$. The latter satisfy the recurrence relation

[TABLE]

This recurrence relation is essentially a formalized Picard iteration of which convergence is assured under the condition of $K,I\in L^{2}$ Tricomi (1985), which for our purposes will always be trivially satisfied. In diagrammatic terms, eqn. (12) corresponds to the sum over all chain diagrams with $\rho$-circles and $K$-bonds on the bounded interval $[0,r]$. The resolvent kernel satisfies the integral equation

[TABLE]

which depends exclusively on the integral kernel $K$. Volterra equations can thus be utilized to compute chain diagrams, which are essential to one-dimensional systems.

III Nearest-Neighbor Interactions

In this section, we develop our approach for nearest-neighbor interactions and present exact results for pair-connectedness functions. There are so far only two systems for which an exact analytical expression for the pair connectedness has been found, namely the one-dimensional ideal gas Domb (1947) and one-dimensional hard “spheres” (i.e. impenetrable line segments) Vericat et al. (1987); Drory (1997). The solutions to these two cases employ vastly different techniques, however, they share the convenient property that interactions (if present at all) are restricted to nearest neighbors.

Consider a system of only pairwise interacting identical particles in one dimension. In the absence of any external field, the one-particle density 111$\rho^{(1)}(r)$ is defined as the grand-canonical ensemble average over all configurations featuring a particle at $r$. $\rho^{(1)}(r)$ equals the number density $\rho$ of the interacting particles. We set coordinates such that one particle is fixed at the origin; without loss of generality we restrict our considerations to $r>0$. The pair-distribution function $g^{(2)}(0,r)=g(|r|)=g(r)$ is a measure of the average density of particles at a distance $r$ from the origin, given that there is a particle at the origin. In contrast to that, the pair connectedness $P$ describes the average density of particles that additionally belong to the same cluster as the particle at the origin. Following eqn. (2), $P$ can be factorized

[TABLE]

where $p$ is an actual probability, i.e., $p\in[0,1]$, the probability that two particles a distance $r$ apart belong to the same connected component. Naturally, within the connectivity shell of a particle $p(r<d)\equiv 1$, hence, $P(r)=g(r)$. Beyond $d$, the connection must be established through at least one mediating particle within the connectivity shell of the first particle. The average density $\omega(\tau)$ of such particles at $\tau\in(0,d)$ can be expressed by the pair-distribution function

[TABLE]

Yet, there might be more than one single particle in the connectivity range of the first particle at the origin. Indeed, if there were an additional particle at $\tau^{\prime}\in(0,\tau)$, this in general would impact $\omega$ and require integrating over all possible configurations of particles in $(0,\tau)$. In order to avoid the difficulties connected to this problem, we instead only look for the particle closest to the origin (in positive direction), ruling out the existence of obstructing correlations due to intermediate particles by definition. Thus, we are interested in the distribution of nearest neighbors. The probability distribution $\omega^{\prime}$ for finding such a nearest neighbor to the first particle at $\tau$, can be decomposed as

[TABLE]

where $\mathcal{P}$ denotes a conditional probability. Torquato et al. wrote down the reverse decomposition of $\omega^{\prime}$ Torquato et al. (1990), i.e., the product of the gap probability and the conditional existence of the particle at $\tau$. Clearly, these descriptions are equivalent through Bayes theorem. However, as we strive to devise a scheme that takes $g(r)$ as input, eqn. (18) bears the advantage that at least $\omega$ is already known. As shown in ref. Torquato et al. (1990), the nearest-neighbor distribution can, in general, not be inferred from $g(r)$ alone, because an exact treatment would require knowledge of the entire hierarchy of density distributions. Yet, in one-dimensional systems with only nearest-neighbor interactions, the decomposition

[TABLE]

is exact for $r_{1}<r_{2}<...<r_{n+1}$ Salsburg et al. (1953) and with that, the entire hierarchy of higher order distribution functions factorizes into products of $g^{(2)}$. Thus, as long as particles only interact with their nearest neighbors, $g(0,r)$ contains all properties of the equilibrium system, for instance the nearest-neighbor distribution $\omega^{\prime}(r)$, which can be constructed in the following way:

The probability of finding the nearest neighbor at a position $r$ is given by the difference between $g(0,r)$ (i.e. the probability of finding a particle at $r$ at all) and the probability of finding at least one particle in between [math] and $r$. It might now be tempting to compute the latter probability by integrating over $g^{(3)}(0,r_{1},r)$

[TABLE]

However, this expression overcounts configurations. (The reader can easily check this statement for the case of the ideal gas, where $g^{(2)}(r,r^{\prime})=1\longrightarrow I_{1}(r)=r\rho$, but $\omega^{\prime}(r)\neq 1-r\rho$, see eqn. (33).) Indeed, imagine there are exactly two particles in between 0 and $r$ placed at $r_{a}$ and $r_{b}$, respectively. The integral in eq. (20) counts this configuration twice - once if $r_{1}=r_{a}$ and again for $r_{1}=r_{b}$ although the configuration are indistinguishable. To account for this, we need to explicitly add the configuration with particles at $0,r_{a},r_{b}$ and $r$ with their corresponding weight given by the four point correlation function $g^{(4)}(0,r_{a},r_{b},r)$. Accounting for up to two mediating particles in general, we thus have to subtract

[TABLE]

However, this term now overcounts configurations with three and more particles in between [math] and $r$. Continued alternating addition and subtraction of terms constructed in this way finally yields the correct $g$-bond expansion of the nearest-neighbor distribution, see fig. 2.

The major advantage of this expansion is that it reveals the correspondence to the following integral equation

[TABLE]

which is a Volterra equation of the second kind, and therefore numerically dealt with. Thus, for systems of only nearest-neighbor interactions, the nearest-neighbor distribution can be computed straightforwardly from the pair-distribution function and vice versa as eqn. (22) is easily inverted. Perhaps, the inverted form

[TABLE]

is even more intuitive as any $g(r)$ is naturally the result of a sequence of nearest neighbors.

If we aim for two particles at [math] and $r$ to be connected, they are either already directly connected ($r<d$), or a nearest neighbor of the particle located at [math] exists within $(0,d)$ which is connected (through an arbitrary number of other particles) with the particle located at $r$. The corresponding integral equation for the pair connectedness reads

[TABLE]

where $\Theta$ denotes the Heaviside step function. It is important to notice that this equation works for nearest-neighbor interactions, because the particle at the origin has no influence on the particle distribution beyond its nearest neighbor. Indeed, the equation suggests that the nearest neighbor can be chosen as a new origin from which to connect to $r$ in a shifted system of coordinates. Equation (24) is trivial for $r<d$. In absence of symmetry breaking external fields, the equation can be recast in the standard Volterra type for $r>d$ via:

[TABLE]

where we introduced the kernel $K(x):=\Theta(d-x)\omega^{\prime}(x)$, accounting for the fact that only the nearest-neighbor distribution within the initial connectivity shell has an influence on the connectivity properties of the system. The inhomogeneity $I(r)$ is given by

[TABLE]

for $d<r<2d$, implementing the initial condition $P\equiv g$ for $r<d$ with $I(r\geq 2d)=0$. Moreover, as presented above, $I$ and $K$ can both be constructed from the pair-distribution function as long as eqn. (19) holds:

[TABLE]

Equations (25)-(28) suffice to recover the analytically known solutions for fully penetrable and impenetrable rods Domb (1947); Vericat et al. (1987); Drory (1997) (we will show this in Subsections III.1 and III.2). While the latter has been derived in its complete form in ref. Drory (1997) through a sophisticated mapping to a specific lattice model, the approach presented here is straight forward and, in particular, generally applicable to any kind of one-dimensional system with nearest-neighbor interaction.

It remains to show that eqn. (25) indeed reproduces the diagrammatic representation of eqn. (4). Volterra equations generate ordered chain diagrams, i.e. in our case diagrams of the structure displayed in fig. 3

with the additional condition that $0<r_{1}<r_{2}<...<r$. We can now replace the nearest-neighbor distribution by the corresponding $g$-bond expansion, presented in fig 4.

Note that the Heaviside bond renders the additional constraint on each individual bond obsolete so that we can use $g$-bonds instead of $g^{\dagger}$-bonds. In order to relate the $g$-bond expansion to the unordered $f$-bond expansion of Coniglio Coniglio et al. (1977), we can identify the diagrams of a specific order in $\rho$ for both approaches. In Coniglio’s expansion (see fig. 1 and eqn. (8)), diagrams are naturally ordered by powers of $\rho$, i.e. all diagrams featuring two black circles form the part of the solution that is quadratic in $\rho$. In contrast to that, the nearest-neighbor distribution as well as $g$ itself already contain contributions to arbitrary order in $\rho$. Instead, the number of dotted black circles in a diagram of our expansion defines the lower limit of nodal circles in the corresponding $f$-bond representation.

The expansions apparently coincide for $r<d$, yielding simply $g^{\dagger}$. We can henceforth ignore any diagram that contains an $f$-bond between the two labeled circles. To first order in $\rho$, there remains only one diagram in Coniglio’s expansion (see fig. 1 top panel). Considering

[math]r$$g^{\dagger}$$\omega^{\prime\dagger}

and taking the zeroth order in $\rho$ for both bonds, we reconstruct

[math]r$$f^{\dagger}$$f^{\dagger}

which corresponds to the first order in Coniglio’s expression except with a dotted circle instead of a proper black one. We can demonstrate that both circles are indeed equivalent by first showing more generally that the region of integration can always be reduced to the interval $[0,r]$ and further that contributions of unordered configurations cancel each other. To this end it is useful to write down $f^{\dagger}$ and $f^{\ast}$ explicitly for the specific conditions that nearest-neighbor interactions provide. Distinguishing between bonds between nearest neighbors (NN) and non-nearest neighbors as well as whether the connectivity shells of the associated particles overlap, we find

[TABLE]

[TABLE]

The integral $\beta^{\dagger}_{1}$, i.e.

[math]r$$f^{\dagger}$$f^{\dagger}

does not contribute to $P$ if $r>d$, as the preceding $e^{\ast}$-bond vanishes in that case. Thus, if the position $r_{1}$ of the intermediate particle in the diagram above is not located within $[0,r]$, either $|r-r_{1}|$ or $|r_{1}-0|$ is larger than $d$ thanks to the one-dimensional nature of the system. This, however, implies that one of the $f^{\dagger}$-bonds and hence the entire diagram vanishes. This argument can be applied to any diagram of the expansion but in a slightly adapted form. Consider the diagram

0rr_{1}$$r_{2}$$f^{\dagger}$$f^{\dagger}$$f^{\dagger}

for $r>d$. The integrand of the integral corresponding to this diagram does not necessarily vanish if one of the mediating particles lies outside of $[0,r]$, for example $r_{1}<0$ while $r_{2}\in[0,r]$. That means that the bond between $r_{1}$ and $r_{2}$ connects non-nearest neighbors, so that the associated $f^{\dagger}$-bond becomes unity. But then the expansion also contains the diagram

0rr_{1}$$r_{2}$$f^{\dagger}$$f^{\dagger}$$f^{\dagger}

where the dashed line represents an $f^{\ast}$-bond. Since $|r_{2}-r_{1}|<d$ (otherwise $f^{\dagger}(r_{1},r_{2})=0$ anyway) also $|r_{2}-0|<d$ so that $f^{\ast}(0,r_{2})=-1$. Therefore, the two integrals considered differ only by sign and thus annihilate each other. For any configuration in which the connection of $f^{\dagger}$-bonds between the white circles features particles not within $[0,r]$, we can repeat this procedure. We link two particles, which the “outlying particle” shares $f^{\dagger}$-bonds with, by an $f^{\ast}$-bond and drop that configuration from the expansion.

However, there is one exception if the particles one would like to link by an $f^{\ast}$-bond are already linked by an $f^{\dagger}$-bond. In this case, the path through the outlying particle is obsolete. One can replace the $f^{\dagger}$-bonds it is connected to by $f^{\ast}$-bonds to obtain different diagrams of the same expansion. One might therefore replace the obsolete bonds immediately by the corresponding $f$-bonds of which there is at least one that connects non-nearest neighbors and hence vanishes. This way it becomes apparent on the diagrammatic level, that the configuration of particles outside of the interval we are interested to bridge, does not influence the connectivity properties within that interval as long as we deal with nearest-neighbor interactions.

Using the same line of reasoning one can show that all configurations that contain an $f^{\dagger}$-bond between non-nearest neighbors will be canceled. Yet, this argument does not directly restrict the $f^{\ast}$-bonds. It should be noted that there cannot be $f^{\ast}$-bonds between nearest neighbors, because a continuous path of $f^{\dagger}$-bonds between the white circles would then require at least one $f^{\dagger}$-bond between non-nearest neighbors, which we ruled out. Moreover, $f^{\ast}(r_{1},r_{2})$ does also vanish for non-nearest neighbors if $|r_{1}-r_{2}|>d$. Therefore, all appearing bonds are in fact short ranged.

In summary, for nearest-neighbor interactions in one dimension, the expansion by Coniglio, eqn. (8), contains all diagrams with $f^{\dagger}$-bonds only between nearest neighbors, $f^{\ast}$-bonds only between non-nearest neighbors, and an $e^{\ast}$-bond between the white circles which are free of connecting circles.

Now we take a second look at eqn. (30) and notice

[TABLE]

for non-nearest neighbor interaction. That at hand we can rewrite our $g^{\dagger}$-bond expansion of fig. 4 by replacing the $\Theta$-bonds by $-f^{\ast}$ bonds

0r$\omega^{\dagger}$ $=$ 0r$+$ 0r$f^{\ast}$ $+$ 0r$f^{\ast}$ $+$ $...\quad,$

thereby eliminating the alternating sign. Then we can insert the expansion of $g^{\dagger}$ and exploit the same arguments as before to find that both expansions can indeed be brought in perfect unison.

However, much more straightforwardly, we can simply put the equation up to the practical test by applying it to problems for which the exact solution is known or at least easily obtained through simulations.

III.1 Fully Penetrable Rods

A one-dimensional ideal gas, i.e., non-interacting fully penetrable connectivity shells, can be solved purely by stochastic tools. However, since the presented framework is straightforward to apply, the ideal gas makes up for a nice test-case system. The integral kernel follows immediately from eqn. (27):

[TABLE]

recovering the well-known exponential distribution of ‘gap lengths’ Zernike and Prins (1927). With $\omega^{\prime}(r)$ known, the inhomogeneity $I$ is readily obtained using eqn. (28):

[TABLE]

The integral equation for the pair connectedness function of the one-dimensional ideal gas thus reads

[TABLE]

The analytical solution via the resolvent kernel $H$ of eqn. (12) is intricate as the ensuing Heaviside integrals are not straightforward to compute. Yet, the Heaviside functions can be eliminated by restricting the integral equation to intervals $[nd,(n+1)d]$ and, progressively, solving it for $n\in\mathbb{N}$. For $n=1$, eqn. (35) is simplified to

[TABLE]

which is solved by the linear function

[TABLE]

In general, the solution can be found by assuming a polynomial of $n$-th degree and comparing the coefficients of all but the leading order in $r$ as well as the coefficient of the $e^{-\rho r}$ term. The equation for the coefficient of leading order is always trivially satisfied. Thus there are $n+1$ equations for $n+1$ coefficients, granting the unique solution. The procedure can be generalized to yield the solution for all $n$ in form of the following recurrence relation:

[TABLE]

Once cast in a closed form, eqn. (37) recovers the known solution obtained before by Domb and others Domb (1947); Torquato (2013). The solution is shown in fig. 5.

III.2 Impenetrable Rods

One of the few non-trivial systems for which the pair-distribution function has been found exactly in an analytic closed form, is the system of one-dimensional impenetrable identical hard rods Vericat et al. (1987); Drory (1997). For instance from classical density functional theory it is known that the corresponding pair-distribution function reads Percus (1976)

[TABLE]

where $\sigma$ denotes the length of the rods. The nearest-neighbor distribution can be computed from eqn. (22). However, the distance distribution to a nearest neighbor is also equivalent to the gap length distribution. This, in turn, can be understood as randomly (i.e. in a uniformly distributed manner) placing $N$ points on a line of a length that corresponds to the free volume. The corresponding probability distribution has already been formulated by Zernike Zernike and Prins (1927) in the form

[TABLE]

which is simply the zeroth order term in eqn. (38). Hence the integral kernel is yet again an exponential, which implies that the resolvent kernel is similar in structure as well. Indeed, the same procedure that worked for the ideal gas also works for impenetrable spheres. As a result, the solution previously reported by Drory Drory (1997) (and partially before also in ref. Vericat et al. (1987)) as

[TABLE]

is found to be the unique solution to eqn. (24), with $g$ and $\omega$ defined by eqn. (38) and eqn. (39), respectively. The agreement of the theory with simulations is shown in fig. 6.

Note, that we used the known solution for the thermal problem for simplicity. This is not required here, as the thermal problem can also be mapped onto a Volterra equation. For impenetrable rods $g$ can be obtained by solving

[TABLE]

Note further, that on the diagrammatic level, the hard core repulsion is not a nearest-neighbor interaction but rather short ranged, as the $e$-bond between non-nearest neighbors is not necessarily unity. Integral equation (25) remains perfectly valid, but the line of reasoning in the comparison to the more general expansion eqn. (8) has to be slightly modified.

III.3 Numerical Examples

So far we were able to reproduce known results straightforwardly because the associated thermal distribution functions were available in a closed analytical form. However, one major virtue of the proposed scheme lies in the fact that it does not require analytic input to work. We can simply sample the pair-distribution function in the first connectivity shell for an arbitrary nearest-neighbor interaction and compute the pair connectedness. Thus we can, for instance, perform a Monte Carlo simulation to extract $g(r)$, numerically solve the ensuing Volterra equation on an equidistant grid to obtain predictions for the pair connectedness and compare them to the simulations. To demonstrate this, we consider the purely repulsive pair potential

[TABLE]

Solving eqn. (22) for this interaction results in the nearest-neighbor distribution depicted in figure 7. As the potential acts only on nearest neighbors, the next-nearest-neighbor distribution $\omega^{(2)}$ is simply the convolution of $\omega^{(1)}:=\omega^{\prime}$ with itself

[TABLE]

which is also shown in fig. 7. Once the hierarchy of nearest-neighbor distributions and therefore the kernel of our integral equation is known, the problem becomes trivial. We solve eqn. (25) numerically yielding the solid line in fig. 8 which as expected is in perfect agreement with the pair connectedness determined by simulations. Notice, that the process is even invertible, i.e. from the pair connectedness the kernel can be reconstructed, yielding the nearest-neighbor distribution which will give you the radial distribution function. That means, for one-dimensional nearest-neighbor interacting systems the pair-distribution functions contains the same information as the pair connectedness.

For short-ranged but not necessarily nearest-neighbor interactions we cannot expect eqn. (25) to be exact. However, it serves as a good approximation if the interaction energy resulting from beyond nearest neighbors is small, i.e. if the potential decays sufficiently fast. As an example, for particles interacting through the Lennard-Jones potential

[TABLE]

theory and simulation results for the pair connectedness cannot be distinguished by eye (see fig. 9).

If you turn your attention to the nearest-neighbor distribution in fig. 10) at large $r$, a slight discrepancy is visible between the simulation data and the theory. Yet, the deviation appears in a regime where the nearest-neighbor distribution is already small whereas the main peak is properly depicted. Moreover, the inhomogeneity

[TABLE]

as the above convolution, weighs the nearest-neighbor distribution around $d$ with the short range $g$ which for any close to hard-core-interaction should be extremely small. The inaccuracy therefore hardly propagates.

While this level of agreement seems to be a lucky coincidence, the following section illustrates how long-ranged interactions can be treated in a more systematic way. (It shall be noted though that the solution of Volterra equations is typically stable against noise in the input.)

IV Generalizations

Arguably, one-dimensional classical systems with nearest-neighbor interactions do not occur in real life on a regular basis. We therefore strive to generalize the depicted procedure to more realistic conditions and discuss the impact of the added complexity.

IV.1 External Fields

If we stay in one dimension for a start, an external field $\phi$ destroys the homogeneity of the system such that the single-particle density $\rho^{(1)}$ is not a constant throughout the system anymore. With that, all distribution functions become explicitly dependent on the positions they are evaluated for. However, eqn. (24) formally already accounted for a potential dependence on two points in space. Thus, the only modification required is to replace the number density by the space-dependent single-particle density:

[TABLE]

assuming $r_{1}<r_{2}$ for simplicity. Supplemented by the initial condition

[TABLE]

this equation is also exact for nearest neighbor-interactions, however, it requires the density profile $\rho^{(1)}(\tau)$ as additional input. Moreover, $g,\omega^{\prime}$ and most importantly $P$ depend on two points in space, such that the equation becomes technically a two-dimensional Volterra equation which is numerically more challenging than its one-dimensional counterpart. The important observation is that on tagging a particle the system still entirely splits, in that there is no coupling between left and right hand side of the particle, respectively. The external field adds a local weight to the integral kernel but that is all there is to it.

Cast in the standard form, the equations that need to be solved read

[TABLE]

Thus, external fields do not add any complexity to the connectivity problem because the diagrammatical structure remains chain-like and can thus be expressed as a Volterra equation. This unfortunately does not apply for the subject of the next section.

IV.2 Long-ranged Interactions

In contrast to external fields, for long-range interactions, we cannot ignore three-particle correlations anymore. Thus, there is no straightforward way to determine the nearest-neighbor distribution from just the radial distribution function. But we can at least attempt to characterize the discrepancy. In one-dimension we can order the particles $r_{1}<r_{2}<...<r_{n}$, so that we know beforehand which particle is neighboring another particle. The system geometry now demands that if particles at $r_{1}$ and $r_{3}$ are connected, the same has to apply for $r_{1}$ and $r_{2}$ as well as $r_{2}$ and $r_{3}$. Therefore all diagrams contributing to $P$ share the same carcass, i.e.

0f^{\dagger}$$f^{\dagger}$$f^{\dagger}$$r-****…$f^{\dagger}$.

This diagram is naturally part of eqn. (8), but every addition of an $f^{\ast}$-bond will result in another diagram of that expansion. Most importantly, all diagrams contributing to $P$ can be constructed in that way. Recall, that our scheme for nearest-neighbor interactions generates chains of the type depicted in fig. 3.

0g^{\dagger}$$r_{1}$$\omega^{\prime\dagger}$$r_{2}$$\omega^{\prime\dagger}$$r-****…$\omega^{\prime\dagger}$.

At this point we hit two obstacles: On the one hand, for long-ranged interactions $\omega^{\prime}$, as obtained from eqn. (22), is not exact, as the pair-distribution hierarchy does not factorize anymore. However, for highly repulsive potentials, the configurations that feature more than one particle within the connectivity shell (to the right) are strongly suppressed energetically. Thus, the additional long-ranged interaction energy hardly alters the short-scale alignment. Figure 11 illustrates this observation for the inverse square potential and varying interaction ranges.

Moreover, it shows that treating the long-ranged interaction pair distribution as a nearest-neighbor interacting system might render $\omega^{\prime}$ negative - hence $\omega^{\prime}$ cannot be safely interpreted as a nearest-neighbor distribution for long range interactions. However, since we need input for the thermal distribution anyway, we might as well use the real, i.e. numerically sampled, nearest-neighbor distribution as input.

On the other hand, the hierarchy of neighbor-distribution functions cannot be generated by iteratively convolving $\omega^{\prime}$. Once the nearest-neighbor is found, we cannot shift the system to its position and expect the same distribution of nearest neighbors to apply for that particle due to the correlation to the previous origin. Essentially, the $P(\tau,r)$ on the right hand side of eqn. (24) is not the regular pair connectedness anymore but rather the pair connectedness under the constraint that there is a particle already placed at the origin. At this point it is useful to formulate eqn. (24) in terms of probabilities to be able to invoke Bayes’ theorem

[TABLE]

Since $p$ is a proper probability, we can treat the constrained probability in the integral according to Bayes theorem

[TABLE]

Thus, we can formally write down a Volterra equation for the probability of particles at 0 and $r$ belonging to the same cluster

[TABLE]

Accordingly, the conditional probability can be absorbed into the kernel. This is hardly surprising as we expect an equation structurally similar to eqn. (9). Notice also that we need the complete pair-distribution function as input where previously only the $g(r)$ within the first connectivity was required. It turns out, that if we assume statistical independence of events in eqn. (49) we can already reproduce the pair connectedness for long-range interacting systems to surprisingly high precision (see fig. 12).

Equation (50) serves as an excellent approximation even with the crudest assumption for the constrained connectivity probability. In this approximation, we essentially compute $P$ as before for nearest-neighbor interactions but also normalize it by the $g$ that is caused by $\omega^{\prime}$ if we assume only nearest-neighbor interactions, i.e. eqn. (23), to get the probability $p$. The substitute $g$ is also shown in figure 12. By multiplying the probability with the sampled $g$, we get a rescaled pair connectedness, which is in excellent agreement with the simulations results.

IV.3 Higher Dimensions

Finally we sketch the application of the proposed scheme to a three dimensional system. The fundamental difference to the one-dimensional case is that there is more than one path that can lead to a connected configuration. To fight at one front at a time, we will only consider the three-dimensional ideal gas, so that we do not need to worry about correlations of particle positions. As the factorization of the hierarchy of density distributions (eqn. (19)) holds, the nearest-neighbor distribution can still be determined with an analogue of eqn. (22). The notion of a nearest-neighbor remains valid and the corresponding distribution function depends exclusively on the distance to a chosen particle. However, since the system does not allow for global order, the idea of acquiring the next-nearest-neighbor distribution through a three-dimensional convolution of nearest-neighbor distributions does not work anymore. We can however switch to eqn. (23) as the nearest-neighbor distribution in that expression appears only with respect to the origin. As a consequence, we find

[TABLE]

which is readily solved to yield

[TABLE]

This expression is still defined on $\mathbb{R}^{3}$ and normalized accordingly, thus the result coincides with an the expression derived, for instance, by Torquato Torquato et al. (1990) once the angular dependencies are integrated out.

The three-dimensional analogue to eqn. (48) for the three-dimensional ideal gas reads

[TABLE]

In order to make further progress, we need to find a way to treat the conditional probability. First we make use of Bayes theorem again to convert the equation into a standard type integral equation.

[TABLE]

Capitalizing on the homogeneity of the system, we choose ${\boldsymbol{r}}=(0,0,r)^{t}$ and parameterize the position of the nearest neighbor by ${\boldsymbol{\tau}}=\tau(\cos(\phi)\sin(\vartheta),\sin(\phi)\sin(\vartheta),\cos(\vartheta))^{t}$. On top of that, the unconditional pair connectedness depends exclusively on the distance between the two arguments so that we obtain

[TABLE]

Introducing $\eta:=\frac{\pi}{6}\rho d^{3}$, the volume fraction of connectivity shells, and substituting $t:=\frac{\tau}{d}$ results in

[TABLE]

Finally, we replace the angular integral by an integration over $u:=|\boldsymbol{\tau}-\boldsymbol{r}|$ which importantly depends monotonically on $\theta$ retrieving the familiar form

[TABLE]

This is the moment we have to leave the realms of exactness to make further progress as we do not know the analytic structure of the correlation

[TABLE]

In order to find sensible approximations we need to understand the function appearing in the enumerator of (58). It is the probability that there is a particle at the origin with a spherical cavity of radius $\tau$ around it devoid of particles given particles at $\boldsymbol{\tau}$ and $\boldsymbol{r}$ which belong to the same cluster. Normalized by the probability of just finding a particle at the origin with no other particles in the open ball $B_{\tau}(\boldsymbol{0})$ we end up with $c(r,t,u)$. However, thanks to the ideal gas, the particle positions are entirely uncorrelated so that the only meaningful information in the condition is the fact that $\boldsymbol{\tau}$ and $\boldsymbol{r}$ belong to the same cluster. Now, if $|\boldsymbol{\tau}-\boldsymbol{r}|<d$, even this piece of information is obsolete, meaning

[TABLE]

The only configurations making a difference are those where $\boldsymbol{\tau}$ and $\boldsymbol{r}$ are necessarily connected through the particle at the origin. In other words, $c$ differs from 1 due to all configurations for which the nearest neighbor turns out to be a deadlock on the path to $\boldsymbol{r}$. However, $c$ is not always larger than 1 due to the obstructing void that comes with the nearest-neighbor condition. On the contrary, as $t$ approaches unity we expect $c$ to be smaller than 1 because in the limit of $t\rightarrow 1$, the probability of the origin having another neighbor vanishes. As a consequence, the origin only blocks the intersect with the particle at $\boldsymbol{\tau}$ for other other potential connectors. Thus,

[TABLE]

On the other hand, if $t$ approaches zero, i.e. the nearest neighbor becomes the origin, the origin again is obsolete and hence

[TABLE]

At this point, we need to close the integral equation by specifying a functional form of $c$. In contrast to the closures typically used in liquid state theory, equation (57) allows for a purely geometrical treatment.

We assume, the probability of a structure to be part of the same cluster as a distant particle is proportional to the volume it provides for other particles to attach, i.e.

[TABLE]

for $u\geq d$. In order for eq. (61) to apply, the proportionality constant even has to be unity. Thus, the full function reads:

[TABLE]

Apparently, condition (59) is satisfied by this choice for $c$ as well. Naturally, the above expression is an approximation as the actual $c$ has to depend on $r$ and $u$ in a non-trivial way. Nevertheless, plugging (63) into equation (57) already yields good agreement with simulation results.

The above kernel has the primary advantage that due to the simple $u$-dependence, the integrals over $u$ and $t$ can be exchanged in order to yield a $t$-integral which can be performed analytically. Thus, eq. (57) can be reduced to a simpler form containing only a single integral which can be considered as a Fredholm integral equation as long as $r$ decays fast enough that $p(r>\xi)=0$ for a finite $\xi$ is a sensible approximation. As typical for Fredholm equations we can make use of the Picard iteration to solve (57) with the specified kernel, the results are illustrated in figures 14 and 15. Keeping in mind the crudeness of the employed approximation, the agreement to our simulation results is remarkably good especially for small volume fractions. In the limit of infinite dilution, the probability of a particle having two neighbors is already heavily suppressed, so that the nearest-neighbor being a deadlock in a connected configuration is effectively impossible. As a consequence $c(r,t,u)=1$ becomes exact for $\eta\rightarrow 0$. For finite volume fractions, the closure asymptotically generates exponential decays which systematically slightly undershoot the simulation results. However, for volume fractions $\eta\gtrsim 0.3$ numerical stability breaks down which is to be expected from Fredholm equations if the integral norm of the kernel becomes too large Tricomi (1985). This behavior is independent of the percolation transition as even exactly at the percolation threshold, the designated power-law solution would still belong to the class $L^{2}$. Additionally, equation (57) can easily be closed to yield solutions decaying like fractals which are perfectly stable numerically. The instability is hence simply an artifact of the closure. There are multiple ways to improve on our choice of $c$ for instance by implementing the geometrical nuances in a less crude fashion. Moreover, since the Picard iteration is computationally not expensive, even brute force methods come to mind. Ultimately, the constrained probability can also be determined by simulations which is as expensive as measuring $p(r)$ right away, but might still be useful as a starting point for more elaborate approximate schemes. However, as this section was supposed to just illustrate, that the range of validity of our derived integral equation exceeds the one-dimensional, further extensions shall be discussed elsewhere.

V Conclusion

We have presented a general method to solve the connectivity problem for one-dimensional systems with an arbitrary nearest-neighbor interaction exactly, given the correct pair-distribution function. For these systems, we showed that the derived integral equation is equivalent to the connectivity Ornstein-Zernike equation, however, it substantially simplifies the derivation of the known exact solutions to continuum percolation models. Moreover, the connectivity properties can be inferred completely from thermal distribution functions. This relation would be of immense practical value if generalized to higher dimensions. For higher-dimensional systems and long-ranged interactions the analogous integral equation still holds, however, it features a constrained probability. In this case, similar to the connectivity Ornstein-Zernike equation, a closure relation is required. Hence one might argue that compared to the standard approach little is gained as approximations are required after all only for a different function. Yet, in contrast to the direct connectivity, the conditional probability appearing in the kernel is an observable of the system. The integral kernel can simply be sampled by Monte Carlo simulation, and especially in view of universality it might reveal useful insight for future approaches. And, as shown for the ideal gas, even approximations based on simple geometrical considerations already lead to decent results.

Acknowledgements.

We acknowledge funding by the German Research Foundation in the project SCHI 853/4-1.