# Equilibrium States Corresponding to Targeted Hyperuniform Nonequilibrium
Pair Statistics

Haina Wang  _Department of Chemistry, Princeton University_ , Princeton, New
Jersey, 08544, USA  Salvatore Torquato  Email<EMAIL_ADDRESS>Department of Chemistry, Department of Physics, Princeton Institute of
Materials, and Program in Applied and Computational Mathematics, Princeton
University, Princeton, New Jersey 08544, USA  School of Natural Sciences,
Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA

###### Abstract

The Zhang-Torquato conjecture [Phys. Rev. E 101, 032124 (2020)] states that
any realizable pair correlation function $g_{2}({\bf r})$ or structure factor
$S({\bf k})$ of a translationally invariant nonequilibrium system can be
attained by an equilibrium ensemble involving only (up to) effective two-body
interactions. To further test and study this conjecture, we consider two
singular nonequilibrium models of recent interest that also have the exotic
hyperuniformity property: a 2D “perfect glass” and a 3D critical absorbing-
state model. We find that each nonequilibrium target can be achieved
accurately by equilibrium states with effective one- and two-body potentials,
lending further support to the conjecture. To characterize the structural
degeneracy of such a nonequilibrium-equilibrium correspondence, we compute
higher-order statistics for both models, as well as those for a hyperuniform
3D uniformly randomized lattice (URL), whose higher-order statistics can be
very precisely ascertained. Interestingly, we find that the differences in the
higher-order statistics between nonequilibrium and equilibrium systems with
matching pair statistics, as measured by the “hole” probability distribution,
provides measures of the degree to which a system is out of equilibrium. We
show that all three systems studied possess the bounded-hole property, and
that holes near the maximum hole size in the URL are much rarer than those in
the underlying simple cubic lattice. Remarkably, upon quenching, the effective
potentials for all three systems possess local energy minima (i.e., inherent
structures) with stronger forms of hyperuniformity compared to their target
counterparts. Our methods are expected to facilitate the self-assembly of
tunable hyperuniform soft-matter systems.

##  I Introduction

Probing and characterizing structural properties of many-body systems in and
out of equilibrium is a crucial task in the understanding of a large variety
of physical, chemical and biological phenomena. [1, 2, 3, 4, 5] An outstanding
challenge is the determination of effective interactions in many-body systems
that accurately yield equilibrium states with prescribed pair statistics.
Solving such inverse problems is a powerful way to tackle the unsolved problem
concerning the realizability of prescribed functional forms of pair statistics
by many-body systems.[6, 7, 8, 9, 10, 11] Such investigations also enable one
to probe systems with identical pair statistics but different higher-body
statistics, which is expected to shed light on the well-known degeneracy
problem of statistical mechanics. [12, 13, 14] Moreover, such effective
potentials can be used to model macromolecules and solutions,[15, 16] and to
design nanoparticles that self-assemble into desired structures, thereby
facilitating material discovery.[17, 18, 19, 20]

Recently, Zhang and Torquato conjectured that any realizable pair correlation
function $g_{2}(\mathbf{r})$ or structure factor $S(\mathbf{k})$ corresponding
to a translationally invariant nonequilibrium system can be attained by an
equilibrium ensemble involving only one- and two-body effective interactions
at positive temperatures.[11] Testing the conjecture requires the precise
determination of the effective interactions for a spectrum of target systems,
including those with the exotic hyperuniform property. [21, 22] Disordered
hyperuniform many-body systems, which can be solid and fluid states, are
unusual amorphous states of matter that lie between a crystal and liquid. They
are like perfect crystals in the way they suppress large-scale density
fluctuations, and yet are like liquids or glasses in that they are
statistically isotropic with no Bragg peaks; [21, 22] see Sec. II.2 for
precise definitions. Disordered hyperuniform states play vital roles in a
variety of different contexts. For example, disordered stealthy ground states
have been discovered [23, 24, 25] corresponding to under soft, long-ranged
interactions, which are highly degenerate and have the “bounded-hole”
property, which is a singular characteristic for a disordered system. [26, 27]
Network structures derived from disordered stealthy point patterns can achieve
complete photonic band gaps and have novel optical properties that were
previously not thought to be possible. [28, 29] Disordered hyperuniform states
also arise in the eigenvalues of random matrices (such as the Gaussian unitary
ensemble), [30, 31, 32] ground states of free fermions [31], as well as one-
component plasmas at positive temperatures, [33, 34, 35] all of which are
exotic fluid states in which the particles interact Coulombically. [36] Other
examples include including glass formation,[37, 38, 22] jamming,[39, 40, 41,
42, 43] rigidity,[38, 44] biology,[45, 46] localization of waves and
excitations,[47, 48, 49, 50] antenna or laser array designs,[51] self-
organization,[52, 53, 54] fluid dynamics,[55, 56, 57] quantum systems,[58, 59,
31, 60, 61] and pure mathematics.[62, 63, 64, 65, 66] Because disordered
hyperuniform states combine the advantages of statistical isotropy and the
suppression of density fluctuations on large scales, they can be endowed with
novel physical properties.[67, 28, 29, 37, 45, 52, 46, 38, 53, 22, 54]

While Zhang and Torquato introduced an algorithm to draw equilibrium classical
particle configurations from canonical ensembles with one- and two-body
interactions that correspond to targeted functional forms for $g_{2}({\bf r})$
or $S({\bf k})$, the algorithm does not generate explicit forms of the
potentials. [11] Very recently, Torquato and Wang developed an inverse
methodology that determines effective interactions with unprecedented
accuracy. [68] Using this procedure, they demonstrated the realizability of
$g_{2}(\mathbf{r})$ for all $r$ and $S(\mathbf{k})$ for all $k$ for two
different nonequilibrium models, including a two-dimensional (2D)
nonhyperuniform random sequential addition process and a 3D hyperuniform
“cloaked” uniformly randomized lattice (URL). [68] However, the Zhang-Torquato
conjecture remains largely untested.

In this paper, we utilize this precise inverse methodology [68] to further
test and study the Zhang-Torquato conjecture for unusual nonequilibrium
hyperuniform systems. Hyperuniform targets are particularly challenging
because at positive temperature $T$ they require a long-ranged pair
interaction $v(r)$ that must be balanced by one-body potentials to stabilize
the equilibrium system. [22]

Prior to the development of our inverse methodology,[68] predictor-corrector
methods,[69, 70, 71, 72] such as Iterative Boltzmann inversion (IBI) [71] and
iterative hypernetted chain inversion (IHNCI),[73, 72] were regarded to be the
most accurate inverse procedures. Both IBI and IHNCI begin with an initial
discretized (binned) approximation of a trial pair potential. The trial pair
potential at each binned distance is iteratively updated to attempt to reduce
the difference between the target and trial pair statistics. However, IBI and
IHNCI cannot treat long-ranged pair interactions required for hyperuniform
targets, nor do they consider one-body interactions that stabilize
hyperuniform equilibrium states [68]; see Sec. II.2 for details. These
algorithms also accumulate random errors in the binned potentials due to
simulation errors in the trial pair statistics, and thus do not achieve the
precision required to probe realizability problems. Moreover, because all
previous methods do not optimize a pair-statistic “distance” functional, they
are unable to detect poor agreement between the target and trial pair
statistics that may arise as the simulation evolves, leading to increasingly
inaccurate corresponding trial potentials, as demonstrated in Ref. 68.

The inverse methodology presented in Ref. 68 improves on previous procedures
in several significant ways. It utilizes a parameterized family of pointwise
basis functions for the potential function at $T>0$, whose initial form is
informed by small- and large-distance behaviors dictated by statistical-
mechanical theory. Pointwise potential functions do not suffer from the
accumulation of random errors during a simulation, resulting in more accurate
interactions.[68] Since it has recently been established [74] that inverse
methods that target only $g_{2}({\bf r})$ or only $S({\bf k})$ for a limited
range of $\mathbf{r}$ or $\mathbf{k}$ may generate effective potentials that
are distinctly different from the unique potential dictated by Henderson’s
theorem, [75] our methodology [68] minimizes an objective function that
incorporates both the target pair correlation function $g_{2}({\bf r})$ and
structure factor $S({\bf k})$ so that both the small- and large-distance
correlations are very accurately captured. For hyperuniform targets, our
methodology is able to optimize the required long-ranged pair potential [76]
as well as the neutralizing background one-body potential;[68] see Sec. IV for
details.

To assess the accuracy of inverse methodologies to target pair statistics, we
introduced [68] the following dimensionless $L_{2}$-norm error:

${\cal E}=\sqrt{D_{g_{2}}+D_{S}},$ (1)

where $D_{g_{2}}$ and $D_{S}$ are $L_{2}$ functions, given by

$D_{g_{2}}=\rho\int_{\mathbb{R}^{d}}[g_{2,T}(\mathbf{r})-g_{2,F}(\mathbf{r};\mathbf{a})]^{2}d\mathbf{r},$
(2)
$D_{S}=\frac{1}{\rho(2\pi)^{d}}\int_{\mathbb{R}^{d}}[S_{T}(\mathbf{k})-S_{F}(\mathbf{k};\mathbf{a})]^{2}d\mathbf{k},$
(3)

where $g_{2,F}(\mathbf{r};\mathbf{a})$ and $S_{F}(\mathbf{k};\mathbf{a})$
represent the final pair statistics at the end of the optimization, which
depend on the vector of potential parameters $\bf a$. We have previously shown
that our method is able to treat challenging near-critical and hyperuniform
targets, [68] which previous methods cannot do. Thus, it is the only available
method to determine effective interactions for nonequilibrium hyperuniform
pair statistics. Moreover, in cases where IBI and IHNCI are able to achieve
optimized potentials, e.g., for equilibrium target pair statistics without
long-range interactions, our inverse methodology generally yields $L_{2}$-norm
errors (1) that are an order of magnitude smaller than those via previous
methods, and reaches the precision required to recover the unique potential
dictated by Henderson’s theorem.[75]

We study two models of recent interest from hitherto unexplored
hyperuniformity classes, including a 2D perfect glass [38] and a 3D critical
absorbing-state model; [2, 52] see Sec. III for detailed descriptions of the
models and Sec. II.2 for the definition of hyperuniformity classes. We show
that the pair statistics of both systems can be achieved by effective
potentials, which lends further support to the Zhang-Torquato conjecture. We
define a nonequilibrium-equilibrium correspondence to be a nonequilibrium and
an equilibrium system with identical number density and pair statistics. Such
correspondences have important consequences, including a capacity to explore
the thermodynamic and dynamic properties of the effectively equivalent
equilibrium systems, such as phase behaviors, ground states and inherent
structures, i.e. local energy minima. [77] Inherent structures are properties
of the energy landscape that are used to estimate various thermodynamic and
dynamic properties, including heat capacity, melting point and glass
transition temperature. [78, 79] Furthermore, structural properties of the
equilibrium states, such as nearest-neighbor probability distribution
functions and percolation threshold, enable one to infer these nontrivial
attributes of the nonequilibrium states, which are crucial in determining
mechanical and electronic properties of materials. [80]

The aforementioned nonequilibrium-equilibrium correspondence also enables one
to probe the degeneracy of structures with the same pair statistics. [12, 13,
14] It is known that for a homogeneous many-body system, one- and two-body
correlations are insufficient to uniquely determine the higher-body
correlation functions $g_{3},g_{4},...$. [9] Equilibrium and nonequilibrium
systems that possess matching pair statistics must have different higher-body
statistics, and such differences reflect the specific dynamics of the
nonequilibrium state. [74]

To study the structural degeneracy of the nonequilibrium-equilibrium
correspondence, we compute higher-order statistics for the perfect glass and
critical absorbing-state models, as well as those for a hyperuniform 3D
cloaked URL, [81] whose higher-order statistics can be determined with much
higher precision than those of typical fluids. The effective potential for the
URL has been previously determined, but the physical properties of this
nonequilibrium-equilibrium pair have not yet been thoroughly studied. [68] We
consider the readily computable “hole” probability density functions
$E_{V}(r)$ and $G_{V}(r)$, which are related to the probability of finding
holes of radius $r$ void of particles and reflect information about all
$n$-particle correlation functions $g_{2},g_{3},g_{4},\dots$; [82] see Sec.
II.3 for precise definitions. We show that each model has the “bounded-hole”
property, [27, 44] i.e., the maximum hole size is bounded, and differences of
hole probability functions between the nonequilibrium and equilibrium systems
provides a useful nonequilibrium index. To study the behavior of the hole
probability functions on approach to the maximum hole size, we introduce a
precise algorithm to compute these functions for the target URL. We apply this
algorithm to show that holes near the maximum hole size in the target URL are
much rarer than those in the underlying simple cubic (SC) lattice.

To study the effect of quenching on the effectively equivalent equilibrium
systems, we compute the inherent structures [77] of the effective potentials
for all three models. We find that all effective potentials yield inherent
structures that are of stronger forms of hyperuniformity compared to those of
their corresponding target systems. Our findings are expected to facilitate
the self-assembly of tunable hyperuniform soft-matter systems.

We begin by providing basic definitions and background in Sec. II. In Sec.
III, we describe the model nonequilibrium hyperuniform systems. Section IV
provides a sketch our inverse methodology. [68] Section V presents results for
the equilibrium systems corresponding to the nonequilibrium pair statistics,
including the optimized potential and configurations (Sec. V.1), higher-order
statistics (Sec. V.2) and inherent structures (Sec. V.3). We provide
concluding remarks in Sec. VI.

##  II Preliminaries and definitions

###  II.1 Pair statistics

We consider many-particle systems in $d$-dimensional Euclidean space
$\mathbb{R}^{d}$ that are completely statistically characterized by the
$n$-particle probability density functions
$\rho_{n}(\mathbf{r}_{1},...,\mathbf{r}_{n})$ for all $n\geq 1$. [83] In the
case of statistically homogeneous systems, $\rho_{1}(\mathbf{r}_{1})=\rho$ and
$\rho_{2}(\mathbf{r}_{1},\mathbf{r}_{2})=\rho^{2}g_{2}(\mathbf{r})$, $\rho$ is
the number density in the thermodynamic limit, $g_{2}(\mathbf{r})$ is the pair
correlation function, and $\mathbf{r}=\mathbf{r}_{2}-\mathbf{r}_{1}$. If the
system is also statistically isotropic, then $g_{2}(\mathbf{r})$ is the radial
function $g_{2}(r)$, where $r=|\mathbf{r}|$. The ensemble-averaged structure
factor $S(\mathbf{k})$ is defined as

$S(\mathbf{k})=1+\rho\tilde{h}(\mathbf{k}),$ (4)

where $h(\mathbf{r})=g_{2}(\mathbf{r})-1$ is the total correlation function,
and $\tilde{h}(\mathbf{k})$ is the Fourier transform of $h(\mathbf{r})$.

For a single periodic configuration containing number $N$ point particles at
positions ${\bf r}_{1},{\bf r}_{2},\ldots,{\bf r}_{N}$ within a fundamental
cell $F$ of a lattice $\Lambda$, the scattering intensity
$\mathcal{I}(\mathbf{k})$ is defined as

$\mathcal{I}(\mathbf{k})=\frac{\left|\sum_{i=1}^{N}e^{-i\mathbf{k}\cdot\mathbf{r}_{i}}\right|^{2}}{N}.$
(5)

For an ensemble of periodic configurations of $N$ particles within the
fundamental cell $F$, the ensemble average of the scattering intensity in the
infinite-volume limit is directly related to structure factor $S(\mathbf{k})$
by

$\lim_{N,V_{F}\rightarrow\infty}\langle\mathcal{I}(\mathbf{k})\rangle=(2\pi)^{d}\rho\delta(\mathbf{k})+S(\mathbf{k}),$
(6)

where $V_{F}$ is the volume of the fundamental cell and $\delta$ is the Dirac
delta function. [22] In simulations of many-body systems with finite $N$ under
periodic boundary conditions, Eq. (5) is used to compute $S(\mathbf{k})$
directly by averaging over configurations.

###  II.2 Hyperuniformity

A hyperuniform point configuration in $d$-dimensional Euclidean space
$\mathbb{R}^{d}$ possesses a structure factor $S(\mathbf{k})$ that goes to
zero as the wave number $\mathbf{k}$ vanishes, i.e.,
$\lim_{\left|\mathbf{k}\right|\rightarrow 0}S(\mathbf{k})=0$, which
corresponds to a local number variance $\sigma^{2}_{N}(R)$ in a spherical
window of radius $R$ that grows slower than $R^{d}$. [21, 22] For hyperuniform
systems whose structure factor is given by a power-law in the vicinity of the
origin, i.e.,

$S(\mathbf{k})\sim|\mathbf{k}|^{\alpha},\quad|\mathbf{k}|\rightarrow 0.$ (7)

The value of the exponent $\alpha>0$ determines three different “classes” of
hyperuniformity, [21, 84, 22] i.e.,

$\sigma^{2}(R)\sim\begin{cases}R^{d-1},\quad\alpha>1\text{ (class I)}\\\
R^{d-1}\ln R,\quad\alpha=1\text{ (class II)}\\\ R^{d-\alpha},\quad
0<\alpha<1\text{ (class III).}\\\ \end{cases}$ (8)

Stealthy hyperuniform systems, which include all perfect crystals and unusual
disordered states, [85, 26] are defined to be those that possess zero
intensity of the structure factor for a set of wavevectors around the origin,
[85] i.e.,

$S({\bf k})=0,\qquad\mbox{for}\;0\leq|{\bf k}|\leq K,$ (9)

where $K>0$. They are class I hyperuniform states that can be roughly regarded
to be those in which $\alpha\rightarrow\infty$.

To achieve equilibrium hyperuniform systems at positive $T$, which are
thermodynamically incompressible, one requires the following long-ranged pair
interaction in the large-$r$ limit: [22]

$v(r)\sim\begin{cases}r^{-(d-\alpha)},\quad d\neq\alpha\\\ -\ln(r),\quad
d=\alpha.\end{cases}$ (10)

A notable example of which is the one-component plasma (OCP), in which
identical point charges interacting via the Coulomb potential in
$\mathbb{R}^{d}$ are immersed in a rigid, uniform background of opposite
charge to ensure overall charge neutrality. [35] Such long-ranged potentials
can be precisely determined via our inverse methodology; see Sec. IV.

###  II.3 Higher-order statistics

Due to the complexity of computing and storing full pointwise information of
all $n$-particle correlation functions $g_{3},g_{4},g_{5},\dots$, this work
considers higher-order statistics that are readily computable, including the
conditional “hole” probability density function $G_{V}(r)$ and $g_{3}$ for
special triangles.

Given that a spherical region $\Omega_{V}(r)$ of radius $r$ is empty of
particles, the quantity $\rho s_{1}(r)G_{V}(r)dr$ is the probability of
finding particles in a spherical shell of volume $s_{1}(r)dr$, where
$s_{1}(r)$ is the surface area of a $d$-dimensional sphere of radius $r$. The
function $G_{V}(r)$ can be expressed in terms of integrals over all the
$n$-particle correlation functions $g_{2},g_{3},g_{4},\dots$. [82]
Importantly, $G_{V}(r)$ is directly related to the void-exclusion probability
function or “hole” probability function $E_{V}(r)$, which gives the
probability of finding a randomly located spherical region of radius $r$ empty
of particles, [80] via the relation

$G_{V}(r)=-\frac{E_{V}^{\prime}(r)}{\rho s_{1}(r)E_{V}(r)}.$ (11)

A many-body system possesses the bounded-hole property if $E_{V}(r)$ has
compact support, i.e., if the maximum hole radius $r_{c}$ is bounded. By
contrast, typical liquids possess holes of arbitrarily large size. [86] The
bounded-hole property characterizes disordered “stealthy” hyperuniform states
[27, 44, 87] and random sequential addition at saturation, [88, 89] but has
been hitherto unexplored for perfect glasses and absorbing-state models.

We also study the three-body statistics for small triangles, and especially
the distribution of bond angles $\theta$. Thus, we express $g_{3}$ in terms of
$\theta$, i.e.,

$g_{3}(r_{1},r_{2},\theta)=\frac{\rho_{3}\left(r_{1},r_{2},\sqrt{r_{1}^{2}+r_{2}^{2}-2ab\cos(\theta)}\right)}{\rho^{3}},$
(12)

where $\rho_{3}(r_{1},r_{2},r_{3})$ is probability density of finding three
particles that form a triangle with side lengths $r_{1},r_{2}$ and $r_{3}$.

##  III Nonequilibrium hyperuniform models

Here, we describe the three nonequilibrium target models we consider, i.e.,
perfect glasses, critical absorbing-state models and URL, as well as the
method of generating their corresponding configurations. They are hyperuniform
systems of recent interest and have important applications in photonics
engineering, [90] packing problems, [91] active matter [92] and geoscience.
[93]

###  III.1 Perfect glasses

Perfect glasses are exotic amorphous states of matter with positive bulk and
shear moduli that banish any crystalline or quasicrystalline phases and form
unique (nondegenerate) disordered states up to trivial symmetries. [38, 94]
These states can be regarded as prototypical glasses since they are out of
equilibrium, maximally disordered, hyperuniform, mechanically rigid with
infinite bulk and shear moduli, and remarkably prohibit the formation of
crystals and quasicrystals from the ground-state manifold due to
configuration-space trapping. The pair statistics of certain perfect glasses
can by realized by equilibrium canonical ensembles, [11] but the explicit
forms of the underlying potentials have heretofore remained unknown.

A perfect glass is created by cooling a many-body system from high to zero $T$
with a total potential energy $\Phi_{N}({\bf r}^{N})$: [38]

$\Phi(\mathbf{r}^{N})=\sum_{0<|\mathbf{k}|<K}\tilde{w}(\mathbf{k})\left[S(\mathbf{k})-S_{0}(\mathbf{k})\right]^{2},$
(13)

where $\mathbf{r}^{N}$ denotes the positions of the $N$ particles in
$\mathbb{R}^{d}$ that are subject to optimization, $K$ is the magnitude of the
largest constrained wavevector, $\tilde{w}(\mathbf{k})$ is a weight function,
and $S_{0}(\mathbf{k})$ is the desired small-$k$ behavior of the perfect-glass
structure factor. The number of independently constrained wavevectors divided
by the total number of degrees of freedom, $d(N-1)$, is a parameter $\chi$
that measures how constrained is the system. Here, we study a class II
hyperuniform perfect glass with the exponent $\alpha=1$. Following Ref. 38, we
choose the parameters $N=2500$,
$\tilde{w}(\mathbf{k})=\left(K/|\mathbf{k}|-1\right)^{3}$,
$S_{0}(\mathbf{k})=|\mathbf{k}|/K$, $\chi=5.1$ and $K=10/a$, where $a$ is a
length scale taken to be unity. These parameters are chosen so that the
position of the first peak in the perfect-glass $g_{2}(r)$ is close to $r=1$.
We generate perfect-glass configurations by finding local minima of
$\Phi(\mathbf{r}^{N})$ using the low-storage BFGS algorithm, [95] starting
from Poisson initial configurations; see Ref. 38 for further details.

###  III.2 Critical absorbing states

Random organization models reveal how chaotically-driven nonequilibrium many-
body systems can self-organize. [2, 96, 52] Such models have been applied to
study the dynamic phase transition of periodically sheared particles at the
onset of irreversibility. [97] Hexner and Levine showed that such critical
absorbing states are class III disordered hyperuniform with $\alpha=0.25$ in
three dimensions. [52]

Following Ref. 52, we generate critical-absorbing states as follows: Starting
from a Poisson initial configuration of $N$ spherical particles with unit
diameter $D$, a particle is deemed “active” if it overlaps with another
particle. In each iteration, all active particles are given a randomly
oriented displacement, whose magnitude is uniformly distributed on $[0,\sigma
D/2]$, where $\sigma$ is a displacement-size parameter. The process is
repeated until an absorbing state is reached where no particles overlap,
resulting in a packing of packing fraction $\phi=\rho\pi D^{3}/6$, or until
$10^{6}$ iterations are performed, implying that the system is an active
state. The critical state with packing fraction $\phi_{c}$ is determined such
that half of the initial configurations with $\rho=6\phi_{c}/(\pi D^{3})$
reach absorbing states, while the rest are active states. In this work, we
choose $N=10^{5}$ and $\sigma=0.1$, and estimate $\phi_{c}$ to be $0.205\pm
0.0005$.

###  III.3 Uniformly randomized lattices

Perturbed lattices serve as important models in cosmology, crystallography and
probability theory. [98, 99] Uniformly randomized lattices are simple
perturbed lattices that are class I hyperuniform with $\alpha=2$. [81] For 3D
URL models, each lattice point in a lattice, here taken to be the SC lattice
$\mathbb{Z}^{3}$, is displaced by a random vector that is uniformly
distributed on $[-b/2,b/2)^{3}$, where the scalar factor $b>0$ is the
perturbation strength. By definition, the lattice constant is set to be unity.
It has been shown that the structure factor for the URL point process contains
Bragg peaks that coincide with the unperturbed lattice as well as a diffuse
part such that $\lim_{k\rightarrow 0}S(k)\sim k^{2}$. [100, 101, 81]
Remarkably, Klatt et al. showed that the Bragg peaks in the structure factors
vanish completely, or become “cloaked”, when $b$ takes integer values. [81]

In this work, we study the 3D cloaked URL with $b=1$. whose effective
potential has a Coulombic asymptotic form, i.e., $v(r;\mathbf{a})\sim 1/r$, as
we have previously determined. [68] Due to the independence of the particles,
URLs are ideal models to study statistical structural descriptors, as they are
more readily computable compared to models with correlated particles. As we
will show in Sec. A Precise numerical algorithm to compute hole probabilities
for URLs, $E_{V}(r)$ [and thus $G_{V}(r)$] for this model can be numerically
evaluated to very high precision via our Monte-Carlo integration technique.

##  IV Inverse methodology

Here we sketch the methodology that we have recently introduced [68] to
determine (up to) pair interactions that yield canonical ensembles that very
accurately match both target $g_{2}(r)$ for all $r$ and target $S(k)$ for all
$k$. Importantly, for a statistically homogeneous system equilibrated under up
to two-body interactions, this technique is able to extract the unique target-
generating pair potential dictated by Henderson’s theorem [75]. The reader is
referred to Ref. 68 for a comprehensive description of the inverse
methodology.

The methodology optimizes a parameterized isotropic potential function
$v(r;{\bf a})$ that can be written as a sum of $n$ smooth pointwise basis
functions, i.e.,

$v(r;{\bf a})=\varepsilon\sum_{j=1}^{n}f_{j}(r/D;a_{j}),$ (14)

where $f_{j}(r/\sigma;a_{j})$ is the $j$th basis function, $a_{j}$ is a vector
of parameters (generally consisting of multiple parameters), ${\bf
a}=(a_{1},a_{2},\ldots,a_{n})$ is the “supervector” parameter whose components
are a collection of all components of all $a_{j}$’s, $\varepsilon$ sets the
energy scale and $D$ is a characteristic length scale, which is taken to be
unity. Examples of the basis functions include the hard core as well as
superexponential-, exponential-, Yukawa- and power-law-damped oscillatory
functions. The components of $a_{j}$ are of four types: dimensionless energy
scales $\varepsilon_{j}$, dimensionless distance scales $\sigma_{j}$,
dimensionless phases $\theta_{j}$, as well as dimensionless exponents $p_{j}$.

The initial form of $v(r;\mathbf{a})$ is informed by the small- and large-
distance behaviors of the target pair statistics $g_{2,T}(r)$ and $S_{T}(k)$,
as dictated by statistical-mechanical theory. [22] For hyperuniform targets,
the large-$r$ behavior of the potential is determined by Eq. (10). To obtain
an initial form of the small- and intermediate-$r$ behavior of
$v(r;\mathbf{a})$, we numerically fit the hypernetted chain (HNC)
approximation [83] for the target pair statistics using the aforementioned
forms of basis functions. The HNC approximation is given by

$\beta v_{\text{HNC}}(r)=h_{T}(r)-c_{T}(r)-\ln[g_{2,T}(r)],$ (15)

where $h_{T}(r)=g_{2,T}(r)-1$ and $c_{T}(r)$ is the target direct correlation
function, whose Fourier transform is given by the Ornstein-Zernike integral
equation [102]

${\tilde{c}}_{T}({\bf k})=\frac{{\tilde{h}}_{T}({\bf k})}{S_{T}({\bf k})}.$
(16)

Next, a nonlinear optimization procedure [95] is used to minimize an objective
function $\Psi({\bf a})$ based on the distance between target and trial pair
statistics in both direct and Fourier spaces:

$\begin{split}\Psi(\mathbf{a})&=\rho\int_{\mathbb{R}^{d}}w_{g_{2}}({\bf
r})\left(g_{2,T}({\bf r})-g_{2}({\bf r};\mathbf{a})\right)^{2}d{\bf r}\\\
&+\frac{1}{\rho(2\pi)^{d}}\int_{\mathbb{R}^{d}}w_{S}({\bf k})\left(S_{T}({\bf
k})-S(k;\mathbf{a})\right)^{2}d\mathbf{k},\end{split}$ (17)

where $g_{2,T}(\mathbf{r})$ and $S_{T}(\mathbf{k})$ are target pair
statistics, $w_{g_{2}}({\bf r})$ and $w_{S}({\bf k})$ are weight functions,
and $g_{2}({\bf r};\mathbf{a})$ and $S(\mathbf{k};\mathbf{a})$ correspond to
an equilibrated $N$-particle system under $v(r;\mathbf{a})$ at a dimensionless
temperature $k_{B}T/\varepsilon=1$, which can be obtained from Monte-Carlo
(MC) (used here with $N=500$) or molecular dynamics simulations under periodic
boundary conditions. The optimization procedure ends when $\Psi(\mathbf{a})$
is smaller than some small tolerance $\epsilon$. If convergence is not
achieved, then a different set of basis functions is chosen and the
optimization process is repeated. If convergence is achieved, we check whether
the effective potential robustly generates the target pair statistics for
systems larger than those used during the optimization process. We performed
MC simulations under the optimized potentials with $N=2500$ for 2D systems and
$N=9261$ for 3D systems, and utilized the dimensionless $L_{2}$ functions [(2)
and (3)] and the dimensionless total $L_{2}$-norm error (1) between target and
trial pair statistics to assess how close these quantities match.

For hyperuniform targets, $v(r;\mathbf{a})$ has the long-ranged asymptotic
form given by Eq. (10), which can be regarded as a generalized Coulombic
interaction of “like-charged” particles. [22] Thus, one requires a
neutralizing background one-body potential to maintain stability. [22, 33, 34,
35] 111Such background terms have been employed to study numerically the one-
component plasma [33, 34] and the Dyson log gas. [35] Importantly, to perform
the MC simulations, the total potential energy involving the long-ranged one-
and two-body potentials is efficiently evaluated using the Ewald summation
technique. [104] It is noteworthy that the optimized potentials via our
methodology generally yields pair statistics that accurately match their
corresponding targets with total $L_{2}$-norm errors that are an order of
magnitude smaller than that of previous methods, in cases where previous
methods are applicable. [68]

##  V Results for the nonequilibrium-equilibrium correspondence

###  V.1 Optimized equilibrium states

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: (a) A 2,500-particle configuration of the 2D perfect glass. (b) A
2,500-particle configuration of the optimized equilibrium state corresponding
to the perfect glass. (c) Short-ranged part of the optimized pair potential.
(d) Targeted and optimized pair correlation functions with $N=2500$. Here we
find that the $L_{2}$ norm function is $D_{g_{2}}=0.0030$. (e) Targeted and
optimized structure factors with $N=2500$. Here we find that the $L_{2}$ norm
function is $D_{S}=0.0030$. The $L_{2}$-norm error is $\mathcal{E}=0.077$. (f)
Log-log plot of the targeted and optimized structure factors, showing their
$k^{1}$ scaling behavior at small $k$.

Here, we presents the effective potentials and the pair statistics for the
equilibrium states corresponding to the target nonequilibrium states. Appendix
B gives the explicit functional forms of the optimized potentials.

Figure 1 shows the configurations of the target and optimized systems for the
perfect glass as well as the corresponding effective pair potential and pair
statistics. As shown in Fig. 1(a) and (b), the nonequilibrium and equilibrium
configurations are visually very similar. The optimized effective pair
potential has the expected asymptotic behavior $v(r)\sim 1/r$, as dictated by
Eq. (10). The short-ranged part of the potential [Fig. 1(c)], i.e.,
$v(r)-4.572/r$, contains two local minima at pair distances 1.20 and 2.10,
respectively. Figures 1(d) and (e) show $g_{2}(r)$ and $S(k)$, respectively,
for the target and optimized systems with $N=2500$. Figure 1(f) depicts the
structure factors on a log-log scale, showing that they are linear in $k$ at
small $k$. We find that the $L_{2}$ functions [Eqs. (2)–(3)] are
$D_{g_{2}}=0.0030$, $D_{S}=0.0030$ and the $L_{2}$-norm error (1) is
$\mathcal{E}=0.077$, showing that the pair statistics of the perfect glass are
in excellent agreement with those of the optimized equilibrium system in both
direct and Fourier spaces. Note that these errors are an order of magnitude
smaller than the errors obtained via IHNCI for typical equilibrium dense
liquids. [68] Remarkably, these results imply that one can reproduce
structures that arise from two-, three and four-body interactions via
effective pair interactions.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: (a) A 1,000-particle configuration of the 3D critical absorbing-
state model. (b) A 1,000-particle configuration of the optimized equilibrium
state corresponding to the critical-absorbing state. (c) Short-ranged part of
the optimized pair potential. (d) Targeted and optimized pair correlation
functions with $N=9261$. Here we find that the $L_{2}$ function is
$D_{g_{2}}=0.0020$. (e) Targeted and optimized structure factors with
$N=9261$. Here we find that the $L_{2}$ function is $D_{S}=0.0022$. The
$L_{2}$ norm error is $\mathcal{E}=0.066$. (f) Log-log plot of the targeted
and optimized structure factors, showing their $k^{0.25}$ scaling behavior at
small $k$.

Figure 2 presents the configurations, effective pair potential and pair
statistics of the target and optimized systems for the absorbing-state model.
The target and optimized configurations closely resemble each other [Fig. 2(a)
and (b)]. The effective pair potential [Fig. 2(c)] contains a hard core and a
sharp minimum at the sphere diameter $r=1$. It decays asymptotically as
$v(r)\sim r^{-2.75}$, which yields the correct the $k^{0.25}$ behavior of
$S(k)$ at small $k$. We see that both $g_{2}(r;\mathbf{a})$ and
$S(k;\mathbf{a})$ with $N=9261$ are in excellent agreement with those of the
target system (Fig. 2(d)–(f)), as manifested by the small values of the
$L_{2}$ functions [Eqs. (2)–(3)] and the $L_{2}$-norm error (1), given by
$D_{g_{2}}=0.0020$, $D_{S}=0.0022$ and $\mathcal{E}=0.066$.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: (a) A portion of a representative 3D 2,744-particle configuration of
a cloaked URL system. Only 512 particles are displayed. (b) A portion of a 3D
configuration of a 2,744-particle system that is equilibrated under the
optimized effective one- and two-body potential for the target 3D cloaked URL.
Only 512 particles are displayed. (c) Optimized pair potential minus its long-
ranged repulsive part $0.940/r$. (d) Targeted and optimized pair correlation
functions with $N=9261$. Here we find that the $L_{2}$ function is
$D_{g_{2}}=4.8\times 10^{-4}$. (e) Targeted and optimized structure factors
with $N=9261$. Here we find the $L_{2}$ function is $D_{g_{2}}=7.5\times
10^{-4}$. The $L_{2}$-norm error is $\mathcal{E}=0.035$. (f) Log-log plot of
the targeted and optimized structure factors, showing their $k^{2}$ scaling
behavior at small $k$. Subfigures (a)–(e) are reproduced from the ones first
presented in Ref. 68.

Figure 3 shows the configurations, effective pair potential and pair
statistics of the target and optimized systems for the 3D URL model, which
have been previously determined and presented in Ref. 68. Figures 3(a) and
3(b) that the target and the optimized systems configurations are visually
indistinguishable. Figure 3(c) shows the short-ranged part of the effective
potential, i.e. $v(r;\mathbf{a})-0.940/r$ against $r$. Figure 3(d)–(f) show
$g_{2}(r)$ and $S(k)$ for the target and optimized systems with $N=9261$. The
pair statistics of the target cloaked URL is in excellent agreement with those
of the optimized equilibrium system in both direct and Fourier space. The
$L_{2}$ functions are $D_{g_{2}}=4.8\times 10^{-4}$ and $D_{S}=7.5\times
10^{-4}$ and the $L_{2}$-norm error is $\mathcal{E}=0.035$, all of which are
remarkably small. For all models, the effective potentials accurately yield
the desired target pair statistics for system sizes much larger than $N=500$
used in the optimization procedure. Thus, the nonequilibrium-equilibrium
correspondences found in this study are robust to further increases in the
system size.

###  V.2 Higher-order statistics

Since the targeted and optimized systems have essentially the same pair
statistics, we expect that due to the aforementioned structural degeneracy,
the dynamics leading to the nonequilibrium states is reflected in the
differences in their higher-order statistics compared to the corresponding
equilibrium states. Here, we present the higher-order statistics of the
nonequilibrium-equilibrium pairs, including $G_{V}(r)$ for all models and
$g_{3}$ for the perfect glass and the critical absorbing state at specific
small triangles. Note that $g_{3}$ and $g_{4}$ for the URL have been discussed
in Ref. 81.

(a)

(b)

(c)

Figure 4: Conditional “hole” probability density functions for three target
and optimized systems. Note the different scales in the axes. (a) 2D perfect
glass. The inset shows a magnified portion of the plots. (b) 3D critical
absorbing-state model. (c) 3D cloaked URL.

In Appendix A, we describe a new and precise numerical algorithm to compute
the hole probabilities for the nonequilibrium cloaked URLs in any space
dimension $d$, which exploits the independence of particles in the URLs to
achieve much higher accuracy than the standard method. [81] Figure 4(a) and
(b) show $G_{V}(r)$ of the nonequilibrium-equilibrium pairs for the 2D perfect
glass, the 3D critical absorbing state and the 3D cloaked URL, respectively.
For all three models, while the target and optimized $G_{V}(r)$ agree closely
in the range where $r$ is below half of the characteristic length scale, they
become increasingly different at larger $r$, revealing the distinctly
different higher-order statistics between the nonequilibrium and equilibrium
systems. Importantly, both target and optimized $G_{V}(r)$ for all systems
studied exhibit steep rises at sufficiently large $r$, indicating that they
possess the bounded-hole property. [27, 44]

The forms of the effective potentials for the perfect glass and critical-
absorbing state allow us to conclude that these models have the bounded-hole
property, since the potentials contain hard cores (an effective hard core in
the perfect-glass case), small but positive repulsions at intermediate $r$,
and long-ranged repulsions at large $r$ that ensure hyperuniformity. Thus,
particles of the resulting packings are well-separated from one another at all
length scales. It is noteworthy that the bounded hole property applies for
perfect glasses with any $\alpha>0$, because of the nature of the long-ranged
two-, three- and four-body interactions.[38] Indeed, we have also computed the
hole probabilities for a different perfect glass with $\alpha=2$ and confirmed
its bounded-hole property.

Using the technique described in Appendix A, we accurately determine that
$E_{V}(r)$ for the target URL in the limit $r\rightarrow r_{c}^{-}$ exhibits a
form first derived in Ref. 87, namely, $E_{V}(r)\sim(r_{c}-r)^{\gamma}$, where
$\gamma=32$. Thus, Eq. (11) gives

$G_{V}(r)\sim\frac{\gamma}{\rho s_{1}(r)(r_{c}-r)}\quad r\rightarrow
r_{c}^{-}$ (18)

i.e., $G_{V}(r)$ has a pole of order one, as derived in Ref. 87. Compared to
the SC lattice, for which $\gamma=d=3$, [87] the very large $\gamma$ value for
the URL indicates that holes close to the critical-hole size are rarer in
disordered systems, as noted in Ref. 87. Also note that unlike the other two
models, $G_{V}(r)$ for the target URL is higher than the optimized equilibrium
$G_{V}(r)$, suggesting that the target URL is locally more homogeneous than
its equilibrium counterpart due to the underlying lattice.

We propose the following scalar nonequilibrium index, $\Gamma$, that reflects
differences in the higher-order correlation functions $g_{3},g_{4},\ldots$ for
a nonequilibrium-equilibrium pair with matching pair statistics:

$\Gamma=\rho\int_{|\mathbf{r}|<r_{\text{max}}}\left[G_{V}(|\mathbf{r}|)-G_{V,T}(|\mathbf{r}|)\right]^{2}d\mathbf{r},$
(19)

where $G_{V}$ and $G_{V,T}$ are conditional “hole” probability density
distributions for the optimized and target systems, respectively, and
$r_{\text{max}}$ is a cutoff radius corresponding to the largest hole radius
detected in an ensemble of finite-size configurations.222If the target and
optimized systems have different largest detected hole radii, $r_{\text{max}}$
is chosen to be the smaller of the two. The aforementioned difference is
within $0.1\rho^{1/d}$ for all three models. By definition, $G_{V}(r)$ in (19)
must be computed at $kT/\varepsilon=1$, because it is only at this
dimensionless temperature that the nonequilibrium and equilibrium pair
statistics match. Note that $\Gamma$ is a purely static nonequilibrium index,
i.e., its computation does not require dynamic information. Importantly, for a
statistically homogeneous equilibrium target system under up to two-body
interactions, one has $\Gamma=0$, because the optimized pair potential
obtained via the inverse methodology (Sec. IV) must agree with the unique
target-generating potential, yielding $G_{V}(r)=G_{V,T}(r)$ for all $r$.

Except for the target URL, whose $G_{V}(r)$ can be evaluated up to arbitrary
precision, we determine $r_{\text{max}}$ and compute $G_{V}(r)$ using $500$
configurations with $N=2500$ for 2D systems and $N=9261$ for 3D systems. Table
1 shows the value of $\Gamma$ for the three models shown in Fig. 4. The URL
has a much larger $\Gamma$ value compared to those of the the other two
systems. This high degree of “nonequilibriumness” is expected for the cloaked
URL, since $g_{4}$ for the target system exhibits long-range order
characteristic of the underlying lattice. [81] An equilibrium fluid cannot
possess any long-range order.

Table 1:  Values of the metric $\Gamma$ [Eq. (19)] for the systems shown in Fig. 4 System | $\Gamma$
---|---
2D perfect glass | 4.9
3D critical-absorbing state | 9.2
3D cloaked URL | 54

We also study the three-body statistics for the perfect glass and the critical
absorbing state at small triangles. In the case of the perfect glass, we
compute the integral

$f(\theta)=\int_{|\mathbf{r}_{1}|=0}^{1.15}\int_{|\mathbf{r}_{2}|=0}^{1.15}g_{3}(r_{1},r_{2},\theta)d\mathbf{r}_{1}d\mathbf{r}_{2},$
(20)

where $g_{3}(r_{1},r_{2},\theta)$ is given by (12). We evaluate $f(\theta)$
instead of $g_{3}$ at individual triangles, because much higher accuracy can
be achieved for the former. Figure 5(a) depicts $f(\theta)$ for the target and
optimized perfect glass systems. We find that compared to the equilibrium
system, the target perfect glass contains 50% more nearly linear small
triangles with $\theta>170^{\circ}$. This difference is likely due to the
actual three and four-body potential in the target perfect-glass system. [38]

(a)

(b)

Figure 5: (a) Plot of
$f(\theta)=\int_{|\mathbf{r}_{1}|=0}^{1.15}\int_{|\mathbf{r}|_{2}=0}^{1.15}g_{3}(\mathbf{r}_{1},\mathbf{r}_{2},\theta)d\mathbf{r}_{1}d\mathbf{r}_{2}$
against $\theta$ for the target and optimized perfect glass systems. (b) Plot
of $g_{3}(1.05,1.05,\theta)$ against $\theta$ for the target and optimized
critical absorbing state systems.

In the case of the critical-absorbing-state model, we compute $g_{3}$ at
individual triangles, as it can be evaluated with high accuracy. Figure 5(b)
shows $g_{3}(r_{1},r_{2},\theta)$ for small isosceles triangles with
$r_{1}=r_{2}=1.05$ for the target and optimized critical absorbing state, from
which it is clear that the target system significantly suppresses the
formation of small triangles that are nearly equilateral. Indeed, their
difference of $g_{3}$ for small equilateral triangles is 90%. This is due to
the specific dynamics during the random organization process, where particles
are displaced to avoid collisions. [2] Dense clusters that contain many small
equilateral triangles are less likely to be found in the nonequilibrium state,
since collisions are more frequent within such clusters. On the other hand,
the equilibrium state forms more dense clusters due to the sharp well of
$v(r)$ at the hard-sphere diameter.

###  V.3 Inherent structures

To study the effect of quenching on the hyperuniform equilibrium systems, we
ascertain the inherent structures of the perfect-glass and critical-absorbing-
state potentials by finding deep local energy minima via the low-storage BFGS
algorithm, [95] starting from equilibrium initial configurations with $N=1000$
for $d=2$ and $N=1728$ for $d=3$, and averaging over 200 configurations. We
find that the energy per particle $\overline{E}_{\text{IS}}$ of these
structures are narrowly distributed above the ground-state energies; see
Appendix C for details. We compute the inherent-structure pair statistics by
averaging over configurations whose energies are in within two standard
deviations around the mean value of $\overline{E}_{\text{IS}}$.

(a)

(b)

(c)

(d)

Figure 6: (a) A 1,000-particle configuration of an inherent structure for the
effective potential for the perfect glass. (b) Pair correlation function of
the inherent structure. (c) Structure factor of the inherent structure. (d)
Log-log plot of the structure factor for the inherent structure, showing the
$k^{2}$ scaling at small $k$.

Figure 6(a) shows a disordered inherent structure for the perfect-glass
effective potential, which contains chains of particles separated by
“channels” with a well-defined width. The first three peaks of the associated
$g_{2}(r)$ [Fig. 6(b)] correspond to the nearest-neighbor distance, the
channel width, and the second-nearest neighbor distance along the chains,
respectively. The first and second peaks occur at $r=1.27$ and $r=2.20$,
respectively, which are close to the locations of the first and second minima
of the short-ranged part of $v(r)$. Figure 6(c) depicts $S(k)$ for the
inherent structures.

To extract the small-$k$ behaviors of the structure factors, we use the
concept of diffusion spreadability introduced in Ref. 64, defined as follows.
Consider the time-dependent problem of mass transfer of a solute between two
phases and assume that the solute is initially distributed in one phase (phase
2) and absent from the other (phase 1). The spreadability $\mathcal{S}(t)$ is
the fraction of total solute present in phase 1 as a function of time. For
sphere packings, which include the perfect glass and the critical absorbing
state, $\mathcal{S}(t)$ can be computed from simulated $S(k)$ via the Fourier-
space expression derived in Ref. 64. Recently, Wang and Torquato [93]
introduced an algorithm that efficiently and accurately extracts the exponent
$\alpha$ from numerical data of $\mathcal{S}(t)$. The algorithm extracts
$\alpha$ as well as a “set-in” time $t_{S}$ of the asymptotic behavior of
$\mathcal{S}(t)$ by solving the following relations via a predictor-corrector
procedure

$\begin{split}|\mathcal{S}(t)-\mathcal{S}_{l}(t;\alpha)|=\epsilon,\quad
t=t_{S},\\\ |\mathcal{S}(t)-\mathcal{S}_{l}(t;\alpha)|<\epsilon,\quad
t>t_{S}\end{split}$ (21)

where $\epsilon>0$ is a convergence criterion and $\mathcal{S}_{l}(t;\alpha)$
is the large-$t$ spreadability approximant. For nonstealthy media, one has
$t^{-(d+\alpha)/2}$. [64] It has been shown that the aforementioned algorithm
is more robust to simulation noise than a direct numerical fit of $S(k)$ at
small $k$, and accurately extracts $\alpha$ with errors less than 1% for a
wide range of models. [93] We find that the inherent structures are
hyperuniform with $\alpha=2.0$. Figure 6(d) depicts $S(k)$ for the inherent
structures on a log-log scale, showing clearly the $k^{2}$ behavior at small
$k$. Thus, the inherent structures of a stronger form of hyperuniformity
(class I) than that of the target structure. Figure 7 shows a disordered
inherent structure for the critical-absorbing-state potential, as well as the
inherent-structure pair statistics, indicating that it is again hyperuniform
with $\alpha=2$. By contrast, $\alpha$ for the URL potential [68] dramatically
increases from 2 at unit temperature to infinity upon quenching, as the
inherent structure is the SC lattice. Thus, in all three cases, the inherent
structures are of higher forms of hyperuniformity than the target models.

To verify that the increase of $\alpha$ upon quenching is robust to the system
size, we also computed the inherent structures for the critical-absorbing
state with $N=512,729$ and $1000$, in addition to $N=1728$. We observed that
while the structure factors for these system sizes are slightly different due
to finite-size effect, they all exhibit $k^{2}$ scaling behaviors at small
$k$, indicating that the inherent structures possess a large hyperuniformity
scaling regime [43] that grows with $N$.

(a)

(b)

(c)

(d)

Figure 7: (a) A 1,000-particle configuration of an inherent structure for the
effective potential for the 3D critical-absorbing state. (b) Pair correlation
function of the inherent structure with $N=1728$. (c) Structure factor of the
inherent structure with $N=1728$. (d) Log-log plot of the structure factor for
the inherent structure, showing the $k^{2}$ scaling at small $k$.

##  VI Conclusions and discussion

We have determined equilibrium systems with effective one- and two-body
interactions that realize the pair statistics of three hyperuniform
nonequilibrium models of recent interest: a 2D perfect glass, a 3D critical-
absorbing state and a 3D cloaked URL, which lends further support to the
Zhang-Torquato conjecture. In the case of the perfect glass, it is noteworthy
that we are able to mimic systems with up to four-body interactions through
effective one-body and isotropic two-body interactions. This finding offers
promise to determining effective one- and two-body interactions that mimic
systems with directional interactions, such as amorphous graphene [106] and
amorphous silicon [107], governed by up to three- and higher-body
interactions. [108, 109] We have shown that all three models considered in
this study have the bounded-hole property.

The volume integral of the squared differences of $G_{V}(r)$ between such
nonequilibrium and equilibrium systems with identical pair statistics enable
us to define a nonequilibrium index $\Gamma$ [Eq. (19)], which can be applied
to measure the “nonequilibriumness” in many systems of practical interest,
including supercooled liquids, [37] defects, [110] self-propelling particles
[5, 111] and biological tissues. [112] We remark that a different static
nonequilibrium index $X$ has been proposed previously, which is based on the
deviation of $S(0)$ from $\rho\kappa_{T}k_{B}T$, where $\kappa_{T}$ is the
isothermal compressibility. [37] In comparison, $\Gamma$ extracts the
nonequilibriumness from structural information alone, and can be applied in
situations where $\kappa_{T}$ (or the pressure) is not readily available. We
also introduced a precise numerical algorithm to compute $E_{V}(r)$ for the
cloaked URL, and applied it to show that holes near the maximum hole size in
the URL are much rarer than those in the underlying SC lattice.

Remarkably, our effective potentials yield hyperuniform deep local energy
minima that are of stronger forms of hyperuniformity (measured by the class or
larger $\alpha$ exponents) relative to the equilibrium states at unit
temperature. This behavior is expected for pair potentials characterized by
steep short-range repulsions, since the local minima can be regarded to be
states with effectively low temperatures relative to the initial higher-
temperature fluid [113] and hence the structure factor $S(k)$ of the quench at
low wavenumbers reflects stronger hyperuniformity. By contrast, the same
reasoning leads to the extended proposition that a quench of an initially
nonhyperuniform fluid state with such interactions to deep local energy minima
have a $S(0)$ that drops but generally not to zero, i.e., it strictly remains
nonhyperuniform, as shown for three well-known models in Appendix C.

A promising avenue for future research is the determination of thermodynamic
and dynamic properties associated with the effective potentials, including
ground states, entropies, free energies, phase diagrams, [114, 115] and their
glass formations. For instance, it has been shown that perfect glasses at
$T=0$ are remarkably zero-entropy states, [94] in contradistinction to normal
glasses that are metastable with respect to crystals, and to disordered
stealthy ground states that possess large extensive entropies. [25] On the
other hand, the equivalent equilibrium system for a perfect glass has positive
entropy, and an effective pair potential [Fig. 1(c)] yields a crystalline
ground state, as shown in Appendix C. Furthermore, we note that the equivalent
nonequilibrium and equilibrium systems possess the same two-body contribution
to the excess entropy per particle due to their identical pair statistics.
[83] However, this excess two-body entropy can be significantly different from
the real entropy even for simple equilibrium liquids. [116, 117] Thus, it is a
fascinating problem to study their higher-order correlation contributions to
the entropy, which is expected to reveal crucial dynamical information.

Finally, we stress that the effective potentials yielding hyperuniform states
at positive $T$ could enable one to produce tunable hyperuniform materials.
Due to the long-ranged nature of these potentials (10), one can use $\rho$ and
$T$ as tuning parameters to generate equilibrium hyperuniform structures whose
exponents $\alpha$ are dictated by (10), or to generate stronger hyperuniform
forms via their inherent structures. While it is challenging to achieve long-
ranged interactions in the laboratory, one could experimentally reproduce the
effective potentials over some finite but large range of $r$ to fabricate
effectively hyperuniform states, i.e., states with very small but nonvanishing
$S(0)$. Subsequently, the deviation of such systems from perfect
hyperuniformity can be characterized via the various quantitative measures
described in Ref. 43.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

The authors gratefully acknowledge the support of the National Science
Foundation under Award No. CBET-2133179. S.T. thanks the Institute for
Advanced Study for their hospitality during his sabbatical leave there.

## Appendix

## A Precise numerical algorithm to compute hole probabilities for URLs

Here, we describe the standard method of computing the void-exclusion
probability or “hole” probability function $E_{V}(r)$ for a many-body system.
For systems with the bounded hole property, sampling holes with radii close to
$r_{c}$ is computationally challenging. To precisely determine the behaviors
of $E_{V}(r)$ and $G_{V}(r)$ as $r\rightarrow r_{c}^{-}$, we introduce an
improved algorithm to compute $E_{V}(r)$ for URLs in all dimensions with any
perturbation strength $b$ that is much more accurate than the standard method.
Although this improved method only applies to URLs, we will see in Sec. V.2
that the results for the hole probabilities of this model confirms the general
trend that holes close to the critical-hole size are rarer in disordered
systems than in ordered systems. [87] We first describe the algorithm specific
to the 3D cloaked URL (i.e., for $d=3$ and $b=1$), and then generalize it to
other $d$ and $b$.

In the standard method, $N_{t}$ test particles are randomly placed in a
$N$-particle configuration
$\mathbf{r}^{N}=\\{\mathbf{r}_{1},\dots\mathbf{r}_{N}\\}$ under periodic
boundary conditions. For each test particle $\mathbf{x}$, one computes the
distance $d_{\mathbf{x}}$ from $\mathbf{x}$ to its nearest real particle
$\mathbf{r}_{i}$. The function $E_{V}(r)$ is estimated to be

$E_{V}(r)=\frac{|\\{\mathbf{x}:d_{\mathbf{x}}>r\\}|}{N_{t}}.$ (22)

The result is then averaged over $n$ configurations. Using $nN=5\times
10^{7}$, $N_{t}/N=10^{5}$, this method yields $E_{V}(r)$ for the target URL in
the range $r\rho^{1/3}<1.3$ with errors $<5\%$. However, computing $E_{V}(r)$
at larger $r$ requires sampling rare events with probabilities smaller than
$10^{-15}$, which is computationally challenging. [87] Thus, we exploit the
independence of the perturbed lattice points to devise a more accurate method
to compute $E_{V}(r)$.

As described in Sec. III.3, the cloaked URL is obtained by perturbing each
lattice point in the SC lattice $\mathbb{Z}^{3}$ by a random vector uniformly
distributed on $[-1/2,1/2)^{3}$. For $i,j,k\in\mathbb{Z}$, let
$C_{ijk}=[i-1/2,i+1/2)\times[j-1/2,j+1/2)\times[k-1/2,k+1/2)$, i.e., the cubic
region accessible to the perturbed lattice point originally at $(i,j,k)$. Let
$\mathbf{x}$ be the position of a test particle and
$\mathcal{B}(\mathbf{x};r)$ be the spherical region of radius $r$ centered at
$\mathbf{x}$.

Since the lattice points are perturbed independently, the probability
$p(\mathbf{x},r)$ that no perturbed lattice point is found in
$\mathcal{B}(\mathbf{x};r)$ is given by the product of the probabilities that
each perturbed point is not in $\mathcal{B}(\mathbf{x};r)$, i.e.,

$p(\mathbf{x},r)=\prod_{ijk}|C_{ijk}\setminus\mathcal{B}(\mathbf{x};r)|,$ (23)

where $|C_{ijk}\setminus\mathcal{B}(\mathbf{x};r)|$ is the volume of the set
difference between $C_{ijk}$ and $\mathcal{B}(\mathbf{x};r)$. Because the
probability distribution of the position of the test particle is uniform, we
have

$E_{V}(r)=\langle
p(\mathbf{x},r)\rangle_{\mathbf{x}\in\mathbb{R}^{3}}=\int_{C_{000}}p(\mathbf{x},r)d\mathbf{x},$
(24)

where the second equality follows from the fact that all $C_{ijk}$ are
equivalent in a URL. In what follows, we let $\mathbf{x}\in C_{000}$.

If $r\geq\sqrt{3}$, the region $\mathcal{B}(\mathbf{x};r)$ covers $C_{000}$,
i.e., $|C_{000}\setminus\mathcal{B}(\mathbf{x};r)|=0$ for any $\mathbf{x}\in
C_{000}$; It follows from Eqs. (23) and (24) that $E_{V}(r)=0$ for
$r\geq\sqrt{3}$. On the other hand, if $0\leq r<\sqrt{3}$, then
$p(\mathbf{x},r)\neq 0$ for $\mathbf{x}=(-1/2,-1/2,-1/2)^{T}$, and thus
$E_{V}(r)$ is nonzero. Therefore, the maximum hole radius in the 3D cloaked
URL is given by $r_{c}=\sqrt{3}$.

For $r<r_{c}$, Eq. (23) reduces to

$p(\mathbf{x},r)=\prod_{|i|\leq 2,|j|\leq 2,|k|\leq
2}|C_{ijk}\setminus\mathcal{B}(\mathbf{x};r)|.$ (25)

The factors in Eq. (23) that are not included Eq. (25) are equal to unity
because $\mathcal{B}(\mathbf{x};r)$ are disjoint from their corresponding
$C_{ijk}$. The 125 factors in Eq. (25) can be evaluated via Monte Carlo
integration technique. [118] We randomly place $N_{s}=10^{6}$ sample points in
each $C_{ijk}$ region and estimate
$|C_{ijk}\setminus\mathcal{B}(\mathbf{x};r)|$ to be the fraction of the number
of sample points not found in $\mathcal{B}(\mathbf{x};r)$. We compute
$p(\mathbf{x},r)$ at all $\mathbf{x}$ vectors on a $M\times M\times M$ regular
mesh on $C_{000}$, where $M$ is a positive integer. The mesh is chosen to be
finer (i.e., $M$ larger) with increasing $r$ to ensure that there are at least
1,600 $\mathbf{x}$ vectors on the support of $p(\mathbf{x},r)$. Subsequently,
we numerically evaluate $E_{V}(r)$ [Eq. (24)] using a Gaussian quadrature of
order 3. [118] This method yields errors of $E_{V}(r)$ on the order of 2% for
$1.6\leq r<r_{c}$, whereas the errors using the standard method is on the
order of 100% in this range. Note that the errors can be further reduced by
increasing $N_{s}$.

To generalize the algorithm to other $d$ and $b$, one simply replaces
$C_{ijk}$ by
$C_{i_{1},...,i_{d}}=[i_{1}-b/2,i_{1}+b/2)\times\dots\times[i_{d}-b/2,i_{d}+b/2)$,
which is the hypercubic region accessible to each perturbed lattice point
originally at $(i_{1},\dots,i_{d})\in\mathbb{Z}^{d}$ in the $d$-dimensional
hypercubic lattice. Equations (24) and (23) will then be replaced by

$E_{V}(r)=\langle
p(\mathbf{x},r)\rangle_{\mathbf{x}\in\mathbb{R}^{d}}=\int_{C_{0,\dots,0}}p(\mathbf{x},r)d\mathbf{x}$
(26)

and

$p(\mathbf{x},r)=\prod_{|i_{m}|\leq\sqrt{d}+1;m=1,\dots,d}|C_{i_{1},\dots,i_{d}}\setminus\mathcal{B}(\mathbf{x};r)|$
(27)

respectively.

## B Optimized effective potentials for the target models

For the 2D perfect glass, the $\alpha=1$ behavior of $S(k)$ at small $k$
implies that $v(r;\mathbf{a})\sim 1/r$ as $r\rightarrow\infty$. The fact that
$g_{2}(r)=0$ for $r\leq 0.88$ implies that there is an effective hard core in
this range. The intermediate-$r$ behavior of $v(r;\mathbf{a})$ was determined
by fitting $v_{\text{HNC}}(r)$ in the range $0.88<r\leq 4$. No re-selection of
basis function was needed. The optimized $v(r)$ is given by

$v(r;\mathbf{a})=\begin{cases}\infty\qquad 0\leq r\leq 0.88\\\
\frac{\varepsilon_{1}}{r}+\varepsilon_{2}\exp\left(-\frac{r}{\sigma_{2}^{(1)}}\right)\cos\left(\frac{r}{\sigma_{2}^{(2)}}+\theta_{2}\right)\\\
+\varepsilon_{3}\exp\left(-(\frac{r}{\sigma_{3}})^{8}\right)\qquad
r>0.88.\end{cases}$ (28)

The optimized parameters are listed in Table 2. The $L_{2}$-norm error is
$\mathcal{E}=0.077$.

For the 3D critical-absorbing-state model, the fact that $\alpha=0.25$ implies
that $v(r;\mathbf{a})\sim 1/r^{2.75}$ at large $r$. Since the target system is
a packing of spheres with unit diameter, $v(r;\mathbf{a})$ has a hard core for
$r\leq 1$. The intermediate-$r$ behavior of $v(r;\mathbf{a})$ was obtained by
fitting $v_{\text{HNC}}(r)$ in the range $1<r\leq 4$. No re-selection of basis
function was needed. The optimized $v(r)$ is given by

$v(r;\mathbf{a})=\begin{cases}\infty\qquad 0\leq r\leq 1\\\
\frac{\varepsilon_{1}}{r^{2.75}}+\varepsilon_{2}\exp\left(-(\frac{r}{\sigma_{2}})^{5}\right)+\varepsilon_{3}\exp\left(-(\frac{r}{\sigma_{3}})^{8}\right)\\\
+\varepsilon_{4}\exp\left(-(\frac{r}{\sigma_{4}})^{10}\right)\qquad
r>1.\end{cases}$ (29)

The optimized parameters are listed in Table 3. The $L_{2}$-norm error is
$\mathcal{E}=0.066$.

The optimized potential for the 3D cloaked URL is given in Ref. 68.

Table 2: Optimized parameters of the effective pair potential for the 2D perfect-glass model. $\varepsilon_{1}$ | 4.572 | $\theta_{2}$ | 1.032
---|---|---|---
$\varepsilon_{2}$ | 12.26 | $\varepsilon_{3}$ | 5.803
$\sigma_{2}^{(1)}$ | 0.3921 | $\sigma_{3}$ | 0.9766
$\sigma_{2}^{(2)}$ | 0.1505 |  |
Table 3: Optimized parameters of the effective pair potential for the 3D critical-absorbing-state model. $\varepsilon_{1}$ | 0.3400 | $\sigma_{3}$ | 1.582
---|---|---|---
$\varepsilon_{2}$ | -0.8120 | $\varepsilon_{4}$ | -49.13
$\sigma_{2}$ | 0.9949 | $\sigma_{4}$ | 0.8827
$\varepsilon_{3}$ | 0.1068 |  |

