Precise algorithms to compute surface correlation functions of two-phase heterogeneous media and their applications
The quantitative characterization of the microstructure of random heterogeneous media in -dimensional Euclidean space via a variety of -point correlation functions is of great importance, since the respective infinite set determines the effective physical properties of the media. In particular, surface-surface and surface-void correlation functions (obtainable from radiation scattering experiments) contain crucial interfacial information that enables one to estimate transport properties of the media (e.g., the mean survival time and fluid permeability) and complements the information content of the conventional two-point correlation function. However, the current technical difficulty involved in sampling surface correlation functions has been a stumbling block in their widespread use. We first present a concise derivation of the small- behaviors of these functions, which are linked to the mean curvature of the system. Then we demonstrate that one can reduce the computational complexity of the problem, without sacrificing accuracy, by extracting the necessary interfacial information from a cut of the -dimensional statistically homogeneous and isotropic system with an infinitely long line. Accordingly, we devise algorithms based on this idea and test them for two-phase media in continuous and discrete spaces. Specifically for the exact benchmark model of overlapping spheres, we find excellent agreement between numerical and exact results. We compute surface correlation functions and corresponding local surface-area variances for a variety of other model microstructures, including hard spheres in equilibrium, decorated “stealthy” patterns, as well as snapshots of evolving pattern formation processes (e.g., spinodal decomposition). It is demonstrated that the precise determination of surface correlation functions provides a powerful means to characterize a wide class of complex multiphase microstructures.
Random heterogeneous media are ubiquitous and arise in many applications
in physics, materials science, biology, and geophysics. Examples of such media include composites Torquato (2013), porous materials Zohdi and Wriggers (2008); Sahimi (2003); Vasseur and Wadsworth (2017), biological tissues Gibson and Ashby (1999); Gevertz and Torquato (2008), and even cosmological structures Peebles (1993). The quantitative characterization of the structure via higher-order correlation functions of these complex media is of importance in many
fields Torquato (2013); Chiu et al. (2013). In general, an infinite set of correlation functions are required to exactly determine the effective physical properties of the media Torquato (2013); Zohdi and Wriggers (2008); Sahimi (2003). However, such complete structural information about the medium is generally not available and hence one must settle for reduced information in the form of lower-order correlation functions. The study of these descriptors has proved fruitful and new applications involving these descriptors are constantly coming up, including in reconstructions using state-of-the-art techniques such as neural networks Mosser et al. (2017); Lubbers et al. (2017).
There are a variety of two-point structural descriptors, including the two-point correlation function Jiao et al. (2007, 2008), two-point cluster function Jiao et al. (2009), surface-surface correlation function Torquato (1986), and surface-void correlation function Torquato (1986), as well as the pore-size density function Torquato (2013), that are practically accessible via computer simulations or imaging techniques. Among them, the most well-known descriptor is the standard two-point correlation function , which can be obtained from scattering experiments Debye and Bueche (1949); Debye et al. (1957). This quantity has been employed to characterize the microstructure and physical properties of heterogeneous materials Torquato (2013), reconstruct the microstructure of heterogeneous materials Jiao et al. (2007, 2008), and recently, to quantify the hyperuniformity of two-phase systems Torquato (2016a, b); Ma and Torquato (2017). Although knowledge of the two-point correlation function has proved to be extremely useful, the corresponding correlation functions that characterize the interface of two-phase media such as the specific surface , surface-surface correlation function , and surface-void correlation function , which contain crucial structural information, have received considerably less attention, especially and . This is due partly to the fact that these surface correlation functions are not as easy to sample as , which we remedy in this paper, as described in Sec. IV.
While the two-point correlation function contains important structural
information, it is usually insufficient to determine both the structure and physical properties of heterogeneous media Yeong and Torquato (1998); Jiao et al. (2009); Torquato (2013). It has been shown that supplementing with surface correlation functions can lead to improved reconstructions
of two-phase media Torquato (2013); Zachary and Torquato (2011); Guo et al. (2014).
These surface correlation functions determine rigorous upper bounds on the fluid permeability of porous media and mean survival time associated with diffusion-controlled reactions among traps. These two-point “interfacial-surface” bounds have been shown to be much sharper than the so-called two-point “void” bound involving alone, reflecting the importance of surface correlation functions. For isotropic media in three-dimensional Euclidean space, these involve the following two key integrals:
where is just another way to write for the void phase if we focus on porous media, i.e., , and is the volume fraction of the void phase. A schematic plot of the correlation functions is shown in Fig. 1(a). For statistically isotropic media, the two-point “interfacial-surface” upper bound for the fluid permeability is given by Torquato (2013)
while the two-point “void” bound is
where is the volume fraction of the solid phase. Similarly, the analogous bounds on the mean survival time are given by Torquato (2013)
where is the diffusion coefficient of the reactant. The fact that the key integrals in these bounds on the fluid permeability are the same as those for the mean survival time is more than a
coincidence. Indeed, is rigorously bounded from above in terms of for general media Torquato (1990).
Moreover, two-point correlation functions determine local volume-fraction and local surface-area fluctuations as measured by the relevant variances. These variances enable one to generalize the concept of hyperuniformity Torquato (2016a), which was originally conceived in the context of point configurations, namely, it refers to the anomalous suppression of density fluctuations on large length scales Torquato and Stillinger (2003); Zachary and Torquato (2009). Notably, it is proven that sphere packings will inherit the hyperuniformity of the underlying point pattern Torquato (2016a). The local volume-fraction variance within a spherical observation window of radius in -dimensional Euclidean space is given by Lu and Torquato (1990)
is the autocovariance function associated with , and is the volume of a -dimensional sphere of radius , and is the scaled intersection volume, the ratio of the intersection volume of two spherical windows of radius whose centers are separated by a distance to the volume of a spherical window. A two-phase system is hyperuniform with respect to volume-fraction variances if decreases more rapidly than for large Torquato (2016a), or equivalently
where is the Fourier transform of . Similarly, the local surface-area variance has been defined by Lu and Torquato (1990)
is the autocovariance function associated with . A two-phase system is hyperuniform with respect to surface-area variances if decreases more rapidly than for large Torquato (2016a), or equivalently
where is the Fourier transform of . It has been suggested that surface-area fluctuations are more sensitive microstructural measures for heterogeneous media than corresponding volume-fraction fluctuations in some cases Torquato (2016b, a). Results obtained in this paper further support this conclusion.
Similar to the two-point correlation function , surface correlation functions can be related to and obtained from the scattering intensity as well Strey (1994); Dietrich and Haase (1995). In the most general case that involves scattering from both the bulk and the surface [see Fig. 1(b) for a schematic plot], the scattering intensity can be written as
where are certain coefficients. When the scattering from the surface is comparable to the bulk, one must consider all these three terms to determine the scattering intensity, while if only bulk or surface scattering is dominant, then one should only care about the corresponding correlation function, this interpretation can potentially provide a general way to understand hyperuniformity in two-phase media.
The rest of the paper is organized as follows: in Sec. II, we provide necessary definitions and background. In Sec. III, we present a concise and simple derivation of the small- behavior of the two-point surface correlation function, which involves the mean curvature of the entire system. In Sec. IV, we introduce and describe a general algorithm that enables the efficient computation of and . We verify the accuracy of our algorithm by applying it to overlapping spheres for which we have exact results Torquato (2013). In Sec. V, we show how to apply the algorithm to treat digitized two-phase media, which is of practical importance. Using Gaussian random fields as an example, we will demonstrate that the image resolution and some drop-out in sampling are crucial in order to obtain reliable results. In Sec. VI we explicitly show results of overlapping spheres, hard spheres in equilibrium and decorated stealthy point patterns. In Sec. VII we explicitly show results of patterns from spinodal decomposition and patterns from the Swift-Hohenberg equation. Using these examples, we demonstrate how surface correlation functions will be very useful for microstructural characterization and can be superior to in certain cases. Finally, in Sec. VIII, we make concluding remarks and discuss the implications of our findings.
Ii Background and Definitions
A two-phase random medium is a domain of space that is partitioned into two disjoint regions that make up : a phase 1 region of volume fraction and a phase 2 region of volume fraction Torquato (2013). The phase indicator function for a given realization is defined as
For statistically homogeneous media, the volume fraction for phase
is a constant. The two-point correlation function is defined as
For homogeneous media, this quantity only depends on the relative displacement vector . The two-point correlation function simplifies as . If the system is also statistically isotropic, then depends only on the radial distance .
The interface indicator function is defined as Torquato (2013)
The specific surface is the expected area of the interface per unit volume, and for homogeneous media is simply the ensemble average of the interface indicator function, i.e.,
The surface-surface correlation function measures the correlation between two points on the interface, and for homogeneous media is defined as
The surface-void correlation function measures the correlation between one point on the interface and the other in the void phase, and for homogeneous media is defined as
Higher-order surface correlation functions are similarly defined Torquato (2013), but the focus in this paper will be the two-point varieties.
Closed-form expressions for the two-point surface correlation functions are very limited. The most notable one is for the model of overlapping spheres Torquato (2013); Chiu et al. (2013), which is generated by circumscribing spheres of radius around each point in a Poisson point process with density . The space interior to the spheres is the solid phase and the space exterior is the void phase Torquato (2013). For statistically homogeneous overlapping spheres in three dimensions, we have
where is a radial distance, is a reduced density and is the Heaviside step function. Here,
is the two-point correlation function for the “void” phase. These relations were first given by Doi Doi (1976).
Using the canonical function Torquato (2013), Torquato derived the following expressions for surface correlation functions for hard spheres:
where is the specific surface, is the radial Dirac delta function, and is the sphere indicator function. The quantity is the total correlation function defined as , where is the pair correlation function, and denotes the convolution of two functions.
The impenetrability constraint alone is not sufficient to specify the hard-sphere model; a hard-sphere system can be in equilibrium or be derived from an infinite number of nonequilibrium ensembles Torquato (2013). The pair correlation function is generally not known for nontrivial hard-sphere models for all densities, an exception being the “ghost” random sequential addition packing model Torquato and Stillinger (2006). For , Torquato used the Percus-Yevick approximation and the Verlet-Weis correction to evaluate these functions for statistically isotropic systems of hard spheres in equilibrium Torquato (1986). These are useful benchmark results that will be used in Sec VI. To date, numerical evaluations of the surface correlation functions have been limited to hard spheres in equilibrium Torquato (1986); Seaton and Glandt (1986) and maximally random jammed sphere packings Klatt and Torquato (2016).
Iii Some Theoretical Remarks on Surface Correlation Functions
iii.1 The small- behavior of and in general
Debye and co-workers Debye and Bueche (1949); Debye et al. (1957) showed that the slope of at the origin () is directly proportional to the specific surface , which enables people to obtain the surface area of the whole system by measuring the tail of a scattering profile. The small- behavior of the two-point surface correlation functions have been derived previously by taking higher-order derivatives of of a dilated interface and then letting the thickness go to zero Teubner (1990); Ciccariello et al. (1981). Here we present a much simpler derivation based on a probabilistic interpretation of the surface correlation functions.
We restrict ourselves to the discussion of systems with interfaces that are differentiable everywhere. This assumption enables us to approximate the vicinity of a point on the interface with planes or spheres in the following discussion.
The small- behavior of the surface-surface correlation function is straightforward to obtain. First, randomly pick a reference point on the interface (with specific surface for the entire system). Second, consider a concentric shell with radius to around the reference point, then a local specific surface of the shell can be defined as . The quantity is then the product of the specific surface of the system and the average local specific surface over the interfaces (the local specific surface is defined at every point on interfaces, thus can be integrated to compute the average). We present a schematic plot that elucidates the derivation in Fig. 2(a) in three dimensions. When is very small, the vicinity of is basically flat for the zeroth-order approximation [see the quadrangle in Fig. 2(a)] and there is no other interface intersecting with the shell. As shown in Fig. 2(a), the area of interface contained in the shell is , and the volume of the shell is , so
For , following the same method we have
Since the zeroth-order approximation of is divergent as , we will not discuss higher-order finite correction terms here Teubner (1990).
The determination of the small- behavior of the surface-void correlation function is more involved. Again, we randomly pick a reference point on the interface. Next consider a “test” sphere of radius centered at the reference point. We denote by the following conditional probability: given a point on the interface, the probability that a uniformly and randomly placed vector emanating from lands in the void phase. The quantity is then the product of the specific surface of the system and the average of over the interfaces, i.e., . Since it is more convenient to illustrate the basic idea behind the computation in two dimensions and the result can be easily generalized to three dimensions, a schematic plot that elucidates our approach is illustrated in Fig. 2(b) in two dimensions, where the shaded area is the part of the small “test” sphere contained in the solid phase. In , under aforementioned assumptions, we can approximate the interface with a circular arc of radius of curvature . Then we can work out the probability up to the first-order approximation with respect to , which writes as
Average out this quantity on interfaces, we get the final result
where is the average of on interfaces.
Following the same procedure in three dimensions, we find
Thus we have
However, in three dimensions, the curvature varies when the normal plane rotates, and hence here is to be interpreted to be the mean curvature at the point. One should also notice that is a signed quantity in general, although we only illustrate the positive situation in Fig. 2(b) for the sake of simplicity. Note our simple approach can be easily extended to derive the small- behavior in higher dimensions. In any dimension, we find
where is the beta function. Using this approach, the connection between the small- behavior of and
mean-curvature interfacial growth problems Von Neumann () is intuitively clear.
Remarks: Note that when the “test” sphere of intersects with the nondifferentiable singularities, such as edges or corners, the derivation above breaks down. Thus, Eqs. (29) and (31) do not hold in general for interfaces that have singularities, even though the integrated mean curvature may still be defined and computed in these systems Mecke (1998). Using the same approach, we obtain in Appendix A some results for certain systems in which the interfaces have singularities. A discussion of these issues can be found in Ref. Ciccariello and Benedetti (1982).
iii.2 Phase-interchange relations for
Here we remark on phase-interchange relations involving the surface-void correlation function . For a two-phase medium, since the sum of indicator functions for phase 1 and phase 2 is unity everywhere, we have
implying that the sum of the two surface-void correlation functions for phases 1 and 2 equals the specific surface, i.e.,
Furthermore, if the two phases are statistically the same, these surface-void correlation functions are constants Teubner (1990), namely,
which is a remarkable relation given that it applies to complex microstructures with such symmetries. This will be verified in Sec. VII.
Iv Precise Algorithms to Compute Both and
Despite the fact that surface correlation functions contain crucial microstructural information, the technical difficulty involved in computing them has been a stumbling block in their widespread use. Methods have been devised to compute the surface correlation functions for dispersions of spheres that rely on dilating the interfaces Seaton and Glandt (1986); Klatt and Torquato (2016). A schematic illustration of how the algorithm works for is presented in Fig. 3, where is the dilation thickness. The algorithm simply measures the two-point probability function of the dilated phase and then one takes the appropriate limit of . Surface-surface correlation functions have also been used as input information to reconstruct two-phase digitized materials by Jiao, Stillinger, and Torquato Jiao et al. (2009). Since reconstruction algorithms require numerous evaluations of evolving microstructures, the surface correlation functions were approximated to improve computational speed.
iv.1 Algorithmic Details
Here we describe efficient general algorithms that enable the precise determination of the surface-surface correlation function and the surface-void correlation function for most situations that one may encounter in simulations and experiments. We consider -dimensional statistically homogeneous and isotropic two-phase systems within a cubic fundamental simulation cell of side length under periodic boundary conditions. We also assume that the interfaces are differentiable almost everywhere with exceptions for corners and edges only.
The idea behind the algorithm is that one can reduce the complexity of the problem by extracting information from a cut of the -dimensional statistically homogeneous and isotropic system with a -dimensional subspace () Torquato (2013). For example, the fully three-dimensional two-point correlation function of such a two-phase system is the same as the one-dimensional two-phase system formed from the cut of the original system with an infinitely long line. Similar ideas can be exploited to compute surface correlation functions as well. A straight line intersects with the interface in and leaves infinitely many intersection points, in principle. We can recover the fully three-dimensional surface correlation functions of the system by analyzing these intersections, but here we need to weight the points in accordance with the fact that the line cuts through the interface at different angles at each intersection point. In particular, because the interface projects to the line differently, each intersection point carries the weight , where is the acute angle between the straight line and the normal vector at the intersection point. From a “dilation” point of view, the straight line will cut through the dilated phase and leave line segments with lengths , then will reduce to the pair correlation function of intersection points with weights in the limit of .
Using this simple observation, the calculation of the surface-surface correlation function consists of the following steps:
1. Generate a straight line parallel to one of the edges of the box at a random position.
2. Find all the intersection points (, ,…) with interfaces of the system along this straight line. Store their positions , ….
3. Find the normal vectors at each intersection point and the angles (the acute one) between the straight line and these norm vectors , …. Compute , ….
4. Bin the distance between every pair of intersection points (suppose the size of each bin is ). Add to the corresponding bin.
5. Normalize the value in each bin by dividing .
6. Repeat the process from the beginning.
7. Compute the average of the results.
The calculation of the surface-void correlation function consists of the following steps:
1. Generate a straight line parallel to one of the edges of the box at a random position.
2. Find all the intersection points (, ,…) with interfaces of the system alone this straight line. Store their positions , ….
3. Find the normal vectors at each intersection point and the angles (the acute one) between the straight line and these norm vectors , …. Compute , ….
4. Generate random points along the straight line. Determine whether each point is in the void phase or not. Suppose , ,… are the points in the void phase, store their positions , ….
5. Bin the distance between every pair of and (suppose the size of each bin is ). Add to the corresponding bin.
6. Normalize the value in each bin by dividing .
7. Repeat the process from the beginning.
8. Compute the average of the results.
A schematic plot that elucidates our new algorithm is shown in Fig. 4. For systems with hard-wall boundary conditions (the usual case for experimental images), the value in the th bin should be multiplied by a factor due to the fact that fewer pairs can be formed near both ends of the boundaries. One can also easily generalize the algorithm to anisotropic media and to higher-order correlation functions such as and Torquato (2013).
For a -dimensional system consisting of particles (or voxels), the complexity for generating a single sampling line is . Computing each pair of intersection points on the line requires . The number of sampling lines is usually a preset number, and thus the overall complexity for computing surface-surface correlation function is . By a similar analysis, we know the complexity for computing surface-void correlation function is . Note that our algorithms are as efficient as the approximation method used in reconstructions Jiao et al. (2009), but with much better accuracy. Actually, as increases, fewer lines are needed, since each line contains more intersection points. If the total number of pairs we want to sample is fixed, both algorithms can give constant time complexity.
iv.2 Testing Against the Benchmark of Overlapping Spheres
Three-dimensional overlapping sphere systems provide an excellent benchmark to test our algorithm, since the surface correlation functions are known exactly; see Eq. (21) and Eq. (22). We generate a single but large configuration consisting of 250,000 overlapping spheres with a reduced density and particle-phase volume fraction . We generate one million straight lines at random locations and on each line we generate 1000 random points () in the case of computing . As we can see from Fig. 5, the theoretical and simulation results for the surface correlation functions are in excellent agreement with one another, even at the nondifferentiable point , indicating that the algorithm works remarkably well. As we discussed in Sec. III, the surface-surface correlation function diverges at the origin. We also ran our algorithm at other particle-phase volume fractions and again find excellent agreement with the corresponding theoretical results.
V Computing Surface Correlation Functions for Digitized Two-Phase Media
Unlike continuous-space microstructures (e.g., overlapping spheres), where we know surfaces exactly, images of heterogeneous materials are necessarily digitized, which presents algorithmic challenges to identify surfaces and normal vectors. We devote this section to the discussion of how to apply the aforementioned algorithm to this practical setting. Considering that experimental images are generally gray scale, we first discuss the case in which the two-phase medium is obtained from a level cut of a digitized scalar field in Torquato (2013). This common way to produce a two-phase medium enables us to identify interface normal vectors by the gradient of the scalar field. We then apply this idea to black and white images by first converting the given two-phase medium to a scalar field.
v.1 Two-phase Media Obtained From Level Cuts of Scalar Fields
Suppose we set a threshold to convert a scalar field to a two-phase medium: regions that satisfy constitute phase 1, and regions that satisfy constitute phase 2. The phase indicator function for phase 1 is given by
where is the Heaviside step function. The interface between two phases is simply the contour defined by . For any point on the contour, the normal vector is defined by the gradient of the scalar field, i.e., .
The algorithm can be implemented in essentially the same way as discussed in Sec. IV, but must be specialized to digitized two-phase media. In order to locate points of intersection of the line with interfaces, we need to find where changes sign along a straight line, and then interpolate the position of the point. The gradient at the point can be computed approximately by the finite differences of its neighboring pixels. However, the most significant difference between dealing with continuous models and digitized media is that the number of sampling straight lines one can afford is bounded by the resolution in the later case. Indeed, for an image, one can only sample at most times if the sampling straight lines are lined up with the grid. Thus the resolution of the image is crucial to obtain reliable results.
Here we use Gaussian random field Adler and Taylor (2009) as an example to demonstrate the importance of resolution. The field is generated by a superposition of 10000 plane waves, as we employed elsewhere Ma and Torquato (2017), to give a rather disordered structure. The results for the surface-surface correlation function are summarized in Fig. 6. Here we considered the field within a fixed square region but with different resolutions , , and . We also include a continuum result which is calculated by directly solving the contour and computing the gradient analytically. It can be seen that as the resolution increases, the numerical results rapidly converge to the continuum result.
v.2 Significance of the Threshold
Figure 6 shows that the computed fluctuates widely when the resolution is low. We discuss the origin of this behavior and how to deal with it in this subsection.
To begin, consider the simple situation illustrated in Fig. 7, which involves a straight line sampling the boundary of a unit circle. The discussion of this case is instructive because the vicinity of the intersection point can be approximated by sphere surfaces in most cases, and when is large enough, is proportional to . So for simplicity, the aim here is to estimate for the lower left quarter of the circle. Suppose the straight line samples from to uniformly along the direction that is perpendicular to itself; then we have
The integral correctly gives the surface area of the lower left quarter of the circle. Notice that although the integrand is divergent at , it is still integrable because the probability of hitting the vicinity of the singularity is proportionally infinitesimally small. However, it is easy to see that the variance of is divergent since the integral of diverges, which implies that large deviations can result when estimating the mean of . However, the simulation results suggest that increasing the number of sampling lines still reduces the fluctuations from the expected value, and a large sampling number yields good estimates, as one can see in Fig. 5. This suggests that the probability of getting a large deviation diminishes when the sampling number is increased. This is indeed the case, as we show in Appendix B.
However, in the case of digitized media, one cannot increase the sampling number arbitrarily. On the other hand, the probability of hitting the vicinity where can be rounded to a relatively large fraction due to the finite resolution. For example, a curved interface can align parallel to the sampling line after the digitization. The consequence is that we are more likely to encounter large deviations, as one can see in Fig. 6, where the abnormal peaks [as well as the universal trend of overestimating ] are due to certain very large values of encountered in the sampling. Although both problems can be alleviated by simply increasing the resolution, it is generally not known a prior that what resolution is required. Furthermore, obtaining high-resolution representations can also be computationally or economically costly, or simply beyond access due to the limitation of experimental techniques or available memory for a simulation. These restrictions force us to come up with a more efficient way to bypass the problems of digitized media. A straightforward way to remove this effect is to simply discard samples when they are larger than a certain threshold , i.e., . The bias induced by this method is usually small and insignificant, but with this small compromise, one can significantly reduce fluctuations (see a detailed analysis in Appendix C). To demonstrate the effect of applying thresholds to digitized media, we take the lowest resolution representation () of the Gaussian random field in the last subsection and recompute with a threshold . The result is shown in Fig. 8. Note that after applying a threshold, the fluctuations are dramatically suppressed and the result is much closer to the continuum result, even comparable to the ones with much higher resolutions in Fig. 6.
v.3 Converting Digitized Two-phase Media into Scalar Fields
We complete our discussion of two-phase media by discussing the case in which all of the information provided about the system is a binary digitized medium. Due to the jagged interface geometry, the transition from one phase to another is sharp and there is no easy way to estimate the direction of the norm vector as we did in the case of Gaussian random fields by computing the gradient of the scalar field.
Here we propose a straightforward method to deal with this situation. We first convert the two-phase medium to a coarse-grained scalar field, then convert it back to a two-phase medium by thresholding. In this way we can again use the algorithm introduced in Sec. V. A. We follow the procedure described in Refs. Torquato (2013) and Blumenfeld and Torquato (1993). By taking pixels in phase as source points, we can convert the two-phase medium into a scalar field by convolving the indicator function with a kernel or filter . Then, the scalar field is
where represents the parameters of the kernel. One of the most common choices of kernels is the Gaussian filter,
where is a length parameter that controls the size of the “influence” region of the filter. By taking a level cut of the scalar field at a threshold , we can then convert the scalar field back into a two-phase medium. The threshold is chosen to retain the original volume fraction of phase .
We include an example of overlapping spheres in two dimensions processed by this method. We prepare a digitized realization of 10,000 overlapping disks under periodic boundary conditions at a particle-phase volume fraction . The resolution is chosen to be that the side length of a pixel is of the diameter of the disk. We apply the Gaussian filter mentioned above for different values of and compute the corresponding of the converted scalar fields. A comparison of a portion of the system before and after applying the filter () is shown in Fig. 9, one can see that the structure of the system is maintained while there is a transition region between two phases. The comparison of computed with different filters is shown in Fig. 10 along with the exact result computed from the continuum model. It is noteworthy that although all the surface-surface correlation functions computed capture the shape of the exact one, they all tend to underestimate the actual function. The possible explanation is that the digitized version loses detailed interfacial information and hence the interface appears to be less curved, which leads to smaller surface areas. However, as decreases and the filter becomes more localized, the difference between the computed and the exact one monotonically diminishes. The smallest value of shown in Fig. 10 is three times of the pixel width, one may expect that when the resolution is high enough, the curve of the digitized version will finally converge to the exact one of the underlying pattern. Although images obtained in experiments may not necessarily be of high resolution, generally they are gray-scale images, which means one can simply use the algorithm described for scalar fields directly.
Vi Results for overlapping and nonoverlapping Sphere packings
In this section, we compute the surface correlation functions of several particle systems, including overlapping spheres, hard-spheres in equilibrium and decorated “stealthy” patterns. Given our abilities to compute the surface-surface correlation function, we can calculate
local surface-area variances through Eq. (10), and compare them with local volume-fraction variances in these systems.
We start by analyzing overlapping spheres in three dimensions. The volume-fraction variance and surface-area variance for the same system studied in Sec. IV are presented in Fig. 11(a), where is the radius of the spherical window and is the radius of particles. These quantities are computed by numerically computing the integrals in Eqs. (7) and (10) as well as through Monte Carlo simulations. In the later method, we generate windows at random positions and calculate the volume-fraction and surface-area variances directly. Since in this model spheres can form very complex clusters, we evaluate the volume fraction and surface area inside each window by generating random points uniformly in the window or on the surface of spheres and counting their fractions inside or on the surface of the clusters correspondingly. The theoretical prediction and simulation results agree very well, as one can see in Fig. 11(a). It is also noteworthy that in Fig. 11(a) the surface-area variance is much larger compared to the volume-fraction variance. However, one can see that in Fig. 11(b), after multiplied by (in order to show the large- behavior of fluctuations), it is clear that there is a crossover of function values around . Further numerical experiments show that the crossover only happens when the particle-phase volume fraction is between 0.57 and 0.7, outside this interval the surface-area variance is always larger than the volume-fraction variance, suggesting that the surface-area variance is a more sensitive descriptor.
We further compare the surface-area variance with the volume-fraction variance of hard spheres in equilibrium in three dimensions at different packing fractions. The variances are again computed using Eqs. (7) and (10); however, the autocovariance functions are not known analytically in this case. As mentioned previously, we use the results included in Ref. Torquato (1986), which was computed using the Percus-Yevick approximation and the Verlet-Weis correction. We include our results in Fig. 12. Note that the surface-area variance is always larger than the volume-fraction variance across a large span of packing fractions. In Fig. 13, the surface-area variance and the volume-fraction variance are compared respectively at different packing fractions. As the packing fraction increases, the hard-sphere system becomes more short-range ordered Torquato et al. (2000), thus the variances are expected to drop. The overall trend of and is consistent with this intuition. However, the volume-fraction variances experience another crossover at small , while the surface-area variances drop monotonically and larger gaps can be seen between the curves for different packing fractions. These results strongly suggest that the surface-area variance is a more sensitive measure of microstructures of the system compared to the volume-fraction variance.
Finally, we compare surface-area variances of hard spheres in equilibrium and overlapping spheres at different volume fractions in Fig. 14. The fact that the hard-sphere systems always suppress surface-area fluctuations variances to a greater degree than those of overlapping spheres, which reflects the stronger pair correlations in the former system.
Besides using to compute surface-area variances, it can itself be used as a “fingerprint” to detect important structural information, such as short-range order or the hyperuniformity of the system. Here we consider a special hyperuniform point patterns called “stealthy” point patterns that were studied in a recent paper Zhang et al. (2016), and we follow the procedure of circumscribing each point with a sphere to make the system a two-phase medium. The “stealthy” patterns are generated in a simulation box with basis vectors , and . We compare two sets of stealthy point patterns, with the parameter and 0.46 (In general, the system with larger will have more short-range order). Each point is decorated with a variable-sized sphere. We include our results for in Fig. 15. When the decorated spheres are very small, such as the case in Fig. 15(a), they do not overlap with one another, and thus should reveal structural features of the underlying point pattern. Clearly, the curve corresponding to in Fig. 15(a) exhibits stronger features, which is consistent with the fact that the pattern is more short-range ordered than that for . As stated in Ref. Zhang et al. (2016), the system loses its hyperuniformity when spheres begin to touch each other. From Fig. 15, it is seen that as the diameter increases, the correlation function begins to lose its features, and ultimately these two “stealthy” cases becomes indistinguishable from each other as well as the corresponding correlation function for overlapping spheres, shown in Fig. 5(a). The dramatic decrease around corresponds to the fact that the correlation between any two points on the same sphere cannot contribute to the function value beyond , and thus reveals a characteristic length scale of the system. Moreover, although the two systems start almost at the same specific surface, the gap between two curves continues to increase, and in the end the system with has a much larger specific surface. This suggests that the system with has greater short-range order that keeps the decorated spheres from overlapping with one another, and thus leads to a larger specific surface.
Vii Results for Snapshots of Evolving Spatial Patterns
In this section, we go beyond the analysis of well-known sphere models and extend the application of our algorithm to other important disordered patterns encountered in the physical and biological sciences. Specifically, we focus on time-dependent pattern formation processes that are governed by the Cahn-Hilliard equation and the Swift-Hohenberg equation. These patterns have recently been shown to be hyperuniform and could have important applications in material science Ma and Torquato (2017). We determine correlation functions of snapshots of these patterns here.
vii.1 Spinodal decomposition patterns from the Cahn-Hilliard equation
The Cahn-Hilliard equation was introduced to describe phase separation by spinodal decomposition Cahn and Hilliard (1958) and has been applied to model alloys Rundman and Hilliard (1967), polymer blends Smolders et al. (1971), and even pattern formations in ecology Liu et al. (2013). In Fig. 16, we show two typical patterns generated by this equation. The left one is at critical quench, in which case the volume-fraction ratio for two phases is 1:1, while the right one is off critical quench and has volume-fraction ratio 2:8. The interface between the two phases is highlighted.
One important feature of the Cahn-Hilliard equation is that the system will enter a “scaling regime” after some time, and the system will remain statistically the same after scaled by a growing characteristic length. This provides an indirect way to check our algorithm on digitized media. We can compute the surface correlation functions at different times and then an appropriate scaling enables them to collapse onto a single curve.
The rescaled surface-surface and surface-void correlation functions at different times are shown in Figs. 17 and 18 for critical and off-critical quenches. The curves for different times do collapse onto each other, as expected, further justifying the accuracy of our algorithm. Note that although the two systems shown in Fig. 16 appear to be structurally different, the corresponding standard autocovariance functions in Fig. 19 are similar to one another. The inability to distinguish the structures of these two systems is easily overcome by complementing with the information content of and , as they differ greatly for these two systems. Specifically, one can see that in Fig. 18, the surface-void correlation function for the critical quench has a flat slope at the origin [as predicted by Eq. (35)], while the one for the off-critical quench has a downward slope at the origin. This can be well explained by the small- behavior of that was derived in Sec. III. A.. From Eq. (29), we know that the slope of at origin is proportional to the mean curvature of the system. In the case of critical quench, the surface consists of both concave and convex parts, whose contributions cancel each other out, and thus the mean curvature is zero. In the case of off critical quench, where the matrix is the solid phase and the droplets are taken to be the void phase, the mean curvature is apparently negative, and thus in Fig. 18(b) we see the curve slopes down initially. This example again demonstrates the value of surface correlation functions in characterizing complex patterns.
vii.2 Patterns from the Swift-Hohenberg equation
The Swift-Hohenberg equation was developed to study Rayleigh-Bénard (RB) convection in hydrodynamics and later it became a subject of interest on its own in pattern formations Cross and Greenside (2009). The pattern produced by this equation is usually labyrinth-like, and the width of the “channel” is determined by a pre-selected wave number . It has been shown that the patterns can have different degrees of hyperuniformity Ma and Torquato (2017) when some tuning parameters are changed, although they may appear to be structurally alike.
Here we compute and compare two surface-surface correlation functions for two patterns generated under different , namely and in the same way in the authors’ previous paper Ma and Torquato (2017). It has been shown that the later one is more long-range ordered, which is also justified in our plot of in Fig. 20. It is evident that the for is much more long-ranged than the one for . Both curves have sharp spikes when is integer times of , which corresponds to the fact that the underlying patterns consist of stripes with width of , leaving roughly parallel interfaces with the same spacing at short scales. Note that spikes in also occur in sphere systems but only at the single location (see Fig. 5 and Fig. 15 for examples), while the corresponding for these systems are smooth functions without sharp transitions (see Refs. Torquato (2013) and Yeong and Torquato (1998) for plots). This again shows that surface-surface correlations can be superior in detecting short-scale microstructural features compared to that of the standard two-point correlation function .
We also evaluate the local surface-area variances in these systems using and Eq. (10). The results are shown in Fig. 21. Note that the surface-area variance for scales like , implying hyperuniformity Torquato (2016a). However, the variance for scales even slower than . The explanation is that in the case of , the corresponding wavelength is too small compared to the pixel size, which makes the numerical integration in Eq. (10) unreliable.
Viii Conclusions and Discussion
In this paper, we developed efficient general algorithms that enable the sampling
of the surface-surface correlation function and the surface-void correlation function with heretofore unattained precision. Our algorithms have advantages over the traditional “dilation” method Seaton and Glandt (1986); Klatt and Torquato (2016). First, the dilation method can only be easily implemented when the interfaces are relatively smooth and easy to be parameterized (e.g., packings of spheres and ellipsoids). However, our algorithms can be easily adapted to treat general complex interfaces. Second, the dilation method is difficult to implement for digitized media, which greatly limits its application to experimental data. By contrast, we have shown that our algorithms can be straightforwardly applied to digitized media. Third, as the dilation thickness approaches to zero, the probability of hitting the dilated phase will proportionally decrease, which requires a large number of samplings to ensure the accuracy, and hence greater computational
time. However, in the extreme situation that the information of a large but single system is available, our algorithms can yield accurate results from a single sample, since it is possible for the straight line to penetrate the interface a sufficiently large number of times. Moreover, our algorithms can be generalized to compute three-point surface correlation functions Torquato (2013) straightforwardly. Application of our algorithms to a variety of model disordered microstructures reveals that surface-surface correlation function is a sensitive descriptor of small-scale structural features, especially compared to the information content of the standard two-point correlation function .
We also showed that the extracted surface correlation functions can be used to compute accurately the surface-area variance, a quantity that can be a more sensitive measure of microstructural fluctuations compared to the volume-fraction variance. Through examples of spinodal decomposition patterns, we showed that surface correlation functions contain information that supplements that of , and the small- behavior of , which is determined by the mean curvature of the system. In two dimensions, the total curvature of a closed simple curve is a constant , implying that when the system approaches a percolation threshold, the absolute value of the mean curvature will drop dramatically due the formation of large clusters. This observation suggests that the surface-void correlation function may aid in detecting the onset of continuum percolation, which is an interesting topic for future exploration. We also showed how surface-surface correlation functions can be used to determine the hyperuniformity of two-phase media using patterns generated by the Swift-Hohenberg equation as examples.
Lower-order correlation functions have been successfully used to infer the physical properties of random media as well as to reconstruct them. This bodes well for their use in machine learning in the area of material optimization Stenzel et al. (2017); Röding et al. (2017). We expect that the algorithms to compute precisely the surface-surface correlation function and the surface-void correlation function presented in this paper will equip the community with powerful computational tools to characterize the structure and physical properties of multiphase media, especially with respect to those physical processes that are intimately linked to the interfaces. In particular, our algorithms can be adapted in reconstruction algorithms Yeong and Torquato (1998); Jiao et al. (2009) with heretofore unattained accuracy without sacrificing computational speed.
A sample Matlab program that enables one to compute the correlation functions , and for three dimensional digitized media can be downloaded at Ref. lin ().
Acknowledgements.The authors are grateful to Michael Klatt, Ge Zhang, JaeUk Kim, Timothy Middlemas, and Duyu Chen for helpful discussions. This work was supported in part by the National Science Foundation under Award No. CBET-1701843.
Appendix A The small- behavior of surface correlation functions of systems with singularities
Since the surface-surface correlation function discussed in Sec. III is only approximated to the zeroth-order, we focus our attention here to the surface-void correlation function of systems with singularities. We include results for certain specific cases, namely, two-dimensional systems that all singularities are corners and three-dimensional overlapping spheres.
For two-dimensional systems that all singularities are corners, suppose the angle formed by a corner to the solid phase is , one can show that the for small can be written as
where is the set of all the points on the interfaces that are differentiable, and is the density of corners in the system. One interesting implication of this formula is that the expression of the surface-void correlation function for packings of equilateral polygons is the same as the expression for packings of their inscribed circles.
Singularities in three dimensions are much more complex to analyze. Here we only take overlapping spheres as an example. Surprisingly, this is a nontrivial model of a heterogeneous material, since the lack of spatial correlation implies that the particles may overlap to form complex clusters, and leave many nondifferentiable edges in the system. Using the same geometric approach in Sec. III, by naively plugging in as in Eq. (31), we should have
The extra negative term implies that we have overestimated the probability of falling into the void phase by neglecting the fact that another sphere (see the upper left sphere in Fig. 22, which we henceforth call the “invading” sphere) can approach to and intersect with our “test” sphere (see the “dotted” sphere in Fig. 22) and reduce its fraction of surface area covered in the void phase. Indeed, any sphere whose centroid lies in the concentric shell with radius to around the reference point will intersect with the “test” sphere (we do not consider spheres that are closer than since then the reference point
would no longer be on the interface). Here we evaluate the reduced fraction of surface area in the void phase of the “test” sphere.
When the “invading” sphere overlaps with both the “test” sphere and the interface it can be difficult to evaluate the extra surface area of the “test” sphere covered by the “invading” sphere. Luckily, since the volume of the shell is , which is already first order, we can approximate the interface around the reference point as a flat plane that divides the “test” sphere into two hemispheres. By symmetry we know the extra surface area covered on average is just of the total surface area covered by the “invading” sphere on average. Suppose the distance between the center of the “invading” sphere and the reference point is , then the surface area of the spherical crown that is covered is
Letting , the fraction of surface area of the “test” sphere that is covered by the “invading” sphere is
Then on average the total fraction that is covered by the “invading” sphere is
Finally, by symmetry, we know the correction term to is or
, as in Eq. (42). This correction term will disappear in the dilute limit, since spheres will not overlap with one another.
One can carry out the same analysis for impenetrable spheres, but the calculation will be much more involved. In this case, the other sphere can only approach from the void phase, and the nonoverlapping condition will restrict its direction to a small range, which will in the end make the correction term of the order as long as is not a delta function, where is the diameter of a sphere. Thus our general formula will apply to hard spheres in equilibrium or random sequential addition (RSA) packings Zhang and Torquato (2013).
Appendix B The probability of getting an abnormal peak
It is instructive to estimate the probability , where is the number of sampling and stands for a given error. To do so, first we can sort such that for . Define
which is the smallest number of elements needed to make the inequality hold. Then we can decompose the probability by conditioning on , i.e.,
We can write down each term explicitly. When is large enough, we have
The reason why Eq. (B) has this scaling behavior is because is effectively zero when is larger than . One can continue this process and it is easy to see that the major contribution to the summation in Eq. (47) comes from the first a few terms. Thus one can see that by increasing one can significantly reduce the chance of getting an abnormal peak in the simulation, which is why the results in Fig. 5 are very smooth, especially compared to the “” case in Fig. 23(a), which has a much smaller sampling number.
Appendix C The effect of setting a threshold
We still take Fig. 7 as an example. By setting a threshold , the new expected value of is
The relative error is then . When , it is approximately . We can easily control this error by using a moderate threshold. For example, by setting , the error is already below . However, with this little compromise we can have a finite second moment, i.e.,
Thus, one can reduce the fluctuation to any level simply by adding more samples. When , the variance is approximated by . In the case of , one can reduce the relative standard deviation to by using around 12000 samples. It is noteworthy that the mean and variance have different asymptotic behaviors, the error diminishes as , while the variance only grows as , which gives us a great flexibility to choose .
The suppression of abnormal peaks by using a threshold can also be deduced from Eq. (47). Since now is bounded, it requires at least terms to have the inequality in Eq. (46), which means the leading terms in Eq. (47) vanish, leaving the probability significantly smaller than the case without a threshold (actually by Hoeffding’s inequality the probability will decrease exponentially).
We test this idea on two models: overlapping spheres and Gaussian random fields. In the former case, we consider the same system in Sec. IV, while restricting ourselves to a relatively small amount of sampling lines (10,000) and compare the computed with different thresholds , 100, 1000, and (no threshold). The result is shown in Fig. 23(a). Clearly, in all these cases, fluctuates around a common curve, while as the threshold is tightened, fluctuations are suppressed. To show this point quantitatively, we focus on the behavior of in the interval . As noted in Sec. IV, when the surface-surface correlation function is a constant in theory and should be a flat line in the plot; while in the simulations , fluctuates around a baseline. To compare the simulation and theoretical results, we compute the average of surface-surface correlation function and the largest deviation from in this interval. As shown in Fig. 23(b), the average of is always very close to the expected value (within error for all cases). However, by applying a more restrictive threshold, there is a trend to reduce the abnormal peaks significantly. We indeed get much smoother curves by making a very minor sacrifice of accuracy.
- Torquato (2013) S. Torquato, Random heterogeneous materials: microstructure and macroscopic properties, Vol. 16 (Springer Science & Business Media, 2013).
- Zohdi and Wriggers (2008) T. I. Zohdi and P. Wriggers, An introduction to computational micromechanics (Springer Science & Business Media, 2008).
- Sahimi (2003) M. Sahimi, Heterogeneous Materials I: Linear transport and optical properties, Vol. 22 (Springer Science & Business Media, 2003).
- Vasseur and Wadsworth (2017) J. Vasseur and F. B. Wadsworth, Bull. Volcanol 79, 77 (2017).
- Gibson and Ashby (1999) L. J. Gibson and M. F. Ashby, Cellular solids: structure and properties (Cambridge University Press, 1999).
- Gevertz and Torquato (2008) J. L. Gevertz and S. Torquato, PLoS Comput. Biol 4, e1000152 (2008).
- Peebles (1993) P. J. E. Peebles, Principles of physical cosmology (Princeton University Press, 1993).
- Chiu et al. (2013) S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke, Stochastic geometry and its applications (John Wiley & Sons, 2013).
- Mosser et al. (2017) L. Mosser, O. Dubrule, and M. J. Blunt, Phys. Rev. E 96, 043309 (2017).
- Lubbers et al. (2017) N. Lubbers, T. Lookman, and K. Barros, Phys. Rev. E 96, 052111 (2017).
- Jiao et al. (2007) Y. Jiao, F. H. Stillinger, and S. Torquato, Phys. Rev. E 76, 031110 (2007).
- Jiao et al. (2008) Y. Jiao, F. H. Stillinger, and S. Torquato, Phys. Rev. E 77, 031135 (2008).
- Jiao et al. (2009) Y. Jiao, F. Stillinger, and S. Torquato, Proc. Natl. Acad. Sci. USA 106, 17634 (2009).
- Torquato (1986) S. Torquato, J. Chem. Phys. 85, 4622 (1986).
- Debye and Bueche (1949) P. Debye and A. Bueche, J. Appl. Phys. 20, 518 (1949).
- Debye et al. (1957) P. Debye, H. Anderson Jr, and H. Brumberger, J. Appl. Phys. 28, 679 (1957).
- Torquato (2016a) S. Torquato, Phys. Rev. E 94, 022122 (2016a).
- Torquato (2016b) S. Torquato, J. Phys.: Condens. Matter 28, 414012 (2016b).
- Ma and Torquato (2017) Z. Ma and S. Torquato, J. Appl. Phys. 121, 244904 (2017).
- Yeong and Torquato (1998) C. L. Y. Yeong and S. Torquato, Phys. Rev. E 57, 495 (1998).
- Zachary and Torquato (2011) C. E. Zachary and S. Torquato, Phys. Rev. E 84, 056102 (2011).
- Guo et al. (2014) E.-Y. Guo, N. Chawla, T. Jing, S. Torquato, and Y. Jiao, Mater Charact 89, 33 (2014).
- Torquato (1990) S. Torquato, Phys. Rev. Lett 64, 2644 (1990).
- Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003).
- Zachary and Torquato (2009) C. E. Zachary and S. Torquato, J. Stat. Mech.: Theory Exp. 2009, P12015 (2009).
- Lu and Torquato (1990) B. Lu and S. Torquato, J. Chem. Phys. 93, 3452 (1990).
- Strey (1994) R. Strey, Colloid Polym. Sci. 272, 1005 (1994).
- Dietrich and Haase (1995) S. Dietrich and A. Haase, Phys. Rep. 260, 1 (1995).
- Doi (1976) M. Doi, J. Phys. Soc. Jpn 40, 567 (1976).
- Torquato and Stillinger (2006) S. Torquato and F. H. Stillinger, Phys. Rev. E 73, 031106 (2006).
- Seaton and Glandt (1986) N. Seaton and E. Glandt, J. Chem. Phys. 85, 5262 (1986).
- Klatt and Torquato (2016) M. A. Klatt and S. Torquato, Phys. Rev. E 94, 022152 (2016).
- Teubner (1990) M. Teubner, J. Chem. Phys. 92, 4501 (1990).
- Ciccariello et al. (1981) S. Ciccariello, G. Cocco, A. Benedetti, and S. Enzo, Phys. Rev. B 23, 6474 (1981).
- (35) J. Von Neumann, Am. Soc. for Metals, Cleveland, 1952 , p. 108.
- Mecke (1998) K. R. Mecke, Int. J. Mod. Phys. B 12, 861 (1998).
- Ciccariello and Benedetti (1982) S. Ciccariello and A. Benedetti, Phys. Rev. B 26, 6384 (1982).
- Adler and Taylor (2009) R. J. Adler and J. E. Taylor, Random fields and geometry (Springer Science & Business Media, 2009).
- Blumenfeld and Torquato (1993) R. Blumenfeld and S. Torquato, Phys. Rev. E 48, 4492 (1993).
- Torquato et al. (2000) S. Torquato, T. M. Truskett, and P. G. Debenedetti, Phys. Rev. Lett. 84, 2064 (2000).
- Zhang et al. (2016) G. Zhang, F. H. Stillinger, and S. Torquato, J. Chem. Phys. 145, 244109 (2016).
- Cahn and Hilliard (1958) J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28, 258 (1958).
- Rundman and Hilliard (1967) K. Rundman and J. Hilliard, Acta Metall 15, 1025 (1967).
- Smolders et al. (1971) C. Smolders, J. Van Aartsen, and A. Steenbergen, Kollold. Z. Z. Polym. 243, 14 (1971).
- Liu et al. (2013) Q.-X. Liu, A. Doelman, V. Rottschäfer, M. de Jager, P. M. Herman, M. Rietkerk, and J. van de Koppel, Proc. Natl. Acad. Sci. USA 110, 11905 (2013).
- Cross and Greenside (2009) M. Cross and H. Greenside, Pattern formation and dynamics in nonequilibrium systems (Cambridge University Press, 2009).
- Stenzel et al. (2017) O. Stenzel, O. Pecho, L. Holzer, M. Neumann, and V. Schmidt, AIChE J. 63, 4224 (2017).
- Röding et al. (2017) M. Röding, P. Svensson, and N. Lorén, Comput. Mater. Sci 134, 126 (2017).
- (49) http://chemlabs.princeton.edu/torquato/links-and-codes.
- Zhang and Torquato (2013) G. Zhang and S. Torquato, Phys. Rev. E 88, 053312 (2013).