# Topological Defects, Inherent Structures, and Hyperuniformity

Duyu Chen<EMAIL_ADDRESS>Materials Research Laboratory, University
of California, Santa Barbara, California 93106, United States  Yu Zheng
Department of Physics, Arizona State University, Tempe, AZ 85287  Yang Jiao
<EMAIL_ADDRESS>Materials Science and Engineering, Arizona State
University, Tempe, AZ 85287  Department of Physics, Arizona State University,
Tempe, AZ 85287

###### Abstract

Disordered hyperuniform systems are exotic states of matter that completely
suppress large-scale density fluctuations like crystals, and yet possess no
Bragg peaks similar to liquids or glasses. Such systems have been discovered
in a variety of equilibrium and non-equilibrium physical and biological
systems, and are often endowed with novel physical properties. While it is
well known that long-range interactions are necessary to sustain
hyperuniformity in thermal equilibrium at positive temperatures, such
condition is not required for the realization of disordered hyperuniformity in
systems out of equilibrium. However, the mechanisms associated with the
emergence of disordered hyperuniformity in nonequilibrium systems, in
particular inherent structures (i.e., local potential-energy minima) are often
not well understood, which we will address from a topological perspective in
this work. Specifically, we consider a representative class of disordered
inherent structures which are constructed by continuously introducing randomly
distributed topological defects (dislocations and disclinations) often seen in
colloidal systems and atomic-scale two-dimensional materials. We demonstrate
that these inherent structures can be viewed as topological variants of
ordered hyperuniform states (such as crystals) linked through continuous
topological transformation pathways, which remarkably preserve
hyperuniformity. Moreover, we develop a continuum theory to demonstrate that
the large-scale density fluctuations in these inherent structures are mainly
dominated by the elastic displacement fields resulted from the topological
defects, which at low defect concentrations can be approximated as
superposition of the displacement fields associated with each individual
defect (strain source). We find that hyperuniformity is preserved as long as
the displacement fields generated by each individual defect decay sufficiently
fast from the source (i.e., the volume integrals of the displacements and
squared displacements caused by individual defect are finite) and the
displacement-displacement correlation matrix of the system is diagonalized and
isotropic. Our results also highlight the importance of decoupling the
positional degrees of freedom from the vibrational degrees of freedom when
looking for disordered hyperuniformity, since the hyperuniformity property is
often cloaked by thermal fluctuations (i.e., vibrational degrees of freedom).

##  I Introduction

Disordered hyperuniform (DHU) systems are exotic states of matter Torquato and
Stillinger (2003); Torquato (2018) that lie between a perfect crystal and
liquid. These systems are similar to liquids or glasses in that they are
statistically isotropic and possess no Bragg peaks, and hence lack any
conventional long-range order, and yet they completely suppress large-scale
density fluctuations like crystals and in this sense possess a hidden long-
range order Torquato and Stillinger (2003); Zachary and Torquato (2009);
Torquato (2018). Specifically, the static structure factor $S(k)$, which is
directly proportional to the scattering intensity measured in scattering
experiments, vanishes for DHU systems in the infinite-wavelength (or zero-
wavenumber) limit, i.e., $\lim_{k\rightarrow 0}S(k)=0$, where $k$ is the
wavenumber Torquato and Stillinger (2003); Torquato (2018). Here $S(k)$ is
defined as $S(k)\equiv 1+\rho\tilde{h}(k)$, where $\tilde{h}(k)$ is the
Fourier transform of the total correlation function $h(r)=g_{2}(r)-1$,
$g_{2}(r)$ is the pair correlation function, and $\rho$ is the number density
of the system. Note that this definition implies that the forward scattering
contribution to the diffraction pattern is omitted. Equivalently, the local
number variance $\sigma_{N}^{2}(R)\equiv\langle N^{2}(R)\rangle-\langle
N(R)\rangle^{2}$ associated with a spherical observation window of radius $R$
grows more slowly than the window volume (i.e., a scaling of $R^{d}$ in
$d$-dimensional Euclidean space) for DHU systems in the large-$R$ limit
Torquato and Stillinger (2003); Torquato (2018), where $N(R)$ is the number of
particles in a spherical window with radius $R$ randomly placed into the
system. The small-$k$ scaling behavior of $S(k)\sim k^{\alpha}$ dictates the
large-$R$ asymptotic behavior of $\sigma_{N}^{2}(R)$, based on which all DHU
systems can be categorized into three classes: $\sigma_{N}^{2}(R)\sim R^{d-1}$
for $\alpha>1$ (class I); $\sigma_{N}^{2}(R)\sim R^{d-1}\ln(R)$ for $\alpha=1$
(class II); and $\sigma_{N}^{2}(R)\sim R^{d-\alpha}$ for $0<\alpha<1$ (class
III) Torquato (2018). It is also noteworthy that the direct correlation
function $c(r)$, defined via the Ornstein–Zernike relation $h(r)=c(r)+\rho
c(r)\bigotimes h(r)$ (where $\bigotimes$ denotes convolution) Ornstein and
Zernike (1914), becomes long-ranged in the sense that it has an unbounded
volume integral. This is in diametric contrast to standard thermal critical
points in which $h(r)$ is long-ranged, and hence a system at a hyperuniform
state is considered an “inverted” critical point Torquato and Stillinger
(2003).

DHU states have been discovered in a wide spectrum of equilibrium and non-
equilibrium physical and biological systems Gabrielli et al. (2002); Donev et
al. (2005); Zachary et al. (2011a); Jiao and Torquato (2011); Chen et al.
(2014); Zachary and Torquato (2011); Torquato et al. (2015); Uche et al.
(2004); Batten et al. (2008, 2009); Lebowitz (1983); Zhang et al. (2015a, b);
Kurita and Weeks (2011); Hunter and Weeks (2012); Dreyfus et al. (2015);
Hexner and Levine (2015); Jack et al. (2015); Weijs et al. (2015); Torquato et
al. (2008); Feynman and Cohen (1956); Jiao et al. (2014); Mayer et al. (2015);
Hejna et al. (2013); Klatt et al. (2019); Lei et al. (2019); Lei and Ni
(2019); Chremos and Douglas (2018); Rumi et al. (2019); Sánchez et al. (2019,
2020); Huang et al. (2021); Torquato (2021); see Ref. Torquato (2018) for a
thorough overview. The exotic structural features of DHU systems appear to
endow such systems with novel physical properties. For example, disordered
hyperuniform dielectric networks were found to possess complete photonic band
gaps comparable in size to photonic crystals, while at the same time
maintaining statistical isotropy, enabling waveguide geometries not compatible
with photonic crystals Florescu et al. (2009); Man et al. (2013). Moreover,
certain disordered hyperuniform patterns have superior color-sensing
capabilities, as demonstrated by avian photoreceptors Jiao et al. (2014).
Recent evidences also suggest that adding disorder into crystalline low-
dimensional materials in a hyperuniform manner through the introduction of
topological defects may enhance electronic transport in such materials Zheng
et al. (2020); Chen et al. (2021); Zheng et al. (2021), which is complementary
to the conventional wisdom of the landmark “Anderson localization” Anderson
(1958) that disorder generally diminishes electronic transport.

