Methodology to Construct Large Realizations of Perfectly Hyperuniform Disordered Packings

TL;DR

This paper introduces a general theoretical method to construct perfectly hyperuniform disordered packings in any dimension, overcoming previous limitations and enabling the design of materials with exotic physical properties.

Contribution

The authors develop a novel methodology to generate perfectly hyperuniform packings from arbitrary tessellations, independent of particle shape or placement, applicable in multiple dimensions.

Findings

Successfully converts large nonhyperuniform packings into hyperuniform ones in 2D and 3D.

Proves hyperuniformity of Hashin-Shtrikman coated-spheres structures.

Demonstrates tunability and potential for discovering new hyperuniform materials.

Abstract

Disordered hyperuniform packings are unusual amorphous states of two-phase materials that are endowed with exotic physical properties. Such hyperuniform systems are characterized by an anomalous suppression of volume-fraction fluctuations at infinitely long-wavelengths, compared to ordinary disordered materials. While there has been growing interest in such singular states of amorphous matter, a major obstacle has been an inability to produce large hyperuniform samples due to practical limitations of conventional methods. To overcome these limitations, we introduce a general theoretical methodology to construct perfectly hyperuniform packings in $d$-dimensional Euclidean space. Specifically, beginning with an initial general tessellation of space that meets a "bounded-cell" condition, hard particles are placed inside each cell such that the packing fraction within the cell becomes identical to the global packing fraction. We prove that the constructed packings are perfectly hyperuniform in the infinite-sample-size limit and their hyperuniformity is independent of particle shapes, positions, and numbers per cell. We employ two distinct types of initial tessellations: Voronoi as well as sphere tessellations. Beginning with Voronoi tessellations, we show that our algorithm can remarkably convert extremely large nonhyperuniform packings into hyperuniform ones in $\mathbb{R}^2$ and $\mathbb{R}^3$. Implementing our theoretical methodology on sphere tessellations, we establish the hyperuniformity of the classical Hashin-Shtrikman multiscale coated-spheres structures that possess optimal effective transport and elastic properties. In addition, the tunability of our methodology offers promise for discovery of novel disordered hyperuniform two-phase materials.

Figures

Tables

Table 1. Table 1: Extrapolated values of χ ~ V ( 0 ) \mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(0\right)}} and hyperuniform metric H 𝐻 H of the constructed packings. While values in the second column are initial packing fractions ϕ init subscript italic-ϕ init \phi_{\mathrm{init}} for RSA and equilibrium hard-sphere liquids (HSL), those represent the vacancy concentration c 𝑐 c for imperfect ℤ d superscript ℤ 𝑑 \mathbb{Z}^{d} . For the progenitor packings, H ∼ 0.1 − 0.01 similar-to 𝐻 0.1 0.01 H\sim 0.1-0.01 . We do not compute H 𝐻 H for the imperfect lattices because of their Bragg peaks. See the Supplemental Material Sup for additional data.

Progenitors/Parameters		${\tilde{\chi }}_{{}_V}\left(0\right)$		$H$
2D RSA $N={10}^7$	${\varphi }_{\mathrm{sat}}\approx 0.5471$	2.29(54)	$\times {10}^{-18}$	1.36(32)	$\times {10}^{-14}$
	0.41025	4.7(15)	$\times {10}^{-18}$	3.4(11)	$\times {10}^{-14}$
	0.27350	1.15(60)	$\times {10}^{-17}$	1.01(52)	$\times {10}^{-13}$
	0.13675	3.1(18)	$\times {10}^{-17}$	3.0(17)	$\times {10}^{-13}$
3D RSA $N={10}^6$	${\varphi }_{\mathrm{sat}}\approx 0.3841$	1.97(31)	$\times {10}^{-12}$	1.51(24)	$\times {10}^{-8}$
	0.288	1.39(31)	$\times {10}^{-12}$	1.27(28)	$\times {10}^{-8}$
	0.192	3.77(85)	$\times {10}^{-12}$	3.90(88)	$\times {10}^{-8}$
	0.096	1.44(32)	$\times {10}^{-11}$	1.66(36)	$\times {10}^{-7}$
2D HSL $N={10}^5$	0.65	2.03(32)	$\times {10}^{-14}$	6.03(98)	$\times {10}^{-11}$
	0.40	6.7(19)	$\times {10}^{-14}$	4.5(13)	$\times {10}^{-10}$
	0.20	2.48(76)	$\times {10}^{-13}$	2.20(67)	$\times {10}^{-9}$
3D HSL $N={10}^6$	0.45	-1.4(32)	$\times {10}^{-13}$	7.0(162)	$\times {10}^{-10}$
	0.30	9.2(47)	$\times {10}^{-13}$	7.6(39)	$\times {10}^{-9}$
	0.20	7.7(32)	$\times {10}^{-13}$	8.0(33)	$\times {10}^{-9}$
Imperfect ${\mathbb{Z}}^2$ $N_s={10}^8$	0.10	6.9(50)	$\times {10}^{-20}$	-
	0.20	2.1(15)	$\times {10}^{-19}$	-
	0.40	2.49(95)	$\times {10}^{-18}$	-
Imperfect ${\mathbb{Z}}^3$ $N_s=27\times {10}^6$	0.10	4.1(71)	$\times {10}^{-15}$	-
	0.20	5.0(20)	$\times {10}^{-14}$	-
	0.40	2.42(76)	$\times {10}^{-13}$	-

Equations

$$\underset{R\to \infty }{\lim }{v}_1\left(R\right){\sigma }_V^2\left(R\right)=0.$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)\sim {\left|\mathbf{k}\right|}^{\alpha }(\alpha >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},$$

$${\mathcal{I}}^{(i)}\left(\mathbf{r}\right)=\{\begin{array}{ll}1,& \mathbf{r}\in \mathrm{phase}i\\ 0,& \mathrm{otherwise}.\end{array}$$

$${\chi }_{{}_V}\left(\mathbf{r}\right)\equiv ⟨{\mathcal{J}}^{(i)}\left({\mathbf{r}}^{\prime }\right){\mathcal{J}}^{(i)}\left({\mathbf{r}}^{\prime }+\mathbf{r}\right)⟩,$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)=\frac{1}{\left|{\mathcal{V}}_F\right|}⟨{\left|{\tilde{\mathcal{J}}}^{(i)}\left(\mathbf{k}\right)\right|}^2⟩,$$

$$\tilde{\mathcal{J}}\left(\mathbf{k}\right)$$

$$=\sum _{j=1}^N\tilde{m}\left(\mathbf{k};{\mathbf{P}}_j\right)e^{-i\mathbf{k}\cdot {\mathbf{r}}_j}-\varphi \left|{\mathcal{V}}_F\right|{\delta }_{\mathbf{k},0},$$

$$\tilde{m}\left(\mathbf{k};a\right)={\left(2\pi a/k\right)}^{d/2}{J}_{d/2}\left(ka\right),$$

$$\tilde{m}\left(\mathbf{k};L\right)=L^d\prod _{l=1}^d\mathrm{sinc}\left(k_{l}L/2\right),$$

$$\mathrm{sinc}\left(x\right)\equiv \{\begin{array}{ll}\frac{\sin x}{x},& x\ne 0\\ 1,& x=0.\end{array}$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)=\rho {\left|\tilde{m}\left(\mathbf{k};a\right)\right|}^{2}S\left(\mathbf{k}\right),$$

$$S\left(\mathbf{k}\right)\equiv \frac{1}{N}⟨{\left|\sum _{j=1}^N{e}^{-i\mathbf{k}\cdot {\mathbf{r}}_j}-N{\delta }_{\mathbf{k},0}\right|}^2⟩.$$

$${\tilde{\chi }}_{{}_V}\left(0\right)={\varphi }^{2}S\left(0\right)/\rho ,$$

$${\tilde{\chi }}_{{}_V}\left(0\right)={\varphi }^2{\kappa }_T{k}_{B}T,$$

$${\sigma }_V^2\left(R\right)\equiv ⟨{\tau }^2\left(\mathbf{x};R\right)⟩-{\varphi }^2,$$

$$H\equiv \frac{{\tilde{\chi }}_{{}_V}\left(\mathbf{k}\to 0\right)}{{\tilde{\chi }}_{{}_V}\left({\mathbf{k}}_{\mathrm{peak}}\right)},$$

$${\ell }_{\max }\ll L.$$

$$\varphi \le {\varphi }_{\max }\equiv \underset{j=1}{\overset{N}{\min }}\left\{\frac{{\left|{\mathbf{P}}_j\right|}_{\max }}{\left|{\mathbf{C}}_j\right|}\right\},$$

$$\left|{\mathcal{V}}_F\right|{\delta }_{\mathbf{k},0}$$

$$=\sum _{j=1}^N{e}^{-i\mathbf{k}\cdot {\mathbf{x}}_j}\tilde{m}\left(\mathbf{k};{\mathbf{C}}_j\right),$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)=\frac{1}{\left|{\mathcal{V}}_F\right|}{\left|\tilde{\mathcal{J}}\left(\mathbf{k}\right)\right|}^2$$

$$=$$

$$\tilde{m}\left(\mathbf{k};\mathcal{V}\right)=$$

$$\displaystyle+\frac{ik_{\alpha}k_{\beta}k_{\gamma}}{6}\mathinner{\mathcal{M}_{\alpha\beta\gamma}\mathopen{\left(\mathcal{V}\right)}}\Bigg{]}+\mathinner{\mathcal{O}\mathopen{\left(k^{4}\right)}},$$

$${\mathcal{M}}_{{\alpha }_1{\alpha }_2\cdots {\alpha }_n}\left(\mathcal{V}\right)\equiv \frac{1}{\left|\mathcal{V}\right|}{\int }_{\mathcal{V}}d\mathbf{r}{r}_{{\alpha }_1}{r}_{{\alpha }_2}\cdots r_{{\alpha }_n},$$

$$\tilde{\mathcal{J}}\left(\mathbf{k}\right)$$

$${\tilde{\mathcal{J}}}_{(1)}\left(\mathbf{k}\right)=$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)\sim {\varphi }^2{\left|\mathbf{k}\right|}^2,$$

$${\tilde{\chi }}_{{}_V}\left(\mathbf{k}\right)\sim {\varphi }^2{\left|\mathbf{k}\right|}^4.$$

$$\tilde{\mathcal{J}}\left(\mathbf{k}\right)$$

$$+$$

$${\tilde{\chi }}_{{}_V}\left(k\right)\sim {\varphi }^2(1-{\varphi }^{2/d}{)}^2{\left|\mathbf{k}\right|}^4.$$

