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

J. Kim  S. Torquato<EMAIL_ADDRESS>[ Department of Physics,
Princeton University, Princeton, New Jersey 08544, USA  Department of
Chemistry, Princeton University, Princeton, New Jersey 08544, USA  Princeton
Institute for the Science and Technology of Materials, Princeton University,
Princeton, New Jersey 08544, USA  Program in Applied and Computational
Mathematics, Princeton University, Princeton, New Jersey 08544, USA

###### Abstract

Disordered hyperuniform dispersions are exotic amorphous two-phase materials
characterized by an anomalous suppression of long-wavelength volume-fraction
fluctuations, which endows 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 and simple construction procedure that guarantees perfect
hyperuniformity for very large sample sizes. This methodology 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.

###### keywords:

Disordered hyperuniformity, Heterogeneous materials, Hashin-Shtrikman coated-
spheres model

††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

$\mathinner{\chi_{{}_{V}}\mathopen{\left(\mathbf{r}\right)}}\equiv\mathinner{S_{2}^{(i)}\mathopen{\left(\mathbf{r}\right)}}-{\phi_{i}}^{2},$
(1)

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]

$\displaystyle\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}=\frac{1}{\left|\mathcal{V}_{d}\right|}\times$
$\displaystyle\left|\sum_{j=1}^{N}\mathinner{\tilde{m}\mathopen{\left(\left|\mathbf{k}\right|;R_{j}\right)}}e^{-i\mathbf{k}\cdot\mathbf{r}_{j}}-\phi\int_{\mathcal{V}_{d}}\mathinner{d\mathbf{r}^{\prime}}e^{-i\mathbf{k}\cdot\mathbf{r}^{\prime}}\right|^{2},$
(2)

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:

$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}=\frac{1}{\left|\mathcal{V}_{d}\right|}\left|\sum_{j=1}^{N}\mathinner{\tilde{m}\mathopen{\left(\left|\mathbf{k}\right|;R_{j}\right)}}e^{-i\mathbf{k}\cdot\mathbf{r}_{j}}\right|^{2},$
(3)

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:

$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(\mathbf{k}\right)}}\sim\left|\mathbf{k}\right|^{\alpha}~{}~{}~{}(\left|\mathbf{k}\right|\to
0).$ (4)

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]:

$\mathinner{\sigma_{V}^{2}\mathopen{\left(R\right)}}\sim\begin{cases}R^{-(d+1)},&\alpha>1~{}(\mathrm{Class~{}I})\\\
R^{-(d+1)}\ln{R},&\alpha=1~{}(\mathrm{Class~{}II})\\\
R^{-(d+\alpha)},&\alpha<1~{}(\mathrm{Class~{}III})\end{cases},$

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:

$\displaystyle N$
$\displaystyle\equiv\lceil\left|\mathcal{V}_{d}\right|/\left(v_{\max}\sum_{m=1}^{\infty}\mathinner{g\mathopen{\left(m\right)}}\right)\rceil,~{}\mathrm{and}$
(5) $\displaystyle\left|\mathcal{V}_{d}\right|$
$\displaystyle=Nv^{(1)}\sum_{i=1}^{\infty}\mathinner{g\mathopen{\left(i\right)}},$
(6)

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:

  1. 1.

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.

Figure 1: Schematics illustrating the tessellation-based procedure to
construct disordered hyperuniform dispersions. (a) One first tessellates the
space with disjoint cells whose maximal lengths are shorter than a certain
length scale $D(\ll L)$, where $\left|C_{j}\right|$ represents the volume of
cell $j$. (b) For a specified packing fraction $0<\phi<1$, within each cell
$C_{j}$ one places a single particle of general shape of volume
$\phi\left|C_{j}\right|$. This ensures that the local-cell packing fraction
$\phi$ is identical to the global packing fraction of the resulting
dispersion. (c) A special case of constructed dispersions with circular disks
and a circular observation window of large radius $R(\gg D)$. The volume-
fraction fluctuations arise only inside a narrow yellow-shaded region whose
effective thickness is thinner than $D$, and thus the resulting fluctuations
are on the order of $R^{d-1}$, i.e., class I hyperuniformity. (d) The
corresponding hyperuniform dispersion of randomly oriented squares. The same
rationale explains the hyperuniformity of the dispersion of squares (d) as
well as the Hashin-Shtrikman multiscale coated-disks structures shown in Fig.
4(b).