While it is well known that effective long-ranged interactions are required to
drive an equilibrium many-particle system to a hyperuniform state, this
condition is not necessary to achieve hyperuniformity in systems out of
equilibrium Torquato (2018). Among the wide spectrum of hyperuniform
nonequilibrium systems discovered previously, many fall into the category of
inherent structures, i.e., local potential-energy minima associated with
certain forms of interactions Torquato (2018). For instance, a variety of
maximally-random-jammed (MRJ) hard-particle packings Zachary et al. (2011a);
Zachary and Torquato (2011); Zachary et al. (2011b, c); Atkinson et al.
(2013); Chen et al. (2014) are demonstrated to be hyperuniform; since in
athermal systems increasing the density plays the same role as decreasing
temperature of a molecular liquid and MRJ packings are local density maxima,
these MRJ packings are considered inherent structures. The amorphous inherent
structures in the quantizer problem also possess a high degree of
hyperuniformity Klatt et al. (2019). Interestingly, avian photoreceptor
patterns are inherent structures associated with isotropic short-range hard-
core repulsions between any pair of cells and isotropic long-range soft-core
repulsions between pairs of cells of the subtype, and they are shown to be
multihyperuniform, i.e., the photoreceptor patterns of both the total
population and the individual cell types are simultaneously hyperuniform Jiao
et al. (2014). Another examples are the inherent structures associated with
the $k$-space overlap potentials, which are shown to be hyperuniform Batten et
al. (2011). It is also noteworthy that not all inherent structures are found
to be hyperuniform Batten et al. (2011). For example, the inherent structures
associated with the Lennard-Jones and steeply repulsive potentials are in
general not hyperuniform due to the dominance of grain boundaries and vacancy
defects Weber and Stillinger (1985); Batten et al. (2011).

Despite the ubiquitous nature of disordered hyperuniform inherent structures,
the mechanisms associated with the emergence of disordered hyperuniformity in
many such systems are still not well understood. In this work, we provide a
topological perspective to shed lights on this issue. In particular, we
consider a representative class of disordered inherent structures in two-
dimensional Euclidean space $\mathbb{R}^{2}$ which can be viewed as defected
states of perfect triangular lattice crystal Halperin and Nelson (1978);
Nelson (1983) obtained by continuously introducing topological defects such as
bound dislocations, free dislocations, and disclinations that are the key
elements in the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) two-stage
melting theory in two dimensions Kosterlitz and Thouless (1973); Young (1979);
Nelson and Halperin (1979). These defects are also commonly seen in 2D
colloidal systems Zahn et al. (1999); Wierschem and Manousakis (2011) and 2D
semiconductors Zheng et al. (2020); Chen et al. (2021) and play an important
role in determining the physical properties of such materials.

Using various structural descriptors, we demonstrate that these inherent
structures preserve the class-I hyperuniformity of the original triangular
lattice crystal. We also show that disclinations result in the strongest
“degradation” of the translational and orientational order of the crystal,
followed by free dislocations and bound dislocations at comparable defect
concentrations. The bond-orientational correlations in these structures
rapidly decay to their long-range values over a short length scale, regardless
of the defect types and concentrations. These behaviors are in stark contrast
to those observed in thermally equilibrium configurations during the 2D
melting process, which shows a two-step change in their translational and
bond-orientational order correlations as temperature increases, corresponding
to the two Kosterlitz-Thouless (KT) type transitions (i.e., solid-hexatic, and
hexatic-liquid). In addition, the structures sampled from this 2D melting
process are typically nonhyperuniform. These results highlight the importance
of decoupling the positional degrees of freedom from the vibrational degrees
of freedom and investigate inherent structures that correspond to local energy
minima of the systems (i.e., positional degrees of freedom) when looking for
disordered hyperuniformity, since the hyperuniformity property is often
cloaked by thermal fluctuations (i.e., vibrational degrees of freedom).

Moreover, we derive a continuum theory to explain the hyperuniformity-
preserving nature of the topological transformations that link the disordered
inherent structures and the original hyperuniform crystals at low defect
concentrations. We demonstrate that the large-scale density fluctuations in
these inherent structures are mainly dominated by the elastic displacement
fields resulted from the topological defects, which at low defect
concentrations can be approximated as superposition of the displacement fields
associated with each individual defect (strain source). Remarkably, the
class-I hyperuniformity of the original crystal is preserved as long as the
displacement fields of individual defects decay sufficiently fast from the
source (i.e., the volume integrals of the displacements and squared
displacements caused by individual defect are finite) and the displacement-
displacement correlation matrix of the system is diagonalized and isotropic.
In addition, the structure factor approaches zero with a universal quadratic
scaling at small wavenumbers, regardless of the types and exact concentrations
of topological defects. Our numerical results and theoretical analysis uncover
the mechanisms underlying the emergence of disordered hyperuniformity in a
wide spectrum of disordered structures, and provide insights to the discovery,
design, and generation of novel disordered hyperuniform materials.

The rest of the paper is organized as follows: in Sec II, we describe the
procedures to generate disordered inherent structures via continuous
topological transformations from the reference triangular lattice crystal
state in $\mathbb{R}^{2}$. In Sec. III, we employ various statistical
descriptors to characterize the large-scale structural features, in particular
hyperuniformity of the resulting inherent structures. In Sec. IV, we derive
continuum theory to explain the class-I hyperuniformity of the inherent
structures. In Sec. V, we provide concluding remarks.

##  II Realizations of disordered inherent structures containing randomly
distributed topological defects

###  II.1 Dislocations and disclinations induced by topological
transformations

$\begin{array}[]{c}\\\
\includegraphics[width=411.93767pt]{fig1.png}\end{array}$

Figure 1: (Color online) Illustration of the formation of inherent structures
containing bound dislocations (top panel), free dislocations (middle panel),
and disclinations with associated dislocations (bottom panel) through series
of topological transformations (i.e., rearrangement of bonding network) and
subsequent structural relaxation in a triangular lattice. Vertices with seven
bonds are highlighted with yellow circles, and vertices with five bonds
highlighted with green circles.

To introduce bound dislocations (i.e., a pair of dislocations that are next to
each other) into the triangular lattice, we first randomly pick a bond in the
lattice. Note that any bond in the triangular lattice is also the short
diagonal of a rhombus. Next, we break this chosen bond and connect the two
vertices associated with the long diagonal of the corresponding rhombus with a
new bond, resulting in a pair of dislocations next to one another. If the
vertices associated with the old and new bonds all possess six bonds before
the transformation, then the transformation would lead to two five-coordinated
vertices and two seven-coordinated vertices; otherwise, we would obtain
higher-order defected structures. Here we impose the constraint that the
vertices after the transformation should each possess at leave five bonds to
ensure local structural stability. The process of introducing a single pair of
bound dislocations is illustrated in the top panel of Fig. 1. We quantify the
amount of bound dislocations by the defect concentration defined as $p\equiv
N_{op}/N_{b}$, where $N_{op}$ is the number of successful topological
transformations, and $N_{b}$ is the number of bonds in the triangular lattice.
Note that a single topological transformation described in this paragraph
would introduce a pair of dislocations, in the absence of other topological
defects.