Full text

I Introduction

A hyperuniform state of matter is characterized by an anomalous suppression of density or volume-fraction fluctuations at infinitely long wavelengths relative to those in typical disordered systems, such as liquids and structural glasses Zachary and Torquato (2009); Torquato and Stillinger (2003); Torquato (2018a). Such hyperuniform states encompass all perfect crystals, many quasicrystals, as well as some exotic disordered systems. Disordered hyperuniform states of matter have been the subject of intense interest across a variety of fields, including physics Scardicchio et al. (2009); Harrison (1970); Pietronero et al. (2002); Lomba et al. (2017); Jancovici (1981); Jiao and Torquato (2011); Feynman and Cohen (1956); Hexner and Levine (2015); Weijs et al. (2015); Chremos and Douglas (2018), materials science Man et al. (2013); Florescu et al. (2013); Leseur et al. (2016); Froufe-Pérez et al. (2017); Le Thien et al. (2017); Chen and Torquato (2018); Zhang et al. (2016), chemistry Rundman and Hilliard (1967); Hashimoto et al. (1986); Ma and Torquato (2017), biology Noh et al. (2010); Jiao et al. (2014); Mayer et al. (2015); Kwon and Kim (2017), and mathematics Dyson (1962); Montgomery (1973); Torquato et al. (2018). The notion of hyperuniformity was first defined in the context of point-particle systems Torquato and Stillinger (2003) and then extended to two-phase heterogeneous systems Zachary and Torquato (2009) and random scalar/vector fields Torquato (2016). General two-phase systems abound in natural and artificial materials, including colloidal suspensions, particulate composites, and concrete Brinker and Scherer (2013); Weeks (2017); Torquato (2002); Patel and Zohdi (2016); Neville (1995); Sahimi (2003); Milton (2002); Mickel et al. (2013). Packings (or dispersions), which are of central concern in this paper, comprise a class of two-phase systems in which nonoverlapping particles are spatially distributed throughout a connected “void” (matrix) phase.

A hyperuniform two-phase system in $d$-dimensional Euclidean space $\mathbb{R}^{d}$ 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}$ in the large-$R$ limit Zachary and Torquato (2009); Torquato (2018a):

[TABLE]

Equivalently, its associated spectral density $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\bm{k}\right)}}$ vanishes as the wavenumber $\left|\bm{k}\right|$ tends to zero Zachary and Torquato (2009); Torquato (2018a). For disordered hyperuniform two-phase systems, $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\bm{k}\right)}}$ typically exhibits the power-law scaling:

[TABLE]

This exponent $\alpha$ is directly related to three distinct classes of hyperuniformity that are categorized based on the large-$R$ scalings of $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}$ Torquato (2018a):

[TABLE]

where classes I and III represent the strongest and weakest forms of hyperuniformity, respectively. Class I systems include all crystals, some quasicrystals, and stealthy hyperuniform systems in which $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}=0$ for $0<k<K$ Batten et al. (2008); Zhang et al. (2015a, b); Torquato et al. (2015).

Disordered hyperuniform two-phase systems Zachary and Torquato (2009); Zhang et al. (2016) are attracting considerable attention due to their unusual physical properties, such as complete isotropic photonic and phononic bandgaps Florescu et al. (2009); Froufe-Pérez et al. (2017); Gkantzounis et al. (2017); López (2018), nearly optimal transport properties Zhang et al. (2016); Chen and Torquato (2018); Le Thien et al. (2017), superior metamaterial designs Wu et al. (2017), dense but transparent materials Leseur et al. (2016), and low-density materials of blackbody-like absorption Bigourdan et al. (2019). Similar to ordinary two-phase systems, their effective properties are also tunable by engineering the phase properties and volume fractions as well as the spatial arrangements Florescu et al. (2009); Froufe-Pérez et al. (2017); Rechtsman and Torquato (2008); Chen and Torquato (2018); Torquato and Chen (2018a); Zhang et al. (2016); Wu et al. (2017). An important class of two-phase systems are the Hashin-Shtrikman structures Hashin (1962); Torquato (2002); Torquato and Pham (2004); Milton (1986, 2002) that are optimal for effective elastic moduli Hashin (1962); Torquato (2002); Milton (2002), thermal (electrical) conductivity Milton (1986), trappiong constant Torquato and Pham (2004), and fluid permeability Torquato and Pham (2004) for given phase volumes and phase properties. Remarkably, certain disordered hyperuniform systems possess nearly optimal transport and elastic properties Zhang et al. (2016); Chen and Torquato (2018); Torquato and Chen (2018a, b).

Theoretical Scardicchio et al. (2009); Hexner and Levine (2015); Tjhung and Berthier (2015); Lomba et al. (2017); Chen et al. (2018), numerical Zachary et al. (2011); Torquato and Jiao (2010); Lubachevsky and Stillinger (1990); Lubachevsky et al. (1991); Atkinson et al. (2016); Uche et al. (2004); Zhang et al. (2016); Chen and Torquato (2018); DiStasio et al. (2018); Wang et al. (2018), and experimental methods Weijs et al. (2015); Kurita and Weeks (2011); Ricouvier et al. (2017); Zito et al. (2015) have been developed to generate disordered hyperuniform packings (dispersions). In practice, however, these methods are system-size limited due to computational cost or imperfections. For instance, the collective-coordinate optimization technique Uche et al. (2004); Zhang et al. (2016); Chen and Torquato (2018); DiStasio et al. (2018) and equilibrium plasma Jancovici (1981); Lomba et al. (2017); Chen et al. (2018) can achieve perfect hyperuniformity, but their long-range interactions lead the computation cost to grow rapidly with system size. Random organization models Hexner and Levine (2015); Tjhung and Berthier (2015) yield disordered hyperuniform packings at critical absorbing states but are limited in producing perfect hyperuniformity because of critical slowing-down phenomena Atkinson et al. (2016); Torquato (2018a). Determinantal point processes are perfectly hyperuniform in the thermodynamic limit, but the current numerical algorithm hardly can generate a realization of more than 100 particles due to accumulated numerical error Scardicchio et al. (2009). Stealthy designs via the superposition procedure DiStasio et al. (2018) provides an efficient means to construct exactly stealthy hyperuniform digitized two-phase systems. However, this construction scheme requires one to prepare many different small systems as building blocks, which is computationally demanding as the system size increases.

Since hyperuniformity is a global property of an infinitely large system, limited sample sizes often make it difficult to ascertain whether such systems are truly hyperuniform or effectively hyperuniform. Furthermore, disordered hyperuniform systems often can include imperfections, such as point vacancies, stochastic displacements Gabrielli (2004), thermally excited phonon modes Kim and Torquato (2018), and rattlers Atkinson et al. (2016) in maximally random jammed (MRJ) packings. Such imperfections can either destroy or degrade the hyperuniformity (even if by a small amount) of otherwise perfectly hyperuniform systems Kim and Torquato (2018); Torquato (2018a). Hence, there is a great need to devise exact and efficient procedures to construct extremely large realizations of disordered hyperuniform two-phase systems. Such capability could then be combined with state-of-the-art photolithographic and 3D printing techniques to fabricate large disordered hyperuniform systems.

In this paper, we introduce a tessellation-based procedure that enables one to generate disordered packings in $\mathbb{R}^{d}$ that are perfectly hyperuniform in the infinite-sample-size limit. Based on this theoretical methodology, we formulate an efficient algorithm to generate extremely large realizations of hyperuniform packings. Specifically, one first tessellates the space into certain disjoint cells that meet a bounded-cell condition, i.e., the maximal length of each cell should be much shorter than side length of the entire tessellation. This mild restriction allows one to employ a wide class of tessellations, as discussed below and in Sec. V. Then, one places hard particles in each cell such that the local packing fractions associated with the cells become identical to the global packing fraction $\phi$; see Fig. 1. Note that for “periodic” tessellations of equal cell volumes (e.g., Voronoi tessellations of Bravais lattices and Weaire-Phelan foam Weaire and Phelan (1994); Conway and Torquato (2006)), our procedure yields periodic packings of identical particles that are also stealthy hyperuniform 111Generally speaking, application of our procedure to any periodic tessellation will result in a periodic packing that is stealthy hyperuniform.. However, our major concerns in this paper are disordered tessellations whose cells have a variability in sizes and shapes. Thus, the application of our procedure to such tessellations generates “disordered” hyperuniform packings in which the particles have a polydispersity in size. We present a detailed theoretical analysis of the small-wavenumber scalings of the spectral density $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}$ for the constructed packings in the case of arbitrary particle shapes. We thereby prove that whenever the initial tessellations meet the bounded-cell condition, these systems are strongly hyperuniform (class I) for any particle shape in the infinite-sample-size limit. In this limit, the system size tends to be infinitely large with other intensive parameters (e.g., number density and packing fractions) held fixed. This procedure is a packing protocol to generate packings of polydisperse particles that is uniquely different from previously known methods Torquato (2018b).

As a proof-of-concept, we verify our theoretical results by numerically constructing packings from certain initial tessellations, and by ascertaining their hyperuniformity from the spectral densities. Our procedure allows for the use of any initial tessellation, including Voronoi tessellations Torquato (2002) and their generalizations Gellatly and Finney (1982); Imai et al. (1985), sphere tessellations, disordered isoradial graphs Boutillier and de Tilière (2012), dissected tessellations Gabbrielli et al. (2012), Delaunay triangulations, and “Delaunay-centroidal” tessellations Florescu et al. (2009); Torquato and Chen (2018b), as long as it meets the bounded-cell condition. For concreteness and simplicity, two types of initial tessellations are considered in this work: Voronoi tessellations (Sec. V) and sphere tessellations, which are necessarily multiscale divisions of space (Sec. VI). Employing Voronoi tessellations of general disordered point patterns in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$, we demonstrate that our methodology enables a remarkable mapping that converts extremely large nonhyperuniform packings (as large as $10^{8}$ particles) into hyperuniform ones. To carry out our simulations, we employ various statistically homogeneous progenitor systems. Based on the same idea, we establish the hyperuniformity of the aforementioned optimal Hashin-Shtrikman multiscale coated-spheres structures Hashin (1962); Torquato (2002); Torquato and Pham (2004); Milton (1986) [see Fig. 1(d)]. Here, we provide a detailed derivation of $\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(k\right)}}$ in the $m$th stage and numerical simulations for two distinct types of cell-volume distributions. It is noteworthy that our methodology only involves calculating cell volumes, which is exactly performed and easy to parallelize. Furthermore, we demonstrate that large samples of many of our designs can be readily fabricated via modern photolithographic and 3D printing techniques Wong and Hernandez (2012); Tumbleston et al. (2015); Shirazi et al. (2015); Zhao and Mason (2018) (Sec. VII). While some of the major results were announced in a brief communication Kim and Torquato (2019), there we focused on applications to sphere packings without detailed derivations. In this work, we treat a broader class of sphere packings as well as packings of nonspherical particles and provide detailed mathematical derivations. We also report associated simulation results that are not contained in Ref. Kim and Torquato (2019).