###  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:

  1. 1.

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.,

$\phi\leq\phi_{\max}\equiv\min_{j=1}^{N}\left\\{\frac{\mathinner{v_{1}\mathopen{\left(R^{\mathrm{max}}_{j}\right)}}}{\left|C_{j}\right|}\right\\},$
(7)

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

Figure 2: Implementation of the tessellation-based methodology via Voronoi
tessellations. (a) A portion of Voronoi tessellation of a progenitor packing
(i.e., a 2D saturated RSA packing) of white circles. (b) A hyperuniform
dispersion (blue disks) of packing fraction $\phi$ that is constructed by our
methodology from (a) with the particle centers fixed. (c) Schematic of two
types of maximal packing fractions $\phi_{\max}^{(1)}$ and $\phi_{\max}$. In
this cell, while at $\phi=\phi_{\max}$ the particle can be as large as the
largest inscribed circle (green circle), at $\phi=\phi_{\max}^{(1)}$ the
particle size can be as large as the largest inscribed circle (black) that is
concentric to the initial particle (white circle).

Figure 3: Simulation results of dispersions constructed from saturated RSA
packings in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$. (a) The probability density
function of Voronoi cell volume $\left|C_{j}\right|$ of 2D saturated RSA
packings on a semi-log scale (larger panel) and a linear scale (inset). Here,
$\overline{\left|C\right|}$ represents the average cell volume. Error bars
represent standard deviations (see Sec. 2.1). (b) Spectral densities of 2D
saturated RSA packings and the constructed dispersions on a log-log (larger
panel) and a linear scale (inset), where $a\equiv\rho^{-1/d}$, $d$ is space
dimension, and $\rho$ is the number density of particle centers. We note that
both progenitor packings and the constructed dispersions are scaled to the
packing fraction $\phi=0.01$. (c) Log-log plot of normalized spectral
densities
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\right)}}/\phi^{2}$ of the
constructed dispersions for various local-cell packing fractions $\phi$. (d-f)
Corresponding results of constructed dispersions from 3D saturated RSA
packings.

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)).

Figure 4: Illustration of the construction of the coated-spheres model via
sphere-tessellations in two dimensions. (a) An initial multiscale-disk tiling
(400th stage) and (b) a dispersion constructed via the tessellation-based
procedure from (a), which is the coated-disks model.

##  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.

Figure 5: Structural characteristics of the simulated coated-disks models
(Sec. 2.5). (a) Covering fraction $\eta_{m}$ of simulated initial tilings as
functions of stage number $m$. (b) Spectral densities for the simulated
multiscale coated-disks model at various $m$’s at the local-cell packing
fraction $\phi=0.5$. Error bars represent standard deviations (see Sec. 2.1).
(c) Spectral density at the minimal wavenumber
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(k\to 0\right)}}$ as a
function of fraction of uncovered space $1-\eta_{m}$.

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]:

$\displaystyle\mathinner{\tilde{\chi}_{{}_{V}}^{(m)}\mathopen{\left(\mathbf{k}\right)}}\leq$
$\displaystyle\left(\frac{\phi(1-\phi^{2/d})}{2(2+d)}\right)^{2}\frac{2}{\left|\mathcal{V}_{d}\right|}\left|\sum_{j=1}^{\mathcal{N}_{m}}\mathinner{v_{1}\mathopen{\left(\mathcal{R}_{j}\right)}}(k\mathcal{R}_{j})^{2}e^{-i\mathbf{k}\cdot\mathbf{r}_{j}}\right|^{2}$
$\displaystyle+2\left|\mathcal{V}_{d}\right|\phi^{2}(1-\eta_{m})^{2},$ (8)

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.

Figure 6: Fabrication of designed hyperuniform dispersions and matrices in
three dimensions. (a) A portion of designed dispersion of spherical particles
(left) and the corresponding matrix (right) at $\phi=0.23$, which are
constructed from a 3D saturated RSA packing. (a) A portion of designed
dispersion of cubical particles (left) and the corresponding matrix (right).
The hyperuniform matrices with spherical or cubical pores can be fabricated
using 3D printing techniques.

##  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.