To generate free dislocations, we start from bound dislocations with the
additional constraint that the initial bound dislocations should consist of
two five-coordinated vertices and two seven-coordinated vertices, and randomly
pick a five-coordinated vertex and a seven-coordinated vertex that are part of
the bound dislocations. We then let these two defected vertices “glide” in the
lattice by continuously breaking existing bonds and forming new bonds. Note
that the direction that the defected vertices can “glide” is fixed once the
two vertices are picked given the local bonding constraints. To form free
dislocations and minimize the spatial correlations of the free dislocations,
we let the defects glide for at least one step; beyond that, the defects have
$1/2$ of probability to stop and $1/2$ of probability to continue gliding at
each lattice site. If the gliding defects stop before hitting any “road
block”, i.e., vertices that are not six coordinated, then we count this as one
successful topological transformation in the context of free dislocations. The
process of introducing a single pair of free dislocations is illustrated in
the middle panel of Fig. 1. We also experiment with other stopping rules, and
find that the details of different stopping rules do not affect the large-
scale structural features of the resulting structures, which is the focus of
this work. We quantify the amount of free dislocations by the defect
concentration defined as $p\equiv N_{op}/N_{b}$, where $N_{op}$ is the number
of successful topological transformations described in this paragraph, and
$N_{b}$ is the number of bonds in the triangular lattice. Note that similar to
the case of bound dislocations, a single topological transformation described
in this paragraph would introduce a pair of dislocations (each consisting of a
5-coordinated vertex and a 7-coordinated vertex), in the absence of other
topological defects.

To generate disclinations, we start from a free dislocation and break the bond
between the seven-coordinated vertex and one of its six-coordinated neighbors
and connect the long diagonal of the corresponding rhombus with a new bond.
This six-coordinated neighbor should have two six-coordinated neighbors that
are not neighbor of the five-coordinated vertex. This bond-breaking and bond-
forming process creates an isolated 5-coordinated vertex, and another
5-coordinated vertex surrounded by two 7-coordinated vertices. We then let one
of the two 7-coordinated vertices and its neighboring 5-coordinated vertex
glide away in the same way as that in the case of free dislocations. If these
steps can be completed, then we count this as a successful topological
transformation in the context of disclinations, which would create an isolated
five-fold disclination, an isolated seven-fold disclination, and two free
dislocations (each consisting of a 5-coordinated vertex and a 7-coordinated
vertex), in the absence of other topological defects. Such a topological
transformation is illustrated in the bottom panel of Fig. 1. Note that the
disclinations are accompanied by free dislocations, which is consistent with
previous observations Zahn et al. (1999); Wierschem and Manousakis (2011) that
disclinations typically arise with free dislocations. We quantify the amount
of disclinations by the defect concentration defined as $q\equiv
N_{op}/N_{b}$, where $N_{op}$ is the number of successful topological
transformations described in this paragraph, and $N_{b}$ is the number of
bonds in the triangular lattice. It is noteworthy that structures containing
disclinations at $q$ should be compared to structures containing bound and
free dislocations at $p=\frac{3}{2}q$ for a fair comparison at the same
effective defect concentration, given the number of 5-coordinated and
7-coordinated vertices that each case generates in the absence of other
topological defects.

###  II.2 Inherent structures

To obtain the inherent structures, we allow the transformed structures to
undergo elastic relaxation by perturbing the positions of the vertices in a
way that drive the bond lengths in the network towards values associated with
the triangular lattice. In particular, this involves local minimization of the
energy function $E$ defined as follows:

$E=\sum_{\mathrm{bonds}}k_{b}(b_{i}-b_{0})^{2},$ (1)

where $b_{i}$ is the bond length associated with bond $i$, and $b_{0}=1$ is
the side length of a triangle in a triangular lattice. Since we are looking at
local energy minima, the choice of the spring constant $k_{b}$ does not affect
the obtained structure, and without loss of generality, we set $k_{b}$ to
unity.

$\begin{array}[]{c}\\\
\includegraphics[width=411.93767pt]{fig2.png}\end{array}$

Figure 2: (Color online) (Top section) Representative inherent structures
containing primarily bound dislocations at defect concentration $p$=0.02,
0.06, 0.10, and 0.17, respectively. (Bottom left section) Representative
inherent structures containing primarily free dislocations at defect
concentration $p$=0.02, and 0.04, respectively. (Bottom right section)
Representative inherent structures containing disclinations at defect
concentration $q$=0.01, and 0.015, respectively. Vertices with more than six
bonds are highlighted in yellow, and vertices with less than six bonds in
green.

We investigate inherent structures containing the aforementioned three types
of topological defects for a wide range of $p$ (or $q$). In the cases of free
dislocations and disclinations, we generate structures up to values close to
saturation, i.e., more topological defects can no longer be inserted into the
system after a sufficiently large number of attempts (e.g., $10N$ attempts,
where $N$ is the number of vertices in the lattice). In the case of bound
dislocations, we stop at $p=0.17$ since increasing $p$ beyond that sometimes
leads to unphysical local bonding networks. In the top section of Fig. 2 we
show representative inherent structures containing primarily bound
dislocations with $N=1200$ particles, in the bottom left section of Fig. 2
representative inherent structures containing primarily free dislocations with
$N=1200$ particles, and in the bottom right section of Fig. 2 representative
inherent structures containing disclinations with $N=1200$ particles.

##  III Structural characterization and hyperuniformity

To characterize the inherent structures with the aforementioned topological
defects, in particular at large length scales, we generate configurations with
$N=10,800$ particles at different $p$ (or $q$) and look at various statistics
including pair statistics such $g_{2}(r)$, $S(k)$ and $\sigma_{N}^{2}(R)$
Torquato and Stillinger (2003); Torquato (2018), and bond-orientational
statistics such as the bond-orientational order metric $Q_{6}$ and correlation
function $C_{6}(r)$ that have been routinely used to study the 2D melting
process Zahn et al. (1999); Wierschem and Manousakis (2011); Li and Ciamarra
(2018). To compute the statistics accurately, we average over 10
configurations at each $p$ (or $q$).

Specifically, the pair correlation function $g_{2}(r)$ is proportional to the
probability density function of finding two centers separated by distance $r$
Torquato (2002), and in practice is computed via the relation

$g_{2}(r)=\frac{\langle N(r)\rangle}{\rho 2\pi r\Delta r},$ (2)

where $\langle N(r)\rangle$ is the average number of particle centers that
fall into the circular ring at distance $r$ from a central particle center
(arbitrarily selected and averaged over all particle centers in the system),
$2\pi r\Delta r$ is the area of the circular ring, and $\rho$ is the number
density of the system Torquato (2002); Atkinson et al. (2013). The static
structure factor $S(k)$ is the Fourier counterpart, and for computational
purposes, $S({k})$ is the angular-averaged version of $S({\bf k})$, which can
be obtained directly from the particle positions ${\bf r}_{j}$, i.e.,

$S({\bf k})=\frac{1}{N}\left|{\sum_{j=1}^{N}\exp(i{\bf k}\cdot{\bf
r}_{j})}\right|^{2}\quad({\bf k}\neq{\bf 0}),$ (3)

where $N$ is the total number of points in the system Zachary and Torquato
(2009); Atkinson et al. (2013); Chen et al. (2014). The trivial forward
scattering contribution (${\bf k}=0$) in Eq. 3 is omitted, which makes Eq. 3
completely consistent with the aforementioned definition of $S(k)$ in the
ergodic infinite-system limit Torquato (2018). To compute $\sigma_{N}^{2}(R)$,
we randomly place circular observation windows with radius $R$ in the system,
and count the number of particles $N(R)$ that fall into the observation
window, which is a random variable. The variance associated with $N(R)$ is
denoted by $\sigma_{N}^{2}(R)\equiv\langle N(R)^{2}\rangle-\langle
N(R)\rangle^{2}$, which measure density fluctuations of particles within a
window of radius $R$. In this work we sample 100,000 windows at each window
radius $R$ to obtain $\sigma_{N}^{2}(R)$.