We present basic mathematical definitions and concepts in Sec. II. Then, we precisely describe the tessellation-based procedure in Sec. III. In Sec. IV we derive the small-$k$ scalings of the spectral densities for the constructed packings. Subsequently, we verify our theoretical results by numerical simulations using Voronoi tessellations and sphere tessellations in Secs. V and VI, respectively. Then, we discuss the feasibility of fabricating our designs in modern technologies in Sec. VII. Finally, we provide concluding remarks in Sec. VIII.

II Background and Definitions

The microstructure of a two-phase system can be described by the phase indicator function associated with phase $i=1,2$ Torquato (2002):

[TABLE]

If the system is statistically homogeneous, then its one-point correlation function is independent of position $\bm{r}$, and identical to the phase-volume fraction $\phi_{i}$, i.e., $\left\langle\mathinner{\mathcal{I}^{(i)}\mathopen{\left(\bm{r}\right)}}\right\rangle=\phi_{i}$, where $\left\langle\cdot\right\rangle$ represents an ensemble average. The autocovariance function can be defined in terms of the mean-zero fluctuating indicator function, $\mathinner{\mathcal{J}^{(i)}\mathopen{\left(\bm{r}\right)}}\equiv\mathinner{\mathcal{I}^{(i)}\mathopen{\left(\bm{r}\right)}}-\phi_{i}$, as follows:

[TABLE]

which is identical for each phase and tends to zero as $r$ increases if the system does not have long-range order. Its Fourier transform $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\bm{k}\right)}}\equiv\int_{\mathbb{R}^{d}}\mathinner{d\bm{y}}e^{-i\bm{k}\cdot\bm{y}}\mathinner{\chi_{{}_{V}}\mathopen{\left(\bm{y}\right)}}$, called the spectral density, is a nonnegative real-valued function of a wavevector $\bm{k}$. In experiments, the spectral densities are directly obtainable from elastic scattering intensities Debye and Bueche (1949) when the wavelength of radiation is larger than atomic distance, but shorter than the length scale of domains. In numerical simulations, the spectral densities are calculated from realizations of the media under the periodic boundary conditions as follows:

[TABLE]

where $\left|\mathcal{V}_{F}\right|$ is the volume of the simulation box, a wavevector $\bm{k}$ is a reciprocal lattice vector of the simulation box, and $\left\langle\cdot\right\rangle$ represents an ensemble average, where $\mathinner{\tilde{\mathcal{J}}^{(i)}\mathopen{\left(\bm{k}\right)}}$ is the Fourier transform of $\mathinner{\mathcal{J}^{(i)}\mathopen{\left(\bm{r}\right)}}$ Torquato (1999).

In the context of two-phase media, a packing can be regarded as domains of a “particle” phase ($N$ generally shaped particles $\mathbf{P}_{1},~{}\mathbf{P}_{2},\cdots,~{}\mathbf{P}_{N}$) that are dispersed throughout a continuous “matrix” (void) phase. For such a packing in a periodic fundamental cell $\mathcal{V}_{F}$, the random variable $\mathinner{\tilde{\mathcal{J}}\mathopen{\left(\bm{k}\right)}}$ associated with the particle phase (dropping the superscript) can be expressed as follows:

[TABLE]

where $\bm{r}_{j}$ is the centroid of $\mathbf{P}_{j}$, $\mathinner{\tilde{m}\mathopen{\left(\bm{k};\mathbf{P}_{j}\right)}}$ is its form factor (i.e., the Fourier transform of the “particle indicator function” $\mathinner{m\mathopen{\left(\bm{r};\mathbf{P}_{j}\right)}}$ with respect to $\bm{r}_{j}$), and $\delta_{\bm{k},\bm{0}}$ represents a Kronecker delta symbol. Since the forward scattering term $\phi\left|\mathcal{V}_{F}\right|\delta_{\bm{k},\bm{0}}$ in Eq. (8) always vanishes at nonzero reciprocal lattice vectors $\bm{k}$’s, this term is often ignored in numerical calculations.

For some special particle shapes, closed-form expressions of $\mathinner{\tilde{m}\mathopen{\left(\bm{k};\mathbf{P}_{j}\right)}}$ are known. For a spherical particle of radius $a$, the associated form factor is

[TABLE]

where $k\equiv\left|\bm{k}\right|$ and $\mathinner{J_{n}\mathopen{\left(x\right)}}$ is the Bessel function of order $n$. For a cubic particle of side length $L$, the corresponding function is

[TABLE]

where $k_{l}$ represents $l$th component of a wavevector $\bm{k}$ and

[TABLE]

We note that Eqs. (9) and (10) have their global maxima at the origin, whose values are identical to their particle volumes.

From Eqs. (6) and (7), one can straightforwardly derive an expression of the spectral density for sphere packings of the identical particle radius $a$ Torquato (2002, 2016):

[TABLE]

where $\rho$ is number density and $\mathinner{S\mathopen{\left(\bm{k}\right)}}$ is the structure factor defined as

[TABLE]

Equation (12) gives

[TABLE]

and thus, one obtains a hyperuniform packing by decorating a hyperuniform point pattern with spheres of an equal size Torquato (2016). In this paper, however, we will not discuss such hyperuniform constructions.

In the case of a one-component many-particle system in thermal equilibrium in which the particles do not overlap, the fluctuation-compressibility relation $\mathinner{S\mathopen{\left(\bm{0}\right)}}=\rho\kappa_{T}k_{B}T$ Torquato (2002); Hansen and McDonald (1990) and Eq. (14) yield

[TABLE]

where $\kappa_{T}\equiv\rho^{-1}(\partial\rho/\partial p)_{T}$ is the isothermal compressibility, $p$ is the pressure, $k_{B}$ is the Boltzmann constant, and $T$ is temperature. Equation (15) implies that any compressible ($\kappa_{T}>0$) one-component system in thermal equilibrium cannot be hyperuniform at a positive temperature Torquato et al. (2015); Torquato (2018a).

The local volume-fraction variance $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}$ associated with spherical windows of radius $R$ is defined as Torquato (2018a); Zachary and Torquato (2009)

[TABLE]

where $\mathinner{\tau\mathopen{\left(\bm{x};R\right)}}$ denotes the local volume fraction of the particle phase inside the spherical window of radius $R$ centered at position $\bm{x}$.

From numerical simulations alone, it is difficult to ascertain whether a system is perfectly hyperuniform because the infinite-sample-size limit is never achievable and $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}$ usually has large relative statistical uncertainties at small wavenumbers. For these reasons, it is desirable to employ alternative criteria to determine whether a system is effectively hyperuniform. A useful empirical criterion to deem a system to be hyperuniform is that the hyperuniformity metric $H$ is less than $10^{-2}$ or $10^{-3}$ Atkinson et al. (2016); Torquato (2018a); Chen et al. (2018), where $H$ is defined by

[TABLE]

where $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\bm{k}_{\mathrm{peak}}\right)}}$ is the spectral density at the first dominant (non-Bragg) peak. Note that this criterion is different from its counterpart for point patterns because of the presence of the form factor in the spectral density Chen et al. (2018); see Eq. (12).

III Tessellation-Based Procedure

Here, we precisely describe the tessellation-based procedure in $d$-dimensional Euclidean space $\mathbb{R}^{d}$. For a periodic cubic fundamental cell $\mathcal{V}_{F}$ of side length $L$ in $\mathbb{R}^{d}$, the procedure is performed as follows:

Divide the simulation box with $N$ disjoint cells $\mathbf{C}_{1}$, $\cdots$, $\mathbf{C}_{N}$ [Figs. 1(a) and 1(b)] in which the maximal characteristic linear cell size $\ell_{\max}\equiv\max_{j=1}^{N}\left\{\max_{\bm{r},\bm{r}^{\prime}\in\mathbf{C}_{j}}\{\left|\bm{r}-\bm{r}^{\prime}\right|\}\right\}$ is much smaller than $L$, i.e.,

[TABLE]

We call this the “bounded-cell” condition. 2. 2.

For a specified local-cell packing fraction $0<\phi<1$, place hard particles of arbitrary shapes of total volume $\phi\left|\mathbf{C}_{j}\right|$ within the $j$th cell $\mathbf{C}_{j}$, and then repeat the same process over all cells [see Figs. 1(c) and 1(d) for illustrative examples with disks].

Application of this procedure results in a packing in which the local-cell packing fraction $\phi$ is identical to global packing fraction. Given an initial tessellation, this 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|\mathbf{C}_{j}\right|$ and $\left|\mathbf{P}_{j}\right|_{\max}$ represent volumes of the $j$th cell and the largest particle of a certain shape which is inscribed in this cell, respectively.

Roughly speaking, the maximal packing fraction $\phi_{\max}$ becomes larger when each cell tends to fully enclose a larger particle. Clearly, the largest particle that a cell can fully circumscribe should be congruent to the cell, and thus the maximal packing fraction will take its largest value ($\phi_{\max}=1$) only if the shape of each particle is similar to its circumscribing cell. For this reason, when only considering spherical particles, packing fraction of the multiscale coated-disks model [Fig. 1(d)] can span up to unity, i.e., $\phi_{\max}=1$.

The rationale behind our methodology can be intuitively understood by considering how the volume of the particle phase within a spherical window of radius $R$ fluctuates around its global mean value $\phi\mathinner{v_{1}\mathopen{\left(R\right)}}$. As shown in Fig. 2, our methodology ensures that the cells only in the gray-shaded region can contribute to the fluctuations, which are proportional to volume of this boundary region. Since the bounded-cell condition ensures that thickness of this boundary region is smaller than $\ell_{\max}$, the resulting variance in local phase-volume grows on the order of $R^{d-1}$ for sufficiently large windows ($R\gg\ell_{\max}$). This implies that $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}\sim R^{-d-1}$, meaning that this packing is strongly hyperuniform (class I). This rationale also explains the hyperuniformity of the constructed packings, regardless of shapes and number of particles inside each cell. However, for simplicity, we henceforth focus on the cases in which each cell contains exactly one particle.

IV General Theoretical Analyses

