Transport, Geometrical and Topological Properties of Stealthy Disordered Hyperuniform Two-Phase Systems

# Transport, Geometrical and Topological Properties of Stealthy Disordered Hyperuniform Two-Phase Systems

G. Zhang Department of Chemistry, Princeton University, Princeton NJ 08544    F. H. Stillinger Department of Chemistry, Princeton University, Princeton NJ 08544    S. Torquato Department of Chemistry, Department of Physics, Princeton Institute for the Science and Technology of Materials, and Program in Applied and Computational Mathematics, Princeton University, Princeton NJ 08544
###### Abstract

Disordered hyperuniform many-particle systems have attracted considerable recent attention, since they behave like crystals in the manner in which they suppress large-scale density fluctuations, and yet also resemble statistically isotropic liquids and glasses with no Bragg peaks. One important class of such systems is the classical ground states of “stealthy potentials.” The degree of order of such ground states depends on a tuning parameter . Previous studies have shown that these ground-state point configurations can be counterintuitively disordered, infinitely degenerate, and endowed with novel physical properties (e.g., negative thermal expansion behavior). In this paper, we focus on the disordered regime () in which there is no long-range order, and control the degree of short-range order. We map these stealthy disordered hyperuniform point configurations to two-phase media by circumscribing each point with a possibly overlapping sphere of a common radius : the “particle” and “void” phases are taken to be the space interior and exterior to the spheres, respectively. The hyperuniformity of such two-phase media depend on the sphere sizes: While it was previously analytically proven that the resulting two-phase media maintain hyperuniformity if spheres do not overlap, here we show numerically that they lose hyperuniformity whenever the spheres overlap. We study certain transport properties of these systems, including the effective diffusion coefficient of point particles diffusing in the void phase as well as static and time-dependent characteristics associated with diffusion-controlled reactions. Besides these effective transport properties, we also investigate several related structural properties, including pore-size functions, quantizer error, an order metric, and percolation thresholds. We show that these transport, geometrical and topological properties of our two-phase media derived from decorated stealthy ground states are distinctly different from those of equilibrium hard-sphere systems and spatially uncorrelated overlapping spheres. As the extent of short-range order increases, stealthy disordered two-phase media can attain nearly maximal effective diffusion coefficients over a broad range of volume fractions while also maintaining isotropy, and therefore may have practical applications in situations where ease of transport is desirable. We also show that the percolation threshold and the order metric are positively correlated with each other, while both of them are negatively correlated with the quantizer error. In the highly disordered regime (), stealthy point-particle configurations are weakly-perturbed ideal gases. Nevertheless, reactants of diffusion-controlled reactions decay much faster in our two-phase media than in equilibrium hard-sphere systems of similar degrees of order, and hence indicate that the formation of large holes is strongly suppressed in the former systems.

## I Introduction

A hyperuniform many-particle system is one in which the structure factor approaches zero in the infinite-wavelength limit.Torquato and Stillinger (2003) In such systems, density fluctuations (measured by the variance of number of particles inside a spherical window) are anomalously suppressed at very large lengths scales, a “hidden” order that imposes strong global structural constraints.Torquato and Stillinger (2003); Torquato (2016a) All structurally perfect crystals and quasicrystals are hyperuniform,Torquato and Stillinger (2003); Zachary and Torquato (2009) but typical disordered many-particle systems, including gases, liquids, and glasses, are not. Disordered hyperuniform many-particle systems are exotic states of amorphous matter that have attracted considerable recent attention.Torquato and Stillinger (2003); Zachary and Torquato (2009); Donev, Stillinger, and Torquato (2005); Zachary, Jiao, and Torquato (2011); Jiao and Torquato (2011); atkinson2016critical (); Kurita and Weeks (2011); Dreyfus et al. (2015); Lesanovsky and Garrahan (2014); Hexner and Levine (2015); Jack, Thompson, and Sollich (2015); De Rosa et al. (2015); Degl’Innocenti et al. (2015); Xie et al. (2013); Muller et al. (2014); Florescu, Torquato, and Steinhardt (2009); Torquato, Zhang, and Stillinger (2015); Uche, Stillinger, and Torquato (2004); Zhang, Stillinger, and Torquato (2015a, b); Batten, Stillinger, and Torquato (2009); Novikov et al. (2014); Torquato (2016a); Xu, Douglas, and Freed (2016) Materials that are simultaneously disordered and hyperuniform can be regarded to be exotic states of matter that lie between a crystal and a liquid; they behave more like crystals in the manner in which they suppress large-scale density fluctuations, and yet they also resemble typical statistically isotropic liquids and glasses with no Bragg peaks.Torquato, Zhang, and Stillinger (2015)

An important class of disordered hyperuniform many-particle systems is comprised of the classical ground states of “stealthy potentials,”Uche, Stillinger, and Torquato (2004); Torquato, Zhang, and Stillinger (2015); Zhang, Stillinger, and Torquato (2015a, b) which are bounded, long-range, pairwise additive potentials designed in Fourier space. These classical ground states are of particular fundamental interest because they can be degenerate and noncrystalline. A nonnegative parameter inversely proportional to the number density, , controls the degree of order of such ground states. For , the ground states are overwhelmingly highly degenerate and disordered. As increases above 0.5, long-range translational and rotational order begins to emerge and eventually the system crystallizes. We have previously studied these disordered ground states, and computed their pair correlation functions,Uche, Stillinger, and Torquato (2004); Uche, Torquato, and Stillinger (2006); Batten, Stillinger, and Torquato (2008); Zhang, Stillinger, and Torquato (2015a); Torquato, Zhang, and Stillinger (2015) structure factors,Uche, Stillinger, and Torquato (2004); Uche, Torquato, and Stillinger (2006); Batten, Stillinger, and Torquato (2008); Zhang, Stillinger, and Torquato (2015a); Torquato, Zhang, and Stillinger (2015) Voronoi cell volume distribution,Uche, Stillinger, and Torquato (2004); Zhang, Stillinger, and Torquato (2015a) and particle-exclusion probabilities.Torquato, Zhang, and Stillinger (2015)

Some initial studies have demonstrated that stealthy hyperuniform systems are endowed with novel thermodynamic and physical properties. For example, their low-temperature excited states are characterized by negative thermal expansion behavior.Batten, Stillinger, and Torquato (2009) It has also been shown that dielectric networks derived from stealthy disordered hyperuniform point configurations possess complete photonic band gaps comparable in size to those of a photonic crystal, while at the same time maintain statistical isotropy, enabling waveguide geometries not possible with photonic crystals as well as high-density disordered transparent materials.Florescu, Torquato, and Steinhardt (2009); Florescu, Steinhardt, and Torquato (2013); Man et al. (2013); Leseur, Pierrat, and Carminati (2016) However, the determination of physical/chemical properties of stealthy disordered hyperuniform materials is generally an unexplored area of research.