On the other hand, the order metric $Q_{6}$ is defined as

$Q_{6}\equiv|\langle\Psi_{6}\rangle|,$ (4)

where

$\Psi_{6}({\bf r}_{i})=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}e^{6\theta_{ij}},$ (5)

and $\langle\cdots\rangle$ denotes ensemble average, $n_{i}$ is the number of
neighbors of vertex $i$ located at ${\bf r}_{i}$, and $\theta_{ij}$ is the
polar angle associated with the vector from vertex $i$ to the $j$-th bonded
neighbor of vertex $i$. The bond-orientational correlation function $C_{6}(r)$
is defined as

$C_{6}(r)\equiv\langle\Psi_{6}({\bf r}_{i})\Psi^{*}_{6}({\bf
r}_{j})\rangle\mid r=|{\bf r}_{i}-{\bf r}_{j}|,$ (6)

where $\Psi^{*}_{6}$ is the complex conjugate of $\Psi_{6}$. In practice, to
calculate $C_{6}(r)$, for each pair of particles located at ${\bf r}_{i}$ and
${\bf r}_{j}$, we compute $\Psi_{6}({\bf r}_{i})\Psi^{*}_{6}({\bf r}_{j})$,
and bin the results according to the distance $r=|{\bf r}_{i}-{\bf r}_{j}|$.
We note that $Q_{6}=1$ and $C_{6}(r)=1$ for a perfect triangular network;
while for isotropic fluid phase, $Q_{6}\approx 0$ and $C_{6}(r)$ decays with
an exponential envelop at large $r$ Zahn et al. (1999); Wierschem and
Manousakis (2011).

###  III.1 Bound dislocations

We first present the results of the inherent structures with primarily bound
dislocations at different $p$. In particular, as shown in Fig. 3(a)-(c), the
structure factor $S(k)$ decreases to essentially zero as $k$ approaches zero
and local number variance $\sigma_{N}^{2}(R)$ grows roughly linearly as $R$
increases at large $R$, indicating the hyperuniformity of these inherent
structures. The pair correlation function $g_{2}(r)$ decays to its long-range
value over a short range of $r$, and the long-range value of $|g_{2}(r)-1|$
and the magnitudes of the Bragg peaks in $S(k)$ also decrease significantly as
$p$ increases, indicating the possible loss of translational order in these
systems as dislocations are introduced into the system. However, we note that
the absence of Bragg peaks alone does not guarantee that the underlying
structure is truly amorphous, since long-range order can be hidden at the two-
point level, but still can be present at higher-point levels, as explicitly
demonstrated by Klatt et al. in the context of random, uncorrelated
displacements of particles on a lattice Klatt et al. (2020). There are clear
wiggles in $S(k)$ at large $k$ and significant oscillations in $g_{2}(r)$ as
well, which are manifestations of the remaining short-range structures in the
defected networks. We further analyze the small-wavenumber behavior of $S(k)$,
and find that the exponent $\alpha$ in $S(k)\sim k^{\alpha}$ oscillates around
2, as shown in Fig. 3(d), demonstrating that bound dislocations preserve the
class-I hyperuniformity of the triangular lattice.

We also analyze the bond-orientional order of the inherent structures. The
results of $Q_{6}$ and $C_{6}(r)$ are shown in Fig. 3(e) and 3(f),
respectively. It can be clearly seen that $Q_{6}$ decreases rapidly as $p$
increases, indicating the loss of the global preferred orientation of the
lattice. On the other hand, $C_{6}(r)$ decays to its long-range value rapidly
over a short length scale regardless of $p$, which can be attributed to the
fact that the bound dislocations are randomly introduced in the system, and
the spatial correlations of defect positions are minimized. The long-range
value of $C_{6}(r)$ also decreases as $p$ increases, indicating the loss of
large-scale orientational correlation as bound dislocations are introduced.

$\begin{array}[]{c}\\\
\includegraphics[width=411.93767pt]{fig3.png}\end{array}$

Figure 3: (Color online) Statistics associated with inherent structures
containing primarily bound dislocations at different defect concentrations $p$
with N=10,800 particles. (a) Structure factor $S(k)$. (b) Local number
variance $\sigma_{N}^{2}(R)$. (c) Log-log plot of $|g_{2}(r)-1|$. (d) Small-
wavenumber scaling exponent $\alpha$ of $S(k)$. (e) Bond-orientational order
metric $Q_{6}$. (f) Bond-orientational order correlation function $C_{6}(r)$.

###  III.2 Free dislocations

We employ similar procedures to investigate inherent structures containing
primarily free dislocations, and the computed statistics are shown in Fig. 4.
Interestingly, there are many structural similarities that these inherent
structures share with inherent structures containing primarily bound
dislocations. For example, these inherent structures also preserve the class-I
hyperuniformity of the triangular lattice, as manifested by the fact that
$S(k)$ essentially decreases to zero with an approximately quadratic scaling
as $k$ approaches zero and $\sigma_{N}^{2}(R)$ increases linearly as $R$
increases at large $R$. Both $g_{2}(r)$ and $C_{6}(r)$ decay to their
respective long-range values over a short length scale, and $Q_{6}$ decreases
rapidly as $p$ increases, indicating the loss of large-scale structural order
as free dislocations are introduced into the system. However, we note that at
the same $p$, free dislocations degrade the translational and orientational
order of the triangular lattice much more than bound dislocations, as
evidenced by $g_{2}(r)$, $Q_{6}$, and $C_{6}(r)$. This is not surprising that
the impact of bound dislocations are much more localized than that of free
dislocations.

It is noteworthy that in colloidal systems during 2D melting, as free
dislocations begin to emerge, the systems start to enter the hexatic phase
regime, and the $h(r)$ and $C_{6}(r)$ typically show an exponential and an
algebraic decay, respectively Zahn et al. (1999); Wierschem and Manousakis
(2011). These behaviors are distinctly different from those of our inherent
structures containing primarily free dislocations, where $h(r)$ and $C_{6}(r)$
decay to their respective long-range values over a short range of $r$ and
oscillate around certain constants afterwards. These differences may be
attributed to the fact that in our systems, the free dislocations are
introduced in a mostly uncorrelated manner, while in those colloidal systems
during 2D melting, the free dislocations arise as a result of thermal
excitation and possess certain degrees of spatial correlation. Our results
suggest that not only the types of topological defects, but also the spatial
correlation of topological defects affect the structural behaviors of the
defected lattices.

$\begin{array}[]{c}\\\
\includegraphics[width=411.93767pt]{fig4.png}\end{array}$

Figure 4: (Color online) Statistics associated with inherent structures
containing primarily free dislocations at different defect concentrations $p$
with N=10,800 particles. (a) Structure factor $S(k)$. (b) Local number
variance $\sigma_{N}^{2}(R)$. (c) Log-log plot of $|g_{2}(r)-1|$. (d) Small-
wavenumber scaling exponent $\alpha$ of $S(k)$. (e) Bond-orientational order
metric $Q_{6}$. (f) Bond-orientational order correlation function $C_{6}(r)$.

###  III.3 Disclinations

Next, we investigate inherent structures containing primarily isolated
disclinations, and the results are shown in Fig. 5. Clearly, the inherent
structures preserve class-I hyperuniformity of the triangular lattice, and the
translational and orientional order of the system are greatly degraded by the
introduced disclinations. Remarkably, both $h(r)$ and $C_{6}(r)$ decay to
their respective long-range values over a short range of $r$ and oscillate
around certain constants afterwards. These results are surprising since
previously isolated disclinations caused by thermal excitation as temperature
increases in colloidal systems were known to induce large-scale structural
distortions, and lead the systems to transition into isotropic liquids, which
essentially lose translational and orientational order, and are generally not
hyperuniform Zahn et al. (1999); Wierschem and Manousakis (2011). In
particular, in those systems $h(r)$ and $C_{6}(r)$ both exhibit an exponential
decay as $r$ increases Zahn et al. (1999); Wierschem and Manousakis (2011).
These different behaviors can be attributed to the fact that in our systems
the disclinations are randomly placed into the systems, which do not affect
large-scale density fluctuations or orientational correlations. Nonetheless,
we find that in our systems disclinations degrade the translational and
orientational order much more than bound and free dislocations at comparable
defect concentration, which is not surprising given that disclinations cause
larger-scale structural distortions than dislocations.

$\begin{array}[]{c}\\\
\includegraphics[width=411.93767pt]{fig5.png}\end{array}$

Figure 5: (Color online) Statistics associated with inherent structures
containing disclinations at different defect concentrations $p$ with N=10,800
particles. (a) Structure factor $S(k)$. (b) Local number variance
$\sigma_{N}^{2}(R)$. (c) Log-log plot of $|g_{2}(r)-1|$. (d) Small-wavenumber
scaling exponent $\alpha$ of $S(k)$. (e) Bond-orientational order metric
$Q_{6}$. (f) Bond-orientational order correlation function $C_{6}(r)$.

##  IV Continuum theory of hyperuniformity in inherent structures containing
topological defects

In this section, we devise a continuum theory to explain our observations from
Sec. III of the impact of the topological defects on hyperuniformity, i.e.,
how the topological transformations involving dislocations and disclinations
preserve the class-I hyperuniformity of the original triangular-lattice
crystal? We note that the introduction of topological defects preserves the
total number of particles in the system, i.e., no particles were removed or
added. Therefore, the impacts on the local number density fluctuations are
resulted from the perturbation of particle positions at the core of the
defects and the associated elastic displacement field. Specifically, we assume
that the particle displacement (at low defect concentrations) at position
$\mathbf{x}$ is the linear superposition of the displacements introduced by
different topological defects at $\mathbf{r}_{1}$, $\cdots$, $\mathbf{r}_{M}$,
where $M$ is the number of topological defects, i.e.,

$\mathbf{u}(\mathbf{x})=\sum_{i=1}^{M}\mathbf{f}(\mathbf{x}-\mathbf{r}_{i}).$
(7)

Therefore, the average displacement field
$\langle\mathbf{u}(\mathbf{x})\rangle$ is given by

$\begin{split}\langle\mathbf{u}(\mathbf{x})\rangle&=\int\sum_{i=1}^{M}\mathbf{f}(\mathbf{x}-\mathbf{r}_{i})P_{M}(\mathbf{r}^{M})d\mathbf{r}^{M}\\\
&=\int\mathbf{f}(\mathbf{x}-\mathbf{r}_{1})\rho_{1s}(\mathbf{r}_{1})d\mathbf{r}_{1}\\\
&=\rho_{s}\int\mathbf{f}(\mathbf{r})d\mathbf{r},\end{split}$ (8)

where $P_{M}(\mathbf{r}^{M})$ is the probability density function Torquato and
Stillinger (2003) associated with finding defects $1,2,\cdots,M$ at position
$\mathbf{r}_{1}$, $\mathbf{r}_{2}$, $\cdots$, $\mathbf{r}_{M}$, and
$\rho_{ms}(\mathbf{r}^{m})(m<M)$ is the reduced generic density function
Torquato and Stillinger (2003) of the defects defined as

$\rho_{ms}(\mathbf{r}^{m})=\frac{M!}{(M-m)!}\int\cdots\int
P_{M}(\mathbf{r}^{M})d\mathbf{r}^{M-m},$ (9)

and because of statistical homogeneity, the one-point density function
$\rho_{1s}(\mathbf{r}_{1})$ is equal to the average defect density $\rho_{s}$
in the system. Similarly, the different components of the displacement-
displacement correlation
$\mathbf{\Psi}_{\mu\nu}(\mathbf{r}=\mathbf{y}-\mathbf{x})\equiv\langle
u_{\mu}(\mathbf{x})u_{\nu}(\mathbf{y})\rangle-\langle
u_{\mu}(\mathbf{x})\rangle\langle u_{\nu}(\mathbf{y})\rangle$ are given by

$\begin{split}\Psi_{\mu\nu}(\mathbf{r})&=\int\sum_{i=1}^{M}\sum_{j\neq
i}^{M}f_{\mu}(\mathbf{x}-\mathbf{r}_{i})f_{\nu}(\mathbf{y}-\mathbf{r}_{j})P_{M}(\mathbf{r}^{M})d\mathbf{r}^{M}\\\
&+\int\sum_{i=1}^{M}f_{\mu}(\mathbf{x}-\mathbf{r}_{i})f_{\nu}(\mathbf{y}-\mathbf{r}_{i})P_{M}(\mathbf{r}^{M})d\mathbf{r}^{M}\\\
&-\int\rho_{s}^{2}f_{\mu}(\mathbf{x}-\mathbf{r}_{1})f_{\nu}(\mathbf{y}-\mathbf{r}_{2})d\mathbf{r}_{1}d\mathbf{r}_{2}\\\
&=\int\rho_{s}^{2}h_{s}(\mathbf{r}_{2}-\mathbf{r}_{1})f_{\mu}(\mathbf{x}-\mathbf{r}_{1})f_{\nu}(\mathbf{y}-\mathbf{r}_{2})d\mathbf{r}_{1}d\mathbf{r}_{2}\\\
&+\int\rho_{s}f_{\mu}(\mathbf{x}-\mathbf{r}_{1})f_{\nu}(\mathbf{y}-\mathbf{r}_{1})d\mathbf{r}_{1},\end{split}$
(10)

where $h_{s}(\mathbf{r})\equiv
g_{2s}(\mathbf{r})-1=[\rho_{2s}(\mathbf{r})-\rho_{s}^{2}]/\rho_{s}^{2}$ is the
total correlation function of the topological defects. If the topological
defects are randomly introduced into the system, then $h_{s}(\mathbf{r})=0$,
which gives

$\begin{split}\Psi_{\mu\nu}(\mathbf{r})&=\int\rho_{s}f_{\mu}(\mathbf{x}-\mathbf{r}_{1})f_{\nu}(\mathbf{y}-\mathbf{r}_{1})d\mathbf{r}_{1}\\\
&=\int\rho_{s}f_{\mu}(\mathbf{r}_{1})f_{\nu}(\mathbf{r}_{1}+\mathbf{r})d\mathbf{r}_{1}\end{split}$
(11)

In the Fourier space, this corresponds to

$\tilde{\Psi}_{\mu\nu}(\mathbf{k})=\rho_{s}\tilde{f}_{\mu}(\mathbf{k})\tilde{f}_{\nu}^{*}(\mathbf{k})=\rho_{s}\tilde{f}_{\mu}(\mathbf{k})\tilde{f}_{\nu}(-\mathbf{k}),$
(12)

where $\tilde{\mathbf{\Psi}}$, $\tilde{\mathbf{f}}$ are the Fourier transforms
of $\mathbf{\Psi}$ and $\mathbf{f}$, respectively, and
$\tilde{\mathbf{f}}^{*}$ is the complex conjugate of $\tilde{\mathbf{f}}$.

Next, we derive the expression for the structure factor $S(k)$ of triangular
lattice affected by displacement fields $\mathbf{u}$. Previously, it was shown
Kim and Torquato (2018) in general that when the displacement field has a
finite variance $\langle|\mathbf{u}|^{2}\rangle$ and the displacement-
displacement correlation matrix is isotropic and diagonalized, i.e.,
$\Psi_{\mu\nu}(\mathbf{r})=\delta_{\mu\nu}\Psi(\mathbf{r})$, the structure
factor $S(\mathbf{k})$ of a displaced hyperuniform point pattern at small
$|\mathbf{k}|$ is approximated by

$\begin{split}S(\mathbf{k})&\approx[|\mathbf{k}|^{2}\Psi(\mathbf{0})+(1-|\mathbf{k}|^{2}\Psi(\mathbf{0}))S_{0}(\mathbf{k})]\\\
&+\rho_{0}|\mathbf{k}|^{2}(\tilde{\Psi}(\mathbf{k})+\int
d(\mathbf{r})h_{0}(\mathbf{r})\Psi(\mathbf{r})e^{-i\mathbf{k}\cdot\mathbf{r}})\end{split}$
(13)

where $S_{0}(\mathbf{k})$ and $h_{0}(\mathbf{r})$ are the structure factor and
total correlation function of the original point patterns. In the cases where
the original point patterns are crystals, we can use the properties of
crystals to simplify Eq. 13. In particular, the structure factor
$S_{0}(\mathbf{k})=0$ holds for $|\mathbf{k}|<K$, where $K$ is the wavenumber
associated with the first Bragg peaks, and pair correlation function
$g_{20}(\mathbf{r})=h_{0}(\mathbf{r})+1$ is simply a collection of
$\delta$-functions at lattice sites and zero otherwise. By Taylor-expanding
the second line of Eq. 13 at small $k$ and invoking the continuum
approximation, we obtain the following expression:

$\begin{split}S(k)&\approx k^{2}\Psi(0)+\rho_{0}k^{2}\tilde{\Psi}(0)\\\
&=\rho_{s}k^{2}\int
f_{1}^{2}(\mathbf{r})d\mathbf{r}+\rho_{0}\rho_{s}k^{2}|\tilde{f}_{1}(0)|^{2}\\\
&=\rho_{s}k^{2}\int
f_{2}^{2}(\mathbf{r})d\mathbf{r}+\rho_{0}\rho_{s}k^{2}|\tilde{f}_{2}(0)|^{2}\end{split}$
(14)

Interestingly, Eq. 14 suggests that as long as the volume integrals of
$\mathbf{f}(\mathbf{r})$ and $|\mathbf{f}(\mathbf{r})|^{2}$ are finite, the
structure factor $S(k)$ of the triangular lattice affected by the displacement
fields generated by the collection of randomly distributed source functions
$\mathbf{f}$ scales as $S(k)\sim k^{2}$, which indicates that such
displacement fields preserve the class-I hyperuniformity of the original
crystals.

Subsequently, we test our continuum theory against numerical examples of
inherent structures investigated in Sec. III. We first check whether the
assumptions of our theory are satisfied in these cases. In Fig. 6(a) we
visualize the magnitude of the displacement field $\mathbf{u}$ in an inherent
structure containing a single pair of bound dislocations. Clearly, the
displacement field concentrates around the center of the topological defect,
i.e., the center of the old broken bond (which is the same as the center of
the new formed bond). We further compute the decay of $|\mathbf{u}(r)|$ as a
function of the distance $r$ from the core of the topological defect, which
appears to decay exponentially as shown in Fig. 6(b), although we note that
$|\mathbf{u}(\mathbf{r})|$ appears to be anisotropic around the core as shown
in Fig. 6(a). The exponential decay suggests that the volume integrals of the
source $\mathbf{f}$ and $|\mathbf{f}|^{2}$ should be finite.

We then look at the cases where a substantial amount of topological defects
are affecting the structures. Specifically, in Fig. 7 we show the spatial
distribution of the vector displacement field $\mathbf{u}(\mathbf{r})$ for
three representative examples: an inherent structure containing primarily
bound dislocations at $p=0.17$, an inherent structure containing primarily
free dislocations at $p=0.04$, and an inherent structure containing
disclinations at $q=0.015$. All of the three fields in Fig. 7 appear to be
approximately isotropic, suggesting a finite $\Psi(0)$. We also compute the
different components of displacement-displacement correlation matrix $\Psi(r)$
for all these three cases, and the results are shown in Fig. 8. The orthogonal
components $\Psi_{12}$ appear to vanish, and the diagonal components are
roughly the same, i.e., $\Psi_{11}(r)\approx\Psi_{22}(r)$, which satisfies the
condition $\Psi_{\mu\nu}(\mathbf{r})=\delta_{\mu\nu}\Psi(\mathbf{r})$ in our
theory. In addition, the diagonal components $\Psi_{11}$ and $\Psi_{22}$ in
Fig. 8 show relatively short-ranged correlations, which are consistent with
the visualizations in Fig. 7.

With the assumptions in our theory largely satisfied by our numerical examples
in Sec. III, we proceed to investigate whether the numerically determined
small-$k$ behavior of $S(k)$ matches the prediction by our theory. Indeed, as
shown in Figs. 3-5, at low and intermediate defect concentrations $p$ and $q$,
the scaling exponent $\alpha$ in $S(k)\sim k^{\alpha}$ at small $k$ oscillates
around 2, matching the quadratic scaling predicted by our theory; however, at
very small or large $p$ and $q$ (relative to saturation), the exponent
$\alpha$ slightly deviates from 2. The deviation at very small $p$ or $q$ can
be attributed to the fact that at these defect concentrations the systems are
not entirely homogeneous, which degrades the accuracy of our continuum-theory
prediction; on the other hand, at large $p$ or $q$, the topological defects
begin to interact with each other and modify the inherent structures
accordingly, which is not taken into account in our continuum theory.

$\begin{array}[]{c}\\\
\includegraphics[width=195.12767pt]{fig6.png}\end{array}$

Figure 6: (Color online) Displacement field $\mathbf{u}$ in an inherent
structure containing a single pair of bound dislocations. (a) Spatial
distribution of normalized $|u(r)|/|u(r)|_{max}$. (b) Decay of
$|\mathbf{u}(r)|$ as a function of the distance $r$ from the core of the
topological defect, i.e., the center of the old broken bond (the same as the
center of the new formed bond), and the red solid line is an exponential fit
(note that the vertical axis is log-scale).

$\begin{array}[]{c}\\\
\includegraphics[width=346.89731pt]{fig7.png}\end{array}$

Figure 7: (Color online) (a) Spatial distribution of the vector displacement
field $\mathbf{u}(\mathbf{r})$ in an inherent structure containing primarily
bound dislocations at $p=0.17$ with $N=1200$ particles. (b) Spatial
distribution of the vector displacement field $\mathbf{u}(\mathbf{r})$ in an
inherent structure containing primarily free dislocations at $p=0.04$ with
$N=1200$ particles. (c) Spatial distribution of the vector displacement field
$\mathbf{u}(\mathbf{r})$ in an inherent structure containing disclinations at
$q=0.015$ with $N=1200$ particles.

$\begin{array}[]{c}\\\
\includegraphics[width=390.25534pt]{fig8.png}\end{array}$

Figure 8: (Color online) (a) Different components of the displacement-
displacement correlation matrix $\mathbf{\Psi}(r)$ of an inherent structure
containing primarily bound dislocations at $p=0.17$ with $N=10,800$ particles.
(b) Different components of the displacement-displacement correlation matrix
$\mathbf{\Psi}(r)$ of an inherent structure containing primarily free
dislocations at $p=0.04$ with $N=10,800$ particles. (c) Different components
of the displacement-displacement correlation matrix $\mathbf{\Psi}(r)$ of an
inherent structure containing disclinations at $q=0.015$ with $N=10,800$
particles.

##  V Conclusions and Discussion

In this work, we made an attempt to elucidate a possible mechanism for the
observed hyperuniformity in disordered inherent structures of a wide spectrum
of systems. In particular, we considered a representative class of disordered
inherent structures which are linked to an original crystal state via
continuous topological transformations involving dislocations and
disclinations. We show via both numerical simulations and theoretical analysis
that these topological transformations preserve the class-I hyperuniformity of
the triangular lattice, and the structure factor $S(k)$ possesses a universal
quadratic scaling as $k$ decreases at small $k$ at low defect concentrations.

Our continuum theory connects the large-scale density fluctuations in these
inherent structures to the elastic displacement fields resulted from the
topological defects. It indicates that class-I hyperuniformity can be
preserved as long as the displacement fields resulted from individual defects
decay fast enough from the source (i.e., the volume integrals of the
displacements and squared displacements caused by individual defect are
finite) and the displacement-displacement correlation matrix of the system is
diagonalized and isotropic. Conceptually, the introduction of topological
defects into a crystal does not affect the average particle density of the
system (since the total number of particles are conserved), and any change in
density fluctuations of the resulting disordered inherent structures could
only come from the elastic displacement fields caused by the topological
defects. As long as these elastic fields are homogenized and sufficiently
localized, the salient features of large-scale density fluctuations of the
original crystal, i.e., hyperuniformity, should be preserved. These results
suggest promising new venues for the discovery, design, and generation of
novel disordered hyperuniform materials.

Moreover, the inherent structures containing dislocations and disclinations
studied in this work are quite different from the equilibrium structures
containing the same type of topological defects in colloidal systems during 2D
melting Zahn et al. (1999); Wierschem and Manousakis (2011) in terms of
various structural features, in particular the hyperuniformity and the
large-$r$ scaling behavior of $g_{2}(r)$ and $C_{6}(r)$. These differences
suggest that not only the types of defects, but also the spatial correlations
of defects are key to understand the impact of defects on the structural
features of crystalline systems. These results also highlights the cloaking
effect of thermal fluctuation on hyperuniformity, and indicates that when
looking for disordered hyperuniformity, in many cases one should probably look
at potential-energy minima (local or global), which are only functions of the
positional degrees of freedom and not affected by the vibrational degrees of
freedom.

Here we studied the introduction of topological defects into triangular
lattice, but given the duality of triangular lattice and honeycomb lattice,
the dual of the various inherent structures obtained in this work are similar
to the network structures obtained previously Chen et al. (2021) that are
topologically transformed from the honeycomb lattice through the introduction
of Stone-Wales (SW) defects. Interestingly, those defected honeycomb network
structures were found to capture the salient structural features of amorphous
graphene and other 2D materials Chen et al. (2021). However, we note that
previously it was demonstrated that as the SW defect concentration reaches
certain critical value, the corresponding systems changes from class-I
hyperuniformity to class-II hyperuniformity, which can be attributed to the
fact that in those systems the bond angles were also regulated because of the
underlying chemical constraints Chen et al. (2021) and the additional coupling
between defects likely modifies the behavior of large-scale density
fluctuations and leads to the change in the class of hyperuniformity at
sufficiently large defect concentrations.

It is also noteworthy that topological defects such as dislocations and
disclinations appear as point defects in two dimensions, while in three
dimensions they are known to appear as line defects. Given this difference, it
would be interesting to explore the impact of topological defects on the
large-scale structural features of crystals in three dimensions, and the
potential findings may shed light on the emergence of hyperuniformity in
certain disordered inherent structures. For example, there were theories
Nelson (1983) suggesting that MRJ packings (or random close packings as termed
by many experimentalists) might be considered as disordered topological
variants of the tetrahedral particle packings, just like that Frank-Kasper
phases Nelson (1983); Reddy et al. (2018); Barbon et al. (2020) are considered
as ordered topological variants.

###### Acknowledgements.

We are grateful for Dr. Jaeuk Kim for very helpful discussion.