## References

  * [1] S. Torquato, F. H. Stillinger, Local density fluctuations, hyperuniformity, and order metrics, Phys. Rev. E 68 (4) (2003) 041113.  doi:10.1103/PhysRevE.68.041113.
  * [2] C. E. Zachary, S. Torquato, Hyperuniformity in point patterns and two-phase random heterogeneous media, J. Stat. Mech: Theory Exp. 2009 (12) (2009) P12015.  doi:10.1088/1742-5468/2009/12/P12015.
  * [3] S. Torquato, Hyperuniform States of Matter, Phys. Rep. 745 (2018) 1 – 95.  doi:10.1016/j.physrep.2018.03.001.
  * [4] M. Florescu, S. Torquato, P. J. Steinhardt, Designer disordered materials with large, complete photonic band gaps, Proc. Natl. Acad. Sci. U.S.A. 106 (49) (2009) 20658–20663.  doi:10.1073/pnas.0907744106.
  * [5] L. S. Froufe-Pérez, M. Engel, J. J. Sáenz, F. Scheffold, Band gap formation and anderson localization in disordered photonic materials with structural correlations, Proc. Natl. Acad. Sci. U.S.A. 114 (36) (2017) 9570–9574.  doi:10.1073/pnas.1705130114.
  * [6] C. López, The true value of disorder, Adv. Opt. Mater. 6 (16) (2018) 1800439\.  doi:10.1002/adom.201800439.
  * [7] G. Zhang, F. Stillinger, S. Torquato, Transport, geometrical, and topological properties of stealthy disordered hyperuniform two-phase systems, J. Chem. Phys. 145 (24) (2016) 244109.  doi:10.1063/1.4972862.
  * [8] D. Chen, S. Torquato, Designing disordered hyperuniform two-phase materials with novel physical properties, Acta Mater. 142 (2018) 152–161.  doi:10.1016/j.actamat.2017.09.053.
  * [9] Q. Le Thien, D. McDermott, C. J. O. Reichhardt, C. Reichhardt, Enhanced pinning for vortices in hyperuniform pinning arrays and emergent hyperuniform vortex configurations with quenched disorder, Phys. Rev. B 96 (9) (2017) 094516.  doi:10.1103/PhysRevB.96.094516.
  * [10] O. Leseur, R. Pierrat, R. Carminati, High-density hyperuniform materials can be transparent, Optica 3 (7) (2016) 763–767.  doi:10.1364/OPTICA.3.000763.
  * [11] C. J. Brinker, G. W. Scherer, Sol-gel science: The physics and chemistry of sol-gel processing, Academic press, 2013.
  * [12] M. Sahimi, Heterogeneous Materials I: Linear Transport and Optical Properties, Vol. 22, Springer-Verlag, New York, 2003.
  * [13] S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Interdisciplinary Applied Mathematics, Springer Science & Business Media, 2002.
  * [14] B. Patel, T. I. Zohdi, Numerical estimation of effective electromagnetic properties for design of particulate composites, Mater. Des. 94 (2016) 546–553.  doi:10.1016/j.matdes.2016.01.015.
  * [15] W. Mickel, S. C. Kapfer, G. E. Schröder-Turk, K. Mecke, Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter, J. Chem. Phys. 138 (4) (2013) 044501.  doi:10.1063/1.4774084.
  * [16] P. Debye, A. M. Bueche, Scattering by an inhomogeneous solid, J. Appl. Phys. 20 (6) (1949) 518–525.  doi:10.1063/1.1698419.
  * [17] O. U. Uche, F. H. Stillinger, S. Torquato, Constraints on collective density variables: Two dimensions, Phys. Rev. E 70 (4) (2004) 046122.  doi:10.1103/PhysRevE.70.046122.
  * [18] S. Torquato, G. Zhang, F. Stillinger, Ensemble theory for stealthy hyperuniform disordered ground states, Phys. Rev. X 5 (2) (2015) 021020.  doi:10.1103/PhysRevX.5.021020.
  * [19] G. Zhang, F. H. Stillinger, S. Torquato, Ground states of stealthy hyperuniform potentials: I. Entropically favored configurations, Phys. Rev. E 92 (2) (2015) 022119.  doi:10.1103/PhysRevE.92.022119.
  * [20] F. J. Dyson, Statistical theory of the energy levels of complex systems. I, J. Math. Phys. 3 (1) (1962) 140–156.  doi:10.1063/1.1703773.
  * [21] B. Jancovici, Exact results for the two-dimensional one-component plasma, Phys. Rev. Lett. 46 (6) (1981) 386.  doi:10.1103/PhysRevLett.46.386.
  * [22] D. Chen, E. Lomba, S. Torquato, Binary mixtures of charged colloids: A potential route to synthesize disordered hyperuniform materials, Phys. Chem. Chem. Phys. 20 (26) (2018) 17557–17562.  doi:10.1039/c8cp02616e.
  * [23] R. Feynman, M. Cohen, Energy spectrum of the excitations in liquid helium, Phys. Rev. 102 (5) (1956) 1189.  doi:10.1103/PhysRev.102.1189.
  * [24] S. Torquato, A. Scardicchio, C. E. Zachary, Point processes in arbitrary dimension from fermionic gases, random matrix theory, and number theory, J. Stat. Mech: Theory Exp. 2008 (11) (2008) P11019.  doi:10.1088/1742-5468/2008/11/P11019.
  * [25] A. Scardicchio, C. E. Zachary, S. Torquato, Statistical properties of determinantal point processes in high-dimensional euclidean spaces, Phys. Rev. E 79 (4) (2009) 041108.  doi:10.1103/PhysRevE.79.041108.
  * [26] A. Donev, F. H. Stillinger, S. Torquato, Unexpected density fluctuations in jammed disordered sphere packings, Phys. Rev. Lett. 95 (9) (2005) 090604.  doi:10.1103/PhysRevLett.95.090604.
  * [27] R. Kurita, E. R. Weeks, Incompressibility of polydisperse random-close-packed colloidal particles, Phys. Rev. E 84 (R) (2011) 030401.  doi:10.1103/PhysRevE.84.030401.
  * [28] Y. Jiao, S. Torquato, Maximally random jammed packings of Platonic solids: Hyperuniform long-range correlations and isostaticity, Phys. Rev. E 84 (4) (2011) 041309.  doi:10.1103/PhysRevE.84.041309.
  * [29] D. Hexner, D. Levine, Hyperuniformity of Critical Absorbing States, Phys. Rev. Lett. 114 (11) (2015) 110602.  doi:10.1103/PhysRevLett.114.110602.
  * [30] J. H. Weijs, R. Jeanneret, R. Dreyfus, D. Bartolo, Emergent hyperuniformity in periodically driven emulsions, Phys. Rev. Lett. 115 (10) (2015) 108301.  doi:10.1103/PhysRevLett.115.108301.
  * [31] R. Garcia-Millan, G. Pruessner, L. Pickering, K. Christensen, Correlations and hyperuniformity in the avalanche size of the oslo model, Europhys. Lett. 122 (5) (2018) 50003.  doi:10.1209/0295-5075/122/50003.
  * [32] D. Di Battista, D. Ancora, G. Zacharakis, G. Ruocco, M. Leonetti, Hyperuniformity in amorphous speckle patterns, Opt. Express 26 (12) (2018) 15594–15608.  doi:10.1364/OE.26.015594.
  * [33] S. Torquato, G. Zhang, M. de Courcy-Ireland, Uncovering multiscale order in the prime numbers via scattering, J. Stat. Mech: Theory Exp. 2018 (9) (2018) 093401\.  doi:10.1088/1742-5468/aad6be.
  * [34] H. L. Montgomery, The pair correlation of zeros of the zeta function, in: Proc. Symp. Pure Math, Vol. 24, 1973, pp. 181–193.
  * [35] M. Sodin, B. Tsirelson, Random complex zeroes, II. Perturbed lattice, Israel J. Math. 152 (1) (2006) 105–124.  doi:10.1007/Bf02771978.
  * [36] H. Noh, S. F. Liew, V. Saranathan, S. G. J. Mochrie, R. O. Prum, E. R. Dufresne, H. Cao, How noniridescent colors are generated by quasi-ordered structures of bird feathers, Adv. Mater. 22 (26-27) (2010) 2871–2880.  doi:10.1002/adma.200903699.
  * [37] Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, C. C. Joseph, S. Torquato, Avian photoreceptor patterns represent a disordered hyperuniform solution to a multiscale packing problem, Phys. Rev. E 89 (2) (2014) 022721.  doi:10.1103/PhysRevE.89.022721.
  * [38] A. Mayer, V. Balasubramanian, T. Mora, A. M. Walczak, How a well-adapted immune system is organized, Proc. Natl. Acad. Sci. U.S.A. 112 (19) (2015) 5950–5955.  doi:10.1073/pnas.1421827112.
  * [39] S. Kwon, J. M. Kim, Hyperuniformity of initial conditions and critical decay of a diffusive epidemic process belonging to the manna class, Phys. Rev. E 96 (1) (2017) 012146.  doi:10.1103/PhysRevE.96.012146.
  * [40] S. Atkinson, G. Zhang, A. B. Hopkins, S. Torquato, Critical slowing down and hyperuniformity on approach to jamming, Phys. Rev. E 94 (1-1) (2016) 012902.  doi:10.1103/PhysRevE.94.012902.
  * [41] M. Hejna, P. J. Steinhardt, S. Torquato, Nearly hyperuniform network models of amorphous silicon, Phys. Rev. B 87 (24) (2013) 245204.  doi:10.1103/PhysRevB.87.245204.
  * [42] M. Florescu, S. Torquato, P. Steinhardt, Effects of random link removal on the photonic band gaps of honeycomb networks, Appl. Phys. Lett. 97 (20) (2010) 201103\.  doi:10.1063/1.3505322.
  * [43] J. Ricouvier, R. Pierrat, R. Carminati, P. Tabeling, P. Yazhgur, Optimizing hyperuniformity in self-assembled bidisperse emulsions, Phys. Rev. Lett. 119 (20) (2017) 208001.  doi:10.1103/PhysRevLett.119.208001.
  * [44] G. Zito, G. Rusciano, G. Pesce, A. Malafronte, R. Di Girolamo, G. Ausanio, A. Vecchione, A. Sasso, Nanoscale engineering of two-dimensional disordered hyperuniform block-copolymer assemblies, Phys. Rev. E 92 (5) (2015) 050601.  doi:10.1103/PhysRevE.92.050601.
  * [45] C. E. Zachary, Y. Jiao, S. Torquato, Hyperuniform long-range correlations are a signature of disordered jammed hard-particle packings, Phys. Rev. Lett. 106 (17) (2011) 178001.  doi:10.1103/PhysRevLett.106.178001.
  * [46] E. Lomba, J. J. Weis, S. Torquato, Disordered hyperuniformity in two-component nonadditive hard-disk plasmas, Phys. Rev. E 96 (6-1) (2017) 062126.  doi:10.1103/PhysRevE.96.062126.
  * [47] S. Torquato, Y. Jiao, Robust algorithm to generate a diverse class of dense disordered and ordered sphere packings via linear programming, Phys. Rev. E 82 (6) (2010) 061302.  doi:10.1103/PhysRevE.82.061302.
  * [48] B. D. Lubachevsky, F. H. Stillinger, Geometric properties of random disk packings, J. Stat. Phys. 60 (5-6) (1990) 561–583.  doi:10.1007/Bf01025983.
  * [49] J. Kim, S. Torquato, Effect of imperfections on the hyperuniformity of many-body systems, Phys. Rev. B 97 (5) (2018) 054105.  doi:10.1103/PhysRevB.97.054105.
  * [50] S. Torquato, Perspective: Basic understanding of condensed phases of matter via packing models, J. Chem. Phys. 149 (2) (2018) 020901.  doi:10.1063/1.5036657.
  * [51] K. Zhao, T. G. Mason, Assembly of colloidal particles in solution, Rep. Prog. Phys. 80 (2018) 126601.  doi:10.1088/1361-6633/aad1a7.
  * [52] K. V. Wong, A. Hernandez, A review of additive manufacturing, ISRN mech eng 2012 (2012) 1–10.  doi:10.5402/2012/208760.
  * [53] J. R. Tumbleston, D. Shirvanyants, N. Ermoshkin, R. Janusziewicz, A. R. Johnson, D. Kelly, K. Chen, R. Pinschmidt, J. P. Rolland, A. Ermoshkin, E. T. Samulski, J. M. DeSimone, Continuous liquid interface production of 3D objects, Science 347 (6228) (2015) 1349–1352.  doi:10.1126/science.aaa2397.
  * [54] S. F. Shirazi, S. Gharehkhani, M. Mehrali, H. Yarmand, H. S. Metselaar, N. Adib Kadri, N. A. Osman, A review on powder-based additive manufacturing for tissue engineering: Selective laser sintering and inkjet 3D printing, Sci. Technol. Adv. Mater. 16 (3) (2015) 033502.  doi:10.1088/1468-6996/16/3/033502.
  * [55] S. Torquato, Exact conditions on physically realizable correlation functions of random media, J. Chem. Phys. 111 (19) (1999) 8832–8837.  doi:10.1063/1.480255.
  * [56] S. Torquato, Hyperuniformity and its generalizations, Phys. Rev. E 94 (2) (2016) 022122.  doi:10.1103/PhysRevE.94.022122.
  * [57] D. J. Durian, Hyperuniformity disorder length spectroscopy for extended particles, Phys. Rev. E 96 (3) (2017) 032910.  doi:10.1103/PhysRevE.96.032910.
  * [58] B. Widom, Random sequential addition of hard spheres to a volume, J. Chem. Phys. 44 (10) (1966) 3888–3894.  doi:10.1063/1.1726548.
  * [59] G. Zhang, S. Torquato, Precise algorithm to generate random sequential addition of hard hyperspheres at saturation, Phys. Rev. E 88 (5) (2013) 053312.  doi:10.1103/PhysRevE.88.053312.
  * [60] B. J. Gellatly, J. L. Finney, Characterisation of models of multicomponent amorphous metals: The radical alternative to the Voronoi polyhedron, J. Non-Cryst. Solids 50 (3) (1982) 313–329.  doi:10.1016/0022-3093(82)90093-X.
  * [61] H. Imai, M. Iri, K. Murota, Voronoi Diagram in the Laguerre Geometry and Its Applications, SIAM J. Comput. 14 (1) (1985) 93–105.  doi:10.1137/0214006.
  * [62] E. F. Krause, Taxicab geometry: An adventure in non-Euclidean geometry, Courier Corporation, 1986.
  * [63] C. H. Rycroft, VORO++: a three-dimensional Voronoi cell library in C++, Chaos 19 (4) (2009) 041111.  doi:10.1063/1.3215722.
  * [64] J. Shier, P. Bourke, An algorithm for random fractal filling of space, Comput. Graphics Forum 32 (8) (2013) 89–97.  doi:10.1111/cgf.12163.
  * [65] J. Kim, S. Torquato, Methodology to Construct Large Realizations of Perfectly Hyperuniform Disordered Packings (2019) in preparation.  arXiv:1901.10006
  * [66] S. Torquato, M. Avellaneda, Diffusion and reaction in heterogeneous media: Pore size distribution, relaxation times, and mean survival time, J. Chem. Phys. 95 (9) (1991) 6477–6489.  doi:10.1063/1.461519.
  * [67] J. Kampf, E. Spodarev, A functional central limit theorem for integrals of stationary mixing random fields, Theory Probab. Its Appl. 63 (1) (2017) 135–150.  doi:10.1137/S0040585X97T988952.
  * [68] Z. Hashin, The elastic moduli of heterogeneous materials, J. Appl. Mech. 29 (1) (1962) 143–150.  doi:10.1115/1.3636446.
  * [69] S. Torquato, D. Pham, Optimal bounds on the trapping constant and permeability of porous media, Phys. Rev. Lett. 92 (25) (2004) 255505.  doi:10.1103/PhysRevLett.92.255505.
  * [70] S. Torquato, D. Chen, Multifunctional hyperuniform cellular networks: optimality, anisotropy and disorder, Multifunct. Mater. 1 (1) (2018) 015001.  doi:10.1088/2399-7532/aaca91.
  * [71] S. Torquato, D. Chen, Multifunctionality of particulate composites via cross-property maps, Phys. Rev. Mater. 2 (9) (2018) 095603.  doi:10.1103/PhysRevMaterials.2.095603.
  * [72] O. Zerhouni, M.-G. Tarantino, K. Danas, F. Hong, Influence of the internal geometry on the elastic properties of materials using 3D printing of computer-generated random microstructures, Society of Exploration Geophysicists, 2018, pp. 3713–3718.  doi:10.1190/segam2018-2998182.1.
  * [73] B. Bhushan, M. Caspers, An overview of additive manufacturing (3d printing) for microfabrication, Microsyst. Technol. 23 (4) (2017) 1117–1124.  doi:10.1007/s00542-017-3342-8.
  * [74] A. Gabrielli, M. Joyce, S. Torquato, Tilings of space and superhomogeneous point processes, Phys. Rev. E 77 (3) (2008) 031125.  doi:10.1103/PhysRevE.77.031125.
  * [75] C. Boutillier, B. de Tilière, Statistical mechanics on isoradial graphs, in: J.-D. Deuschel, B. Gentz, W. König, M. von Renesse, M. Scheutzow, U. Schmock (Eds.), Probability in Complex Physical Systems, Springer Berlin Heidelberg, 2012, pp. 491–512.