In this paper, we investigate steady-state and time-dependent diffusion properties of certain decorations of stealthy disordered hyperuniform ground-state point configurations in two and three dimensions. In particular, we derive two-phase heterogeneous media from point configurations by decorating the point configurations with spheres (circles); specifically, all points are circumscribed by spheres of radius that generally may overlap with one another. By varying the radius, the fraction of space occupied by the spheres will vary. We study the effective transport properties of these disordered two-phase systems, including the effective diffusion coefficient,Torquato (2013) and static and time-dependent characteristics of diffusion-controlled reactions at the interfaces between the two continuous phases, as well as the trapping rate (or its inverse, the mean survival time) as well as the principal (largest) relaxation time.Reck and Prager (1965); Torquato and Avellaneda (1991) Quantifying the effective diffusion coefficient is of importance not only because it has direct applications (e.g., diffusion of fuel and oxygen in a fuel cell Zamel et al. (2010), diffusion tensor magnetic resonance imaging, Sen and Torquato (1989); Tuch et al. (2001) regulation and metabolism of normal organs, Comper (1996); Gevertz and Torquato (2008) and drug release from porous matrices Lemaire, Belair, and Hildgen (2003)), but also because its determination translates immediately into equivalent results for the effective thermal and electric conductivity, the effective dielectric constant, and the effective magnetic permeability for reasons of mathematical analogy,Torquato (2013) and is therefore related to a host of applications. Diffusion-controlled reactions arise in widely different processes, such as heterogeneous catalysis,Baiker (1999) gas sensor operation,Sakai et al. (2001) cell metabolism,Rohde and Price (2007) crystal growth,Watson and Müller (2009) and nuclear magnetic resonance (NMR).Banavar and Schwartz (1987); Mitra and Sen (1992); Straley et al. (1987)

These transport properties are related to several statistical geometrical and topological characteristics, which we therefore also study. These include the pore-size functions (the distribution of the distance from a randomly chosen location in the void phase to the closest phase boundary),Torquato (2010) the quantizer error (a moment of the pore-size function, which is related to the principal relaxation time),Torquato and Avellaneda (1991); Torquato (2010) the order metric (a measure of the translational order of point configurations),Torquato, Zhang, and Stillinger (2015) and the percolation threshold or the critical radius (the radius of the spheres at which a specific phase becomes connected) of each phase.Rintoul and Torquato (1997); Lorenz and Ziff (2001); Quintanilla, Torquato, and Ziff (2000); Quintanilla and Ziff (2007)

We compare the aforementioned physical and geometrical properties of our two-phase system derived from decorated stealthy ground states, as a function of the tuning parameter , with those of two other two-phase media: (1) equilibrium disordered (fluid) hard-sphere systems and (2) decorated Poisson point processes (ideal-gas configurations). The former has short-range order that is tunable by its volume fraction but no long-range order. The latter has neither short-range order nor long-range order. Through comparison, we find that some of these quantities are dramatically affected by the degree of long-range order, while other quantities are much more sensitive to the degree of short-range order. Because many of these quantities depend on the density, we re-scale all systems to unit number density to ensure a fair comparison.

Among our major findings, we show that these transport, geometrical and topological properties of our two-phase media are generally distinctly different from those of equilibrium hard-sphere systems and spatially uncorrelated overlapping spheres. At high values, the stealthy disordered two-phase media can attain nearly maximal effective diffusion coefficient, while also maintaining isotropy. This novel property could have practical implications, e.g., optimal and isotropic drug release from designed nanoparticles. Stealthy ground states tend to ideal gases configurationally in the limit. Torquato, Zhang, and Stillinger (2015) Nevertheless, we find that even in the low- regime, our two-phase media have much lower principal relaxation time than that of equilibrium hard-sphere systems of similar degrees of order, indicating that the formation of large holes in the stealthy systems is strongly suppressed. Lastly, we also find that the aforementioned geometrical and topological quantities are strongly correlated with each other.

The rest of the paper is organized as follows: In Sec. II, we give precise definitions of the stealthy potential and the aforementioned transport, geometrical and topological quantities. In Sec. III, we present our numerical method to calculate them. We present our results in Sec. IV and conclusions in Sec. V.

## Ii Mathematical definitions and background

### ii.1 Preliminaries

This paper studies properties of point-particle systems as well as two-phase heterogeneous media derived from certain decorations of these point configurations. A point-particle system consists of point particles with a certain probability density function , where , , …, is the particle positions, in a simulation box of volume under periodic boundary conditions in -dimensional Euclidean space , where is 2 or 3. The number density is defined as . The “Poisson point process” (also called “ideal gas”) is produced by the probability density function that does not depend on particle positions . The equilibrium hard-sphere point process of radius is another point process with equal to a positive constant if the distance between every pair of points is larger than and zero otherwise.

A realization of a two-phase medium can be mathematically described as a partition of a domain of space with volume into two separate regions, and . It is characterized by an indicator function, , where is any position in the two-phase medium. The indicator function is one if and zero if . The volume fraction of phase 1 is given by , where denotes an ensemble average. That of the other phase is given by . Let be the interface between and , the specific surface, i.e., the total area of divided by , is given by:

 s=<|∇I(x)|>. (1)

The two-phase media that we consider here are derived from point configurations by decorating the point configurations with spheres (circles); specifically, each point is circumscribed by a sphere of radius that generally may overlap with one another. Therefore, it is composed of a void region (phase 1) and a particle region (phase 2). When such a mapping is applied to a Poisson point process, the decorated system is also called “fully penetrable spheres” Quintanilla and Ziff (2007) or “spatially uncorrelated spheres.” Lorenz and Ziff (2001)

### ii.2 Stealthy potentials and their entropically favored ground states

Consider point processes that are obtained from the canonical ensemble probability distribution function defined by

 P(rN)=exp[−βΦ(rN)]/Z, (2)

where is an interaction potential, is the inverse temperature, and is the partition function. Of particular interest in this paper is the “stealthy” interaction potential:

 Φ(rN)=12vF∑0