## References

  * Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003).
  * Torquato (2018) S. Torquato, Phys. Rep. 745, 1 (2018).
  * Zachary and Torquato (2009) C. E. Zachary and S. Torquato, J. Stat. Mech. Theor. Exp. 2009, P12015 (2009).
  * Ornstein and Zernike (1914) L. S. Ornstein and F. Zernike, Proc. Akad. Sci. 17, 793 (1914).
  * Gabrielli et al. (2002) A. Gabrielli, M. Joyce, and F. S. Labini, Phys. Rev. D 65, 083523 (2002).
  * Donev et al. (2005) A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 95, 090604 (2005).
  * Zachary et al. (2011a) C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. Lett. 106, 178001 (2011a).
  * Jiao and Torquato (2011) Y. Jiao and S. Torquato, Phys. Rev. E 84, 041309 (2011).
  * Chen et al. (2014) D. Chen, Y. Jiao, and S. Torquato, J. Phys. Chem. B 118, 7981 (2014).
  * Zachary and Torquato (2011) C. E. Zachary and S. Torquato, Phys. Rev. E 83, 051133 (2011).
  * Torquato et al. (2015) S. Torquato, G. Zhang, and F. H. Stillinger, Phys. Rev. X 5, 021020 (2015).
  * Uche et al. (2004) O. U. Uche, F. H. Stillinger, and S. Torquato, Phys. Rev. E 70, 046122 (2004).
  * Batten et al. (2008) R. D. Batten, F. H. Stillinger, and S. Torquato, J. Appl. Phys. 104, 033504 (2008).
  * Batten et al. (2009) R. D. Batten, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 103, 050602 (2009).
  * Lebowitz (1983) J. L. Lebowitz, Phys. Rev. A 27, 1491 (1983).
  * Zhang et al. (2015a) G. Zhang, F. H. Stillinger, and S. Torquato, Phys. Rev. E 92, 022119 (2015a).
  * Zhang et al. (2015b) G. Zhang, F. H. Stillinger, and S. Torquato, Phys. Rev. E 92, 022120 (2015b).
  * Kurita and Weeks (2011) R. Kurita and E. R. Weeks, Phys. Rev. E 84, 030401(R) (2011).
  * Hunter and Weeks (2012) G. L. Hunter and E. R. Weeks, Rep. Prog. Phys. 75, 066501 (2012).
  * Dreyfus et al. (2015) R. Dreyfus, Y. Xu, T. Still, L. A. Hough, A. G. Yodh, and S. Torquato, Phys. Rev. E 91, 012302 (2015).
  * Hexner and Levine (2015) D. Hexner and D. Levine, Phys. Rev. Lett. 114, 110602 (2015).
  * Jack et al. (2015) R. L. Jack, I. R. Thompson, and P. Sollich, Phys. Rev. Lett. 114, 060601 (2015).
  * Weijs et al. (2015) J. H. Weijs, R. Jeanneret, R. Dreyfus, and D. Bartolo, Phys. Rev. Lett. 115, 108301 (2015).
  * Torquato et al. (2008) S. Torquato, A. Scardicchio, and C. E. Zachary, J. Stat. Mech. Theor. Exp. 2008, P11019 (2008).
  * Feynman and Cohen (1956) R. P. Feynman and M. Cohen, Phys. Rev. 102, 1189 (1956).
  * Jiao et al. (2014) Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, J. C. Corbo, and S. Torquato, Phys. Rev. E 89, 022721 (2014).
  * Mayer et al. (2015) A. Mayer, V. Balasubramanian, T. Mora, and A. M. Walczak, Proc. Natl. Acad. Sci. USA 112, 5950 (2015).
  * Hejna et al. (2013) M. Hejna, P. J. Steinhardt, and S. Torquato, Phys. Rev. B 87, 245204 (2013).
  * Klatt et al. (2019) M. A. Klatt, J. Lovrić, D. Chen, S. C. Kapfer, F. M. Schaller, P. W. A. Schönhöfer, B. S. Gardiner, A. Smith, G. E. Schröder-Turk, and S. Torquato, Nat. Commun. 10, 1 (2019).
  * Lei et al. (2019) Q.-L. Lei, M. P. Ciamarra, and R. Ni, Sci. Adv. 5, eaau7423 (2019).
  * Lei and Ni (2019) Q.-L. Lei and R. Ni, Proc. Natl. Acad. Sci. U.S.A. 116, 22983 (2019).
  * Chremos and Douglas (2018) A. Chremos and J. F. Douglas, Phys. Rev. Lett. 121, 258002 (2018).
  * Rumi et al. (2019) G. Rumi, J. A. Sánchez, F. Elías, R. C. Maldonado, J. Puig, N. R. C. Bolecek, G. Nieva, M. Konczykowski, Y. Fasano, and A. B. Kolton, Phys. Rev. Res. 1, 033057 (2019).
  * Sánchez et al. (2019) J. A. Sánchez, R. C. Maldonado, N. R. C. Bolecek, G. Rumi, P. Pedrazzini, M. I. Dolz, G. Nieva, C. J. van der Beek, M. Konczykowski, C. D. Dewhurst, et al., Commun. Phys. 2, 1 (2019).
  * Sánchez et al. (2020) J. A. Sánchez, G. Rumi, R. C. Maldonado, N. R. C. Bolecek, J. Puig, P. Pedrazzini, G. Nieva, M. I. Dolz, M. Konczykowski, C. J. van der Beek, et al., Sci. Rep. 10, 1 (2020).
  * Huang et al. (2021) M. Huang, W. Hu, S. Yang, Q.-X. Liu, and H. P. Zhang, Proc. Natl. Acad. Sci. U.S.A. 118 (2021).
  * Torquato (2021) S. Torquato, Proc. Natl. Acad. Sci. U.S.A. 118 (2021).
  * Florescu et al. (2009) M. Florescu, S. Torquato, and P. J. Steinhardt, Proc. Natl. Acad. Sci. U.S.A. 106, 20658 (2009).
  * Man et al. (2013) W. Man, M. Florescu, E. P. Williamson, Y. He, S. R. Hashemizad, B. Y. C. Leung, D. R. Liner, S. Torquato, P. M. Chaikin, and P. J. Steinhardt, Proc. Natl. Acad. Sci. U.S.A. 110, 15886 (2013).
  * Zheng et al. (2020) Y. Zheng, L. Liu, H. Nan, Z.-X. Shen, G. Zhang, D. Chen, L. He, W. Xu, M. Chen, Y. Jiao, et al., Sci. Adv. 6, eaba0826 (2020).
  * Chen et al. (2021) D. Chen, Y. Zheng, L. Liu, G. Zhang, M. Chen, Y. Jiao, and H. Zhuang, Proc. Natl. Acad. Sci. U.S.A. 118, e2016862118 (2021).
  * Zheng et al. (2021) Y. Zheng, D. Chen, L. Liu, Y. Liu, M. Chen, H. Zhuang, and Y. Jiao, Phys. Rev. B 103, 245413 (2021).
  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Zachary et al. (2011b) C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. E 83, 051308 (2011b).
  * Zachary et al. (2011c) C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. E 83, 051309 (2011c).
  * Atkinson et al. (2013) S. Atkinson, F. H. Stillinger, and S. Torquato, Phys. Rev. E 88, 062208 (2013).
  * Batten et al. (2011) R. D. Batten, F. H. Stillinger, and S. Torquato, J. Chem. Phys. 135, 054104 (2011).
  * Weber and Stillinger (1985) T. A. Weber and F. H. Stillinger, Phys. Rev. B 31, 1954 (1985).
  * Halperin and Nelson (1978) B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 (1978).
  * Nelson (1983) D. R. Nelson, Phys. Rev. B 28, 5515 (1983).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 (1973).
  * Young (1979) A. Young, Phys. Rev. B 19, 1855 (1979).
  * Nelson and Halperin (1979) D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 (1979).
  * Zahn et al. (1999) K. Zahn, R. Lenke, and G. Maret, Phys. Rev. Lett. 82, 2721 (1999).
  * Wierschem and Manousakis (2011) K. Wierschem and E. Manousakis, Physical Rev. B 83, 214108 (2011).
  * Li and Ciamarra (2018) Y.-W. Li and M. P. Ciamarra, Phys. Rev. Mater. 2, 045602 (2018).
  * Torquato (2002) S. Torquato, _Random Heterogeneous Materials:
Microstructure and Macroscopic Properties_ (Springer-Verlag, New York, 2002).

  * Klatt et al. (2020) M. A. Klatt, J. Kim, and S. Torquato, Phys. Rev. E 101, 032118 (2020).
  * Kim and Torquato (2018) J. Kim and S. Torquato, Phys. Rev. B 97, 054105 (2018).
  * Reddy et al. (2018) A. Reddy, M. B. Buckley, A. Arora, F. S. Bates, K. D. Dorfman, and G. M. Grason, Proc. Natl. Acad. Sci. U.S.A. 115, 10233 (2018).
  * Barbon et al. (2020) S. M. Barbon, J.-A. Song, D. Chen, C. Zhang, J. Lequieu, K. T. Delaney, A. Anastasaki, M. Rolland, G. H. Fredrickson, M. W. Bates, et al., ACS Macro Lett. 9, 1745 (2020).