Here, we derive an asymptotic expression for the spectral density for the constructed packings in the small-wavenumber limit. Consider an initial tessellation $\{\mathbf{C}_{j}\}_{j=1}^{N}$ of a cubic periodic fundamental cell $\mathcal{V}_{F}$ of side length $L$ in $\mathbb{R}^{d}$. Since a fundamental cell is the union of all cells $\mathbf{C}_{1}$, $\cdots$, $\mathbf{C}_{N}$ of the tessellation, the Fourier transform of $\mathcal{V}_{F}$ can be decomposed as follows:

[TABLE]

where $\delta_{\bm{k},\bm{0}}$ represents the Kronecker delta symbol, $\bm{k}$ is a reciprocal lattice vector of $\mathcal{V}_{F}$, and for the $j$th cell $\mathbf{C}_{j}$, $\bm{x}_{j}$, $\mathinner{m\mathopen{\left(\bm{r};\mathbf{C}_{j}\right)}}$, and $\mathinner{\tilde{m}\mathopen{\left(\bm{k};\mathbf{C}_{j}\right)}}$ represent its centroid, the indicator function with respect to $\bm{x}_{j}$, and the form factor, respectively.

The application of our procedure to this tessellation yields a particle packing that consists of $N$ particles $\mathbf{P}_{1}$, $\cdots$, $\mathbf{P}_{N}$ whose centroids are $\bm{r}_{1}$, $\cdots$, $\bm{r}_{N}$, respectively, with the identical local-cell packing fraction $\phi$. Using the decomposition of Eq. (20) and the spectral density given in Eqs. (6) and (7), we obtain the following general expression:

[TABLE]

where $\Delta\bm{X}_{j}\equiv\bm{r}_{j}-\bm{x}_{j}$, and $\mathinner{\tilde{m}\mathopen{\left(\bm{k};\mathbf{P}_{j}\right)}}$ is the form factor of a particle $\mathbf{P}_{j}$.

Due to the bounded-cell condition, $\left|\bm{r}\right|<\ell_{\max}\ll L$ for every $\bm{r}$ in each cell $\mathbf{C}_{j}$, or equivalently, one can consider small wavevectors satisfying that $\left|\bm{k}\cdot\bm{r}\right|\sim\left|\bm{r}\right|/L\ll 1$. Thus, the form factor of $\mathcal{V}=(\mathbf{C}_{j}\text{ or }\mathbf{P}_{j}$) can be well approximated by its Taylor series about $\bm{k}=\bm{0}$:

[TABLE]

where the Einstein summation convention is employed,

[TABLE]

is the $m$th moment of the mass distribution of $\mathcal{V}(=\mathbf{C}_{j}\text{ or }\mathbf{P}_{j})$ that is normalized by its volume $\left|\mathcal{V}\right|$, and $r_{\alpha_{j}}$ represents the $\alpha_{j}$th Cartesian component of a vector $\bm{r}$. We note that since Eq. (23) refers to the moments with respect to the centroid of $\mathcal{V}$, its first moment is identically zero.

Using Eq. (22), the term $\mathinner{\tilde{\mathcal{J}}\mathopen{\left(\bm{k}\right)}}$ given in Eq. (21) can be written as follows (see Appendix B for details):

[TABLE]

where

[TABLE]

and $\mathinner{\tilde{\mathcal{J}}_{(2)}\mathopen{\left(\bm{k}\right)}}$ and $\mathinner{\tilde{\mathcal{J}}_{(3)}\mathopen{\left(\bm{k}\right)}}$ are defined in Eq. (51). Since the leading order terms of $\mathinner{\tilde{m}\mathopen{\left(\bm{k};\mathbf{P}_{j}\right)}}$ and $\mathinner{\tilde{m}\mathopen{\left(\bm{k};\mathbf{C}_{j}\right)}}$ exactly cancel each other due to the local constraint $\left|\mathbf{P}_{j}\right|=\phi\left|\mathbf{C}_{j}\right|$ for all $j$, these three terms $\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}$, $\mathinner{\tilde{\mathcal{J}}_{(2)}\mathopen{\left(\bm{k}\right)}}$, and $\mathinner{\tilde{\mathcal{J}}_{(3)}\mathopen{\left(\bm{k}\right)}}$ exhibit the power-law scalings in the small-wavenumber limit. Assuming typical tessellations that exhibit $\left|\mathinner{\tilde{\mathcal{J}}_{(2)}\mathopen{\left(\bm{k}\right)}}\right|\sim\left|\bm{k}\right|^{2}$ and $\left|\mathinner{\tilde{\mathcal{J}}_{(3)}\mathopen{\left(\bm{k}\right)}}\right|\sim\left|\bm{k}\right|^{3}$, the first term can have two scalings [either $\left|\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}\right|\sim\left|\bm{k}\right|$ or $\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}\sim\mathinner{\mathcal{O}\mathopen{\left(\left|\bm{k}\right|^{2}\right)}}$], depending on the particle displacements $\Delta\bm{X}_{j}$ with respect to their cell centroids. Therefore, the scaling of $\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}$ determines that of the spectral density (21).

Specifically, whenever $\left|\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}\right|\sim\left|\bm{k}\right|$, which can be achieved when particle displacements $\Delta\bm{X}_{j}$ are uncorrelated with one another 222One such set of structures are perturbed crystal packing that result from stochastically displacing particles within a cell of a periodic Voronoi tessellation Gabrielli (2004); Kim and Torquato (2018)., the spectral density of the constructed packings tends to zero quadratically in $\left|\bm{k}\right|$; specifically,

[TABLE]

which corresponds to class I hyperuniformity. When $\left|\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}\right|\sim\mathinner{\mathcal{O}\mathopen{\left(\left|\bm{k}\right|^{2}\right)}}$, the packings are more strongly hyperuniform with a new scaling given by

[TABLE]

For example, this scaling can be achieved when $\Delta\bm{X}_{j}=\bm{0}$ for all $j$ [i.e., $\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}=0$] or when $\left|\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}\right|\sim\left|\bm{k}\right|^{2}$ due to certain spatial correlations in $\Delta\bm{X}_{j}$ [see Sec. V]. Therefore, the manipulation of particle displacements $\Delta\bm{X}_{j}$ enables us to engineer either quadratic or quartic scalings of the spectral density.

It is noteworthy that the appearance of $\phi^{2}$ factor in Eqs. (26) and (27) is common for small-$\left|\bm{k}\right|$ scalings of the spectral densities of all statistically homogeneous sphere packings: see Eq. (12). Here, we note that the theoretical results (26) and (27) can be straightforwardly generalized to the cases of multiple particles are added in each cell. This is achieved by dividing cells such that each subdivided cell should circumscribe a single particle with an identical local-cell packing fraction.

Now, we consider a special case where all particles are similar to the associated cells in the sense that the particles have identical shapes and orientations with their cells, but have different sizes ($\mathbf{P}_{j}=\phi^{1/d}\mathbf{C}_{j}$ for $j=1,\cdots,N$). Then, the $n$th moments of particles and cells can be related as $\mathinner{\mathcal{M}_{\alpha_{1}\cdots\alpha_{n}}\mathopen{\left(\mathbf{P}_{j}\right)}}=\phi^{n/d}\mathinner{\mathcal{M}_{\alpha_{1}\cdots\alpha_{n}}\mathopen{\left(\mathbf{C}_{j}\right)}}$. Substituting this expression into Eq. (24), we obtain

[TABLE]

Again, whenever $\left|\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}\right|\sim\left|\bm{k}\right|$, the spectral density of the constructed packings shows a scaling $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\bm{k}\right)}}\sim\phi^{2}\left|\bm{k}\right|^{2}$. However, in the special case of $\Delta\bm{X}_{j}=\bm{0}$ for all $j$ (i.e., $\tilde{\mathcal{J}}_{(1)}=0$), the resulting scaling is

[TABLE]

Both hyperuniform cases belong to the class I.

Here, we should note that all theoretically predicted scalings of $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}$, given in Eqs. (26), (27), and (29), are analytic at the origin, i.e., the power exponents are even positive integers. This implies that autocovariance functions $\mathinner{\chi_{{}_{V}}\mathopen{\left(\left|\bm{r}\right|\right)}}$ of the constructed packings must decay to zero exponentially fast (or faster) as $\left|\bm{r}\right|\to\infty$ Torquato (2018a).

V Hyperuniform Packings from Voronoi Tessellations

In this section, we formulate an efficient numerical algorithm that is based on our tessellation-based methodology to generate very large disordered hyperuniform packings in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$. This can be accomplished from various types of tessellations Imai et al. (1985); Gellatly and Finney (1982); Boutillier and de Tilière (2012); Florescu et al. (2009); Torquato and Chen (2018b); Gabbrielli et al. (2012), but we focus on Voronoi tessellations of disordered nonhyperuniform point patterns in this section. For a given point pattern, the Voronoi cell of a point is defined as the region of space closer to this point than any other points, and the Voronoi tessellation is the collection of all Voronoi cells Torquato (2002). The computational time for the Voronoi tessellation of a point pattern of $N$ particles in $\mathbb{R}^{d}$ is at most of the order of $\mathinner{\mathcal{O}\mathopen{\left(N\log N+N^{\lceil(d-1)/2\rceil}\right)}}$ Aurenhammer and Klein (2000). Due to such small computation cost, the implementation of our methodology in the case of Voronoi tessellations enables a remarkably efficient mapping that converts very large nonhyperuniform point patterns or packings into very large disordered hyperuniform packings.

We choose the progenitor point patterns for the Voronoi tessellations not only to meet the bounded-cell condition but so that the resulting packing is easy to fabricate. From a practical viewpoint, it is useful to employ similar particle shapes and sizes to construct the packings. All of these conditions can be readily fulfilled by considering progenitor point patterns derived from dense hard-sphere packings. In what follows, we elaborate the bounded-cell condition.

V.1 The bounded-cell condition

The bounded-cell condition in the initial tessellation is a central requirement to ensure hyperuniformity via the tessellation-based procedure. For Voronoi tessellations, the largest cell size $\ell_{\max}$ is on the order of the largest nearest neighbor distance, which in turn is of the order of the largest radius $r_{\max}$ of holes (i.e., spherical regions that are empty of particle centers); $\ell_{\max}\sim r_{\max}$. Therefore, the bounded-cell condition (18) can be satisfied if the relative size $r_{\max}/L$ of the largest hole to the sample size [or equivalently, their volume ratio $\delta\equiv\mathinner{v_{1}\mathopen{\left(r_{\max}\right)}}/L^{d}$] is much smaller than a certain small value.

The information about largest hole size is contained in the void-exclusion or hole probability function $\mathinner{E_{V}\mathopen{\left(r\right)}}$, which gives the probability for finding a hole of radius $r$ when it is randomly placed in a point-particle system in the thermodynamic limit Torquato (2002). Clearly, if the hole probability has compact support, i.e., $\mathinner{E_{V}\mathopen{\left(r\right)}}=0$ at any $r>\mathcal{D}$ for a certain length $\mathcal{D}$, then Voronoi tessellations of the associated point patterns always meet the bounded-cell condition for relatively small sample sizes (say $10^{d+1}$ in space dimension $d$). Examples of such systems include all crystals, disordered stealthy hyperuniform point patterns Zhang et al. (2017); Ghosh and Lebowitz (2017), and the saturated random sequential addition (RSA) packings Torquato (2002); Zhang and Torquato (2013). Specifically, RSA is a time-dependent process that irreversibly, randomly, and sequentially adds nonoverlapping spheres into space. In the infinite-time limit the resulting packing does not have any available space to add further particles, called saturated.

However, many other disordered point patterns, including Poisson point patterns, equilibrium hard-sphere liquids, and unsaturated RSA packings, can possess arbitrarily large holes in the thermodynamic limit. Finite-size samples of these systems tend to have larger holes as the sample size grows [see Fig. LABEL:sub@fig:MaxGapDist for 1D Poisson point patterns], but the hole size relative to the sample size decreases as the sample size grows [see Fig. LABEL:sub@fig:MaxGapDist_scaled], implying that samples, in fact, tend to meet the bounded-cell condition (18).

We now rigorously show that for sufficiently large statistically homogeneous point patterns, Voronoi tessellations of almost every realization should obey the bounded-cell condition, as long as their hole probabilities are smaller than or equal to that of the Poissonian counterparts for large hole radii. For this purpose, we will show that Voronoi tessellations of almost every Poisson point pattern obey the bounded-cell condition for sufficiently large sample sizes. We begin by considering, for finite-size samples, the hole probability $\mathinner{E_{V}\mathopen{\left(r\right)}}$ which can be interpreted as the probability that a finite sample possesses at least a single hole of radius greater than $r$ (i.e., $\ell_{\max}\sim r_{\max}\geq r$). Thus, $\mathinner{E_{V}\mathopen{\left(L[\delta/\mathinner{v_{1}\mathopen{\left(1\right)}}]^{1/d}\right)}}$ is, in turn, the probability that the Voronoi tessellation of a single sample does not meet the bounded-cell condition (18). For Poisson point patterns of $N$ particles and volume $V=L^{d}$ in $\mathbb{R}^{d}$, its hole probability can be straightforwardly obtained as follows 333 In a Poisson point pattern, the probability that a single particle is placed in the space such that a hole of radius $r$ is empty of a particle is $(V-\mathinner{v_{1}\mathopen{\left(r\right)}})/V$ because all positions are equally probable. Since particle positions are mutually independent, the hole probability $\mathinner{E_{V}\mathopen{\left(r\right)}}$ that $N$ particles are placed in the space such that this hole is empty of particles is simply $[(V-\mathinner{v_{1}\mathopen{\left(r\right)}}/V]^{N}$.:

[TABLE]

where $\delta\equiv\mathinner{v_{1}\mathopen{\left(r\right)}}/V$. This quantity converges to a well-known expression in the thermodynamic limit Torquato (2002) $\mathinner{\exp\mathopen{\left(-\rho\mathinner{v_{1}\mathopen{\left(r\right)}}\right)}}$, where $\rho$ is the number density. Equation (30) decays exponentially fast for a given $\delta<1$ as particle number $N$ increases; see Fig. LABEL:sub@fig:Ev_Poisson. This implies that for sufficiently large sample sizes, almost every realization of a Poisson point pattern should obey the bounded-cell condition. Thus, the constructed packings from Voronoi tessellations of Poisson point patterns should be hyperuniform, which is consistent with a recent study of random fields Kampf and Spodarev (2017).

We will employ correlated and statistically homogeneous point patterns as the progenitor configurations, which are even more likely to obey the bounded-cell condition than Poisson-point counterparts at the same number density. In the ensuing discussion, we show that such nonhyperuniform point patterns can be converted into disordered hyperuniform packings.

V.2 Spherical particles

Here, we numerically implement our procedure by solely rescaling particle volumes without changing particle centers and shapes. It results in disordered sphere packings whose Voronoi tessellations are identical to those of their progenitor packings and we show that they are exactly hyperuniform of class I.

Due to the fixed particle centers, we define two alternative maximal packing fractions $\phi_{\max}^{(1)}$ and $\phi_{\max}^{(2)}$; see Fig. 4. For $\phi_{\max}^{(1)}$, the volume $\left|\mathbf{P}_{j}\right|_{\max}$ of the largest particle of a Voronoi cell in definition (19) is that of the smallest spherical particle inscribed in the cell. For $\phi_{\max}^{(2)}$, the largest particle is the particle in the progenitor packing. By definition, we have the following inequalities: $\phi_{\max}\geq\phi^{(1)}_{\max}\geq\phi^{(2)}_{\max}$, where $\phi^{(1)}_{\max}=\phi^{(2)}_{\max}$ occurs only when every particle in the progenitor has a neighbor in contact, and $\phi_{\max}=\phi^{(1)}_{\max}$ occurs only when every particle in the progenitor is inscribed in its Voronoi cell. Obviously, if a constructed packing has the packing fraction $\phi\leq\phi_{\max}^{(2)}$, none of its particles can be larger than those in the progenitor packing. For 2D and 3D saturated RSA packings, values of $\phi^{(1)}_{\max}$ are around $0.35$ or $0.25$, respectively. These values tend to decrease as the packing fraction of the progenitor packings is smaller, and they become arbitrarily small for Poisson point patterns (i.e., RSA with the zero packing fraction). Additional values of $\phi^{(1)}_{\max}$ and $\phi^{(2)}_{\max}$ are summarized in Sec. IB in the Supplemental Material Sup .

We employ three different types of nonhyperuniform progenitor point patterns in both $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$: RSA packing Torquato (2002); Zhang and Torquato (2013), equilibrium hard-sphere liquids, and lattice-packings with point vacancies. First, we employ RSA packings for various values of initial packing fractions $\phi_{\mathrm{init}}$. These RSA packings are efficiently generated by the voxel-list algorithm Zhang and Torquato (2013). Second, as progenitor packings, we use equilibrium hard-sphere liquid configurations with a range of values of the initial packing fraction $\phi_{\mathrm{init}}$; see the Supplemental Material Sup for employed parameters. In 2D, its $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\bm{0}\right)}}$ can be well approximated from Eq. (15) and a formula for the pressure 444The isothermal compressibility $\kappa_{T}\equiv-V^{-1}\mathop{\left(\frac{\partial V}{\partial p}\right)_{T}}$ of 2D equilibrium hard-sphere liquids with the packing fraction $\phi$ is computed from the equation of states $p/(\rho k_{B}T)=1+2\phi\mathinner{g_{2}\mathopen{\left(D\right)}}$, where $\mathinner{g_{2}\mathopen{\left(D\right)}}$ is the contact value of the radial pair-correlation Torquato (2002). The contact value is $\mathinner{g_{2}\mathopen{\left(D\right)}}=\mathinner{G\mathopen{\left(D\right)}}=(1-0.436\phi)/(1-\phi)^{2}$ for $0\leq\phi\leq\phi_{f}$; see Ref. Torquato1995_2 as follows:

[TABLE]

where $\phi$ is the packing fraction of the decorated spheres.

The last type of progenitor packings are vacancy-riddled square and simple cubic lattices in two and three dimensions, respectively. These progenitor systems are characterized by the number $N_{s}$ of initial lattice sites and the fraction $c$ of point vacancies to $N_{s}$. While these imperfect lattices still have Bragg peaks, they are not hyperuniform Durian (2017); Kim and Torquato (2018) and it is easy to generate extremely large samples ($N\sim 10^{8}$). Furthermore, for high values of $c$, their cell-volume distributions are similar to those for Poisson point patterns. Therefore, this investigation will immediately demonstrate our discussion in Sec. V.1.

We compute the Voronoi tessellations of the three aforementioned types of progenitor packings via VORO++ library Rycroft (2009) (see Appendix A for the implementation), and then perform our methodology without changing particle centers. To get a visual sense of the resulting hyperuniform packings, we show representative but small hyperuniform sphere packings in 2D and 3D derived from RSA initial conditions in Fig. 5.

We compute the local volume-fraction variances $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}$ of certain progenitor and constructed packings by Monte Carlo sampling of windows Torquato (2002); Quintanilla and Torquato (1997) for many values of $R$ up to $L$ (i.e., side length of the simulation box). For illustrative purposes, Fig. 6 plots such variances for progenitor packings that are 2D saturated RSA packings of various system sizes ($N=10^{4}$, $10^{5}$, and $10^{6}$). We note that these results for $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}$ are reliable only up to a window radius $R<L/4$ Dreyfus et al. (2015); Chieco et al. (2017) due to the finite-size effects. The constructed packings exhibit a common scaling $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}\sim R^{-(d+1)}$, as predicted from the heuristic rationale associated with Fig. 2. This rationale also predicts that such a scaling starts from window sizes $R\gtrsim 5\ell_{\max},$ where $\ell_{\max}$ is the maximal cell length defined in step 18 [cf. Eq. (18)], from which the region where volume-fraction fluctuations arise (the gray-shaded region in Fig. 2) can be effectively regarded as window boundary. Since the Voronoi tessellations of 2D saturated RSA that have a small $\ell_{\max}$ ($\lesssim 1.7\rho^{-1/2}$), the class I hyperuniform scaling is well-established in length scales over several orders of magnitude (even beyond the reliable regime) for even relatively small samples ($N=10^{4}$) and the largest ones ($N=10^{6}$); see Fig. 6.

Figures 7, 8, and 9 summarize the simulation results for packings that were converted from RSA, equilibrium hard-sphere liquids, and vacancy-riddled lattices, respectively. By construction, these constructed packings have particle-volume distributions that are identical to those of Voronoi cell volumes in their progenitors; see Figs. LABEL:sub@fig:2DRSA_Voronoicellvolume, LABEL:sub@fig:3DRSA_Voronoicellvolume, LABEL:sub@fig:2D_HSF_volDist, and LABEL:sub@fig:volDist_2DImLat. Furthermore, since particle centers are fixed in the procedure, the progenitor packings and their associated constructed packings have identical Voronoi tessellations and thus obviously have identical local statistics associated with Voronoi cells. This implies that local statistics alone generally may or may not determine hyperuniformity, which we elaborate in Sec. VIII.

Isotropic spectral densities of progenitor and constructed packings are computed by using Eqs. (6) and (9). To compare them, we rescale all packings to a common packing fraction $\phi=0.01$; see Figs. LABEL:sub@fig:chivk_2DRSA, LABEL:sub@fig:3DRSA_chivk, LABEL:sub@fig:2D_HSF_chiv, and LABEL:sub@fig:chiv_2DImLat. For small wavenumbers (i.e., $0<k\overline{a}\ll 1$, where $\overline{a}\sim 0.1\rho^{-1/d}$ is the mean of particle radii), as shown in these figures, the spectral densities of the progenitors are significantly different from those of the constructed packings. While progenitor packings are clearly nonhyperuniform (i.e., $H\sim 0.1-0.01$ for RSA packings and equilibrium hard-sphere liquids), the constructed packings are class I hyperuniform with a common power-law scaling $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\left|\bm{k}\right|\right)}}\sim\left|\bm{k}\right|^{4}$. Similar to the case of $\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}$ shown in Fig. 6, this power-law scaling is particularly well-established over several orders of magnitude, even for relatively small samples [$N=10^{3}$ as shown in Figs. LABEL:sub@fig:2DRSA_chivk_size and LABEL:sub@fig:3DRSA_chivk_size]. This $k^{4}$-scaling results from $\left|\mathinner{\tilde{\mathcal{J}}_{(1)}\mathopen{\left(\bm{k}\right)}}\right|\sim\left|\bm{k}\right|^{2}$, which comes from special correlations in $\Delta\bm{X}_{j}$; see Sec. IV. Some values of $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(0\right)}}$ and $H$ are summarized in Table 1; see, for additional data, Sec. IC in the Supplemental Material Sup .

However, for intermediate wavenumbers (i.e., $1<k\overline{a}<3$), the spectral densities of the constructed packings largely resemble those of their progenitors; see, for example, Fig. LABEL:sub@fig:2D_HSF_chiv2. This result reflects that local statistics of both packings are identical at corresponding length scales. At large wavenumbers (i.e., $k\overline{a}>3$), however, the spectral densities of progenitor and constructed packings again become different from each other because the spectral density at this regime mainly depends on the particle-size distributions; see Fig. LABEL:sub@fig:2D_HSF_chiv2.

V.3 Nonspherical particles

As shown in the derivation in Sec. IV, the tessellation-based procedure can result in class I hyperuniformity even for nonspherical particle shapes, including regular and random polyhedrons. However, since it is generally nontrivial to compute the form factors of those particle shapes, we consider a simple geometric shape, i.e., a $d$-dimensional hypercubic particles; see Eq. (10).

For simplicity, we take advantage of the sphere packings constructed in Sec. V.2 by replacing the spherical particles with the inscribed cubic particles. In this case, the corresponding maximal packing fraction for cubes will be related to that of spheres as follows: $\phi_{\max}^{\mathrm{cube}}=\frac{\mathinner{\Gamma\mathopen{\left(1+d/2\right)}}2^{d-1}}{(d\pi)^{d}}\phi_{\max}^{\mathrm{sphere}}$. We consider two cases: (identically) oriented particles and randomly oriented particles. Figure 10 summarizes the simulation results for 2D and 3D packings converted from saturated RSA packings of $N\approx 10^{6}$. We note that the results for squares or cubes are largely indistinguishable from those for spheres because both particle shapes have comparable length scales and the Fourier transforms of both a cube (10) and a sphere (9) are isotropic near the origin and have similar profiles. However, these spectral densities will exhibit different behaviors at intermediate and large wavenumbers, depending on particles shapes and orientations.

VI Multiscale Coated-Spheres Model

The tessellation-based methodology can be applied to any type of tessellations of space besides the Voronoi tessellations, as long as they obey the bounded-cell condition (18). One interesting example is a tessellation of space that consists of nonoverlapping spherical cells. Since a monodisperse sphere packings in $d>1$ cannot occupy all space, spherical cells in such tessellations should have a polydispersity in size down to the infinitesimally small; see Fig. 1(b). For these (multiscale) sphere tessellations, the bounded-cell condition can be guaranteed by making the ratio of the largest cell-volume to the sample volume, $\delta=v_{\max}/\left|\mathcal{V}_{F}\right|$, sufficiently small. Implementing our methodology from these tessellations should result in hyperuniform packings regardless of particle shapes, positions, and numbers per cell.

A simple but important case of such disordered hyperuniform packings are ones in which each cell includes only a single particle to which the cell is concentric; see Fig. 1(d). Since particle shapes are similar to those for cells, the local-cell packing fraction $\phi$ of those packings can span up to unity, as mentioned in Sec. III. In addition, we can use the expression of the second moment of a sphere of radius $R$, ${\mathcal{M}_{\alpha\beta}}=\frac{1}{(2+d)}R^{2}\delta_{\alpha,\beta}$, $\Delta\bm{X}_{j}=\bm{0}$, and Eq. (28) to obtain the following power-law scaling for the spectral density:

[TABLE]

where radius of cell $j$ is denoted by $\mathcal{R}_{j}$.

It is crucial to observe that cells and the associated particles form composite spheres that fill all space. Each composite sphere is comprised of a spherical core (particle) of one 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 global phase volume fraction, which is guaranteed by our procedure. Thus, this structure is a packing (dispersion) in which spherical particles (blue regions) are “well-separated” from one another in a fully connected (continuous) void (matrix) phase (red region). This is traditionally understood as a key characteristic of optimal two-phase structures Torquato and Chen (2018b, a).

Surprisingly, these hyperuniform multiscale structures are identical to the Hashin-Shtrikman multiscale coated-spheres models Hashin (1962); Torquato (2002); Milton (1986); Torquato and Pham (2004), which are perhaps the most famous results from the theory of heterogeneous materials. These special particle composites are optimal for the effective electric (thermal) conductivity and bulk elastic modulus for a given phase fractions and phase properties. This observation will shed light on the origin of nearly optimal transport properties of disordered stealthy hyperuniform packings Zhang et al. (2016); Chen and Torquato (2018) and cellular networks derived from hyperuniform systems Torquato and Chen (2018b).

The coated-spheres structures are infinitely degenerate with varying degrees of order/disorder. The most ordered structures would be ones, derived from certain deterministic sphere tessellations, such as Apollonian gaskets, or initial lattice packings of identical spheres in which particles are added in a sequential multiscale manner.

VI.1 Simulation model

In what follows, we describe our model to simulate the coated-spheres model and ascertain their hyperuniformity. We first construct very dense disordered sphere packings via a multistage version of the RSA procedure, in which volume of the inserted spheres reduced in every stage. These dense “precursor” packings are later used to simulate ideal sphere tessellations. The coated-spheres model is then simulated by scaling spheres in the precursor packings at a certain ratio.

Specifically, these precursor packings are constructed from an empty simulation box $\mathcal{V}_{F}$ in $\mathbb{R}^{d}$ under the periodic boundary conditions. The construction procedure has two control parameters: the upper bound $v_{\max}$ on cell volumes and a positive decreasing function $\mathinner{g\mathopen{\left(i\right)}}$ of positive integers $i$, which satisfies $\mathinner{g\mathopen{\left(1\right)}}=1$ and $\sum_{j=1}^{\infty}\mathinner{g\mathopen{\left(j\right)}}<\infty$. The infinite-sample-size limit can be achieved as the upper bound $v_{\max}$ tends to be infinitesimally small with the simulation box fixed. Subsequently, one determines the prescribed number $N$ of spheres that will be inserted in every stage and the maximum cell volume $v^{(1)}$ to fill all space:

[TABLE]

where $G\equiv\sum_{j=1}^{\infty}\mathinner{g\mathopen{\left(j\right)}}$ and $\lceil x\rceil$ is the ceiling function. Every sphere in the $m$th stage has volume $v^{(m)}=v^{(1)}\mathinner{g\mathopen{\left(m\right)}}$ and the associated diameter $D_{m}$. Using these parameters, the precursor packings are constructed by the following steps:

In the first stage ($m=1$), irreversibly, randomly, and sequentially add (i.e., via RSA procedure) nonoverlapping spheres of a diameter $D_{1}$. The insertion stops when $N$ spheres are added, unless the packing becomes saturated. 2. 2.

In the $m$th stage ($m>1$), nonoverlapping spheres of a diameter $D_{m}$ ($<D_{m-1}$) are added via the RSA procedure. Make sure that newly inserted spheres do not overlap with the spheres in the previous stages $1,2,\cdots,m-1$. The insertion stops when the packing reaches to a prescribed covering fraction $N/{\left|\mathcal{V}_{F}\right|}\sum_{i=1}^{m}v^{(i)}$, unless the packing becomes saturated. 3. 3.

The procedure stops when it reaches to a prescribed number of stages; see Fig. 1(b).

Note that the aforementioned voxel-list algorithm Zhang and Torquato (2013) is implemented in step 1 and 2. Using the $m$th stage precursor packings, we simulate the coated-spheres model of the inclusion volume fraction $\phi$ by reducing these spheres at a volume ratio $\phi$. Figure 11 shows a representative but small disordered multiscale coated-disks construction.

We note that due to saturation at each stage, the number $N_{m}$ of spheres inserted in the $m$th stage can be different from the prescribed number $N$. Let $\mathcal{N}_{m}\equiv\sum_{i=1}^{m}N_{i}$ and $\eta_{m}=\sum_{i=1}^{m}N_{i}v^{(i)}/\left|\mathcal{V}_{F}\right|$ denote the total number of spheres at the end of the $m$th stage and the fraction of space covered by these spheres, called the covering fraction, respectively. By construction, $\eta_{m}$ cannot exceed the prescribed covering fraction in the $m$th stage:

[TABLE]

Thus, in a finite $m$th stage, the precursor packings will not cover all space but have gaps that can be only covered by smaller spheres in the next stages. As the number of stages tends to be infinite, those gaps are eventually covered by spheres of size down to the infinitesimally small, i.e., $\eta_{m}\to 1$ as $m\to\infty$.

In this work, we consider two different types of volume scalings: a power-law scaling and an exponential scalings. For the power-law scaling, the cell volume in the $m$ stage is determined by

[TABLE]

for a given scaling exponent $p>1$. From this relation, it is straightforward to derive the relation between maximum cell volume $v^{(1)}$ and the prescribed insertion number $N$, and the prescribed covering fraction:

[TABLE]

where $\mathinner{\zeta\mathopen{\left(p\right)}}$ denotes the Riemann zeta function.

For the exponential scaling, cell volume in the $m$th stage is determined by

[TABLE]

for a given scaling base $q>1$. The analogues of Eqs. (37) and (38) are

[TABLE]

VI.2 Theoretical analyses

While in the limit of $m\to\infty$ our coated-spheres model is predicted to be strongly hyperuniform [see Eq. (32)], our model in a finite $m$th stage will be nearly hyperuniform, rather than perfectly hyperuniform, due to the uncovered gaps ($\eta_{m}<1$). To estimate the degree of hyperuniformity of our model in the $m$th stage, we consider the associated spectral density $\mathinner{\tilde{\chi}^{(m)}_{{}_{V}}\mathopen{\left(\bm{k}\right)}}$, given by

[TABLE]

where $\phi$ denotes the local-cell packing fraction. Without any prior knowledge of cell-volume distribution, its rigorous bound can be obtained by using some simple inequalities (see Appendix C.1):

[TABLE]

where

[TABLE]

Here, the constant term $\phi^{2}\left|\mathcal{V}_{F}\right|(1-\eta_{m})^{2}$ is the largest volume-fraction fluctuations that can arise from uncovered gaps in the $m$th stage coated-spheres model. Importantly, the constant term depends on $\phi(1-\eta_{m})$, i.e., the deviation between the local-cell packing fraction $\phi$ and the global packing fraction $\phi\eta_{m}$, which will vanish as the uncovered gaps become completely filled in the limit of $m\to\infty$. Thus, in this limit, the upper bound implies that our model becomes perfectly hyperuniform of class I, consistent with Eq. (32).

Although the upper bound (43) is rigorous, it is a gross overestimation compared to corresponding simulation results. Thus, we obtain a better estimate of Eq. (42) by assuming that the uncovered gaps are spatially uncorrelated, which effectively removes the ensemble average of their cross terms, yielding

[TABLE]

where the summation term becomes dominant as $\left|k\right|\to 0$ because $\mathinner{F\mathopen{\left(\bm{k};\phi\right)}}$ shows a power-law scaling; see Appendix C.2 for details.

Further estimation of the second term in Eq. (45) requires the cell-size distribution. Consider two distinct cell-size scalings discussed in Sec. VI.1: a power-law form (36) and an exponential functional form (39). For the power-law scaling, using Eqs. (37) and (38), the summation in Eq. (45) can be written as

[TABLE]

where $\mathinner{f\mathopen{\left(m\right)}}\equiv\left(\sum_{j=m+1}^{\infty}j^{-2p}\right)/\left(\sum_{j=m+1}^{\infty}j^{-p}\right)^{2}$. Substituting Eq. (46) into Eq. (45) yields an approximation of the residual spectral density in the small-wavenumber limit as follows:

[TABLE]

For the exponential scaling, Eqs. (40) and (41) are used to simplify Eq. (45) as follows:

[TABLE]

which results in the following leading-order term of the spectral density in the $m$th stage:

[TABLE]

VI.3 Simulation results

For the purpose of illustration and simplicity, we specialize to two dimensions to generate multiscale-disk tessellations for the aforementioned cell-size scalings: power-law and exponential functional forms, as defined by Eqs. (36) and (39), respectively. Using the power-law functional forms, the simulations proceed up to the $400$th stage, yielding a scaling exponent $p$ ranging from 1.5 to 1.8, and the ratios of the maximal cell volume to the sample volume $\delta=10^{-3}$ and $10^{-4}$. For the exponential functional forms, the constructions proceed to achieve around covering fraction $\eta_{m}\approx 0.99$ for values of a scaling base $q$=$1.05$, $1.10$, and $1.20$, and the ratio $\delta=10^{-3}$.

For both types of scaling functions, as the scaling parameters increases, the smaller is the number of stages needed to achieve a prescribed covering fraction. Instead, for even larger scaling parameters, the precursor packings are more likely to be saturated Shier and Bourke (2013), which often results in significant computational costs. This is because as the packing approaches to a saturation state, the voxel-list algorithm that we employ subdivides current voxels with an increasingly finer resolution, which requires increasingly larger computer memory and computational times. For this reason, we choose scaling parameters smaller than 2.

To obtain multiscale-disk tessellations that nearly fill space ($\eta_{m}\approx 0.95$), these procedures need to continue up to around a few hundred stages [Fig. LABEL:sub@fig:2DZeta_coveringfraction and LABEL:sub@fig:2DPower_coveringfraction]. For $m\lesssim 10$, cell volumes in the exponential scalings do not change much, compared to those in the power-law scalings. Thus, the tessellations in the exponential scalings can cover a larger fraction of the simulation box with a smaller number of stages than those for the power-law scaling. Instead, for the exponential scaling, the precursor packings more easily achieve saturation from relatively low covering fractions ($\eta_{m}\approx 0.65$). Then, the simulated covering fraction never keeps up with the prescribed covering fraction (41); see Fig. LABEL:sub@fig:2DPower_coveringfraction. However, the number of total cells $\mathcal{N}_{m}$ required to reach to a certain covering fraction will be similar for both cases; see Fig. LABEL:sub@fig:2DZeta_totalNumber and LABEL:sub@fig:2DPower_totalnumber.

From these precursor polydisperse packings, we simulate the coated-disks model with the inclusion volume fraction $\phi=0.5$. To ascertain their hyperuniformity, we compute the associated spectral densities $\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(\bm{k}\right)}}$. As shown in Fig. LABEL:sub@fig:2DZeta_chiv and LABEL:sub@fig:2DPower_chiv, the packings in the finite stages (around $\eta_{m}\approx 0.95$) are effectively hyperuniform. Importantly, this comes from the fact that the precursor packings have uncovered gaps, and thus the local-cell packing fraction $\phi$ and their global packing fraction $\phi\eta_{m}$ have a discrepancy. The consequent volume-fraction fluctuations are estimated by Eqs. (47) and (49). Figures LABEL:sub@fig:2DZeta_chiv0 and LABEL:sub@fig:2DPower_chiv0 show that while the numerical results in $\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(k=k_{\min}\right)}}$ and these theoretical estimations conform to each other for low covering fractions ($\eta_{m}\lesssim 0.8$), they significantly deviate for relatively high covering fractions. This is because the uncovered gaps are no longer spatially uncorrelated in high covering fractions.

However, for high covering fractions, the residual spectral densities $\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(k=k_{\min}\right)}}$ vanish with an error on the order of $(1-\eta_{m})^{2}$. Thus, our upper bound (43), which although is a gross overestimation, correctly predicts the scaling behavior of $\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(k=k_{\min}\right)}}$. Both the theoretical and numerical results show that the coated-spheres model should be strongly hyperuniform in the limit of $m\to\infty$.

It is noteworthy that as the number $m$ of stages increases, the spectral densities in Fig. LABEL:sub@fig:2DZeta_chiv and LABEL:sub@fig:2DPower_chiv tend to resemble one another for intermediate wavenumbers. In addition, the terms $\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(k=k_{\min}\right)}}$, as shown in Fig. LABEL:sub@fig:2DZeta_chiv0 and LABEL:sub@fig:2DPower_chiv0, collapse to a single curve that seemingly only depends on the fraction of uncovered space. From these observations, we surmise that the spectral densities of the multiscale coated-spheres structures ($\eta_{m}\to\infty$) are identical to one another, regardless of the cell-volume distributions of the precursor sphere packings. This universal behavior of spectral densities is apparently related to the fact that the coated-spheres structures have identical effective properties.

VII Fabrication of Our Designs

Disordered hyperuniform structures, constructed by numerical simulations have been produced via modern fabrication technologies Man et al. (2013); Ma et al. (2016). Here, we explicitly discuss how our designed hyperuniform packings (dispersions) can be fabricated via state-of-the-art technologies, such as 2D photolithographic Zhao and Mason (2018) and 3D printing techniques Wong and Hernandez (2012); Shirazi et al. (2015); Tumbleston et al. (2015).

Photolithography is a microfabrication technique that uses light to transfer a designed 2D pattern in the photomask on the photo-sensitive chemicals coated on a flat substrate. Then, after a series of chemical treatments, the desired pattern is engraved on the material, or material is deposed on the pattern. These methods are widely used in industry and research because of their high efficiency and exact control over the shapes and size of the patterns. Instead, diffraction tends to round all sharp corners in the designed pattern, the radius of which is associated with the minimum feature size. State-of-the-art photolithography techniques are capable of creating patterns on a 30-cm-diameter wafer with the minimum feature size down to 25 nm Zhao and Mason (2018). With equipment of 1.5 $\mu$m minimum feature size, one can readily fabricate our 2D hyperuniform disk packings (probably derived from Voronoi tessellations) that include more than one million particles.

Three-dimensional hyperuniform packings can be fabricated using 3D additive manufacturing techniques Wong and Hernandez (2012); Shirazi et al. (2015); Tumbleston et al. (2015). 3D printing refers to various processes that solidify materials layer by layer to create a 3D object (e.g., fused filament fabrication, stereolithography, and selective laser melting). A printed structure must be topologically connected such that it can be mechanically self-supporting after the procedure. Thus, the void (matrix) phase of our 3D hyperuniform packings(dispersions) can be readily printed Zerhouni et al. (2018). Due to recent developments in 3D printing methods, some commercial desktop 3D printers can print a sample with dimensions $125\times 125\times 125~{}\mathrm{cm}^{3}$ in 50 h with around 100 $\mu\mathrm{m}$ $XY$ resolution and 20 $\mu\mathrm{m}$ in $Z$ resolution. Setting the minimum pore size is 300 $\mu$m, such devices can readily fabricate our 3D hyperuniform sphere packings that include up to 50 million pores. If pore sizes are on the order of the resolutions, then spherical pores will be suitable to reduce the effect of thermal deformations during the printing process.

VIII Conclusions and Discussion

In summary, we have introduced the tessellation-based methodology to construct disordered hyperuniform packings (dispersions) in $d$-dimensional Euclidean space $\mathbb{R}^{d}$. This procedure is simple to implement because it only involves computing cell volumes, each of which can be performed exactly and in parallel. Furthermore, it guarantees that the constructed packings are perfectly hyperuniform of class I in the infinite-sample-size limit, as we analytically showed in Eqs. (26), (27), and (29). We have implemented numerically our methodology from two distinct types of disordered tessellations: Voronoi tessellations of nonhyperuniform progenitor packings and sphere tessellations.

In the case of Voronoi tessellations, we demonstrated that our methodology provides a remarkable mapping that converts virtually all samples of any statistically homogeneous point pattern, whether hyperuniform or nonhyperuniform, into perfectly hyperuniform packings. Furthermore, the fact that it is easy to create a Voronoi tessellation of a very large point pattern with the combination of our efficient procedure enables us to construct very large samples (of the order of $10^{8}$ particles) of perfectly hyperuniform disordered packings. Such large system sizes have not been possible for previous numerical methods.

In the case of sphere tessellations, we established that the optimal Hashin-Shtrikman multiscale dispersions are hyperuniform of class I. In addition to the fact that the spherical inclusions in such dispersions are “well-separated” from one another Torquato and Chen (2018b), hyperuniformity is apparently another important structural attribute to attain optimal effective transport and elastic properties Milton (2002); Torquato (2002); Hashin (1962); Milton (1986); Torquato and Pham (2004). In this regard, it is noteworthy that some disordered hyperuniform packings Zhang et al. (2016); Chen and Torquato (2018) and cellular networks derived from hyperuniform systems Torquato and Chen (2018b) have been reported to possess (nearly) optimal effective transport and elastic properties. Thus, it will be interesting to investigate the physical properties of hyperuniform packings constructed by our procedure.

Furthermore, it is important to observe that our procedure allows many distinct types of tessellations as long as they meet the bounded-cell condition (18). Besides Voronoi tessellations and sphere tessellations studied in this work, examples include disordered isoradial graphs Boutillier and de Tilière (2012), Delaunay triangulations, “Delaunay-centroidal” tessellations Florescu et al. (2009); Torquato and Chen (2018b), dissected tessellations Gabbrielli et al. (2012), and various generalizations of Voronoi tessellations, such as Laguerre tessellations Gellatly and Finney (1982); Imai et al. (1985) and tessellations in Manhattan distances Krause (1986). Since constructing progenitor packings is the most time-consuming step in our numerical implementations, one would increase sample sizes by employing other tessellations that can be efficiently generated.

We note that similar to our procedure, a 1D model Gabrielli et al. (2004) and the “equal-volume tessellation” Gabrielli et al. (2008) enable the generation of hyperuniform point mass patterns and hyperuniform point patterns from certain initial tessellations, respectively. Specifically, these two models place a point mass Gabrielli et al. (2004) and point particles Gabrielli et al. (2008) in each cell such that mass densities and number densities in the cells become identical, respectively. Thus, the hyperuniform systems in both models can be regarded to be a zero-$\phi$ limit of the hyperuniform packings in our procedure. Importantly, however, both previous models Gabrielli et al. (2004, 2008) do not consider volume-fraction fluctuations, which is our central concern and accounted for by our procedure. In addition, the mass density of point masses in this 1D model does not change the structural characteristics of the resulting systems, which again is different from our methodology. We also note that the equal-volume tessellation is less versatile in its initial tessellations than our procedure because all cell volumes in an initial tessellation should be integer multiples of a common unit volume.

We have shown that our procedure enables a mapping that converts a nonhyperuniform packing into a hyperuniform one without changing the initial tessellation, i.e., the tessellation is an invariant under the transformation. In other words, packings that have identical tessellations can either be nonhyperuniform or hyperuniform by simply tuning local characteristics. It immediately follows that any tessellation statistic, including the distributions of nearest-neighbor distances and Voronoi-cell volumes, are identical for both the progenitor nonhyperuniform point patterns and their corresponding constructed hyperuniform packings. These results reinforce previous observations that local structural characteristics may or may not determine the hyperuniformity of systems, e.g., substantial local clustering of particles may not be inconsistent with hyperuniformity Zachary and Torquato (2011) and disordered systems with appreciable short-range order are often not hyperuniform Zhang and Torquato (2013); see Ref. Torquato (2018a). Moreover, our results also show that transitions from a nonhyperuniform state to a hyperuniform state can be achieved without correlated movements of particles/mass across all length scales. Such transitions stand in contrast to those in some previous hyperuniform systems in thermal equilibrium Jancovici (1981); Uche et al. (2004); Batten et al. (2008); Zhang et al. (2015a); Lomba et al. (2017); Chen and Torquato (2018) as well as nonequilibrium Hexner and Levine (2015); Tjhung and Berthier (2015); Weijs et al. (2015); Kurita and Weeks (2011); Ricouvier et al. (2017); Torquato and Jiao (2010); Lubachevsky and Stillinger (1990), in which the transitions always involve collective rearrangements of the particles.

It should not go unnoticed that our methodology also allows ones to tune particle shapes, positions, and numbers within each cell with preserving hyperuniformity of the constructed packings. For instance, one can engineer these two-phase systems at large length scales (i.e., $k\ll 1$) by choosing nonspherical particles with various aspect ratios and placing them away from the centroids of the associated cells. However, at intermediate length scales the two-phase media are structurally similar to their progenitor patterns (see Sec. V.2). One can tune local structures of the hyperuniform packings to achieve the “well-separated” condition, which is a necessary requirement to attain the optimal transport and mechanical properties (see Sec. VI). Moreover, the small-$k$ scaling of the spectral density can either be $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}\sim k^{2}$ or $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}\sim k^{4}$ by engineering particle displacements $\Delta\bm{X}_{j}$ with respect to the associated cell centroids (see Sec. IV). Due to this tunability, our methodology allows ones to design an enormous class of hyperuniform packings, including (nearly) optimal structures. Combining our computational designs with the aforementioned 2D and 3D fabrication techniques Kim and Torquato (2019) is expected to accelerate the discovery of novel disordered hyperuniform two-phase materials.

There are several other open questions for future exploration. For example, could other types of initial tessellations lead to scaling behaviors of the spectral density besides quadratic or quartic shown here? To what extent can our procedure be generalized by relaxing the identical local-cell packing faction constraint (e.g., giving certain spatial correlations in local-cell packing fractions)? Can such generalized versions of our procedure construct “disordered” stealthy (i.e., $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}=0$ for $k<K$) Batten et al. (2008); Zhang et al. (2015a, b); Torquato et al. (2015) or class II hyperuniform (i.e., $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}\sim k$, such as MRJ packings Jiao and Torquato (2011); Atkinson et al. (2016))? If possible, then this would certainly allow one to obtain more general scaling behaviors.

Acknowledgements

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

Appendix A Voronoi tessellations

The Voronoi tessellation of a given progenitor packing is computed via VORO++ library Rycroft (2009). To enhance the performance, we divide the particle centroids of the progenitor packing into several domains and compute the Voronoi tessellations for each domain in parallel. To avoid any deformation in Voronoi cells due to domain boundaries, we choose domains in the following steps:

Divide the simulation box into several disjoint subdomains in cubic shape. 2. 2.

Add a marginal region surrounding each subdomain, and these two regions form a domain. Here, the thickness $W$ of the marginal region is chosen as $6\rho^{-1/d}$, where $\rho$ is the number density.

Then, for each domain, we compute the Voronoi cells of points inside “subdomains.” For point patterns with larger holes, the thickness of the marginal region should be increased accordingly.

Since VORO++ is designed for three-dimensional geometries, its 2D implementation is performed in a quasi 3D simulation box whose height (in $z$ component) is unity. To avoid any possible “memory leakage” in this 2D implementation, the number of particles within a domain should be smaller than $10^{5}$.

Appendix B Derivation of Eq. (24)

First, we rewrite the Fourier transform of the particle indicator function $\mathinner{\tilde{\mathcal{J}}\mathopen{\left(\bm{k}\right)}}$, given in Eq. (21):

[TABLE]

Here, we substitute the form factors $\mathinner{\tilde{m}\mathopen{\left(\bm{k};\mathbf{P}_{j}\right)}}$ and $\mathinner{\tilde{m}\mathopen{\left(\bm{k};\mathbf{C}_{j}\right)}}$ with their Taylor expansions (22), which yields

[TABLE]

In Eq. (50), we used the identical local-cell packing fraction condition, i.e., $\left|\mathbf{P}_{j}\right|=\phi\left|\mathbf{C}_{j}\right|$ for $j=1,\cdots,N$. In Eq. (51), $(\Delta\bm{X}_{j})_{\gamma}$ is the $\gamma$th Cartesian component of a vector $\Delta\bm{X}_{j}$, and we apply the Taylor expansion $\left(e^{-i\bm{k}\cdot\Delta\bm{X}_{j}}-1\right)=-ik_{\gamma}(\Delta\bm{X}_{j})_{\gamma}+\mathinner{\mathcal{O}\mathopen{\left(k^{2}\right)}}$, which is a good approximation due to the bounded-cell condition (18).

Appendix C Derivations of the spectral density of the $m$th stage coated-spheres model

Here, we compute upper bounds on the spectral density $\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(\bm{k}\right)}}$ defined in Eq. (42) in three cases of cell-volume distributions: unknown, a power-law scaling, and an exponential scaling.

C.1 Upper bounds of the spectral density

From Eq. (42),

[TABLE]

Now, note the following inequality for two complex numbers $A$ and $B$:

[TABLE]

which is obtained by successively applying the triangle inequality and then the inequality of arithmetic and geometric means, and thus equality occurs if and only if $A=B$. Using this inequality, we obtain a rigorous upper bound on the spectral density described in Eq. (53):

[TABLE]

Subsequently, we derive an upper bound on the second term in Eq. (56) by applying the triangle inequality:

[TABLE]

where the inequality in Eq. (57) comes from the fact that $\left|\mathinner{\tilde{m}\mathopen{\left(\bm{k};R\right)}}\right|\leq\mathinner{v_{1}\mathopen{\left(R\right)}}$; see Eq. (9). We note that the last term in Eq. (57) can be regarded as the largest volume-fraction fluctuations contributed from the uncovered gaps. Combining Eqs. (56) and (57), we obtain a rigorous bound as follows:

[TABLE]

where $\mathinner{F\mathopen{\left(\bm{k};\phi\right)}}$ is defined by Eq. (44). We remark that this bound is derived without any prior knowledge of the cell-volume distribution.

C.2 Derivation of Eq. (45)

We can approximate the spectral density given in Eq. (42) by assuming that its contributions from composite spheres [the first term in Eq. (53)] and uncovered gaps are uncorrelated, which effectively removes an ensemble average of their cross terms:

[TABLE]

Here, the second term can be further simplified by assuming that the uncovered gaps are spatially uncorrelated:

[TABLE]

where the last approximation comes from $\left|\mathinner{\tilde{m}\mathopen{\left(k;R\right)}}\right|\approx\mathinner{v_{1}\mathopen{\left(R\right)}}+\mathinner{\mathcal{O}\mathopen{\left(k^{2}\right)}}$ when $kR\ll 1$; see Eq.(9).