where the sum is over all reciprocal lattice vector ’s of the simulation box such that , ,

 Φ0=[N(N−1)−∑0

is a constant independent of the particle positions , and the second equal sign in Eq. (3) can be proved by Parseval’s theorem. Such potential is interesting not only because it is a pairwise additive potential [as the right side of Eq. (3) shows], but also because it allows one to directly tune the structure factor

 S(k)=|~n(k)|2/N. (5)

The ground state (i.e., or zero-temperature state) of this potential is obtained by constraining for all .Uche, Stillinger, and Torquato (2004); Torquato, Zhang, and Stillinger (2015)

Let be half the number of points in the summation of Eq. (3) 111Since , is the number of independent constraints. ; the parameter

 χ=Md(N−1) (6)

determines the degree to which the ground states are constrained and therefore the degeneracy and disorder of the ground states.Uche, Stillinger, and Torquato (2004) For , the ground states are typically disordered and uncountably infinitely degenerate.Torquato, Zhang, and Stillinger (2015); Zhang, Stillinger, and Torquato (2015a) Therefore, there are multiple ways to assign different weights (i.e., probabilities) to different sets of ground states. One way of particular interest is the zero-temperature () limit of Eq. (2). Ground states drawn from such distribution are called “entropically favored ground states”.Torquato, Zhang, and Stillinger (2015); Zhang, Stillinger, and Torquato (2015a) It is interesting to note that in the and limit, both entropically favored ground states of stealthy potentials and equilibrium hard-sphere point processes tend to Poisson point process geometrically. In the rest of the paper this fact will be frequently used to test our simulation results since many properties of the Poisson point process have been studied previously.

### ii.3 Transport properties

This paper studies the following steady-state and time-dependent diffusion properties in phase 1 (the void phase) of decorated entropically favored ground states of stealthy potentials, and compare them with that of decorated Poisson point process and equilibrium disordered (fluid) hard-sphere system at unit number density.

#### ii.3.1 Effective diffusion coefficient

Consider the steady-state diffusion problem of some species with concentration field in a two-phase medium in which phase 1 is the space in which diffusion occurs and phase 2 are “obstacles” that the diffusing species cannot enter. In phase 1, the flux of the species, , is predicted by Fick’s first law:

 J(x)=D∇c(x), x∈V1 (7)

where is a diffusion coefficient which we set to unity for simplicity. However, Eq. (7) is valid only in phase 1 and has to be paired with the following Neumann boundary condition:

 n⋅J=0, on ∂V, (8)

where is the normal vector of the surface. We see that the inclusion of such obstacles adds a complicated boundary condition and makes the overall diffusion problem difficult. Nevertheless, on a length scale much larger than the characteristic length of the obstacles, the system can be homogenized Torquato (2013) and characterized by an “effective” diffusion coefficient, , defined by the average Fick’s first law:

 =De<∇c(x)>, for any x (9)

where angular brackets denote ensemble averages.

The effective diffusion coefficient of an isotropic two-phase medium must satisfy the Hashin-Shtrikman (HS) upper bound.Hashin and Shtrikman (1963) For our case where phase 1 has unit diffusion coefficient and phase 2 cannot be entered, this bound in dimensions is given by:

 De≤d−1d−1+ϕ2. (10)

This bound is optimal because it is realizable by certain model microstructures, including the “coated-sphere model” described in Ref. Hashin, 1962, and is therefore the best possible bound for isotropic systems given volume-fraction information only.

#### ii.3.2 Diffusion-controlled reactions

Consider the problem of diffusion and reaction among absorbing “traps” in the random medium. Let phase 1 be the region in which diffusion occurs and phase 2 be the trap region, the diffusion process in phase 1 is governed by the same Fick’s first law but with time dependency:

 J(x,t)=D∇c(x,t), in V1. (11)

This equation, combined with the conservation of the diffusing species inside phase 1, , yields Fick’s second law:

 ∂c(x,t)∂t=D△c(x,t), in V1. (12)

If phase 2 are absorbing “traps” (rather than impenetrable obstacles as in the aforementioned effective diffusion problem), the boundary condition has to be changed. In the diffusion-controlled limit, i.e., when the reaction rate at the interface is infinite, we have the following boundary condition:Torquato and Avellaneda (1991)

 c(x,t)=0, on ∂V. (13)

If we also set the initial concentration to be uniform outside of traps:

 c(x,0)=c0, in V1, (14)

then we have the survival problem. The “survival probability,” is equal to the fraction of reactant not yet absorbed at time :Torquato and Avellaneda (1991); Torquato (2013)

 p(t)=∫Rdc(x,t)dx∫Rdc(x,0)dx. (15)

The mean survival time of the reactant is the zeroth moment of :222Note that the commonly used notation for the mean survival time is Torquato and Avellaneda (1991); Torquato (2013). Here, we use to avoid confusion with the order metric .

 Tmean=∫∞0p(t)dt. (16)

The survival probability can be decomposed as a sum of exponential functions:

 p(t)=∞∑n=1Inexp(−t/Tn), (17)

where are coefficients and are relaxation times. The largest relaxation time is called “principal relaxation time” and by convention denoted . These quantities can be measured directly by NMR experiments since in NMR experiment of fluid-saturated porous media, proton magnetization decays mainly on the phase boundary.Straley et al. (1987); Banavar and Schwartz (1987); Mitra and Sen (1992)

It is worth noting that although the above problems involve differential equations, , , and can actually be calculated much more efficiently by simulating Brownian motions using the so-called “first-passage time” technique. (See the Sec. III for details.) The effective diffusion coefficient can be found from the ratio of the mean square displacement of such Brownian particles and the time spent. The survival probability is equal to the probability that a Brownian particle have never reached any trap at time . The mean survival time, , can be calculated by integrating but can also be calculated, more easily, by finding the average time needed for a particle to reach a trap the first time. It is also worth noting that while the effective diffusion coefficient is identically zero as long as the void phase is not percolating, and are both positive until the spheres cover the entire space.Torquato (2013)

### ii.4 Geometrical and topological properties

This paper also studies the following geometrical and topological properties that are intimately related to the aforementioned diffusion characteristics.

#### ii.4.1 Hyperuniformity and stealthiness in many-particle systems and two-phase media

As we have explained earlier, a hyperuniform many-particle system is one in which the structure factor, Eq. (5), approaches zero in the limit. The name “hyperuniform” refers to an anomalous suppression of density fluctuations: Consider random placements of a spherical observation window of radius in a -dimensional many-particle system. The number of points contained in such window, , is a random variable. For a uniform but not hyperuniform many-particle system (e.g., ideal gas without a gravity field), for large scales as . For a hyperuniform system, for large grows more slowly than . It has been proved that the above-mentioned two conditions of hyperuniformity, and for large grows more slowly than , are mathematically equivalent.Torquato and Stillinger (2003)

A similar definition exists for two-phase media.Torquato (2016a) One can compute the volume fraction of either phase inside a spherical observation window of radius and find its variance. For large , this variance scales as for typical (non-hyperuniform) random two-phase media and decreases faster than for hyperuniform two-phase media. An equivalent condition for hyperuniformity is that , where

 χV(k)=1vF|J(k)|2 (18)

is called the “spectral density” and is the Fourier transform of .Torquato (1999)

Stealthy hyperuniform many-particle systems or two-phase media are subsets of hyperuniform many-particle systems or two-phase media in which or is zero for a range of vectors around the origin, i.e.,

 (19)

where is some positive number. For the many-particle systems mentioned in this paper, the ground state of “stealthy” potentials are stealthy and hyperuniform while equilibrium hard-sphere systems and Poisson point process are neither stealthy nor hyperuniform.

#### ii.4.2 Packing and packing fraction

When we decorate a point-particle configuration by replacing points with spheres of radius , the whole collection of spheres is considered a “sphere packing” if each pair of point particles is separated by a distance of at least (i.e., if the spheres do not overlap). The fraction of space occupied by the union of spheres, , is called the packing fraction . Of particular interest in this paper is the maximum packing radius , which is half the minimum separation distance between two particles, and maximum packing fraction , which is the volume fraction of phase 2 when .

Why should we study the maximum packing fraction? One important reason is that when we decorate a point configuration and map it into a two-phase medium, if spheres do not overlap, then the spectral density of the two-phase medium is proportional to the structure factor of the underlying point configuration:Torquato (2016b)

 ~χV(k)=ϕ2v1(a)(2πa|k|)dJ2d/2(|k|a)S(k) (a≤amaxp), (20)

where is the volume of a -dimensional sphere of radius and is the Bessel function of order . Therefore, a decorated stealthy point configuration is a stealthy two-phase medium if . When , however, Eq. (20) no longer holds and we will see in Sec. IV.1 that decorated systems are generally no longer stealthy or hyperuniform.

#### ii.4.3 Nearest-neighbor and pore-size functions

Given a point-particle system, the void-exclusion probability is the probability that a spherical cavity of radius , centered at a random location, is empty of particles. A related quantity is , the probability density function of the distance to the nearest particle from a randomly chosen location. A different interpretation of is that if each point particle is replaced with a sphere of radius , then is the volume fraction of the space outside of the spheres, i.e.,

 EV(a)=ϕ1=1−ϕ2. (21)

Since is the negative derivative of , is the specific surface .Torquato (2013)

Another quantity related to is the scaled dimensionless quantizer error . For a point configuration with positions , a quantizer is a device that takes as an input a position in and outputs the nearest point of the configuration to . Assuming is uniformly distributed, one can define a mean square error, which can be obtained from via the relation:Torquato (2010)

 G=2ρ2dd∫∞0rEV(r)dr. (22)

Finally, two more related quantities can be defined for two-phase media. The pore-size cumulative distribution function, , is defined as the fraction of pore space (i.e., space covered by phase 1) which has a pore radius larger than . The function of our decorated system is trivially related to of the underlying point-particle system:

 F(δ)=EV(δ+a)EV(a). (23)

Moreover, the associated pore-size probability density function is given by . This pore-size function at the origin is related to the specific surface, , by

 P(δ=0)=sϕ1. (24)

It is interesting to note that the moments of are related to the mean survival time and principle relaxation time via the following rigorous lower bounds Torquato and Avellaneda (1991):

 Tmean≥1D(∫∞0F(δ)dδ)2, (25)

and

 T1≥2D∫∞0δF(δ)dδ. (26)

We see that is proportional to the first moment of in the limit and is therefore related to the principal relaxation time.

#### ii.4.4 Order metric τ

We will be studying the above properties for systems of varying degrees of order. Therefore, it is desirable to have a way to quantify such order. Moreover, since the underlying point-configurations we study include both stealthy ground states, which have long-range order, and equilibrium liquid hard-sphere systems, which have short-range order, we desire an order metric that reflects short-range order and long-range order equally well. A suitable choice is the order metric , introduced in Ref. Torquato, Zhang, and Stillinger, 2015 and defined as:

 τ=1Dd∫∞0[g2(r)−1]2dr=1(2π)dDd∫∞0[S(k)−1]2dk, (27)

where is some characteristic length scale, is the pair correlation function,Chandler (1987) is the angular average of , and the second equal sign can be proved by Parseval’s theorem. In this paper, we simply let because we always rescale the configuration to make the number density unity.

#### ii.4.5 Percolation threshold and critical radius

Since the effective diffusion coefficient is trivially zero when the void phase is topologically disconnected, it is important to quantify when the phases are connected. To do this, we will be considering the percolation properties of the systems. As we specified earlier, we map point configurations into two-phase media by replacing each point with a sphere of radius . For phase 2, the critical or percolation radius, , is the minimum such that a connected part of phase 2 becomes infinite in size. The percolation volume fraction, , is the fraction of space occupied by the union of spheres of radius .

We can define similar percolation characteristics of the void phase.Rintoul (2000); Priour Jr (2014) The percolation radius of the void phase, , is defined as the maximum such that there is still an infinite-sized connected part of phase 1. The percolation volume fraction, , is the volume fraction of phase 1 at radius . In two dimensions, it is very rare to have both phases percolating simultaneously (see Ref. Sheng and Kohn, 1982 for such a rare example). In our case, and . In three dimensions, however, both phases can simultaneously percolate, i.e., the two-phase system is bicontinuous. Indeed, this is the case for our 3D systems and hence we must compute and separately.

## Iii simulation details

### iii.1 Generating entropically favored stealthy ground states

We generate entropically favored ground states of stealthy potentials using the same protocol as our previous work.Zhang, Stillinger, and Torquato (2015a) This protocol involves performing molecular dynamics (MD) simulations at a very low temperature ( in 2D and in 3D in dimensionless units), taking snapshots periodically, and performing a local energy minimization starting from each snapshot. Because the MD temperature is sufficiently low, the snapshots before energy minimization are already very close to ground states. Therefore, the ground states produced by the subsequent energy minimization closely follow the canonical distribution in the zero-temperature limit. We generate 20,000 configurations per value, same as Ref. Zhang, Stillinger, and Torquato, 2015a. The only two differences between this work and our previous work Zhang, Stillinger, and Torquato (2015a) are (1) system sizes are different (see Appendix A for our choice of system sizes and the justification), and (2) each configuration is rescaled to unit number density (in order to ensure a fair comparison).

### iii.2 Generating equilibrium disordered hard-sphere systems

We also generate equilibrium disordered hard-sphere systems via standard Monte-Carlo techniques in order to compare their statistics with entropically favored stealthy ground states’ statistics. Depending on the packing fraction , an equilibrium hard-sphere system can be disordered (liquid-like) or crystalline. Disordered equilibrium hard-sphere system exists for in 2D and in 3D.Torquato (2013) Therefore, the packing fraction we used include , 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, and 0.68 in 2D and , 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, and 0.48 in 3D. For each in each dimension, we generate equilibrium hard-sphere systems with , 300, and 500 particles. In each case, the system was first equilibrated with trial moves. After that, we sample a configuration every trial moves until we obtain 20,000 configurations. Similar to the stealthy ground states, we keep the number density . Therefore, we adjust sphere radius to attain a desired packing fraction.

### iii.3 Calculating survival probability, mean survival time, and principal relaxation time

Because the method we used to calculate the effective diffusion coefficient is an extension of the method to calculate survival probability and mean survival time, we will explain the latter method first. The survival probability and mean survival time can be calculated by simulating particles undergoing Brownian motions.

The Brownian motion can be simulated very efficiently using the first-passage-time technique.Torquato and Kim (1989) The key idea of this technique is that for a Brownian particle at a particular location, let be the distance between it and the closest phase boundary. Construct a sphere centered at the particle with radius (which is called a first-passage-time sphere). Let be the time needed for the particle to reach the surface of such sphere for the first time, the distribution of can be calculated analytically. In 3D, the cumulative distribution function (CDF) of isTorquato and Kim (1989)

 F(tR)=1+2∞∑m=1(−1)mexp(−Dm2π2tRR2). (28)

In 2D, Ref. Torquato and Kim, 1989 did not provide the distribution of . Here we find the following explicit 2D expression for :

 F(tR)=1−2∞∑m=1exp(−Dw2mtR/R2)wmJ1(wm), (29)

where is the Bessel function of order , and is the th root of . The mean of , in any dimension, is simply .

Therefore, the Brownian motion inside the first-passage-time sphere does not need to be simulated in detail. One simply moves the particle to a random location on the surface of such sphere, and increase the time by a certain amount, as detailed below. When calculating the mean survival time , the time increment can simply be , the mean of . When calculating , however, the time increment has to be a random number drawn from the distributions given in Eq. (28) or (29). The process of finding , moving the particle, and increasing the time is repeated until the Brownian particle gets very close () to a trap, at which time the Brownian particle is deemed trapped. In our implementation, Eqs. (28)-(29) are pre-computed and tabulated to accelerate the simulation. For each configuration, we simulate 10 Brownian trajectories to calculate and 1000 trajectories to calculate . When calculating , each trajectory is additionally sampled 100 times, with different random time increments drawn from distributions (28)-(29).

After calculating , we calculate the principal relaxation time by fitting in the range to the asymptotic equation

 ln[p(t)]≈c+t/T1, (30)

where and are fitting parameters.

### iii.4 Calculating effective diffusion coefficient

The effective diffusion coefficient can also be calculated using first-passage-time techniques. kim1990determination (); Kim and Torquato (1992a); torquato1999effective () In this case, however, the Brownian particle cannot be deemed trapped when it is sufficiently close to the phase boundary because phase 2 is now non-absorbing obstacles rather than absorbing traps. Instead, we construct a first-passage-time sphere of radius , find a random place on the surface of the first-passage-time sphere that is outside of the obstacle phase, and move the Brownian particle to that random place. Although this first-passage-time sphere contains two phases, the mean time taken for the Brownian particle to reach such surface could still be computed analytically and was given in Ref. Kim and Torquato, 1992a:

 tR=R2(1+v2/v1)2d, (31)

where is the volume of the obstacle phase divided by the volume of the conducting phase inside the first-passage-time sphere and can be found analytically.

The process of constructing a first-passage-time sphere and moving the point particle is repeated to form a Brownian trajectory. In the infinite-time limit, the effective diffusion coefficient is given by:Kim and Torquato (1992a)

 De=limt→∞<|R(t)|2>2dt, (32)

where is the mean-squared displacement of a Brownian particle at time . In practice, in a finite-time simulation, one should only consider the time regime in which the mean square displacement is strongly linear in time, since for sufficiently early times the mean square displacement is either ballistic or grows faster than linear in time.torquato1999effective () We find by fitting versus and extracting the slope of the line after some sufficiently large dimensionless time. We define the unit of time to be

 t∗=1ρ2/dD, (33)

and set both and to be unity. The point in time in which first becomes a strongly linear function occurs when the Brownian particle sufficiently samples the two-phase system such that it can be viewed effectively as Brownian motion in a homogeneous medium. For the microstructures that we considered here, we find that the linear regime occurs in the dimensionless time interval , i.e., we determine from the linear relationship

 <|R(t)|2>=(2dDe)t+c, 40

To get a sense of the possible behaviors of the mean square displacements as a function of time, we show examples in Fig. 1 at several values of for a three-dimensional system at and indicate the linear fit in each case. This fitting procedure works especially well near percolation, which is the most difficult regime to simulate. For this particular system, the void phase stops percolating at obstacle radius . Figure 1 shows that for , is linear with . For , only a fraction of configurations still have a percolating void phase, and our fitting procedure was able to distinguish the initial uprise in (contributions mainly from Brownian particles moving inside a “cage”, i.e., a disconnected part of the void phase) from the steady increase in (contributions from Brownian particles that are in a percolating part of the void phase). For , all Brownian particles are caged, and the fit has a virtually zero slope (and therefore produces a virtually zero ).

We simulate 1 Brownian trajectory per configuration to calculate . In Fig. 2, we compare the computed with the distribution of the void-phase percolation threshold and find that becomes zero right after all configurations stop percolating. The fact that our measured diminishes to zero at the percolation threshold indicates that our choice of the fitting range in Eq. (34) is appropriate.

### iii.5 Calculating percolation thresholds

Generally speaking, the precise calculation of the percolation threshold of disordered systems require very large system sizes. For example, to accurately determine the percolation threshold of 3D fully penetrable spheres, Ref. Lorenz and Ziff, 2001 employed systems of up to particles. The whole system is divided into smaller cubes and the content particles in each cube is generated only when such cube is being probed.

Unfortunately, our protocol of low-temperature MD and a subsequent energy minimization does not allow us to save time by only generating required parts of the configuration. Moreover, in order to accurately follow the canonical-ensemble distribution at zero-temperature limit, the MD temperature has to be so low such that many () time steps are required to produce a sufficiently long trajectory. The requirement of a very large number of time steps forces us to further sacrifice system size. As a result, our system is limited to several hundred particles. Therefore, accurate determination of the percolation threshold is extremely challenging. Thus we experimented with two advanced algorithms to minimize finite-size effect in order to obtain relatively accurate results. We will first explain how to use these two methods to determine the percolation threshold for the particle phase, and then describe the generalizations to the void phase.

One of them, which we call “ maximum method,” is described in Ref. Newman and Ziff, 2000. Starting from a random particle in a configuration, one randomly chooses two of its periodic images in two different directions. The quantity (denoted as in Ref. Newman and Ziff, 2000) is defined as the probability that this particle is connected to one of the chosen periodic images but not the other. At the percolation threshold, attains its maximum. Therefore, one can numerically find as a function of sphere radius and find its maximum in order to find the percolation threshold. In our implementation, we calculate for various ’s starting from , with increment , until develops a peak and then returns to zero. We then select all data points such that , where is the maximum of , and perform a quadratic fit of the selected data points. The maximum of the fitted function gives the percolation radius .

Ref. Newman and Ziff, 2000 measures the percolation radii using several different system sizes and then extrapolates to the infinite-system-size limit. However, when we perform the same fitting procedure using different system sizes, we did not find a clear trend: In each dimension, for some values larger systems produces larger while for other values larger systems produces smaller . Moreover, the extrapolated , as a function of , is not as smooth as the un-extrapolated one. We therefore conclude that random noise is probably more important than finite-size effect in this case and extrapolation is not proper. Thus, we will simply use of our largest system as an estimate of the infinite-system-size .

After finding the percolation radius , we determine the percolation volume fraction . One could have simply read this quantity from a plot of the quantity , since . However, we decide to use a somewhat more accurate method: we divide the whole simulation box into 1200012000 pixels (in 2D) or 120012001200 voxels (in 3D) and find out if the center of each pixel or voxel is inside any sphere of radius . We then count the number of pixels or voxels that are centered inside spheres to find out the volume fraction. From our experience, this procedure gives us a four-significant-figures precision in .

The other method we employed, which we call “ intersection method,” is introduced in Ref. ziff2016percolation, . At a given radius , define to be the size of the largest cluster, (denoted as in Ref. ziff2016percolation, ) is defined as:

 M2=⟨s2max⟩−⟨smax⟩2⟨smax⟩2, (35)

where denotes an ensemble average. As Ref. ziff2016percolation, shows, at the percolation threshold, is the same for different system sizes. Therefore, one can compute as a function of , and find the intersection of for different system sizes to find . Following Ref. ziff2016percolation, , we use three different ’s for each value, and perform an extrapolation to find in the infinite- limit. After that, we use the same procedure discussed in the previous paragraph to calculate from .

We have used both methods to calculate the percolation volume fraction of decorated stealthy ground states at various ’s in 2D and 3D. They are presented in Fig. 3. Figure 3 also presents for decorated Poisson point processes obtained from Refs. Lorenz and Ziff, 2001 and Quintanilla, Torquato, and Ziff, 2000, which are in 2D and in 3D. These results can be used as benchmarks since Poisson point processes are geometrically equivalent to entropically favored stealthy ground states at . We see that in 2D, while both methods give results that approaches the Poisson value very well in the limit, the maximum method produces much smoother results. In 3D, however, although both methods produce relatively smooth results, only results from the intersection method approaches the Poisson value very well in the limit. Therefore, we decide to choose the maximum method in 2D and the intersection method in 3D for the rest of the paper. It is interesting to note that to our knowledge, the maximum method has been demonstrated to work well in 2D Newman and Ziff (2000) but not in 3D, while the intersection method has been demonstrated to work well in 3D ziff2016percolation () but not in 2D. It is possible that these two methods are just more suited to their respective dimensions.

Besides the percolation threshold of the spheres, we also study the percolation threshold of the void phase. In two dimensions, the percolation radius of the void phase, , is equal to the percolation radius of the spheres, . In three dimensions, however, has to be calculated separately. We compute in three dimensions by performing a Voronoi tessellation of each configuration, and then computing of the Voronoi vertices. As in the particle-phase case, the intersection of at different system sizes gives . Similar to the particle-phase case, can then be converted to by digitization, the result of which is presented in Fig. 4. Similar to the particle-phase case, we compare for our systems with that for the decorated Poisson point processes obtained from Ref. Priour Jr, 2014, . Combining the and results, we see that as increases from 0 to 0.46, the range for bicontinuity moves upwards, from to , respectively.

### iii.6 Calculating EV(r), G, and τ

The quantities and are calculated by first computing . For each configuration of point particles, random locations in 2D or random locations in 3D are generated in the simulation box. For each location, the distance from it to its nearest particle is found. These distances are then binned to yield . We then integrate using trapezoidal rule to find . The quantizer error is obtained by another integration of , using trapezoidal rule, as Eq. (22) shows. The numerically obtained always have compact support, and thus the above-mentioned integrations does not need to be truncated.

The order metric can be computed from either or , as Eq. (27) shows. We have tried both approaches. The real-space integration in Eq. (27) is truncated at half the simulation box side length and the reciprocal space integration in Eq. (27) is truncated at , where is the cutoff of the stealthy potential (as detailed in Sec. II.2).

## Iv Results

We present visualizations of our two-phase systems derived from decorated stealthy ground states in Figs. 5 - 7. In three dimensions, we present separate figures for the particle phase and the void phase for clarity. In the rest of the section, we present the above-mentioned properties of our two-phase systems, and compare them with decorated Poisson point process and hard-sphere point process.

### iv.1 Packing fraction and stealthiness

We present the maximum packing fraction of decorated stealthy ground states, , in Fig. 8. In each dimension, as increases, remains to be zero for up to about 0.3 and then start to increase. This indicates that for , particles in entropically favored stealthy ground states can become arbitrarily close to each other. As becomes higher, particles develop an effective hard core that are impenetrable. The development of such hard core was also observed in Ref. Uche, Stillinger, and Torquato, 2004.

When we decorate a stealthy ground state and map it into a two-phase medium, if , then Eq. (20) ensures that the resulting two-phase medium is also stealthy. However, if , will the resulting two-phase medium also be stealthy or hyperuniform? To answer this question, we decorated a two-dimensional stealthy ground state of particles at with several different sphere radii , digitized the resulting two-phase medium into pixels, and calculated the spectral density using Eq. (18). The result is presented in Fig. 9. For this particular system, the maximum packing radius is . We see that for , is zero for . For , however, is positive and does not tend to zero as . Therefore, a decorated stealthy ground state is generally neither stealthy nor hyperuniform if .

### iv.2 Effective diffusion coefficient

We present the calculated effective diffusion coefficient for our two-phase systems derived from decorated stealthy ground states in Fig. 10. It is interesting to note that in two dimensions, the curves of cross over each other for different values of : while for smaller higher produces a higher , for larger higher produces a smaller . An explanation for such phenomenon will be presented in the next paragraph.

It is also useful to plot versus the particle-phase volume fraction, , by mapping to using Eq. (21). We present such plots in Fig. 11. These plots show that higher values (more ordered arrangements of the obstacle phase) always produce higher at the same volume fraction, which is consistent with our intuition: a more ordered arrangement of the obstacles leaves more space between them, and produces a higher . So why did we see the opposite relationship between and in Fig. 10, except for smaller in 2D? It turns out that a lower induces more overlap between the spherical obstacles and thus results in a lower . This in turn produces a higher .

With plotted versus , it is interesting to compare our result with the HS upper bound given in Eq. (10). We make such comparison in Fig. 11. Our result is consistent with the upper bound for any and except for small fluctuations, but the bound is sharp only for smaller .

If it is desired to find structures that maximizes , then any one of the degenerate structures that achieves the HS bound is optimal. We see that our two-phase systems derived from decorated stealthy ground states at very high ’s are very close to being optimal for up to 0.4-0.5. In two dimensions, this range coincides with at high ’s. Since decorated stealthy ground states loses stealthiness as increases beyond , our results suggest that a loss of stealthiness causes to stop being optimal. In three dimensions, however, is less sensitive to structures. Thus, although decorated stealthy ground states (at high ’s) stops being a packing at around , does not deviate from the optimal value until about . Our observation that is less sensitive to structures in 3D than in 2D is consistent with the trends indicated in Ref. Torquato, 1998, which found that in the infinite- limit, is given exactly by the arithmetic average of the diffusion coefficients of the two phases (weighted by their volume fraction), independent of the structure.

In Fig. 11 we also present of decorated lattice structures (i.e., periodic arrays of spherical inclusions). Since these lattice structures are stealthy with even higher values, unsurprisingly, their sticks with the HS upper bound for an even larger range.

Lastly, we would like to mention a difference between the support of as a function of in 2D versus 3D. While for for our two-phase systems in 2D diminishes to zero at , in 3D does not vanish until . This difference emerges from the difference in the topological (connectedness) characteristics of the void phase between these dimensions.

### iv.3 Survival probability and mean survival time

We have computed the mean survival time, , of a diffusing reactant with unit diffusion coefficient, in our two-phase systems derived from decorated stealthy ground states. These results are summarized in Fig. 12. For comparison, the same quantity for equilibrium disordered hard-sphere systems are also included. Clearly, increasing order (increasing for stealthy ground states or increasing for hard spheres) suppresses . However, there is a crossover between the curves for stealthy ground states and that for equilibrium disordered hard spheres. This crossover is expected because as increases, an equilibrium hard-sphere system becomes more ordered, and therefore comparable to a stealthy two-phase medium with a higher . In Fig. 13, we plot versus for and 0.5. We see that in 2D, is somewhat more sensitive to than in 3D.

In Fig. 14 we present the survival probability, , at , for our two-phase systems derived from decorated stealthy ground states and equilibrium disordered hard spheres. The same crossover phenomenon also appears here, suggesting that the long-range order possessed by stealthy ground states suppresses at large more efficiently, while the short-range order possessed by equilibrium disordered hard spheres suppresses at small more efficiently.

In Fig. 15 we present the principal relaxation time for our two-phase systems derived from decorated stealthy ground states. It turns out that is much more sensitive to in 2D than in 3D. More interestingly, one can compare of stealthy ground states and equilibrium disordered hard disks at the same order metric . We present such comparison in Fig. 16. In two dimensions, one can see that at , of equilibrium disordered hard disks is much higher than that of our two-phase systems derived from decorated stealthy ground states with similar ’s. As we explained earlier, is related to the pore-size distribution. Therefore, our results suggest that hyperuniformity suppresses the formation of large holes, even in the very disordered regime. As increases to 0.5, however, the difference between the two systems diminishes. Our finite-sized simulation results suggest that at this value of , even equilibrium hard-sphere systems suppress the formation of large holes very well. 333We should clarify that in the infinite-system-size limit, of equilibrium disordered hard-sphere systems is actually infinite because of a non-zero probability of forming arbitrarily large holes (i.e., is non-zero for arbitrarily large ). Torquato and Avellaneda (1991) For finite-sized systems, however, is finite. For decorated stealthy ground states, it is unclear whether or not would be infinite in the infinite-system-size limit.

### iv.4 Geometrical and Topological Properties

The percolation volume fraction for both phases in 2D and 3D was already presented in Fig. 3 and 4. The void-exclusion probability , quantizer error , and order metric are presented in Figs. 17-19. We see that in each dimension, as increases, increases, at any decreases, decreases, and increases.

The order metric can be computed from either or , as shown in Fig. 19. In 2D, the results from these two approaches have good consistency. However, in 3D, computed from is often slightly lower than computed from . We discovered that this is because is still oscillating around 1 at half the simulation box side length, where the integration in Eq. (27) has to be cut off. Therefore, such a cutoff should make computed from too low. We thus use computed from in the rest of the paper. It is seen that is very sensitive at detecting the rise in short-range and long-range order as increases.

### iv.5 Correlations between geometrical properties, and comparison with equilibrium disordered hard-sphere systems

In Fig. 20, we explore the correlation between , , and , and compare that for entropically favored stealthy ground states with that for equilibrium disordered hard-sphere systems. It is interesting to note that these two systems behave very differently. At the same , stealthy systems give lower values of , indicating is more sensitive to long-range order than to short-range order. At the same or , stealthy systems give higher values of , indicating that is even more sensitive to long-range order than to short-range order.

## V Conclusions and discussion

In this work, we decorated stealthy disordered hyperuniform point configurations of different degrees of order with spheres of various radii, and computed several transport and structural properties of these decorated systems. The transport properties that we studied include effective diffusion coefficient , mean survival time , survival probability , and principal relaxation time . The structural properties examined include hyperuniformity and stealthiness, maximum packing fraction , the void-exclusion probability , the order metric , and the percolation thresholds and . We showed that the order metric is an exquisite detector of both short- and long-range translational order. While all geometrical and topological quantities are strongly correlated (positive correlation between and , and negative correlation between and the former two quantities), the relation between the physical quantities are more complex: While increases as increases or as decreases, and increases as decreases or as decreases. Therefore, there is no simple relationship between and or . Another reason why there is no such relationship is because if phase 1 ceases to percolate, then becomes zero but and are still positive Torquato (2013).

Besides finding correlations between geometrical and topological properties, we find that in the highly disordered () regime, of our two-phase systems derived from decorated stealthy ground states is much lower than that of equilibrium hard-sphere system. Together with the void-exclusion probability, low suggests that the formation of large holes is strongly suppressed, even though the configuration appears completely disordered.

In the higher order () regime, of our disordered isotropic two-phase systems derived from decorated stealthy ground states is very close to the Hashin-Shtrikman upper bound for , where is the maximum packing fraction. Since such decorated systems maintain stealthiness if and only if , our results suggest a connection between stealthiness and the ability to have a nearly optimal (maximal) . The fact that stealthy disordered two-phase media have nearly optimal could have practical implications, e.g., optimal and isotropic drug release from designed nanoparticles. Although nearly optimal can also be achieved by lattice structures (i.e., periodic arrays of inclusions), the latter are always anisotropic. Thus, if one desires isotropic two-phase media with highest possible at a specific volume fraction, disordered stealthy two-phase media could be the best choice.

Disordered stealthy ground states are uncountably infinitely degenerate.Zhang, Stillinger, and Torquato (2015a) The maximum packing fraction varies among configurations. In the future, it would be interesting to design algorithms that sample stealthy ground states with a bias toward configurations with higher values. With such an algorithm, one would be able to design isotropic two-phase media with nearly optimal with very high .

Lastly, we would like to mention that although here we only study the diffusion problem of point Brownian particles, the diffusion problem of finite-sized Brownian particles has also been of interest.Ashworth et al. (2015) It is noteworthy that our results can be trivially extended to the latter case. The diffusion of Brownian particles of radius among obstacles of radius is equivalent to the diffusion of point Brownian particles among obstacles of radius . This mapping was previously exploited to quantify diffusion of finite-sized spheres in various models of porous media.Kim and Torquato (1992b)

Interestingly, one can relate the transport properties computed here (, , and ) to different physical properties of the same systems via cross-property relations, including those that relate them to the elastic moduli, Gibiansky and Torquato (1993, 1996) as well as fluid permeability. Avellaneda and Torquato (1991); Torquato (1990) In future work, we will carry out such analyses.

###### Acknowledgements.
The authors thank Duyu Chen and Jaeuk Kim for their careful reading of the manuscript.

## Appendix A System sizes

As discussed in Sec. III.1, it is nontrivial to choose the system size and parameter , especially because one of our protocol to calculate the percolation threshold requires three different ’s for each . We enumerated all possible choices of ’s and ’s for and picked up some values that allow at least three different choices of ’s. Our choice of and in 2D and 3D are listed in Tables 1 and 2. It is desirable to consider values of higher than 0.4133… in 3D, but our enumeration did not find such a value that satisfies the above condition. Therefore, we chose three more ’s that allow to be very close to but makes differ in the fifth decimal place. See the caption for Table 2 for details. Except for the percolation threshold calculation, we only use the largest for each and .

## Appendix B Properties of stealthy point configurations and decorated systems

In this section we tabulate all of the physical and geometrical properties of stealthy point configurations and decorated systems that we study in this paper (Tables III and IV).

## References

• Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003).
• Torquato (2016a) S. Torquato, Phys. Rev. E 94, 022122 (2016a).
• Zachary and Torquato (2009) C. E. Zachary and S. Torquato, J. Stat. Mech. Theor. Exp. , P12015 (2009).
• Donev, Stillinger, and Torquato (2005) A. Donev, F. H. Stillinger,  and S. Torquato, Phys. Rev. Lett. 95, 090604 (2005).
• Zachary, Jiao, and Torquato (2011) C. E. Zachary, Y. Jiao,  and S. Torquato, Phys. Rev. Lett. 106, 178001 (2011).
• Jiao and Torquato (2011) Y. Jiao and S. Torquato, Phys. Rev. E 84, 041309 (2011).
• (7) S. Atkinson, G. Zhang, A. B. Hopkins, and S. Torquato, Phys. Rev. E 94, 012902 (2016).
• Kurita and Weeks (2011) R. Kurita and E. R. Weeks, Phys. Rev. E 84, 030401 (2011).
• Dreyfus et al. (2015) R. Dreyfus, Y. Xu, T. Still, L. A. Hough, A. G. Yodh,  and S. Torquato, Phys. Rev. E 91, 012302 (2015).
• Lesanovsky and Garrahan (2014) I. Lesanovsky and J. P. Garrahan, Phys. Rev. A 90, 011603 (2014).
• Hexner and Levine (2015) D. Hexner and D. Levine, Phys. Rev. Lett. 114, 110602 (2015).
• Jack, Thompson, and Sollich (2015) R. L. Jack, I. R. Thompson,  and P. Sollich, Phys. Rev. Lett. 114, 060601 (2015).
• De Rosa et al. (2015) C. De Rosa, F. Auriemma, C. Diletto, R. Di Girolamo, A. Malafronte, P. Morvillo, G. Zito, G. Rusciano, G. Pesce,  and A. Sasso, Phys. Chem. Chem. Phys. 17, 8061 (2015).
• Degl’Innocenti et al. (2015) R. Degl’Innocenti, Y. D. Shah, L. Masini, A. Ronzani, A. Pitanti, Y. Ren, D. S. Jessop, A. Tredicucci, H. E. Beere,  and D. A. Ritchie, in SPIE OPTO (International Society for Optics and Photonics, 2015) pp. 93700A–93700A.
• Xie et al. (2013) R. Xie, G. G. Long, S. J. Weigand, S. C. Moss, T. Carvalho, S. Roorda, M. Hejna, S. Torquato,  and P. J. Steinhardt, Proc. Natl. Acad. Sci. 110, 13250 (2013).
• Muller et al. (2014) N. Muller, J. Haberko, C. Marichy,  and F. Scheffold, Adv. Opt. Mater. 2, 115 (2014).
• Florescu, Torquato, and Steinhardt (2009) M. Florescu, S. Torquato,  and P. J. Steinhardt, Proc. Natl. Acad. Sci. 106, 20658 (2009).
• Torquato, Zhang, and Stillinger (2015) S. Torquato, G. Zhang,  and F. H. Stillinger, Phys. Rev. X 5, 021020 (2015).
• Uche, Stillinger, and Torquato (2004) O. U. Uche, F. H. Stillinger,  and S. Torquato, Phys. Rev. E 70, 046122 (2004).
• Zhang, Stillinger, and Torquato (2015a) G. Zhang, F. H. Stillinger,  and S. Torquato, Phys. Rev. E  (2015a).
• Zhang, Stillinger, and Torquato (2015b) G. Zhang, F. H. Stillinger,  and S. Torquato, Phys. Rev. E  (2015b).
• Batten, Stillinger, and Torquato (2009) R. D. Batten, F. H. Stillinger,  and S. Torquato, Phys. Rev. Letters 103, 050602 (2009).
• Novikov et al. (2014) D. S. Novikov, J. H. Jensen, J. A. Helpern,  and E. Fieremans, Proc. Nat. Acad. Sci. 111, 5088 (2014).
• Xu, Douglas, and Freed (2016) W. Xu, J. Douglas,  and K. Freed, Macromolecules 49, 8341 (2016).
• Uche, Torquato, and Stillinger (2006) O. U. Uche, S. Torquato,  and F. H. Stillinger, Phys. Rev. E 74, 031104 (2006).
• Batten, Stillinger, and Torquato (2008) R. D. Batten, F. H. Stillinger,  and S. Torquato, J. Appl. Phys. 104, 033504 (2008).
• Florescu, Steinhardt, and Torquato (2013) M. Florescu, P. J. Steinhardt,  and S. Torquato, Phys. Rev. B 87, 165116 (2013).
• Man et al. (2013) W. Man, M. Florescu, E. P. Williamson, Y. He, S. R. Hashemizad, B. Y. C. Leung, D. R. Liner, S. Torquato, P. M. Chaikin,  and P. J. Steinhardt, Proc. Natl. Acad. Sci. 110, 15886 (2013).
• Leseur, Pierrat, and Carminati (2016) O. Leseur, R. Pierrat,  and R. Carminati, Optica 3, 763 (2016).
• Torquato (2013) S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Interdisciplinary Applied Mathematics (Springer New York, 2013).
• Reck and Prager (1965) R. Reck and S. Prager, J. Chem. Phys. 42, 3027 (1965).
• Torquato and Avellaneda (1991) S. Torquato and M. Avellaneda, J. Chem. Phys. 95, 6477 (1991).
• Zamel et al. (2010) N. Zamel, N. G. Astrath, X. Li, J. Shen, J. Zhou, F. B. Astrath, H. Wang,  and Z.-S. Liu, Chem. Engi. Sci. 65, 931 (2010).
• Sen and Torquato (1989) A. K. Sen and S. Torquato, Phys. Rev. B 39, 4504 (1989).
• Tuch et al. (2001) D. S. Tuch, V. J. Wedeen, A. M. Dale, J. S. George,  and J. W. Belliveau, Proc. Nat. Acad. Sci. 98, 11697 (2001).
• Comper (1996) W. Comper, Extracellular Matrix, Extracellular Matrix No. v. 1 (Taylor & Francis, 1996).
• Gevertz and Torquato (2008) J. L. Gevertz and S. Torquato, PLoS Comput. Biol. 4, e1000152 (2008).
• Lemaire, Belair, and Hildgen (2003) V. Lemaire, J. Belair,  and P. Hildgen, Int. J. Pharm. 258, 95 (2003).
• Baiker (1999) A. Baiker, Chem. Rev. 99, 453 (1999).
• Sakai et al. (2001) G. Sakai, N. Matsunaga, K. Shimanoe,  and N. Yamazoe, Sensor. Actuat. B Chem. 80, 125 (2001).
• Rohde and Price (2007) R. A. Rohde and P. B. Price, Proc. Nat. Acad. Sci. 104, 16592 (2007).
• Watson and Müller (2009) E. B. Watson and T. Müller, Chem. Geol. 267, 111 (2009).
• Banavar and Schwartz (1987) J. R. Banavar and L. M. Schwartz, Phys. Rev. Lett. 58, 1411 (1987).
• Mitra and Sen (1992) P. P. Mitra and P. N. Sen, Phys. Rev. B 45, 143 (1992).
• Straley et al. (1987) C. Straley, A. Matteson, S. Feng, L. M. Schwartz, W. E. Kenyon,  and J. R. Banavar, Appl. Phys. Lett. 51, 1146 (1987).
• Torquato (2010) S. Torquato, Phys. Rev. E 82, 056109 (2010).
• Rintoul and Torquato (1997) M. Rintoul and S. Torquato, J. Phys. A 30, L585 (1997).
• Lorenz and Ziff (2001) C. D. Lorenz and R. M. Ziff, J. Chem. Phys. 114, 3659 (2001).
• Quintanilla, Torquato, and Ziff (2000) J. Quintanilla, S. Torquato,  and R. M. Ziff, J. Phys. A 33, L399 (2000).
• Quintanilla and Ziff (2007) J. A. Quintanilla and R. M. Ziff, Phys. Rev. E 76, 051115 (2007).
• (51) Since , is the number of independent constraints.
• Hashin and Shtrikman (1963) Z. Hashin and S. Shtrikman, J. Mech. Phys. Solids 11, 127 (1963).
• Hashin (1962) Z. Hashin, J. Appl. Mech. 29, 143 (1962).
• (54) Note that the commonly used notation for the mean survival time is Torquato and Avellaneda (1991); Torquato (2013). Here, we use to avoid confusion with the order metric .
• Torquato (1999) S. Torquato, J. Chem. Phys. 111, 8832 (1999).
• Torquato (2016b) S. Torquato, J. Phys.: Condens. Matter 28, 414012 (2016b).
• Chandler (1987) D. Chandler, Introduction to Modern Statistical Mechanics (Oxford University Press, 1987).
• Rintoul (2000) M. Rintoul, Phys. Rev. E 62, 68 (2000).
• Priour Jr (2014) D. Priour Jr, Phys. Rev. E 89, 012148 (2014).
• Sheng and Kohn (1982) P. Sheng and R. Kohn, Phys. Rev. B 26, 1331 (1982).
• Torquato and Kim (1989) S. Torquato and I. C. Kim, Appl. Phys. Lett. 55, 1847 (1989).
• (62) I. C. Kim, and S. Torquato, J. Appl. Phys. 68, 3892 (1990).
• Kim and Torquato (1992a) I. C. Kim and S. Torquato, J. Appl. Phys. 71, 2727 (1992a).
• (64) S. Torquato, I. C. Kim, and D. Cule, J. Appl. Phys. 85, 1560 (1999).
• Newman and Ziff (2000) M. Newman and R. M. Ziff, Phys. Rev. Lett. 85, 4104 (2000).
• (66) R. Ziff and S. Torquato, https://arxiv.org/abs/1611.00279 (2016) .
• Torquato (1998) S. Torquato, J. Mech. Phys. Solids 46, 1411 (1998).
• (68) We should clarify that in the infinite-system-size limit, of equilibrium disordered hard-sphere systems is actually infinite because of a non-zero probability of forming arbitrarily large holes (i.e., is non-zero for arbitrarily large ). Torquato and Avellaneda (1991) For finite-sized systems, however, is finite. For decorated stealthy ground states, it is unclear whether or not would be infinite in the infinite-system-size limit.
• Ashworth et al. (2015) J. C. Ashworth, M. Mehr, P. G. Buxton, S. M. Best,  and R. E. Cameron, Adv. Healthc. Mater. 4, 1317 (2015).
• Kim and Torquato (1992b) I. C. Kim and S. Torquato, J. Chem. Phys. 96, 1498 (1992b).
• Gibiansky and Torquato (1993) L. Gibiansky and S. Torquato, Phys. Rev. Lett. 71, 2927 (1993).
• Gibiansky and Torquato (1996) L. Gibiansky and S. Torquato, Proc. R. Soc. A,  452, 253 (1996).
• Avellaneda and Torquato (1991) M. Avellaneda and S. Torquato, Phys. Fluids A 3, 2529 (1991).
• Torquato (1990) S. Torquato, Phys. Rev. Lett. 64, 2644 (1990).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters