New Tessellation-Based Procedure to Design Perfectly Hyperuniform Disordered Dispersions for Materials Discovery

TL;DR

This paper introduces a tessellation-based method to create large, perfectly hyperuniform disordered dispersions, enabling the design of novel materials with unique physical properties and overcoming previous computational and experimental limitations.

Contribution

A new simple construction procedure guarantees perfect hyperuniformity in large samples, transforming nonhyperuniform dispersions into hyperuniform ones and establishing their properties analytically and numerically.

Findings

Method guarantees perfect hyperuniformity in large samples.

Numerical demonstrations in 2D and 3D confirm the approach.

Establishes hyperuniformity of Hashin-Shtrikman dispersions.

Abstract

Disordered hyperuniform dispersions are exotic amorphous two-phase materials characterized by an anomalous suppression of long-wavelength volume-fraction fluctuations, endowing them with novel physical properties. While such unusual materials have received considerable attention, a stumbling block has been an inability to create large samples that are truly hyperuniform due to current computational and experimental limitations. To overcome such barriers, we introduce a new, simple construction procedure that guarantees perfect hyperuniformity for very large sample sizes. It involves tessellating space into cells and then inserting a particle into each cell such that the local-cell particle packing fractions are identical to the global packing fraction. We analytically prove that such dispersions are perfectly hyperuniform in the infinite-sample-size limit. Our methodology enables a remarkable mapping that converts a very large nonhyperuniform disordered dispersion into a perfectly hyperuniform one, which we numerically demonstrate in two and three dimensions. A similar analysis also establishes the hyperuniformity of the famous Hashin-Shtrikman multiscale dispersions, which possess optimal transport and elastic properties. Our hyperuniform designs can be readily fabricated using modern photolithographic and 3D printing technologies. The exploration of the enormous class of hyperuniform dispersions that can be designed and tuned by our tessellation-based methodology paves the way for accelerating the discovery of novel hyperuniform materials.

Equations

$${\chi }_{{}_V}\left(\mathbf{r}\right)\equiv S_2^{(i)}\left(\mathbf{r}\right)-{{\varphi }_i}^2,$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)=\frac{1}{\left|{\mathcal{V}}_d\right|}\times$$

$${\left|\sum _{j=1}^N\tilde{m}\left(\left|\mathbf{k}\right|;R_j\right)e^{-i\mathbf{k}\cdot {\mathbf{r}}_j}-\varphi {\int }_{{\mathcal{V}}_d}d{\mathbf{r}}^{\prime }{e}^{-i\mathbf{k}\cdot {\mathbf{r}}^{\prime }}\right|}^2,$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)=\frac{1}{\left|{\mathcal{V}}_d\right|}{\left|\sum _{j=1}^N\tilde{m}\left(\left|\mathbf{k}\right|;R_j\right)e^{-i\mathbf{k}\cdot {\mathbf{r}}_j}\right|}^2,$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)\sim {\left|\mathbf{k}\right|}^{\alpha }(\left|\mathbf{k}\right|\to 0).$$

$${\sigma }_V^2\left(R\right)\sim \{\begin{array}{ll}{R}^{-(d+1)},& \alpha >1(\mathrm{Class}\mathrm{I})\\ R^{-(d+1)}\ln R,& \alpha =1(\mathrm{Class}\mathrm{II})\\ R^{-(d+\alpha )},& \alpha <1(\mathrm{Class}\mathrm{III})\end{array},$$

$$N$$

$$\left|{\mathcal{V}}_d\right|$$

$$\varphi \le {\varphi }_{\max }\equiv \underset{j=1}{\overset{N}{\min }}\left\{\frac{v_1\left(R_j^{\max }\right)}{\left|C_j\right|}\right\},$$

$${\tilde{\chi }}_{{}_V}^{(m)}\left(\mathbf{k}\right)\le$$

$$+2\left|{\mathcal{V}}_d\right|{\varphi }^2(1-{\eta }_m{)}^2,$$

Full text

††journal: Acta Materialia

url]http://chemlabs.princeton.edu/torquato/

1 Introduction

Disordered hyperuniform materials [1, 2, 3] are exotic amorphous states of matter that are like crystals in the manner in which their large-scale density fluctuations are anomalously suppressed and yet behave like liquids or glasses in that they are statistically isotropic without any Bragg peaks. Recent results offer glimpses into the remarkable physical properties that such unusual correlated disordered materials can possess, including complete isotropic photonic/phononic band gaps [4, 5, 6], nearly optimal transport and mechanical properties [7, 8, 9], and dense but transparent materials [10].

The hyperuniformity concept was first introduced in the context of point patterns [1] and later extended to two-phase systems [2, 3], such as composites, colloidal suspensions, and polymer blends [11, 12, 13, 14, 15]. A hyperuniform two-phase material is one in which the local volume-fraction variance $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}$ inside a spherical observation window of radius $R$ decays faster than $R^{-d}$ for large $R$, which is the scaling for typical disordered systems. In the strongest form of hyperuniformity, called “class I,” the large-$R$ scaling of the variance is $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}\sim R^{-(d+1)}$ [2, 3]. Equivalently, its spectral density function $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}$, obtainable from scattering experiments [16], vanishes as the wavenumber $\left|\mathbf{k}\right|$ tends to zero [2, 3], and hence hyperuniform systems encompass all periodic and special disordered systems.

To date, a variety of disordered hyperuniform systems have been identified, including classical equilibrium systems [17, 18, 19, 20, 21, 22], quantum systems [23, 24, 25], maximally random jammed packings [26, 27, 28], non-equilibrium critical states [29, 30], non-equilibrium dynamical systems [31], random speckle patterns [32], number theory [20, 33, 34, 35], and biological systems [36, 37, 38, 39]; see also a recent review [3] and references therein. While some of these systems [17, 18, 19, 20, 21, 22, 23, 24, 25, 34, 35] are proved to be perfectly hyperuniform in the infinite-sample-size limit, others are effectively hyperuniform, i.e., $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(0\right)}}$ is not exactly zero but small compared to the peak value of the spectral density [3, 40]. Although effective hyperuniformity is sufficient to exhibit desirable physical properties, it has been shown that driving systems into perfect hyperuniformity will further improve their performance [41, 42].

Importantly, hyperuniformity is a property of infinitely large systems, which implies that one needs an arbitrarily large sample to ascertain whether the system is perfectly hyperuniform. However, current experimental [27, 30, 43, 44] and computational [22, 40, 45, 46, 47, 48] methods are limited in their capacity to create perfectly hyperuniform materials – increasing sample size comes at the expense of perfect hyperuniformity leading to effective hyperuniformity at best. For instance, while the collective-coordinate optimization technique [7, 8, 17] ensures the generation of perfectly hyperuniform disordered many-particle systems, even at finite wavelengths, its computational cost grows rapidly with sample size. Furthermore, hyperuniformity can be degraded or destroyed (even if by a small amount) due to the inevitable presence of imperfections in otherwise perfectly hyperuniform systems [40, 49].

Hence, there is a great need to devise systematic procedures to construct extremely large realizations of perfectly hyperuniform disordered two-phase systems, which would ensure the best property performance. To achieve this goal, we introduce a new tessellation-based procedure that ensures the creation of perfectly hyperuniform dispersions in two and three dimensions for any sample size, including the infinite-sample-size limit.

We begin by stating preliminary definitions and computational methods (see Sec. 2). Then, we introduce a tessellation-based procedure to construct disordered hyperuniform dispersions, i.e., two-phase systems in which hard (nonoverlapping) “particles” are spatially distributed throughout a connected “matrix” phase. This procedure is a new and special type of packing protocol that is different from previously known techniques [50]. We first describe this procedure and then prove that such constructions produce perfectly hyperuniform dispersions in the infinite-sample-size limit (see Sec. 3). As a proof-of-concept, we numerically verify our analytical results by constructing hyperuniform dispersions from two types of initial tessellations: Voronoi tessellations (see Sec. 4) and multiscale-sphere tessellations (see Sec. 5). The former case provides a remarkable mapping that converts a very large nonhyperuniform packing (at least as large as $10^{7}$ particles) to a hyperuniform dispersion. The latter case enables us to establish for the first time the hyperuniformity of the famous optimal multiscale coated-spheres model. The challenging numerical task of constructing such disordered multiscale dispersions is carried out and their hyperuniformity is verified. Finally, we demonstrate how such a capability can be combined with state-of-the-art 2D photolithographic [51] and 3D printing techniques [52, 53, 54] to design and fabricate very large disordered hyperuniform materials (see 6). We make concluding remarks in Sec. 7.

2 Preliminaries and computational methods

2.1 Spectral density

The microstructure of a two-phase system is often described by two-point correlation function $\mathinner{S_{2}^{(i)}\mathopen{\left(\mathbf{r}\right)}}$ that measures the probability that two points separated by $\mathbf{r}$ are simultaneously located in phase $i$. The autocovariance function is defined as

[TABLE]

which is identical for each phase, and tends to zero as $r$ increases when the system is in the absence of long-range order [13]. Its Fourier transform $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}$ at a wavevector $\mathbf{k}$ is called the spectral density and is a nonnegative real-valued function of $\mathbf{k}$. Importantly, the spectral density is directly measurable from elastic scattering experiments [16].

Consider a sphere-packing in a periodic simulation box $\mathcal{V}_{d}$ in $d$-dimensional Euclidean space $\mathbb{R}^{d}$, which consists of $N$ spheres of different radii $R_{1},R_{2},\cdots,R_{N}$. Its spectral density can be written as [55, 13, 56]

[TABLE]

where $\phi$ is the packing fraction and $\left|\mathcal{V}_{d}\right|$ is the volume of a simulation box. Here, $\mathinner{\tilde{m}\mathopen{\left(k;R\right)}}=\left(2\pi R/k\right)^{d/2}\mathinner{J_{d/2}\mathopen{\left(kR\right)}}$ represents the Fourier transform of a sphere of radius $R$, $\mathinner{J_{\nu}\mathopen{\left(x\right)}}$ is the Bessel function of order $\nu$, $L$ is side length of the unit cell, and a wavevector $\mathbf{k}\equiv\frac{2\pi}{L}(n_{1},n_{2},\cdots,n_{d})$ corresponds to the reciprocal lattice vectors of the unit cell.

In numerical simulations, Eq. (2) can be simplified as follows:

[TABLE]

using the fact that the space integral in Eq. (2) effectively vanishes at every nonzero wavevector. We note that Eqs. (2) and (3) cannot be applied to overlapping spherical particles. For an ensemble of multiple packings, the spectral density is computed by averaging spectral densities of individual packings, and the associated errors are estimated from the sample standard deviations. The angular-averaged spectral density, of central concern in this paper, is computed by averaging the spectral densities at wavevectors whose magnitudes are in the same bin, and error bars in the wavenumbers represent the standard deviations.

2.2 Classes of hyperuniform two-phase systems

For disordered hyperuniform two-phase systems, the spectral density $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}$ frequently exhibits power-law scaling behavior in the small wavenumber limit:

[TABLE]

Values of the positive exponent $\alpha$ define three classes of hyperuniformity associated with the large-$R$ behaviors of local volume-fraction variance [2, 3, 57]:

[TABLE]

where class I is the strongest form of hyperuniformity.

2.3 Random sequential addition (RSA) procedure

RSA is a time-dependent process that irreversibly, randomly, and sequentially places nonoverlapping spheres into space [58, 13]. In its infinite-time limit, the resulting packing becomes saturated (i.e., there is no available space to add other particles) and yet nonhyperuniform [59]. We generate exactly saturated RSA packings in a finite amount of time by using the voxel-list algorithm, developed by Zhang and Torquato [59].

2.4 Voronoi tessellations

For a point pattern, the Voronoi cell associated with a given point is the region of space closer to this point than to any other points, and the Voronoi tessellation is the collection of all Voronoi cells [13]. While Voronoi tessellations can be naturally generalized to Laguerre [60, 61] and Manhattan [62] distances, in this work we focus on Euclidean distances.

We compute the Voronoi tessellation of a given point pattern by the VORO++ library [63]. To improve the computational time, we divide a point pattern into several $d$-dimensional hypercubic domains with overlapping marginal regions whose thickness is $6a$, where $a\equiv\rho^{-1/d}$; see Fig.1 in Supplementary Material. Then, Voronoi tessellations of individual domains are computed in parallel. Here, we note that in 2D implementations, the particle number in each domain should be smaller than $10^{5}$ to avoid any possible memory leakage, which may occur because the VORO++ library is designed for 3D Voronoi tessellations.

2.5 Simulated sphere tessellations

A sphere tessellation is a partition of space ($\mathbb{R}^{d}$) by nonoverlapping spheres with a polydispersity in size down to the infinitesimally small. Such multiscale-sphere tessellations are infinitely degenerate with varying degrees of order/disorder. The most ordered ones would be derived from lattice packings or certain deterministic procedures, such as Apollonian gaskets. In this work, we focus on disordered sphere tessellations that are constructed by implementing a multicomponent version of the aforementioned RSA packing procedure, which is a multi-stage process.

Specifically, this multi-stage RSA process has a control parameter and a control function, namely, an upper bound $v_{\max}$ on cell volumes and a positive decreasing function $\mathinner{g\mathopen{\left(i\right)}}$ for natural numbers $i$, where $\mathinner{g\mathopen{\left(1\right)}}=1$ and its infinite sum exists (i.e., $\sum_{i=1}^{\infty}\mathinner{g\mathopen{\left(i\right)}}<\infty$). Using these factors, we determine the prescribed number $N$ of spheres that are inserted in each stage, and the largest cell volume $v^{(1)}$ to fill all space:

[TABLE]

where $\lceil x\rceil$ is the ceiling function and $\left|\mathcal{V}_{d}\right|$ is the volume of the simulation box. In the $m$th stage, the sphere has volume $v^{(m)}=v^{(1)}\mathinner{g\mathopen{\left(m\right)}}$ and diameter $D_{m}$, and the prescribed covering fraction is $N\sum_{i=1}^{m}v^{(i)}/\left|\mathcal{V}_{d}\right|$.

Using these parameters, sphere tessellations are constructed by the following steps:

Begin with an empty hypercubic simulation box in $\mathbb{R}^{d}$, whose volume is $\left|\mathcal{V}_{d}\right|$, under the periodic boundary conditions. 2. 2.

In the first stage ($m=1$), $N$ nonoverlapping spheres of an identical volume $v^{(1)}$ (i.e., diameter $D_{1}$) are irreversibly, sequentially, and randomly added in the simulation box (i.e., RSA procedure). The insertion in this stage stops only when the packing reaches to a prescribed covering fraction $Nv^{(1)}/\left|\mathcal{V}_{d}\right|$, unless the packing becomes saturated. At the end of the first stage, $N_{1}$($\leq N$) spheres are added and the associated covering fraction is $\eta_{1}\equiv N_{1}v^{(1)}/\left|\mathcal{V}_{d}\right|$. 3. 3.

In the $m$th stage ($m>1$), nonoverlapping spheres of an identical volume $v^{(m)}(<v^{(m-1)})$ are inserted via the RSA procedure into the packing generated up to the $(m-1)$st stage. Here, the spheres in the $m$the stage also should be nonoverlapping to those in the previous stages. The insertion stops when the packing reaches to a prescribed covering fraction $N\sum_{i=1}^{m}v^{(i)}/\left|\mathcal{V}_{d}\right|$, unless the packing becomes saturated. At the end of this stage, $N_{m}$ spheres are additionally inserted, and thus the packing has $\mathcal{N}_{m}(\equiv\sum_{i=1}^{m}N_{i})$ spheres in total and a covering fraction $\eta_{m}=(\sum_{i=1}^{m}N_{i}v^{(i)}/\left|\mathcal{V}_{d}\right|)$. 4. 4.

This procedure is repeated until it reaches to a prescribed number of stages.

In the RSA procedure of steps 2 and 3, we implement the aforementioned voxel-list algorithm [59]. In this work, we consider a power-law scaling $\mathinner{g\mathopen{\left(m\right)}}=1/m^{p}$ whose infinite series converges to the Riemann zeta function $\mathinner{\zeta\mathopen{\left(p\right)}}$, and thus schedule covering fraction in the $m$th stage is $\sum_{n=1}^{m}n^{-p}/\mathinner{\zeta\mathopen{\left(p\right)}}$. We choose scaling exponents $1<p<2$ because at scaling exponent higher than 2 makes it more difficult to insert additional particles as the number $m$ of stages increases [64]; see the Supplementary Material for details of the simulation parameters employed.

3 Theoretical Analyses

In this section, we precisely describe the tessellation-based procedure. We analytically prove that the constructed dispersions are strongly hyperuniform by deriving the small-$\left|\mathbf{k}\right|$ scaling of the spectral densities for spherical particles and subsequently provide an intuitive rationale of the theoretical results.

3.1 Tessellation-based procedure

Consider a hypercubic simulation box of side length $L$ in $d$-dimensional Euclidean space $\mathbb{R}^{d}$ under the periodic boundary conditions. Our procedure, graphically illustrated in Fig. 1, consists of the following steps:

Divide the simulation box into $N$ disjoint cells whose maximal lengths must be smaller than a given length scale $D$, which also should be much smaller than $L$; see Fig. 1(a). 2. 2.

Fill the $j$th cell $C_{j}$ with a single particle or multiple particles of general shape of total volume $\phi\left|C_{j}\right|$; see Fig. 1(b) for the simple instance of a single particle per cell. Then, repeat the same process over all cells. The global packing fraction of the constructed dispersion is identical to the local-cell packing fraction $\phi$.

We call the restriction on the cell size in step 1 the bounded-cell condition. We note that given an initial tessellation, our construction is realizable only when the local-cell packing fraction $\phi$ in step 2 is smaller than or equal to the maximal packing fraction $\phi_{\max}$, i.e.,

[TABLE]

where $\left|C_{j}\right|$ and $R^{\mathrm{max}}_{j}$ represent volume of the $j$th cell and radius of the largest sphere inscribed in this cell, respectively. Roughly speaking, the maximal packing fraction $\phi_{\max}$ becomes larger when the cell shapes become more isotropic (sphere-like).

3.2 General theoretical analyses

Here, we sketch the proof that the constructed dispersions are hyperuniform, the details of which are provided in Ref. [65]. We start from a general expression (2) for the spectral density of a sphere packing of various particle radii [55, 13, 56]. We now decompose $\phi\int_{{}_{\mathcal{V}_{d}}}\mathinner{d\mathbf{r}^{\prime}}e^{-i\mathbf{k}\cdot\mathbf{r}^{\prime}}$ into the Fourier transform of each cell, and rearrange each term with that of the associated particle. Importantly, the bounded-cell condition guarantees that the Fourier transforms of both a cell $C_{j}$ and the associated particle $j$ can be well approximated by their Taylor series about $\mathbf{k}=\mathbf{0}$. It follows that the tessellation-based procedure ensures that the leading order terms in these series exactly cancel one another such that the remaining terms exhibit power-law scalings in the wavenumber, i.e., $\left|\mathbf{k}\right|^{\alpha}$ as $\left|\mathbf{k}\right|\to 0$. A detailed analysis given in Ref. [65] shows that in the small-$\left|\mathbf{k}\right|$ limit, $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}\sim\phi^{2}\left|\mathbf{k}\right|^{2}$ is achievable, in general. However, the behavior $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}\sim\phi^{2}\left|\mathbf{k}\right|^{4}$ is achieved when the special condition $\mathbf{k}\cdot\left[\sum_{j=1}^{N}\Delta\mathbf{r}_{j}\mathinner{v_{1}\mathopen{\left(R_{j}\right)}}e^{-i\mathbf{k}\cdot\mathbf{r}_{j}}\right]=0$ is satisfied, where $\Delta\mathbf{r}_{j}$ represents the displacement between the centroids of cell $j$ and the associated particle. These results clearly show that the constructed dispersions are hyperuniform of class I for $\phi\leq\phi_{\max}$, given in Eq. (7), which also guarantees that the particles remain inside the associated cells. Note that $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}$ is exactly the same for either the particle or matrix (space exterior to particles) phases [13] and hence the matrix phase is also strongly hyperuniform.

To understand intuitively why a dispersion constructed by this procedure is hyperuniform, it is useful to consider the local volume-fraction variance $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}$. Imagine sampling the constructed dispersion with many randomly placed observation windows of radius $R(\gg D)$. Clearly, the volume-fraction fluctuations will be concentrated only in the small region near the boundary of the window, as shown in Fig. 1(c). Consequently, the variance in the phase-volume will be proportional to the window surface area in the large-$R$ limit, i.e., ${\mathinner{v_{1}\mathopen{\left(R\right)}}}^{2}\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}\sim R^{d-1}$, and in turn $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}\sim R^{-(d+1)}$, implying that such dispersions are strongly hyperuniform (class I). Our methodology is valid for any particle shape as shown in Fig. 1(d) and allows for the addition of multiple particles per cell provided that local-cell packing fraction is identical for each cell. However, for concreteness, we mainly analyze dispersions that have only a single spherical particle in each cell.

4 Disordered hyperuniform dispersions derived from Voronoi tessellations

As a proof-of-concept, we now implement the tessellation-based procedure numerically via Voronoi tessellations (see Sec. 2.4) of disordered and “nonhyperuniform” sphere packings; see Fig. 2. Thus, performing our methodology from Voronoi tessellations provides an efficient mapping that remarkably converts very large ($N\sim 10^{7}$) nonhyperuniform progenitor point patterns to perfectly hyperuniform dispersions.

To properly implement our methodology, the Voronoi tessellations of the progenitor point patterns should obey the bounded-cell condition of step 1. Surprisingly, when sample sizes are sufficiently large, this condition is satisfied by Voronoi tessellations of virtually all statistically homogeneous point patterns, including Poisson point patterns. This seems counterintuitive because in the infinite-sample-size limit Poisson point patterns as well as many other disordered ones can possess arbitrarily large holes, or equivalently, arbitrarily large Voronoi cells. Nonetheless, a detailed analysis given in Ref. [65] definitely demonstrates that practically all samples of these systems meet the bounded-cell condition as the sample size grows. To sketch the basic idea, we note that for Poisson point patterns of $500,000$ points, the probability that a single finite-size sample possesses a hole larger than $0.1\%$ of its sample size is around $10^{-22}$, and this probability becomes exponentially lower as sample size grows. Interestingly, while such rare events of large hole formation play a crucial role in many physical phenomena, such as the diffusion-controlled reactions [66], the bounded-cell condition is not governed by these rare events. This has the crucially practical implication that any statistically homogeneous point pattern, whether hyperuniform or not, can be used to numerically implement our mapping, which is consistent with a recent study of random fields generated from Voronoi tessellations for the special case of a Poisson point process [67].

In this paper, we choose to perform the tessellation-based procedure from Voronoi tessellations of saturated RSA packings; see Fig. 2(a) and Sec. 2.4. In contrast to uncorrelated Poisson point patterns, saturated RSA packings unquestionably meet the bounded-cell condition for even relatively small sample sizes (say, $10^{d+1}$ particles in $d$ dimensions). From saturated RSA packings via the voxel-list algorithm [59], for simplicity we construct dispersions (Fig. 2(b)) by solely scaling particle sizes without changing their positions, and hence their maximal packing fractions are defined by $\phi_{\max}^{(1)}$, illustrated in Fig. 2(c).

Figure 3 summarizes simulation results of our methodology from Voronoi tessellations of saturated RSA packings. We have constructed large saturated RSA packings with $N\approx 10^{7}$ and $10^{6}$ particles in $d=2$ and $3$, respectively. For these progenitor packings, the maximal packing fractions $\phi_{\max}^{(1)}$ are computed from Eq. (7) by replacing $R_{j}^{\max}$ with a half of the nearest-neighbor distance of particle $j$. The corresponding values of $\phi_{\max}^{(1)}$ are $0.360(4)$ and $0.252(7)$ in $d=2$ and $d=3$, respectively. Details of the simulation parameters employed, computational times, and machine type are summarized in the Supplementary Material. By construction, our dispersions always have a particle-size distribution that is identical to the cell-size distribution (Figs. 3(a) and 3(d)), which can be well approximated by the Gamma or the log-normal distributions. Thus, from the definition of a saturated RSA packing, we can estimate bounds on the particle volumes at the local-cell packing fraction $\phi$ as follows: $\phi_{\mathrm{sat}}<\frac{\mathinner{v_{1}\mathopen{\left(R_{j}\right)}}}{\phi\overline{\left|C\right|}}<2^{d}\phi_{\mathrm{sat}}$, where $\phi_{\mathrm{sat}}$ represents the saturated packing fraction of RSA in $\mathbb{R}^{d}$ and $\overline{\left|C\right|}$ is the mean cell-volume. In fact, the largest particle volumes are about three times as large as the smallest ones in both $d=2$ and $3$. Finally, the spectral density exhibits the behavior $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\left|\mathbf{k}\right|\right)}}\sim\phi^{2}\left|\mathbf{k}\right|^{4}$ as $\left|\mathbf{k}\right|$ goes to zero (Figs. 3(b)-3(c) and 3(e)-3(f)). Thus, the constructed dispersions in both $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$ are strongly hyperuniform (class I) and their deviation term $\sum_{j=1}^{N}\Delta\mathbf{r}_{j}\mathinner{v_{1}\mathopen{\left(R_{j}\right)}}$ effectively vanishes. Finally, it is noteworthy that the spectral densities of both progenitor and constructed dispersions collapse onto a single curve for intermediate wavenumbers; see insets in Fig. 3(b,e). This behavior is expected, since both dispersions have identical local statistics (see Fig. 3(a) and (d)).

5 Hyperuniformity of the optimal coated-sphere structures

It is noteworthy that our tessellation-based methodology can be applied to other tessellations besides Voronoi tessellations that consist of polyhedral cells. For example, although space cannot be tessellated by identical nonoverlapping spheres, it can be tessellated by nonoverlapping spheres with a polydispersity in size down to the infinitesimally small; see, for example, the multiscale-sphere tessellation in Fig. 5(a) and Sec. 2.5. Application of our tessellation-based procedure to multiscale-sphere tessellations should produce hyperuniform dispersions, regardless of the particle shapes and cell positions with varying degrees of order/disorder. We now show that the famous Hashin-Shtrikman multiscale coated-spheres structures from the theory of heterogeneous materials are a special case derived from such sphere tessellations. This broad class of particulate composites are optimal for the effective thermal (electrical) conductivity and elastic moduli for prescribed phase properties and volume fractions [68, 13, 69]. Such an optimal two-phase material consists of composite spheres that are composed of a spherical core of one phase (dispersed phase) that is surrounded by a concentric spherical shell of the other phase such that the fraction of space occupied by the core phase is equal to its overall phase volume fraction. The composite spheres fill all space, implying that their sizes range down to the infinitesimally small. Thus, this structure is a dispersion, i.e., the matrix (red region) is a fully connected (continuous) phase and the inclusions (blue regions) are “well-separated” from one another [70, 71]; see Fig. 5(b). Here, we, for the first time, prove that these Hashin-Shtrikman multiscale optimal dispersions are in fact hyperuniform, which apparently is related to their optimal transport and elastic properties.

It is crucial to observe that in the coated-spheres model, the composite spheres comprise a multiscale-sphere tessellation (see Fig. 5(a)) and the fraction of space occupied by the core phase is equal to its global phase volume fraction. In other words, these optimal structures can be constructed via our methodology from multiscale-sphere tessellations. Assuming the bounded-cell condition, this implies that the coated-spheres model should be hyperuniform, whether disordered or not. Further analysis shows that for the coated-spheres model of packing fraction $\phi$ in $\mathbb{R}^{d}$, the small-$\left|\mathbf{k}\right|$ scaling of the associated spectral density will exhibit the following behavior: $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\left|\mathbf{k}\right|\right)}}\sim\left(\phi(1-\phi^{2/d})\right)^{2}\left|\mathbf{k}\right|^{4}$, implying class I; see Ref. [65] for details.

We carry out the challenging numerical task of constructing such disordered multiscale dispersions and verify their hyperuniformity. To do so, first we implement a multicomponent version of the aforementioned RSA packing procedure, which amounts to a multi-stage process; see Sec. 2.5. Then, we simulate the coated-spheres model of the inclusion volume fraction $\phi$ by reducing sphere volumes in a precursor packing at a volume ratio of $\phi$ without moving their centers. In a finite $m$th stage, the precursor packing will not cover all space and hence there will be gaps in which smaller spheres can be added in next stages. Thus, as the number of stages increases in the limit of $m\to\infty$, those gaps are eventually covered by spheres of size down to the infinitesimally small, i.e., $\lim_{m\to\infty}\eta_{m}=1$. To achieve a nearly complete tiling of space (say, $\eta_{m}\approx 0.95$) in $d=2$, this process requires around a few hundred stages; see Fig. 5(c).

To estimate the degree of hyperuniformity in the $m$th stage, we employ an upper bound of the spectral density $\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(\mathbf{k}\right)}}$ with any prior knowledge of cell-volume distributions [65]:

[TABLE]

where $\phi$ stands for the local-cell packing fraction, and $\mathcal{R}_{j}$ denotes radius of cell $j$. In short, this upper bound is obtained by assuming that complex Fourier components of all uncovered gaps constructively interfere, and thus the second term in Eq. (8) represents the largest possible volume-fraction fluctuations due to these gaps. This rigorous upper bound (8) indicates that in the limit of $\eta_{m}\to 1$, the spectral density is at least on the order of $k^{4}$ for small $\mathbf{k}$’s, implying that our coated-spheres construction becomes strongly hyperuniform in the limit of $m\to\infty$.

In what follows, we focus on numerical results in two-dimensional cases for simplicity. The simulations proceed up to the 400th stage with the following parameters. Particle volumes in the $m$th stage is determined by a power-law scaling $\mathinner{v_{1}\mathopen{\left(D_{m}/2\right)}}=\mathinner{v_{1}\mathopen{\left(D_{1}/2\right)}}m^{-p}$ for a scaling exponent $p$ that ranges from 1 to 2 concerning computational efficiency. An upper bound $v_{\max}$ on cell volumes is chosen to achieve $v_{1}(D_{1}/2)\lesssim 0.001\left|\mathcal{V}_{2}\right|$; see Sec. 2.5 and the Supplementary Material. Figure 4(a)-4(b) illustrate our multiscale-disk tiling model in the $400$th stage with $p=1.5$, and the resulting coated-disks model at the local-cell packing fraction $\phi=0.25$. Figure 5(a) summarizes that covering fraction of the $m$th stage of our model for various values of a scaling exponent $p$. This shows that in finite stages $m$, our model has the uncovered gaps that contribute to the long-wavelength fluctuations and thus becomes nearly hyperuniform, rather than perfectly hyperuniform. However, Fig. 5(c) clearly demonstrates as those gaps are filled ($\eta_{m}\to 1$) with smaller composite disks, the associated fluctuations $\mathinner{\tilde{\chi}_{V}^{(m)}\mathopen{\left(k\to 0\right)}}$ vanish on the order of $(1-\eta_{m})^{2}$. This scaling behavior is consistent with our upper bound (8) although our prediction is grossly overestimated. Thus, from both our theoretical and numerical results, we conclude that the ideal multiscale coated-spheres model with $\eta_{m}=1$ will be strongly hyperuniform.

6 Fabrication of the computational designs

Importantly, our designed hyperuniform structures obtained from Voronoi tessellations can be easily fabricated via either photolithographic or 3D-printing techniques. State-of-the-art photolithography fabrication methods are highly suitable for mass production, and capable of creating a 2D pattern on a wafer up to $30$ cm in diameter with $25$ nm in minimal feature size [51]. Since any sharp corners in the designed particle shapes are rounded up to the minimal feature size, it is easier to fabricate circular particles than polygonal ones. When the minimal feature size is around $1.5~{}\mu$m, one can easily fabricate our two-dimensional hyperuniform designs with more than one million particles. For three-dimensional designs, modern 3D-printing techniques [54] can be applied. Since a printed structure must be topologically connected phase that is mechanically self-supporting, one needs to print the connected matrix phase corresponding to our hyperuniform dispersions [72] as shown in right panels in Fig. 6. Due to recent advances in 3D-printing techniques, even commercial desktop 3D-printers can achieve around $100~{}\mu$m in XY-resolution, $20~{}\mu$m in Z-resolution [73], and can print a sample with dimensions $125\times 125\times 125$ cm3 in 50 hours. Such devices should be able to print our designs of around 50 million pores when the smallest pore diameter is set around $300\mu$m.

7 Conclusions

In summary, we provide an efficient procedure that is capable of constructing very large disordered dispersions that are perfectly hyperuniform. Unlike many previous methods that have been used to generate disordered hyperuniform materials, our procedure is simple to implement as it only involves constraining the local-cell packing fraction, which is independent of the rest of a system. Furthermore, all computations to determine particle volumes can be exactly performed and easily parallelized. Our methodology shares some similarities with the “equal-volume tessellation” to obtain hyperuniform point patterns [74] in which the local number density within each cell is identical. However, our procedure is more versatile in that it can be applied to Voronoi tessellations of disordered point patterns and multiscale-sphere tessellations that we discussed in this paper.

For simplicity, we have focused in this work on constructed dispersions of spherical particles and sketched the proof of their hyperuniformity. In Ref. [65], we provide detailed derivations of the spectral densities of hyperuniform dispersions consisting of particles of arbitrary shape and show how any statistically homogeneous progenitor packing satisfies the bounded-cell condition. There, we also present additional simulation results for a variety of progenitor sphere packings as well as packings of nonspherical particles. In the case of coated-spheres model, we provide additional mathematical details as well as numerical results using different cell-volume scalings [65].

In this work, we have shown that local statistics alone generally does not determine hyperuniformity, which although reinforces previous observations (see a recent review [3] and references therein), is not commonly understood. While the nonhyperuniform progenitor configurations have no influence on the behavior of the spectral density around the origin, since the constructed dispersion is hyperuniform of class I, it starts to determine its behavior away from the origin at intermediate to large wavenumbers. Indeed, at intermediate wavenumbers the spectral densities of the progenitor and constructed dispersions collapse onto one another [65] because local statistics, such as Voronoi cell volumes and nearest-neighbor distance distributions, are identical for both systems. This outcome is consistent with the fact that hyperuniformity is a global property of a system.

It should not go unnoticed that many other tessellations can be employed in our tessellation-based methodology, including generalizations of Voronoi tessellations, such as Laguerre or radical tessellations [61, 60] and Voronoi tessellations in Manhattan distance [62], Delaunay triangulations [13], “Delaunay-centroidal” tessellations [4, 70], and disordered isoradial graphs [75] whenever they meet the bounded-cell condition. Furthermore, our methodology allows one to tune particle shapes and their numbers within each cell while preserving hyperuniformity. This exploration of different tessellations and particle geometries/numbers in each cell enables one to generate an enormous class of hyperuniform dispersions and hence represents fertile ground for future research.

We have established for the first time that the optimal Hashin-Shtrikman multiscale dispersions are indeed hyperuniform. This finding suggests that hyperuniformity is evidently a crucial characteristic to achieve optimality with respect to effective transport and elastic properties [13, 68, 69]. Thus, not only do the spheres have to be “well-separated” from one another (as is traditionally understood) but the entire dispersion should possess the global property of hyperuniformity to achieve optimal effective physical properties, which heretofore had not been known. While it has begun to be shown that some disordered hyperuniform dispersions have desirable and nearly optimal transport, mechanical and optical properties [10, 7, 8], it will be of great interest to study the physical properties of the enormous class of hyperuniform dispersions that can be designed and tuned by our tessellation-based methodology. These computational designs can subsequently be combined with the aforementioned 2D and 3D fabrication techniques to accelerate the discovery of hyperuniform two-phase materials.

Acknowledgements

We thank G. Zhang, T. G. Mason, A. B. Hopkins, and M. A. Klatt for very helpful discussions. This work was partially supported by the U.S. National Science Foundation under Grant No. DMR-1714722.