## C Inherent structures and ground states

Table 4: Inherent-structure and ground-state energies levels for hyperuniform and nonhyperuniform models with $N=1000$ and the small-wavevector behaviors of the inherent structures. System | Ground-state structure | $\overline{E}_{0}/\varepsilon$ | $\overline{E}_{\text{IS}}/\varepsilon$ | $\overline{E}_{\text{IS}}/\overline{E}_{0}$ | $S_{\text{IS}}(0)$ | $\alpha_{\text{IS}}$
---|---|---|---|---|---|---
2D perfect glass | Fig. 8(a) | -5.973 ††footnotemark:  |  $-5.594\pm 0.0033$ ††footnotemark:  | 94% | 0 | 2
3D cloaked URL | sc | -1.081 ††footnotemark:  | -1.081 ††footnotemark:  | 100.0% | 0 | $\infty$
3D critical-absorbing state | Fig. 8(b) | -7.104 ††footnotemark:  |  $-5.29\pm 0.035$ ††footnotemark:  | 74% | 0 | 2
3D GC | bcc [119] | 2.284297 | $2.284521\pm 4.1\times 10^{-6}$ | 100.010% | $8.9\times 10^{-8}$ | …
3D PL | fcc [120] | 6.066 | $6.971\pm 0.026$ | 115% | 0.0018 | …
3D LJ | hcp [120] | -8.593 | $-7.836\pm 0.018$ | 91% | 0.0047 | …

Here, we present details on the energy levels of the inherent structures for
the hyperuniform models considered in this study, including the 2D perfect
glass, the 3D critical-absorbing state and the 3D cloaked URL. We also
consider inherent structures for several well-studied 3D models that yield
nonhyperuniform states at positive $T$ in three dimensions, including Lennard-
Jones (LJ), inverse power-law (PL) and Gaussian-core (GC) [121] models. For
all models, we ascertain the the inherent structures using the low-storage
BFGS algorithm, [95] staring from equilibrium initial states away from phase
transitions. The inherent-structure energies per particle
$\overline{E}_{\text{IS}}$ are compared with the corresponding ground-state
energies $\overline{E}_{0}$.

22footnotetext: Energy value includes the “neutralizing” one-body potential.

To find the ground states for the hyperuniform models, we first applied the
simulated annealing algorithm [122] on configurations with $N\geq 100$
particles, starting from the equilibrium liquid states that match the target
pair statistics. The ground state for the URL effective potential is found to
be the SC lattice. However, for the perfect glass and the critical-absorbing
state, we found that even with very slow cooling rates, the annealing
procedures resulted in disordered metastable states and were unable to find
the global minima.333With cooling rate $T_{i+1}=0.999T_{i}$, where $i$ is the
iteration index, the annealing procedures still result in metastable states
for the hyperuniform models. Therefore, we searched for ground-state
candidates by optimizing crystalline structures via simulated annealing over
$n$-particle bases, where $n=1,2,\dots,10$. The parameters subject to
optimization are unit cell vectors and coordinates of the $n-1$ particles in
the interior of a unit cell. We observe that $v(r)$ for the perfect glass
gives the same optimized crystal structure for all even $n$, whose energy is
below the optimized energies for odd $n$, and that $v(r)$ for the critical-
absorbing state yields identical optimized structures for all $n$. Thus, it is
highly likely that the ground state for the perfect-glass potential is a two-
particle basis, whereas that for the critical-absorbing-state potential is a
Bravais lattice. Figure 8 shows the ground-state candidates for both models.
The ground-state candidate for the perfect glass contains zig-zag chains of
particles, whereas the one for the critical-absorbing state contains planar
sheets of triangle lattice.

(a)

(b)

Figure 8: Structures of the ground-state candidates associated with $v(r)$ for
(a) the perfect glass and (b) the critical-absorbing state.

The potentials for the nonhyperuniform models are given by

$v_{\text{LJ}}(r)=4\varepsilon\left[(r/\sigma)^{-12}-(r/\sigma)^{-6}\right],$
(30) $v_{\text{PL}}(r)=4\varepsilon\left[(r/\sigma)^{-12}\right],$ (31)
$v_{\text{GC}}(r)=\varepsilon\exp\left[-(r/\sigma)^{2}\right],$ (32)

where $\varepsilon$ and $\sigma$ are energy and distance scales, respectively.
The LJ and PL models are strongly repulsive at small $r$. While
$v_{\text{GC}}(r)$ is a soft interaction, the model possesses an effective
hard core at low temperatures. [124, 125] We compute the inherent structures
for Eqs. (30)–(32) starting from equilibrium configurations at
$\rho\sigma^{3}=1,kT/\varepsilon=1$. These $\rho,T$ values are chosen so that
the initial states are dense liquids away from phase transitions for all three
potentials. [120, 119, 125]

Table 4 lists values of $\overline{E}_{0}$ and $\overline{E}_{\text{IS}}$ for
the hyperuniform and nonhyperuniform models with $N=1000$, as well as the
small-wavevector behaviors extracted using the diffusion spreadability. [93]
For all models except the URL, $\overline{E}_{\text{IS}}$ values are narrowly
distributed above $\overline{E}_{0}$. The distributions are Gaussian and are
insensitive to the system size. For the perfect glass,
$\overline{E}_{\text{IS}}/\overline{E}_{0}=94\%$, and the standard deviation
$\sigma_{\text{IS}}$ of $\overline{E}_{\text{IS}}$ is $5\times
10^{-4}\overline{E}_{0}$. For the critical-absorbing state,
$\overline{E}_{\text{IS}}/\overline{E}_{0}=74\%$ and
$\sigma_{\text{IS}}=0.005\overline{E}_{0}$. The inherent structure for the URL
is the SC lattice, identical to the ground state. The fact that the structure
is crystalline (periodic) means that it is stealthy hyperuniform and hence
$\alpha\to\infty$. [22]

Importantly, the three systems that are hyperuniform at positive $T$ yield
hyperuniform inherent structures with increased values of $\alpha$ compare to
those of the equilibrium states at positive $T$, which is consistent with the
proposition that the quench reflects stronger forms of hyperuniformity, as
reported in Sec. V.3 of the main article. In particular, $\alpha$ for the URL
potential dramatically increases from 2 at unit temperature to infinity upon
quenching. On the other hand, the pair potentials for nonhyperuniform systems
yield nonhyperuniform inherent structures. Note that while the inherent
structure for the GC model is nearly hyperuniform, [22] it is not perfectly
hyperuniform, which again is consistent with the proposition stated in the
main article.

## References

  * Brenner and Shokef [2007] N. Brenner and Y. Shokef.  Nonequilibrium statistical mechanics of dividing cell populations.  _Phys. Rev. Lett._ , 99:138102, Sep 2007.  doi: 10.1103/PhysRevLett.99.138102.  URL https://link.aps.org/doi/10.1103/PhysRevLett.99.138102.
  * Corté et al. [2008] L. Corté, P. M. Chaikin, J. P. Gollub, and D. J. Pine.  Random organization in periodically driven systems.  _Nat. Phys._ , 4:420–424, 2008.  ISSN 17452481.  doi: 10.1038/nphys891.
  * Torquato et al. [2000] S. Torquato, T. M. Truskett, and P. G. Debenedetti.  Is random close packing of spheres well defined?  _Phys. Rev. Lett._ , 84:2064–2067, 2000.  doi: 10.1103/PhysRevLett.84.2064.
  * Spagnolo et al. [2020] B. Spagnolo, A. A. Dubkov, and D. Valenti.  New trends in nonequilibrium statistical mechanics: classical and quantum systems.  _J. Stat. Mech. Theory Exp._ , 2020:094001, 2020.  ISSN 17425468.  doi: 10.1088/1742-5468/abaed0.
  * Huang et al. [2021] M. Huang, W. Hu, S. Yang, Q. X. Liu, and H. P. Zhang.  Circular swimming motility and disordered hyperuniform state in an algae system.  _Proc. Natl. Acad. Sci. U.S.A._ , 118:e2100493118, 2021\.  ISSN 10916490.  doi: 10.1073/pnas.2100493118.
  * Yamada [1961] M. Yamada.  Geometrical Study of the Pair Distribution Function in the Many-Body Problem.  _Prog. Theor. Exp. Phys._ , 25(4):579–594, 04 1961.  ISSN 0033-068X.  doi: 10.1143/PTP.25.579.
  * Costin and Lebowitz [2004] O. Costin and J. Lebowitz.  On the construction of particle distributions with specified single and pair densities.  _J. Phys. Chem. B._ , 108:19614–19618, 2004.  doi: 10.1021/jp047793m.
  * Uche et al. [2006] O. U. Uche, F. H. Stillinger, and S. Torquato.  On the realizability of pair correlation functions.  _Physica A_ , 360:21–36, 2006.  doi: 10.1016/j.physa.2005.03.058.
  * Torquato and Stillinger [2006] S. Torquato and F. H. Stillinger.  New conjectural lower bounds on the optimal density of sphere packings.  _Exp. Math._ , 15:307–331, 2006.  doi: 10.1080/10586458.2006.10128964.
  * Kuna et al. [2007] T. Kuna, J. L. Lebowitz, and E. R. Speer.  Realizability of point processes.  _J. Stat. Phys._ , 129:417–439, 2007.  doi: 10.1007/s10955-007-9393-y.
  * Zhang and Torquato [2020] G. Zhang and S. Torquato.  Realizable hyperuniform and nonhyperuniform particle configurations with targeted spectral functions via effective pair interactions.  _Phys. Rev. E_ , 101:032124, 2020.  doi: 10.1103/PhysRevE.101.032124.
  * Gereben et al. [1994] O. Gereben, L. Pusztai, and A. Baranyai.  Calculation of the three-particle contribution to the configurational entropy for two different models of amorphous Si.  _Phys. Rev. B_ , 49:13251–13254, May 1994.  doi: 10.1103/PhysRevB.49.13251.
  * Jiao et al. [2010] Y. Jiao, F. H. Stillinger, and S. Torquato.  Geometrical ambiguity of pair statistics: Point configurations.  _Phys. Rev. E_ , 81:011105, 2010.  doi: 10.1103/PhysRevE.81.011105.
  * Stillinger and Torquato [2019] F. H. Stillinger and S. Torquato.  Structural degeneracy in pair distance distributions.  _J. Chem. Phys._ , 150:204125, 2019.  doi: 10.1063/1.5096894.
  * Hadley and McCabe [2012] Kevin R. Hadley and Clare McCabe.  Coarse-grained molecular models of water: a review.  _Mol. Simul._ , 38(8-9):671–681, 2012.  doi: 10.1080/08927022.2012.671942.
  * Kmiecik et al. [2016] Sebastian Kmiecik, Dominik Gront, Michal Kolinski, Lukasz Wieteska, Aleksandra Elzbieta Dawid, and Andrzej Kolinski.  Coarse-grained protein models and their applications.  _Chem. Rev._ , 116(14):7898–7936, 2016.  doi: 10.1021/acs.chemrev.6b00163.  URL https://doi.org/10.1021/acs.chemrev.6b00163.
  * Rechtsman et al. [2006] M. C. Rechtsman, F. H. Stillinger, and S. Torquato.  Designed isotropic potentials via inverse methods for self-assembly.  _Phys. Rev. E_ , 73:011406:1–12, 2006.  doi: 10.1103/PhysRevE.73.011406.  Erratum, 75, 019902(E) (2006).
  * Torquato [2009] S. Torquato.  Inverse optimization techniques for targeted self-assembly.  _Soft Matter_ , 5:1157–1173, 2009.  doi: 10.1039/B814211B.
  * Ma et al. [2020] Z. Ma, E. Lomba, and S. Torquato.  Optimized large hyperuniform binary colloidal suspensions in two dimensions.  _Phys. Rev. Lett._ , 125:068002, 2020.  doi: 10.1103/PhysRevLett.125.068002.
  * Sherman et al. [2020] Z. M. Sherman, M. P. Howard, B. A. Lindquist, R. B. Jadrich, and T. M. Truskett.  Inverse methods for design of soft materials.  _J. Chem. Phys._ , 152:140902, 2020.  ISSN 00219606.  doi: 10.1063/1.5145177.
  * Torquato and Stillinger [2003] S. Torquato and F. H. Stillinger.  Local density fluctuations, hyperuniform systems, and order metrics.  _Phys. Rev. E_ , 68:041113, 2003.  doi: 10.1103/PhysRevE.68.041113.
  * Torquato [2018] S. Torquato.  Hyperuniform states of matter.  _Phys. Rep._ , 745:1–95, 2018.  doi: 10.1016/j.physrep.2018.03.001.
  * Uche et al. [2004] O. U. Uche, F. H. Stillinger, and S. Torquato.  Constraints on collective density variables: Two dimensions.  _Phys. Rev. E_ , 70:046122, 2004.
  * Batten et al. [2008] R. D. Batten, F. H. Stillinger, and S. Torquato.  Classical disordered ground states: Super-ideal gases, and stealth and equi-luminous materials.  _J. Appl. Phys._ , 104:033504, 2008.
  * Torquato et al. [2015] S. Torquato, G. Zhang, and F. H. Stillinger.  Ensemble theory for stealthy hyperuniform disordered ground states.  _Phys. Rev. X_ , 5:021020, 2015.  doi: 10.1103/PhysRevX.5.021020.
  * Zhang et al. [2016a] G. Zhang, F. H. Stillinger, and S. Torquato.  Transport, geometrical and topological properties of stealthy disordered hyperuniform two-phase systems.  _J. Chem. Phys_ , 145:244109, 2016a.
  * Zhang et al. [2017] G. Zhang, F. H. Stillinger, and S. Torquato.  Can exotic disordered “stealthy” particle configurations tolerate arbitrarily large holes?  _Soft Matter_ , 13:6197–6207, 2017.  doi: 10.1039/C7SM01028A.
  * Florescu et al. [2009a] M. Florescu, S. Torquato, and P. J. Steinhardt.  Designer disordered materials with large complete photonic band gaps.  _Proc. Nat. Acad. Sci. U. S. A._ , 106:20658–20663, 2009a.  doi: 10.1073/pnas.0907744106.
  * Man et al. [2013] W. Man, M. Florescu, K. Matsuyama, P. Yadak, G. Nahal, S. Hashemizad, E. Williamson, P. Steinhardt, S. Torquato, and P. Chaikin.  Photonic band gap in isotropic hyperuniform disordered solids with low dielectric contrast.  _Opt. Express_ , 21:19972–19981, 2013.  doi: 10.1364/OE.21.019972.
  * Dyson [1970] F. J. Dyson.  Correlations between eigenvalues of a random matrix.  _Comm. Math. Phys._ , 19:235–250, 1970.
  * Torquato et al. [2008] S. Torquato, A. Scardicchio, and C. E. Zachary.  Point processes in arbitrary dimension from Fermionic gases, random matrix theory, and number theory.  _J. Stat. Mech.: Theory Exp._ , 2008:P11019, 2008.
  * Lacroix-A-Chez-Toine et al. [2019] B. Lacroix-A-Chez-Toine, J. A. M. Garzón, C. S. H. Calva, I. P. Castillo, A. Kundu, S. N. Majumdar, and G. Schehr.  Intermediate deviation regime for the full eigenvalue statistics in the complex Ginibre ensemble.  _Phys. Rev. E_ , 100:012137, Jul 2019.  doi: 10.1103/PhysRevE.100.012137.  URL https://link.aps.org/doi/10.1103/PhysRevE.100.012137.
  * Hansen [1973] J. P. Hansen.  Statistical mechanics of dense ionized matter. i. equilibrium properties of the classical one-component plasma.  _Phys. Rev. A_ , 8:3096–3109, Dec 1973.  doi: 10.1103/PhysRevA.8.3096.  URL https://link.aps.org/doi/10.1103/PhysRevA.8.3096.
  * Gann et al. [1979] R. C. Gann, S. Chakravarty, and G. V. Chester.  Monte Carlo simulation of the classical two-dimensional one-component plasma.  _Phys. Rev. B_ , 20:326–344, Jul 1979.  doi: 10.1103/PhysRevB.20.326.  URL https://link.aps.org/doi/10.1103/PhysRevB.20.326.
  * Dyson [1962] F. J. Dyson.  Statistical theory of the energy levels of complex systems. I.  _J. Math. Phys._ , 3:140–156, 1962.  doi: 10.1063/1.1703773.
  * Lewin [2022] Mathieu Lewin.  Coulomb and riesz gases: The known and the unknown.  _Journal of Mathematical Physics_ , 63(6):061101, 2022.  doi: 10.1063/5.0086835.
  * Marcotte et al. [2013a] É. Marcotte, F. H. Stillinger, and S. Torquato.  Nonequilibrium static growing length scales in supercooled liquids on approaching the glass transition.  _J. Chem. Phys._ , 138:12A508, 2013a.  doi: 10.1063/1.4769422.
  * Zhang et al. [2016b] G. Zhang, F. H. Stillinger, and S. Torquato.  The perfect glass paradigm: Disordered hyperuniform glasses down to absolute zero.  _Sci. Rep._ , 6:36963, 2016b.  doi: 10.1038/srep36963.
  * Torquato and Stillinger [2007] S. Torquato and F. H. Stillinger.  Toward the jamming threshold of sphere packings: Tunneled crystals.  _J. Appl. Phys._ , 102:093511, 2007.  Erratum, 103, 129902 (2008).
  * Atkinson et al. [2016] S. Atkinson, G. Zhang, A. B. Hopkins, and S. Torquato.  Critical slowing down and hyperuniformity on approach to jamming.  _Phys. Rev. E_ , 94:012902, 2016.
  * Coniglio et al. [2017] A. Coniglio, M. Pica Ciamarra, and T. Aste.  Universal behaviour of the glass and the jamming transitions in finite dimensions.  _Soft Matter_ , 13:8766–8771, 2017.
  * Rissone et al. [2021] P. Rissone, E. I. Corwin, and G. Parisi.  Long-range anomalous decay of the correlation in jammed packings.  _Phys. Rev. Lett._ , 127:038001, Jul 2021.  doi: 10.1103/PhysRevLett.127.038001.  URL https://link.aps.org/doi/10.1103/PhysRevLett.127.038001.
  * Torquato [2021] S Torquato.  Structural characterization of many-particle systems on approach to hyperuniform states.  _Phys. Rev. E_ , 103:052126, 2021.
  * Ghosh and Lebowitz [2018] S. Ghosh and J. L. Lebowitz.  Generalized stealthy hyperuniform processes: Maximal rigidity and the bounded holes conjecture.  _Comm. Math. Phys._ , 363:97–110, 2018.  doi: 10.1007/s00220-018-3226-5.
  * Jiao et al. [2014] Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, J. C. Corbo, and S. Torquato.  Avian photoreceptor patterns represent a disordered hyperuniform solution to a multiscale packing problem.  _Phys. Rev. E_ , 89:022721, 2014.  doi: 10.1103/PhysRevE.89.022721.
  * Mayer et al. [2015] A. Mayer, V. Balasubramanian, T. Mora, and A. M. Walczak.  How a well-adapted immune system is organized.  _Proc. Nat. Acad. Sci. U. S. A._ , 112:5950–5955, 2015\.  doi: 10.1073/pnas.1421827112.
  * Florescu et al. [2009b] M. Florescu, S. Torquato, and P. J. Steinhardt.  Complete band gaps in two-dimensional photonic quasicrystals.  _Phys. Rev. B_ , 80:155112, 2009b.
  * Marcotte et al. [2013b] É. Marcotte, F. H. Stillinger, and S. Torquato.  Designed diamond ground state via optimized isotropic monotonic pair potentials.  _J. Chem. Phys._ , 138:061101, 2013b.
  * Zhang et al. [2019] H. Zhang, H. Chu, H. Giddens, W. Wu, and Y. Hao.  Experimental demonstration of luneburg lens based on hyperuniform disordered media.  _Appl. Phys. Lett._ , 114:053507, 2019.
  * Sgrignuoli and Dal Negro [2021] F. Sgrignuoli and L. Dal Negro.  Hyperuniformity and wave localization in pinwheel scattering arrays.  _Phys. Rev. B_ , 103:224202, Jun 2021.  doi: 10.1103/PhysRevB.103.224202.  URL https://link.aps.org/doi/10.1103/PhysRevB.103.224202.
  * Christogeorgos et al. [2021] O. Christogeorgos, H. Zhang, Q. Cheng, and Y. Hao.  Extraordinary Directive Emission and Scanning from an Array of Radiation Sources with Hyperuniform Disorder.  _Phys. Rev. Appl._ , 15:014062, 2021.  doi: 10.1103/PhysRevApplied.15.014062.
  * Hexner and Levine [2015] D. Hexner and D. Levine.  Hyperuniformity of critical absorbing states.  _Phys. Rev. Lett._ , 114:110602, 2015.  doi: 10.1103/PhysRevLett.114.110602.
  * Hexner et al. [2017] D. Hexner, P. M. Chaikin, and D. Levine.  Enhanced hyperuniformity from random reorganization.  _Proc. Nat. Acad. Sci. U. S. A._ , 114:4294–4299, 2017\.  doi: 10.1073/pnas.1619260114.
  * Ma and Torquato [2019] Z. Ma and S. Torquato.  Hyperuniformity of generalized random organization models.  _Phys. Rev. E_ , 99:022115, 2019.  doi: 10.1103/PhysRevE.99.022115.
  * Lei et al. [2019] Q.-L. Lei, M. Pica Ciamarra, and R. Ni.  Nonequilibrium strongly hyperuniform fluids of circle active particles with large local density fluctuations.  _Sci. Adv._ , 5(1):eaau7423, 2019.  doi: 10.1126/sciadv.aau7423.  URL https://www.science.org/doi/abs/10.1126/sciadv.aau7423.
  * Ding et al. [2018] Z. Ding, Y. Zheng, Y. Xu, Y. Jiao, and W. Li.  Hyperuniform flow fields resulting from hyperuniform configurations of circular disks.  _Phys. Rev. E_ , 98:063101, 2018.
  * Duerinckx and Gloria [2022] M. Duerinckx and A. Gloria.  Sedimentation of random suspensions and the effect of hyperuniformity.  _Ann. PDE_ , 8(1):1–66, 2022.  doi: 10.1007/s40818-021-00115-0.
  * Reatto and Chester [1967] L. Reatto and G. V. Chester.  Phonons and the properties of a Bose system.  _Phys. Rev._ , 155:88–100, 1967.
  * Feynman [1998] R. P. Feynman.  _Statistical Mechanics_.  Westview Press, Boulder, Colorado, 1998.
  * Abreu et al. [2017] L. D. Abreu, J. M. Pereira, J. L. Romero, and S. Torquato.  The Weyl-Heisenberg ensemble: Hyperuniformity and higher landau levels.  _J. Stat. Mech.: Th. and Exper._ , 2017:043103, 2017.
  * Crowley et al. [2019] P. J. D. Crowley, C. R. Laumann, and S. Gopalakrishnan.  Quantum criticality in Ising chains with random hyperuniform couplings.  _Phys. Rev. B_ , 100:134206, Oct 2019.  doi: 10.1103/PhysRevB.100.134206.  URL https://link.aps.org/doi/10.1103/PhysRevB.100.134206.
  * Sarnak and Strömbergsson [2006] P. Sarnak and A. Strömbergsson.  Minima of Epstein’s zeta function and heights of flat tori.  _Inventiones Math._ , 165:115–151, 2006.
  * Brauchart et al. [2019] J. S. Brauchart, P. J. Grabner, and W. Kusner.  Hyperuniform point sets on the sphere: deterministic aspects.  _Constructive approximation_ , 50:45–61, 2019.
  * Torquato et al. [2018] S. Torquato, G. Zhang, and M. de Courcy-Ireland.  Uncovering multiscale order in the prime numbers via scattering.  _J. Stat. Mech.: Theory & Exper._, 2018:093401, 2018.
  * Torquato et al. [2019] S. Torquato, G. Zhang, and M. de Courcy-Ireland.  Hidden multiscale order in the primes.  _J. Phys. A: Math. & Theoretical_, 52:135002, 2019.
  * Brauchart et al. [2020] J. S. Brauchart, P. J. Grabner, W. Kusner, and J. Ziefle.  Hyperuniform point sets on the sphere: probabilistic aspects.  _Monatsh. Math._ , 192(4):763–781, 2020.  doi: 10.1007/s00605-020-01439-y.
  * Donev et al. [2005] A. Donev, F. H. Stillinger, and S. Torquato.  Unexpected density fluctuations in disordered jammed hard-sphere packings.  _Phys. Rev. Lett._ , 95:090604, 2005.  doi: 10.1103/PhysRevLett.95.090604.
  * Torquato and Wang [2022] S. Torquato and H. Wang.  Precise determination of pair interactions from pair statistics of many-body systems in and out of equilibrium.  _Phys. Rev. E_ , 106:044122, Oct 2022.  doi: 10.1103/PhysRevE.106.044122.  URL https://link.aps.org/doi/10.1103/PhysRevE.106.044122.
  * Levesque et al. [1985] D. Levesque, J. J. Weis, and L. Reatto.  Pair interaction from structural data for dense classical liquids.  _Phys. Rev. Lett._ , 54:451–454, Feb 1985.  doi: 10.1103/PhysRevLett.54.451.  URL https://link.aps.org/doi/10.1103/PhysRevLett.54.451.
  * Lyubartsev and Laaksonen [1995] A. P. Lyubartsev and A. Laaksonen.  Calculation of effective interaction potentials from radial distribution functions: A reverse monte carlo approach.  _Phys. Rev. E_ , 52:3730–3737, Oct 1995.  doi: 10.1103/PhysRevE.52.3730.  URL https://link.aps.org/doi/10.1103/PhysRevE.52.3730.
  * Soper [1996] A. K. Soper.  Empirical potential Monte Carlo simulation of fluid structure.  _Chem. Phys._ , 202:295–306, 1996.
  * Heinen [2018] M. Heinen.  Calculating particle pair potentials from fluid-state pair correlations: Iterative Ornstein–Zernike inversion.  _J. Comput. Chem._ , 39:1531–1543, 2018.
  * Levine and Steinhardt [1984] D. Levine and P. J. Steinhardt.  Quasicrystals: A new class of ordered structures.  _Phys. Rev. Lett._ , 53:2477–2480, 1984.
  * Wang et al. [2020] H. Wang, F. H. Stillinger, and S. Torquato.  Sensitivity of pair statistics on pair potentials in many-body systems.  _J. Chem. Phys._ , 150:124106, 2020.
  * Henderson [1974] R. L. Henderson.  A uniqueness theorem for fluid pair correlation functions.  _Phys. Lett. A_ , 49:197–198, 1974.
  * Torquato and Stillinger [2008] S. Torquato and F. H. Stillinger.  New duality relations for classical ground states,.  _Phys. Rev. Lett._ , 100:020602, 2008.
  * Stillinger and Weber [1982] Frank H. Stillinger and Thomas A. Weber.  Hidden structure in liquids.  _Phys. Rev. A_ , 25:978–989, Feb 1982.  doi: 10.1103/PhysRevA.25.978.  URL https://link.aps.org/doi/10.1103/PhysRevA.25.978.
  * Stillinger and Weber [1985a] F. H. Stillinger and T. A. Weber.  Inherent structure of liquids in the hard-sphere limit.  _J. Chem. Phys._ , 83:4767–775, 1985a.  doi: 10.1063/1.449840.
  * Stillinger [1998] F. H. Stillinger.  Enumeration of isobaric inherent structures for the fragile glass former o-terphenyl.  _J. Phys. Chem. B_ , 102(15):2807–2810, 1998\.  doi: 10.1021/jp973144h.
  * 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.  Cloaking the underlying long-range order of randomly perturbed lattices.  _Phys. Rev. E_ , 101:032118, 2020.  doi: 10.1103/PhysRevE.101.032118.
  * Torquato et al. [1990a] S. Torquato, B. Lu, and J. Rubinstein.  Nearest-neighbor distribution functions in many-body sytems.  _Phys. Rev. A_ , 41:2059–2075, 1990a.  doi: 10.1103/PhysRevA.41.2059.
  * Hansen and McDonald [1986] J. P. Hansen and I. R. McDonald.  _Theory of Simple Liquids_.  Academic Press, New York, 1986.
  * Zachary and Torquato [2009] C. E. Zachary and S. Torquato.  Hyperuniformity in point patterns and two-phase heterogeneous media.  _J. Stat. Mech.: Theory & Exp._, 2009:P12015, 2009.  doi: 10.1088/1742-5468/2009/12/P12015.
  * Torquato [2016] S. Torquato.  Disordered hyperuniform heterogeneous materials.  _J. Phys.: Cond. Mat_ , 28:414012, 2016.
  * Torquato et al. [1990b] S. Torquato, B. Lu, and J. Rubinstein.  Nearest-neighbor distribution function for systems of interacting particles.  _J. Phys. A: Math. Gen._ , 23:L103–L107, 1990b.  doi: 10.1088/0305-4470/23/3/005.
  * Middlemas and Torquato [2020] T. M. Middlemas and S. Torquato.  Nearest-neighbor functions for disordered stealthy hyperuniform many-particle systems.  _J. Stat. Mech.: Theory & Exp._, 2020:103302, 2020.  doi: 10.1088/1742-5468/abb8cb.
  * Torquato et al. [2006] S. Torquato, O. U. Uche, and F. H. Stillinger.  Random sequential addition of hard spheres in high Euclidean dimensions.  _Phys. Rev. E_ , 74:061308, 2006.  doi: 10.1103/PhysRevE.74.061308.
  * Zhang et al. [2013] G. Zhang, F. H. Stillinger, and S. Torquato.  Probing the limitations of isotropic pair potentials to produce ground-state structural extremes via inverse statistical mechanics.  _Phys. Rev. E_ , 88:042309, 2013.  doi: 10.1103/PhysRevE.88.042309.
  * Yu et al. [2021] S. Yu, C.-W. Qiu, Y. Chong, S. Torquato, and N. Park.  Engineered disorder in photonics.  _Nature Rev. Mater._ , 6:226–243, 2021.
  * Wilken et al. [2021] S. Wilken, R. E. Guerra, D. Levine, and P. M. Chaikin.  Random close packing as a dynamical phase transition.  _Phys. Rev. Lett._ , 127:038002, Jul 2021.  doi: 10.1103/PhysRevLett.127.038002.  URL https://link.aps.org/doi/10.1103/PhysRevLett.127.038002.
  * Reichhardt and Reichhardt [2014] C. Reichhardt and C. J. O. Reichhardt.  Absorbing phase transitions and dynamic freezing in running active matter systems.  _Soft Matter_ , 10(38):7502–7510, 2014.  doi: 10.1039/C4SM01273A.
  * Wang and Torquato [2022] H. Wang and S. Torquato.  Dynamic measure of hyperuniformity and nonhyperuniformity in heterogeneous media via the diffusion spreadability.  _Phys. Rev. Appl._ , 17:034022, 2022.  doi: 10.1103/PhysRevApplied.17.034022.
  * Zhang et al. [2017] G. Zhang, F. H. Stillinger, and S. Torquato.  Classical many-particle systems with unique disordered ground states.  _Phys. Rev. E_ , 96:042146, 2017.  doi: 10.1103/PhysRevE.96.042146.
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal.  On the limited memory BFGS method for large scale optimization.  _Math. Program._ , 45:503–528, 1989.  doi: 10.1007/BF01589116.
  * Henkel et al. [2008] M. Henkel, H. Hinrichsen, and S. Lübeck.  _Non-Equilibrium Phase Transitions - Volume 1: Absorbing Phase Transitions_.  Springer, New York, 2008.
  * Pine et al. [2005] D. J. Pine, J. P. Gollub, J. F. Brady, and A. M. Leshansky.  Chaos and threshold for irreversibility in sheared suspensions.  _Nature_ , 438:997–1000, 2005.  doi: 10.1038/nature04380.
  * Welberry et al. [1980] T. R. Welberry, G. H. Miller, and C. E. Carroll.  Paracrystals and growth-disorder models.  _Acta Crystallogr., Sect. A: Found. Crystallogr._ , 36:921–929, 1980.
  * Gabrielli et al. [2002] A. Gabrielli, M. Joyce, and F. Sylos Labini.  Glass-like universe: Real-space correlation properties of standard cosmological models.  _Phys. Rev. D_ , 65:083523, 2002.
  * Gabrielli and Torquato [2004] A. Gabrielli and S. Torquato.  Voronoi and void statistics for superhomogeneous point processes.  _Phys. Rev. E_ , 70:041105, 2004.
  * Kim and Torquato [2018] J. Kim and S. Torquato.  Effect of imperfections on the hyperuniformity of many-body systems.  _Phys. Rev. B_ , 97:054105, 2018.  doi: 10.1103/PhysRevB.97.054105.
  * Ornstein and Zernike [1914] L. S. Ornstein and F. Zernike.  Accidental deviations of density and opalescence at the critical point of a single substance.  _Proc. Akad. Sci. (Amsterdam)_ , 17:793–806, 1914.
  * Note [1] Note1.  Such background terms have been employed to study numerically the one-component plasma [33, 34] and the Dyson log gas. [35].
  * Ewald [1921] P. P. Ewald.  Die berechnung optischer und elektrostatischer gitterpotentiale.  _Ann. Phys._ , 369(3):253–287, 1921.  doi: https://doi.org/10.1002/andp.19213690304.
  * Note [2] Note2.  If the target and optimized systems have different largest detected hole radii, $r_{\text{max}}$ is chosen to be the smaller of the two. The aforementioned difference is within $0.1\rho^{1/d}$ for all three models.
  * Chen et al. [2021] D. Chen, Y. Zheng, L. Liu, G. Zhang, M. Chen, Y. Jiao, and H. Zhuang.  Stone–wales defects preserve hyperuniformity in amorphous two-dimensional networks.  _Proc. Natl. Acad. Sci._ , 118:e2016862118, 2021.
  * Hejna et al. [2013] M. Hejna, P. J. Steinhardt, and S. Torquato.  Nearly hyperuniform network models of amorphous silicon.  _Phys. Rev. B_ , 87:245204, 2013.
  * Lindsay and Broido [2010] L. Lindsay and D. A. Broido.  Optimized tersoff and brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene.  _Phys. Rev. B_ , 81:205441, May 2010.  doi: 10.1103/PhysRevB.81.205441.  URL https://link.aps.org/doi/10.1103/PhysRevB.81.205441.
  * Stillinger and Weber [1985b] F. H. Stillinger and T. A. Weber.  Computer simulation of local order in condensed phases of silicon.  _Phys. Rev. B_ , 31:5262–5271, 1985b.
  * Kamien and Selinger [2001] R. D. Kamien and J. V. Selinger.  Order and frustration in chiral liquid crystals.  _J. Phys. Cond. Matter_ , 13(3):R1, 2001.  doi: 10.1088/0953-8984/13/3/201.
  * Zhang and Snezhko [2022] B. Zhang and A. Snezhko.  Hyperuniform active chiral fluids with tunable internal structure.  _Phys. Rev. Lett._ , 128:218002, May 2022.  doi: 10.1103/PhysRevLett.128.218002.  URL https://link.aps.org/doi/10.1103/PhysRevLett.128.218002.
  * Bi et al. [2016] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning.  Motility-driven glass and jamming transitions in biological tissues.  _Phys. Rev. X_ , 6:021011, Apr 2016.  doi: 10.1103/PhysRevX.6.021011.
  * Zhang et al. [2021] J. Zhang, W. Zheng, S. Zhang, D. Xu, Y. Nie, Z. Jiang, and N. Xu.  Unifying fluctuation-dissipation temperatures of slow-evolving nonequilibrium systems from the perspective of inherent structures.  _Sci. Adv._ , 7(31):eabg6766, 2021.  doi: 10.1126/sciadv.abg6766.
  * Miller and Frenkel [2004] M. A. Miller and D. Frenkel.  Phase diagram of the adhesive hard sphere fluid.  _J. Chem. Phys._ , 121(1):535–545, 2004.  doi: 10.1063/1.1758693.
  * Pàmies et al. [2009] J. C. Pàmies, A. Cacciuto, and D. Frenkel.  Phase diagram of Hertzian spheres.  _J. Chem. Phys._ , 131(4):044514, 2009.  doi: 10.1063/1.3186742.
  * Lomba et al. [2017] E. Lomba, J-J. Weis, and S. Torquato.  Disordered hyperuniformity in two-component non-additive hard disk plasmas.  _Phys. Rev. E_ , 96:062126, 2017.
  * Zu et al. [2020] M. Zu, A. Bupathy, D. Frenkel, and S. Sastry.  Information density, structure and entropy in equilibrium and non-equilibrium systems.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2020(2):023204, feb 2020.  doi: 10.1088/1742-5468/ab684b.  URL https://dx.doi.org/10.1088/1742-5468/ab684b.
  * Press et al. [1986] W H Press, B P Flannery, S Teukolsky, and W T Vetterling.  _Numerical recipes, The art of scientific computing_.  Cambridge University Press, 1 1986.  URL https://www.osti.gov/biblio/5271321.
  * Prestipino et al. [2005] S. Prestipino, F. Saija, and P. V. Giaquinta.  Phase diagram of the gaussian-core model.  _Phys. Rev. E_ , 71:050102, May 2005.  doi: 10.1103/PhysRevE.71.050102.  URL https://link.aps.org/doi/10.1103/PhysRevE.71.050102.
  * Travesset [2014] A. Travesset.  Phase diagram of power law and lennard-jones systems: Crystal phases.  _J. Chem. Phys._ , 141(16):164501, 2014.  doi: 10.1063/1.4898371.
  * Stillinger and Weber [1981] F. H. Stillinger and T. A. Weber.  Gaussian core model in two dimensions. i. melting transition.  _J. Chem. Phys._ , 74(7):4015–4019, 1981.  doi: 10.1063/1.441581.
  * Kirkpatrick et al. [1983] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi.  Optimization by simulated annealing.  _Science_ , 220:671–680, 1983.
  * Note [4] Note4.  With cooling rate $T_{i+1}=0.999T_{i}$, where $i$ is the iteration index, the annealing procedures still result in metastable states for the hyperuniform models.
  * Stillinger [1976] F. H. Stillinger.  Phase transitions in the Gaussian core system.  _J. Chem. Phys._ , 65:3968, 1976.
  * Zachary et al. [2008] C. E. Zachary, F. H. Stillinger, and S. Torquato.  Gaussian-core model phase diagram and pair correlations in high Euclidean dimensions.  _J. Chem. Phys._ , 128:224505, 2008.