# Designer Pair Statistics of Disordered Many-Particle Systems with Novel
Properties

Haina Wang  _Department of Chemistry, Princeton University_ , Princeton, New
Jersey, 08544, USA  Salvatore Torquato  Email<EMAIL_ADDRESS>Department of Chemistry, Department of Physics, Princeton Materials Institute,
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

Knowledge of exact analytical functional forms for the pair correlation
function $g_{2}(r)$ and its corresponding structure factor $S(k)$ of
disordered many-particle systems is limited. For fundamental and practical
reasons, it is highly desirable to add to the existing data base of analytical
functional forms for such pair statistics. Here, we design a plethora of such
pair functions in direct and Fourier spaces across the first three Euclidean
space dimensions that are realizable by diverse many-particle systems with
varying degrees of correlated disorder across length scales, spanning a wide
spectrum of hyperuniform, typical nonhyperuniform and antihyperuniform ones.
This is accomplished by utilizing an efficient inverse algorithm that
determines equilibrium states with up to pair interactions at positive
temperature that precisely match targeted forms for both $g_{2}(r)$ and
$S(k)$. Among other results, we realize an example with the strongest
hyperuniform property among known positive-temperature equilibrium states,
critical-point systems (implying unusual 1D systems with phase transitions)
that are not in the Ising universality class, systems that attain self-similar
pair statistics under Fourier transformation, and an experimentally feasible
polymer model. We show that our pair functions enable one to achieve systems
with a wide range of translational order and self-diffusion coefficients $\cal
D$, which are inversely related to one another. One can design other
realizable pair statistics via linear combinations of our functions or by
applying our inverse procedure to other desirable functional forms. Our
approach facilitates the inverse design of materials with desirable physical
and chemical properties by tuning their pair statistics.

##  I Introduction

The determination of equilibrium and nonequilibrium properties of disordered
many-particle systems in $d$-dimensional Euclidean space $\mathbb{R}^{d}$
based on structural information contained in correlation functions is of
fundamental and practical importance in statistical and condensed-matter
physics, chemistry, mathematics, biological sciences and materials engineering
[1, 2, 3, 4, 5, 6, 7, 8]. The pair correlation function $g_{2}(\textbf{r})$
and the corresponding structure factor $S(\textbf{k})$ are widely used
structural functions due to their computational simplicity and experimental
accessibility through diffraction measurements. They enable one to determine
key equilibrium bulk properties, such as the pressure, the isothermal
compressibility [7, 4] and the two-body excess entropy [9]. Moreover, pair
statistics are crucial inputs for the determination of transport properties of
many-particle systems, including Enskog theories for the self-diffusion
coefficient [10, 11, 12] and the shear and bulk viscosities [13], as well as
effective temperatures of structural glasses [14] and the time-dependent
diffusion spreadability [15].

Despite the importance of such pair functions in direct and Fourier spaces
[i.e., $g_{2}(\textbf{r})$ and $S(\textbf{k})$, respectively], there are only
a relatively small number of disordered many-particle systems for which we
know exact analytical pair functions over their entire parameter regime. These
include determinantal point processes across space dimensions [2, 16, 1, 17,
18, 19, 20, 21], equilibrium one-dimensional (1D) hard rods [22], “ghost”
random sequential addition (RSA) process across dimensions [23], as well as 2D
and 3D unit-step function $g_{2}(\textbf{r})$ at all densities up to the
“terminal packing fraction” $1/2^{d}$ [24]. 1515footnotetext: These include
$d$-dimensional ground-state free fermions [2, 16], eigenvalues for certain
random matrices [1, 17, 18], 1D Lorentzian $g_{2}(r)$ [19] and 2D Ginibre
ensemble [20]. Given our current limited knowledge of analytical pair
functions, it is highly desirable to expand the existing database of radial
functional forms of $g_{2}(r)$ and $S(k)$ that can be achieved by
statistically homogeneous and isotropic disordered many-particle systems,
where $r=|\mathbf{r}|$ and $k=|\mathbf{k}|$. Substantially increasing such a
database would enable one to probe the space of equilibrium and nonequilibrium
properties analytically in comprehensive ways not heretofore possible.
Furthermore, adding to such a database would facilitate the inverse design of
materials with desirable mechanical, optical and chemical properties, which
could be further accelerated by machine learning techniques [25].

A hypothesized functional form for the pair correlation function $g_{2}(r)$ is
not necessarily realized by a many-particle system at number density $\rho$,
even if the corresponding structure factor $S(k)$ is nonnegative for all $k$,
as is should be [26, 19, 27]. This realizability problem has a rich and long
history [28, 29, 26, 19, 27, 30, 31], but it is a notoriously difficult
problem to solve because an uncountable number of necessary and sufficient
conditions are generally required [30, 31], which are not possible to check in
practice. This places great importance on the need to formulate algorithms to
construct particle configurations that realize targeted hypothetical
functional forms of the pair statistics [32, 24].

In this work, we design a plethora of exact functional forms for pair
correlations in direct and Fourier spaces across the first three Euclidean
space dimension that are realizable by diverse many-particle systems with
varying degrees of correlated disorder across length scales with novel
structural and physical properties. As suggested by Zhang and Torquato [32],
we solve this realizability problem here by employing an efficient inverse
algorithm [33] that determines equilibrium states with up to pair interactions
$v(r)$ at positive temperature $T$ that precisely match targeted forms for
both $g_{2}(r)$ and $S(k)$. Our methodology minimizes an objective function
that incorporates pair statistics in both direct and Fourier spaces, so that
both the small-, intermediate- and large-distance correlations are very
accurately captured; see Sec. II for details.

For purposes of illustration, we consider a broad range of functional forms
that span across a wide spectrum of hyperuniform and nonhyperuniform classes.
Disordered hyperuniform systems anomalously suppress large-scale density
fluctuations compared to typical disordered systems, such as liquids [34]. A
point pattern in $\mathbb{R}^{d}$ is hyperuniform if the local number variance
$\sigma^{2}(R)=\langle N^{2}(R)\rangle-\langle N(R)\rangle^{2}$ of the number
of points $N(R)$ found in a spherical window of radius $R$ grows slower than
$R^{d}$ in the large-$R$ limit [34, 35]. Equivalently, a hyperuniform state is
characterized by a structure factor that vanishes in the infinite-wavelength
limit, i.e., $\lim_{k\rightarrow 0}S(k)=0$. Often, the structure factor of
hyperuniform systems takes a radial power-law form in the vicinity of the
origin, i.e., $S(k)\sim k^{\alpha}$, where $\alpha$ is a positive exponent.
The value of $\alpha$ determines different large-$R$ scaling behaviors of the
number variance [34, 36, 35], 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}$ (1)

Disordered hyperuniform systems are exotic states of matter that have been
receiving great attention because they connect a variety of seemingly
unrelated systems that arise in physics, chemistry, materials science,
mathematics, and biology [35]. By contrast, for any nonhyperuniform system,
the local variance has the following large-$R$ scaling behaviors [37]:

$\sigma^{2}(R)\sim\begin{cases}R^{d}\quad\alpha=0\text{ (typical
nonhyperuniform)}\\\ R^{d-\alpha}\quad\alpha<0\text{ (antihyperuniform).}\\\
\end{cases}$ (2)

For a “typical” nonhyperuniform system, $S(0)$ is bounded. In antihyperuniform
systems, $S(0)$ is unbounded, i.e., $\lim_{k\rightarrow 0}S(k)=+\infty$, and
hence are diametrically opposite to hyperuniform systems. Typical
nonhyperuniform states have bounded nonzero value of $S(k)$ as $k\rightarrow
0$ [37]. Antihyperuniform systems include those at thermal critical points
(e.g., liquid-vapor and magnetic critical points) [38, 39, 40, 41, 42], as
well as fractals [43] and certain substitution tilings [44].

Table I lists the target $g_{2}(r)$ and $S(k)$ functions studied in this work.
The hyperuniform pair functions that we target include $d$-dimensional
Gaussians, $d$-dimensional generalization of the one-component plasma (OCP)
and their Fourier duals [32], $d$-dimensional hyperbolic secant function
$g_{2}(r)$, and a 2D Gaussian-damped polynomial. The typical nonhyperuniform
functions include a 1D Hermite-Gaussian function, a hypothetical 3D
hyposurficial state [34], as well as 1D and 2D ghost RSA packings [23]. These
functional forms are chosen for their fundamental significance, their
capability to yield systems in which particles form clusters of various sizes,
and their potential applications in materials design. It is noteworthy that
the Gaussian functions, the 1D hyperbolic secant function, the 2D Gaussian-
damped polynomial and the 1D Hermite-Gaussian function have the property that
they are self-similar under Fourier transformations. Self-similar pair
statistics enable one to exploit duality relations in the phase diagrams of
soft-matter systems. For example, via duality relations, information about
low-density classical ground states under short-ranged interactions can be
mapped to high-density ground states under long-ranged interactions [45, 46].
This knowledge facilitates the determination of ground-state structures over
large density ranges [46], which may provide good solutions to covering and
quantizer problems that arise in discrete geometry and have many applications
[47], including wireless communication network layouts.

Table 1. The analytical target pair statistics studied here and shown to be realizable by statistically homogeneous and isotropic disordered many-particle systems in various spatial dimensions. The small-$k$ scaling exponents $\alpha$ in Eq. (1) are shown. All states are at $\rho=1$. The definitions of the special functions and symbols in the pair statistics are given in the table footnotes. Functional form |  $d$ |  $g_{2}(r)$ |  |  $S(k)$ |  |  $\alpha$
---|---|---|---|---|---|---
Gaussian |  1, 2, 3 |  $1-\exp(-\pi r^{2})$ |  (3) |  $1-\exp(-\frac{k^{2}}{4\pi})$ |  (4) |  2
Generalized OCP[48] |  3 |  $1-\exp(-\frac{4\pi}{3}r^{3})$ |  (5) |  $\begin{aligned} &1-\,_{1}F_{4}\left(1;\frac{1}{3},\frac{2}{3},\frac{5}{6},\frac{7}{6};-\frac{k^{6}}{20736\pi^{2}}\right)\\\ &+\frac{2\pi\text{ber}_{\frac{2}{3}}\left(\frac{\left(k^{2}\right)^{3/4}}{3\sqrt{\pi}}\right)}{3\sqrt{3}}\\\ &+\frac{2\pi\text{ber}_{-\frac{2}{3}}\left(\frac{\left(k^{2}\right)^{3/4}}{3\sqrt{\pi}}\right)}{3\sqrt{3}}\end{aligned}$ |  (6) |  2
Fourier dual of OCP |  1 |  $1-\frac{1}{\pi r^{2}+1}$ |  (7) |  $1-\exp(-k/\pi)$ |  (8) |  1
Fourier dual of OCP[48] |  3 |  $\begin{aligned} &1-\,_{1}F_{4}\left(1;\frac{1}{3},\frac{2}{3},\frac{5}{6},\frac{7}{6};-\frac{1}{324}\pi^{4}r^{6}\right)\\\ &+\frac{2\pi\text{ber}_{\frac{2}{3}}\left(\frac{2}{3}\sqrt{2}\pi\left(r^{2}\right)^{3/4}\right)}{3\sqrt{3}}\\\ &+\frac{2\pi\text{ber}_{-\frac{2}{3}}\left(\frac{2}{3}\sqrt{2}\pi\left(r^{2}\right)^{3/4}\right)}{3\sqrt{3}}\end{aligned}$ |  (9) |  $1-\exp(-\frac{2}{3(2\pi)^{2}}k^{3})$ |  (10) |  3
Hyperbolic secant $g_{2}(r)$ |  1 |  $1-\operatorname{sech}(\pi r)$ |  (11) |  $1-\operatorname{sech}(k/2)$ |  (12) |  2
Hyperbolic secant $g_{2}(r)$[49] |  2 |  $1-\operatorname{sech}(2\sqrt{\pi G}r)$ |  (13) |  No closed analytical form |  |  2
Hyperbolic secant $g_{2}(r)$ |  3 |  $1-\operatorname{sech}(\frac{\pi^{4/3}}{2^{1/3}}r)$ |  (14) |  $\displaystyle 1-\frac{2^{2/3}\sqrt[3]{\pi}\tanh\left(\frac{k}{2^{2/3}\sqrt[3]{\pi}}\right)\operatorname{sech}\left(\frac{k}{2^{2/3}\sqrt[3]{\pi}}\right)}{k}$ |  (15) |  2
Gaussian-damped polynomial |  2 |  $\displaystyle 1-\frac{1}{3}e^{-\pi r^{2}}\left(4\pi^{2}r^{4}-8\pi r^{2}+3\right)$ |  (16) |  $\displaystyle 1-\frac{e^{-\frac{k^{2}}{4\pi}}\left(k^{4}-8\pi k^{2}+12\pi^{2}\right)}{12\pi^{2}}$ |  (17) |  2
Hermite-Gaussian[50] |  1 |  $1-\lambda\left(4r^{4}-12r^{2}+3\right)e^{-\frac{r^{2}}{2}}$ |  (18) |  $1-\lambda\sqrt{2\pi}\left(-4k^{4}+12k^{2}-3\right)e^{-\frac{k^{2}}{2}}$ |  (19) |  0
Hyposurficial[51] |  3 |  $\displaystyle 1+\frac{\exp(-r^{*})}{4\pi}-\frac{\exp(-r^{*})\sin(r^{*})}{r^{*}}$ |  (20) |  $\displaystyle\frac{6k^{*8}+12k^{*6}+19k^{*4}+24k^{*2}+16}{6(k^{*2}+1)^{2}(k^{*2}-2k^{*}+2)(k^{*2}+2k^{*}+2)}$ |  (21) |  0
Ghost RSA[52] |  1 |  $\displaystyle\frac{2\Theta(2r-1)}{2-(1-r)\Theta(2-2r)}$ |  (22) |  $1-4\cos(k)\left[\text{Ci}\left(\frac{3k}{2}\right)+\text{Ci}(2k)\right]-4\sin(k)\left[\text{Si}\left(\frac{3k}{2}\right)+\text{Si}(2k)\right]-\frac{2\sin(k)}{k}$ |  (23) |  0
Ghost RSA[53] |  2 |  $\displaystyle\frac{2\Theta(r^{*}-1)}{2-\frac{2}{\pi}\left[\cos^{-1}\left(\frac{r^{*}}{2}\right)-\frac{r}{2}\sqrt{1-\frac{r^{*2}}{4}}\right]\Theta(2-r^{*})}$ |  (24) |  No closed analytical form |  |  0
Antihyperuniform[54] |  1, 2, 3 |  $\displaystyle\Theta(r/D-1)\left(\frac{A}{(r/D)^{d-\frac{1}{2}}}+1\right)$ |  (25) |  See Appendix |  |  $-\frac{1}{2}$

We also, for the first time, demonstrate the realizability of antihyperuniform
systems across the first three space dimensions corresponding to critical
points that do not belong to the Ising universality class, including a
remarkable 1D critical-point state, which is to be contrasted with typical 1D
systems in which phase transitions cannot exist. Note that the pair functions
in Table I possess analytical forms in both direct and Fourier spaces in all
but two cases, the 2D hyperbolic secant $g_{2}(r)$ and the 2D ghost RSA
packing, for which $g_{2}(r)$ are known analytically but $S(k)$ must be
numerically evaluated to high precision. The reasons for choosing the targets
in Table I are detailed in Sec. III.

The functional forms of our pair statistics and their corresponding potentials
can be used to model soft-core polymers with varying magnitudes of short- and
long-range repulsive forces [Eqs. (I)–(I), (I)–(I)], “sticky” macromolecules
with oscillatory interactions [Eqs. (I)–(I)], as well as polymer-grafted
nanoparticles with both hard- and soft-core interactions [Eqs. (I)–(I)].
Notably, we discover that one of our models possesses a stronger hyperuniform
property (i.e., larger exponent $\alpha$) than any other positive-temperature
equilibrium state known previously. To illustrate that our procedure can be
used to guide experimental polymer design, we also study a case of prescribed
pair statistics whose functional forms are motivated by the polymer reference
interaction site model (PRISM) [55, 56] and show that they can be realized. It
is important to note that one can realize many designed pair functions at
different densities and temperatures other than those explicitly shown in
Table I, either by creating linear combinations of our pair functions or by
applying our inverse procedure to other desirable functional forms. We further
show that our pair functions exhibit distinctly different translational order
metric $\tau$ (defined in Sec. V) [57] and self-diffusion coefficients
$\mathcal{D}$, and that $\tau$ and $\mathcal{D}$ are inversely related for our
3D models.

Section II provides a comprehensive description of our inverse methodology. In
Sec. III, we present our results for the target and optimized pair statistics
and the corresponding effective potentials. In Sec. IV, we probe the
realizability of the prescribed polymer pair statistics motivated by the PRISM
model. In Sec. V, we show that the translational order and the self-diffusion
coefficients are inversely correlated to one another for our 3D designs. In
Sec. VI, we discuss the implications of our work for outstanding problems in
statistical mechanics and for materials design.

##  II Inverse Methodology

To determine effective one- and two-body interactions from prescribed pair
statistics, we employ an accelerated version of an inverse algorithm that we
developed recently [33], which enables one to recover the unique pair
potential corresponding to given pair statistics, as dictated by Henderson’s
theorem [58]. Our current methodology improves upon the original one [33] by
incorporating automatic differentiation [59] and canonical ensemble reweighing
[60] techniques, achieving significantly faster convergence of potentials. A
flowchart of the the algorithm is shown in Fig. 1.

Figure 1: Flowchart of the inverse methodology developed in this work,
assuming homogeneous target system.

The methodology uses a parametrized potential $v(r;\mathbf{a})$, whose initial
functional form is informed by the small- and large-distance behaviors of the
target pair statistics [35]. The scalar components $a_{j}$ of the potential
parameter vector $\mathbf{a}$ are of three types: dimensionless energy scales
$\varepsilon_{j}$, dimensionless distance scales $\sigma_{j}$ and
dimensionless phases $\theta_{j}$ [33]; see Appendix B for the potential
parameters in each case. 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 for the target pair statistics; see
details in Ref. [33]. The large-$r$ behavior of $v(r;\mathbf{a})$ is given by
the target direct correlation function $c(r)$ under mild conditions [4, 33].
For hyperuniform targets, the large-$r$ asymptotic form of $v(r)\sim-c(r)$ is
given by Eq. (28) [35]. That is, the potential is long-ranged in the sense
that its volume integral is unbounded. Thus, one requires a neutralizing
background one-body potential to maintain stability [35, 61, 62, 1], and the
Ewald summation technique [63] is used to compute the total potential energy.
For nonhyperuniform targets, the volume integral of $c(r)$ is bounded and one
only requires two-body interactions.

Once the initial form of $v(r;\mathbf{a})$ is chosen, the low-storage BFGS
algorithm [64] 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}$ (26)

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})$ are
equilibrated pair statistics under $v(r;\mathbf{a})$ obtained from Monte-Carlo
simulations. As done in our previous works[33, 24, 65], we use Gaussian weight
functions $w_{g_{2}}(r)=\exp[-(r/4)^{2}]$ and $w_{S}(k)=\exp[-(k/2)^{2}]$.

The BFGS algorithm requires the gradient of the objective function
$\Psi(\mathbf{a})$. In this work, we replace the finite difference method in
Ref. [33] by automatic differentiation (AD), a technique that quickly
generates precise gradients on input parameters of computer programs by
repeatedly applying the chain rule [59]. AD has the advantage that its
computational complexity for gradients is independent of the number of
parameters [59]. Because standard implementations of AD [66] can only compute
gradient for deterministic functions, but $\Psi(\mathbf{a})$ depends on
stochastically simulated pair statistics, one requires a deterministic
function $\psi(\mathbf{a})$ that approximates $\Psi(\mathbf{a})$ in the
vicinity of given $\mathbf{a}$. In this work, we define $\psi(\mathbf{a})$ via
the canonical ensemble reweighing technique [60] described below.

Let $\mathcal{C}_{i}$ be a set of configurations equilibrated under
$v(r;\mathbf{a}_{i})$ drawn from a simulated canonical ensemble in iteration
$i$, and let $\mathbf{a}^{\prime}$ be a perturbed potential parameter in the
vicinity of $\mathbf{a}_{i}$. The ensemble-averaged pair statistics under
$v(r;\mathbf{a}^{\prime})$ are estimated by reweighing the pair statistics of
each configuration in $\mathcal{C}_{i}$ by a Boltzmann factor involving the
difference in total configurational energies
$\Phi(\mathbf{r}^{N})=\sum_{i<j}^{N}v(|\mathbf{r}_{i}-\mathbf{r}_{j}|)$ under
$v(r;\mathbf{a}_{i})$ and $v(r;\mathbf{a}^{\prime})$ [60], i.e.,

$\varphi(x;\mathbf{a}^{\prime})=\frac{\sum_{\mathbf{r}^{N}\in\mathcal{C}_{i}}\varphi(x;\mathbf{r}^{N})\exp[-(\Phi(\mathbf{r}^{N};\mathbf{a}^{\prime})-\Phi(\mathbf{r}^{N};\mathbf{a}_{i}))/T]}{\sum_{\mathbf{r}^{N}\in\mathcal{C}_{i}}\exp[-(\Phi(\mathbf{r}^{N};\mathbf{a}^{\prime})-\Phi(\mathbf{r}^{N};\mathbf{a}_{i}))/T]},$
(27)

where $\varphi(x;\mathbf{r}^{N})$ is $g_{2}(r)$ or $S(k)$ computed from each
configuration $\mathbf{r}^{N}$ in $\mathcal{C}_{i}$. The approximate objective
function $\psi(\mathbf{a}^{\prime})$ at the perturbed parameters
$\mathbf{a}^{\prime}$ can then be computed by inserting (27) into (26). For a
fixed set of simulated configurations $\mathcal{C}_{i}$,
$\psi(\mathbf{a}^{\prime})$ is a deterministic function of
$\mathbf{a}^{\prime}$. Thus, one can apply BFGS optimization with AD-generated
gradients to find the optimal parameters $\mathbf{a}_{i+1}$ that minimize
$\psi(\mathbf{a}^{\prime})$. In the next iteration, simulations are performed
under $v(r;\mathbf{a}_{i+1})$ to obtain a new set of configurations
$\mathcal{C}_{i+1}$. The iterations repeat until $\Psi(\mathbf{a}_{i})$
computed from $\mathcal{C}_{i}$ is smaller than some given convergence
threshold $\epsilon$, set to be $10^{-3}$ in this study. If convergence is not
achieved, a different set of basis functions is chosen and the optimization
process is repeated. Note that this inverse methodology accelerates our
original algorithm [33] as it requires just one simulation per iteration for
$n$ potential parameters, as opposed to $2n$ simulations using finite
differences. For systems away from the critical point, convergence can be
achieved within 5 hours on an Intel 2.8 GHz CPU with 15 cores and 4GB memory
per core.

For systems near the critical point, common inverse methods based on
predictor-corrector approaches, such as IBI [67], do not perform well due to
critical slowing down and the large number of iterations required for the
trial potential to converge [33]. IBI updates a binned trial potential $v(r)$
based on the difference between the target and simulated $g_{2}(r)$ in each
iteration. By contrast, our inverse method uses an analytical potential form
dictated by statistical-mechanical theory. This form is robust to variations
in the large-$r$ behavior of the simulated $g_{2}(r)$, which can cause
undesirable long-ranged oscillation in $v(r)$ in the IBI approach. Thus, the
potential in our method converges in much fewer iterations. For the 1D
antihyperuniform state, which is the most challenging target among the ones
studied here, convergence was achieved within 24 hours.

##  III Results for Designed Pair Correlation Functions

In this section, we describe the targeted analytical forms for pair statistics
in direct and Fourier spaces in various space dimensions summarized in Table I
and then show their realizability by equilibrium states with up to pair
interactions via the inverse algorithm described in Sec. II. For each target,
we show that the optimized pair statistics precisely match the prescribed
functions. The initial functional forms of the effective potentials are chosen
based on the small- and large-$r$ behaviors of the target pair statistics, as
dictated by statistical mechanical theory [33]. The asymptotic decay rate of
the direct correlation function $c(r)$, which is linked to the total
correlation function $h(r)=g_{2}(r)-1$ via the Ornstein-Zernike equation [68],
determines whether one needs both one-body and two-body interactions (for
hyperuniform states) or only two-body interactions (for nonhyperuniform
states); see Sec. II for details. Note that certain pair interactions for
hyperuniform states are not square integrable, but the total energy is
nevertheless bounded due to the background one-body interaction that maintains
overall charge neutrality [35, 33]. Unless otherwise stated, we take $\rho=1$,
and the lengths are expressed in units of $\rho^{-1/d}$. We also set the
reduced temperature $k_{B}T$ to be unity for all the equilibrium states, where
$k_{B}$ is the Boltzmann constant.

###  III.1 $d$-dimensional Hyperuniform Gaussian Pair Statistics

We consider $d$-dimensional generalizations of Gaussian pair statistics given
by Eqs. (I) and (I) in Table I. These targets are hyperuniform with $\alpha=2$
in Eq. (1). Note that the pair functions are self-similar in under Fourier
transformations up to scaling, and are independent of the dimensionality.
Gaussian pair statistics incorporate soft-core short-ranged repulsion found in
macromolecules, and thus are commonly used to model polymer systems [69, 70,
56, 71]. The Gaussian pair functions represent the low-density, high-
temperature limit of the pair statistics for the popular Gaussian-core model
[72, 73, 74, 12], which is known to exhibit anomalous phase behaviors such as
reentrent melting [75, 74].

(a)

(b)

(c)

(d)

Figure 2: (a)–(b) Target (solid curves) and optimized (dashed curves)
$g_{2}(r)$ and $S(k)$ for Gaussian pair statistics. For clarity, the pair
statistics for $d=2$ and $d=3$ are shifted up by 1 and 2, respectively. (c) A
400-particle configuration of the 2D equilibrium state with Gaussian pair
statistics. (d) A 125-particle configuration of the 3D equilibrium state with
Gaussian pair statistics. The point particles are shown as spheres for
visualization purposes.

Figure 2(a) and (b) show the target and optimized $g_{2}(r)$ and $S(k)$,
respectively, for the hyperuniform Gaussian pair statistics in one, two and
three dimensions. It is clear that the target and optimized pair statistics
agree very well with one another. The precise form of the corresponding
effective potential is given in the Appendix and plotted in Fig. 12(a).
Torquato [35] showed that for hyperuniform targets, the long-ranged asymptotic
form of $v(r)$ is given by

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

which can be regarded as a generalized Coulombic interaction of “like-charged”
particles. For $\alpha=2$, the long-ranged parts of the potentials are given
by the Coulomb interactions $v(r)\sim r^{-(d-2)}$ in the respective
dimensions. Note that for the 2D Gaussian pair statistics, it is known that
the effective pair potential is exactly given by the 2D Coulombic interaction
$v(r)=-2\ln(r)$ [20]. Our inverse procedure has accurately recovered this
potential.

Figure 2(c) and (d) show representative configurations of the 2D and 3D
equilibrium states with Gaussian pair statistics. Note that while particles
can get arbitrarily close, both configurations do not contain large clusters
or large spherical holes that are void of particles. Indeed, the particles
tend to form short chains that are well separated. These observations are
consistent with the soft repulsive interactions in the effective potentials.

###  III.2 $d$-dimensional Hyperuniform Generalized OCP and Their Fourier
Duals

We consider 3D generalizations [32] of the hyperuniform OCP given by Eqs.
(I)–(I) as well as the 1D and 3D Fourier duals of OCP given by Eqs. (I)–(I) in
Table I. A OCP is an equilibrium system of identical point particles of charge
$e$ interacting via the Coulomb potential and immersed in a rigid uniform
background of opposite charge to ensure overall charge neutrality [20]. OCPs
are classical analogs of the quantum jellium model and are used as as a first
approximation to dense plasmas [76]. The pair statistics for the 2D OCP, which
is equivalent to the Ginibre ensemble, are exactly given by the Gaussian
functions (I) and (I) [20].

Zhang and Torquato [32] proposed a generalization of the OCP pair correlation
function to arbitrary dimensions, given by

$g_{2}(r)=1-\exp[-v_{1}(r)],$ (29)

where $v_{1}(r)$ is the volume of a $d$-dimensional sphere of radius $r$. The
corresponding structure factors at $\rho=1$ are hyperuniform with $\alpha=2$.
For $d=1$, the generalized OCP pair statistics are identical to a Lorenztian
target

$g_{2}(r)=1-\exp(-2r),$ (30)

which has been shown to be realizable by a determinantal point process [19].
Thus, we study the 3D case here, whose realizability is unknown.

For the $d$-dimensional Fourier duals of OCP, the structure factors assume the
same form as the OCP pair correlation function (29) up to scaling, i.e.,

$S(k)=1-\exp[-v_{1}(k/2\pi)],$ (31)

for which the exponent in Eq. (1) is given by $\alpha=d$. This implies that
$S(k)$ is nonanalytic at the origin for any odd dimension, and thus the
corresponding $g_{2}(r)$ has power-law asymptotic $r^{-2d}$ [33].

We achieved excellent agreement between target and optimized pair statistics
for both OCP and its Fourier dual. Figure 3(a) and (b) show the target and
optimized pair statistics for the 3D OCP. The corresponding configuration
[Fig. 3(c)] is similar to that for the 3D Gaussian pair statistics, as both
systems have soft-core repulsive interactions. However, the OCP has the larger
force magnitudes as $r\rightarrow 0$, as shown in Fig. S1 in the Supplemental
Information (SI), which plots the forces $-dv(r)/dr$ corresponding to the
effective potentials for our 3D hyperuniform models. The difference in the
small-$r$ repulsion is due to the fact that the OCP $g_{2}(r)$ increases as
$r^{d}$ in the vicinity of the origin, rather than $r^{2}$ for the Gaussian.
Our capacity to design systems with varying magnitudes of short-range
repulsive forces enables coarse-grained modeling of polymers with different
levels of chain rigidity and penetrability [77, 78].

Figure 3(d)–(f) present the results for the Fourier dual of OCP in the first
three spatial dimensions. Note that the exponent $\alpha$ in (1) for the
simulated 3D Fourier dual of OCP, i.e., $\alpha=3$, is larger than that of any
other positive-temperature equilibrium states known previously. The largest
known exponent prior to this work was $\alpha=2$ for systems under Coulombic
interactions [35]. Our result here shows that states with $\alpha>2$ can be
attained by pair potentials that are longer-ranged than the $d$-dimensional
Coulombic interactions. While it is experimentally challenging to create
nanoparticles with such long-ranged molecular interactions, they are
hypothesized to be attainable by topological defects in the gauge fields in
the form of space-time instantons [79, 80]. Hyperuniform systems with large
$\alpha$ are desirable for optical purposes, because they tend to be
effectively transparent for a wide range of wave numbers compared to
nonhyperuniform ones [81]. We note further that despite having strong
large-$r$ repulsion, particles in the 3D Fourier dual of OCP are less
repulsive at small $r$ than those in the 3D OCP, as shown in Fig. S1.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: (a)–(b) Target and optimized $g_{2}(r)$ and $S(k)$ for 3D
generalized OCP. (c) A 125-particle configuration of the 3D equilibrium state
with the generalized OCP pair statistics. (d)–(e) Target (solid curves) and
optimized (dashed curves) for 1D and 3D Fourier dual of OCP. For clarity, the
pair statistics for the 3D case are shifted up by 1. (f) A 125-particle
configuration of the 3D equilibrium state with the Fourier dual OCP pair
statistics. In all 3D configurations, the point particles are shown as spheres
for visualization purposes.

###  III.3 $d$-dimensional Hyperuniform Hyperbolic Secant Function $g_{2}(r)$

We consider the hyperuniform states characterized by the hyperbolic secant
function $g_{2}(r)$ given by Eqs. (I)–(I) in Table I. As in the case of the
Gaussian pair statistics, these are soft-core models and the structure factors
scale as $k^{2}$ at small $k$. Moreover, in the 1D case, $g_{2}(r)$ and $S(k)$
are self-similar up to scaling. We have chosen this target for its capability
to model polymers with softer repulsive interactions than the Gaussian core
[82] and to generate configurations with a higher degree of clustering than
the Gaussian pair statistics. A hyperbolic secant probability distribution is
leptokurtic, i.e., it has smaller probability density around the mean but has
a heavier tail than the corresponding Gaussian distribution with the same mean
and variance [82]. For a hyperuniform state at unit number density, the
hyperbolic secant $g_{2}(r)$ rises more rapidly at small $r$ but decay to
unity more slowly at large $r$ compared to the Gaussian $g_{2}(r)$. Thus, we
propose that hyperbolic secant pair functions and other leptokurtic functions
[82] can be applied to describe polymer systems with “ultrasoft” particles,
such as star polymers [83, 84].

Figure 4 (a) and (b) show the target and optimized pair statistics for the
hyperbolic secant $g_{2}(r)$. Because these pair functions imply softer core
interactions than those for the Gaussian pair statistics, the corresponding 3D
configuration [Fig. 4(c)] contains more clusters with four or more closely
spaced particles with pair distances $r\leq 0.4\rho^{-1/3}$ compared to the
Gaussian case. Clusters in latter system are primarily dimer and trimer; see
Fig. 2(d). Cluster statistics can be quantified by mapping the point
configurations to a two-phase medium generated by decorating each point
particle by a sphere of a certain radius and then using theoretical tools from
percolation theory [85, 7, 86]. Figure S1 shows that the effective potential
for the hyperbolic secant $g_{2}(r)$ is the least repulsive at small $r$ among
the 3D hyperuniform models. As remarked in Sec. III.2, softer interactions are
useful in describing polymers molecules with high penetrability.

(a)

(b)

(c)

Figure 4: (a)–(b) Target (solid curves) and optimized (dashed curves) for
1D–3D hyperbolic secant $g_{2}(r)$. For clarity, the pair statistics for $d=2$
and $d=3$ are shifted up by 1 and 2, respectively. (c) A 125-particle
configuration of the 3D equilibrium state with hyperbolic secant $g_{2}(r)$.
The point particles are shown as spheres for visualization purposes.

###  III.4 2D Hyperuniform Gaussian-Damped Polynomial and 1D Nonhperuniform
Hermite-Gaussian Pair Statistics

We consider the 2D hyperuniform Gaussian-damped polynomial pair statistics
given by Eqs. (I) and (I) in Table I, as well as the 1D nonhyperuniform
Hermite-Gaussian pair statistics given by (I) and (I). These pair functions
mimic those found in dense polymer systems, which often exhibit oscillations
due to the coordination shells [87, 88]. In both cases, the pair functions are
self-similar under Fourier transformations up to scaling. Notably, the forms
of the 2D pair functions are chosen to closely resemble the pair statistics of
polymer chains modeled by freely rotating fused hard spheres [88]. The
oscillatory $g_{2}(r)$ function [Fig. 5(a)] dictates that the corresponding
many-particle system must have significant particle clustering. Also note that
the form of the 1D total correlation function $h(r)$ is a Hermite-Gaussian
function of order 4. Hermite-Gaussian functions have been used to describe
semiflexible polymers such as DNA and actin [89]. We can precisely realize the
pair functions in both cases, as shown in Fig. 5 and 6. The effective
potentials for both cases are oscillatory, and may be realized by “sticky”
macromolecules, such as DNA-coated nanoparticles, as one can achieve precise
control of attraction and repulsion at specific radial distances for these
particles [90, 91]; see Appendix for the precise functional forms. Compared to
the configuration for the 2D Gaussian pair function [Fig. 2(c)], the
configuration for the 2D Gaussian-damped polynomial [Fig. 5(c)] contains
considerably fewer chains. Instead, it contains many clusters with more than
four closely spaced particles due to the peak in $g_{2}(r)$ at
$r\rho^{1/2}\sim 0.6$. Despite this heterogeneity on the small to intermediate
scale, Fig. 5(c) is homogeneous on the large scale as a result of the
prescribed hyperuniformity, i.e., there are no large clusters or spherical
holes void of particles on the order of the system size.

(a)

(b)

(c)

Figure 5: (a)–(b) Target and optimized $g_{2}(r)$ and $S(k)$ for the 2D
Gaussian-damped polynomial pair statistics. (c) A 400-particle configuration
of the 2D equilibrium state with the Gaussian-damped polynomial pair
statistics.

(a)

(b)

Figure 6: Target and optimized (a) $g_{2}(r)$ and (b) $S(k)$ for the 1D
hyposurficial Hermite-Gaussian $g_{2}(r)$ with $\lambda=1/20$.

###  III.5 3D Hyposurficial State

We consider, for the first time, the possible realizability of the
nonhyperuniform hyposurficial state specified by Eqs. (I) and (I) in Table I,
which applies in three dimensions. Hyposurficial systems are a special class
of nonhyperuniform states that can be observed in nonequilibrium phase
transitions in amorphous ices [92], supercritical equilibrium liquids [33] and
certain hard-core systems [34]. They behave like ideal gases in their large-
scale local number variance $N(R)$, because they lack a “surface-area” term
proportional to $R^{d-1}$ in the growth of $N(R)$, i.e., $\sigma^{2}(R)\sim
AR^{d}+o(R^{d-1})$ as $R\rightarrow\infty$, where $A$ is a positive constant.
Hyposurficial pair statistics obey the following sum rule [34]:

$\int_{0}^{\infty}r^{d}h(r)dr=0,$ (32)

which implies that they must generally contain both negative and positive
correlations. Ideal gases are trivially hyposurficial. Figure 7 shows that we
can precisely realize the designed hyposurficial target. The configuration
[Fig. 7(c)] shows spatial heterogeneities, i.e., there exists significant
particle clustering as well as large spherical holes void of particles.

(a)

(b)

(c)

Figure 7: (a)–(b) Target and optimized $g_{2}(r)$ and $S(k)$ for the 3D
hyposurficial state. (c) A 125-particle configuration of the 3D equilibrium
hyposurficial state. The point particles are shown as spheres for
visualization purposes.

###  III.6 1D and 2D Ghost RSA

We consider the nonhyperuniform pair statistics corresponding to the 1D and 2D
ghost RSA packings, given by Eqs. (I)–(I) in Table I. The nonequilibrium ghost
RSA process [23, 27] is a special case of a generalization of the standard RSA
process [93, 94, 7] and is the only known model for which all $n$-particle
correlation functions are exactly solvable [23]. In the ghost RSA process,
spherical “test” particles of diameter $D$ are added continually to
$\mathbb{R}^{d}$ during time $t\geq 0$ according to a translationally
invariant Poisson process of density $\eta$ per unit time, i.e., $\eta$ is the
number of points per unit volume and time, here taken to be unity without loss
of generalization. A test sphere centered at position $\mathbf{r}$ at time $t$
is retained if and only if it does not overlap with any test sphere added in
the time interval $[0,t)$ [23]. The saturation packing fraction for the ghost
RSA process is $\phi(t=+\infty)=1/2^{d}$ [23].

We have previously shown that the 3D ghost RSA pair statistics at the terminal
packing fraction is realizable by an equilibrium state with up to pair
interactions [24]. Here, we investigate whether the same is true for $d=1$ and
2. It is not obvious that the ghost RSA processes are realizable in these
lower dimensions because, generally speaking, it is known that the lower the
space dimension, the more difficult it is to satisfy realizability conditions
[27]. This dimensionality effect arises on account of the decorrelation
principle [27], which states that spatial correlations that exist for a
particular model in lower dimensions diminish as the space dimension becomes
larger.

Figure 8 shows that we can achieve the nonequilibrium 1D and 2D ghost RSA
packings with the effective potentials in Fig. 12(i), which contain both a
hard-core and soft repulsion beyond the diameter. The effective soft repulsion
ensures that particles are well-separated, and corresponds to the exclusion of
real particles by the “ghost” test particles in the nonequilibrium process.
This combination of hard core and soft repulsion in $v(r)$ can be found in
colloidal suspensions of highly charged nanoparticles [95]. Note that the soft
repulsion is more significant for the 1D case than the 2D case, reflecting the
aforementioned decorrelation principle [27]. For $d=3$, the ghost RSA
$g_{2}(r)$ [23] is nearly equal to unity for $D\leq r\leq 2D$, and the soft
repulsion in its effective potential [24] is weaker than that in the 2D case.
In the limit that $d\rightarrow\infty$, the ghost RSA $g_{2}(r)$ approaches
the unit-step $g_{2}(r)$ [27], and the effective potential would consequently
converge to that of the hard-sphere potential. The effective potentials also
exhibit weak attractive wells [Fig. 8(a)], which offset the oscillations in
$g_{2}(r)$ that would occur in pure equilibrium hard-sphere system, ensuring
that $g_{2}(r)$ is flat beyond two diameters.

(a)

(b)

Figure 8: Target (solid curves) and optimized (dashed curves) (a) $g_{2}(r)$
and (b) $S(k)$ for the 1D and 2D ghost RSA packings. For clarity, the
$g_{2}(r)$ curve for $d=2$ is shifted up 2 and the $S(k)$ curve for $d=2$ is
shifted up by 1.

###  III.7 $d$-dimensional Antihyperuniform Pair Statistics

We study the antihyperuniform pair statistics (I) in one- two and three
dimensions. Note that the target $g_{2}(r)$ functions imply a hard core with
diameter $D$. The functional forms of $S(k)$ as well as the packing fraction
$\phi=\rho v_{1}(D/2)$ and the parameter $A$ are given in the Appendix A. The
values of $\phi$ and $A$ chosen ensure the necessary condition that $S(k)$
must be nonnegative at all $k$. As mentioned in the Introduction, $S(k)$ for
an antihyperuniform state has the small-$k$ scaling behavior $S(k)\sim
k^{\alpha}$, where $\alpha$ is negative. This divergence behavior of $S(k)$ at
small $k$ are characteristic of fluids at thermal critical points, which can
have fractal-like structures [38, 40]. For our antihyperuniform targets in any
dimension, we take $\alpha=-1/2$, implying that $h(r)$ has a long-ranged
power-law decay behavior at large $r$ as $h(r)\sim r^{-(d-1/2)}$ [35].
Realizable systems with such pair statistics would correspond to a different
universality class than the Ising class, which gives the large-$r$ behavior
$h(r)\sim r^{-(d-2+\eta)}$, where $\eta=1/4$ for $d=2$ and $\eta=0.036$ [42].

Computing equilibrium properties at critical points using Monte-Carlo or
molecular dynamics simulations is known to be challenging due to the critical
slowing down phenomena, i.e. the relaxation time to the equilibrium state
diverges at the critical point [38, 96]. Our inverse method reduces the time
required to reach the targeted critical-point pair statistics compared to
previous methods, as it uses parametrized analytical pair potentials rather
than binned potentials [33]; see Sec. II for details. We have realized the
antihyperuniform critical-point pair statistics in one, two and three
dimensions, as shown in Fig. 9. The 2D and 3D critical-point states (Fig. 9(c)
and (d)] clearly show large clusters and large empty spherical holes on the
order of the system size, confirming that antihyperuniformity corresponds to
divergent long-range number density fluctuations.

The realizability of the 1D critical state is unusual because it is commonly
thought that phase transitions cannot exist in 1D systems, but this only true
when the interactions are sufficiently short-ranged [97, 98]. Indeed, in the
Appendix, we show that the effective pair potential has a long-ranged
attractive part, i.e., $v(r)\sim-r^{-(d-1/2)}$ as $r\rightarrow\infty$. This
long-range behavior in the interaction explains why our 1D, 2D and 3D
critical-point states do not belong to the Ising university class that are
characterized by short-ranged interactions [39, 40, 41, 42].

(a)

(b)

(c)

(d)

Figure 9: (a)–(b) Target (solid curves) and optimized (dashed curves)
$g_{2}(r)$ and $S(k)$ for $d$-dimensional antihyperuniform pair statistics
(I). For clarity, $g_{2}(r)$ for $d=2$ and $d=3$ are shifted up by 1.5 and 3,
respectively, and $S(k)$ for $d=2$ and $d=3$ are shifted up by 1 and 2,
respectively. (c) A 400-particle configuration of the 3D equilibrium
antihyperuniform state. (d) A 125-particle configuration of the 3D equilibrium
antihyperuniform state. The disk and sphere diameters in the configurations
are the hard-core diameter $D$.

##  IV Realizability of Pair Statistics from a 3D Polymer Model

To illustrate that our procedure can be used to guide experimental polymer
design, we now examine a case of prescribed pair statistics whose functional
forms are motivated by the PRISM model [55, 56]. According to this model, the
total correlation function for polymer chains is given by [56]

$h(r;\rho)=\frac{3}{4}\sqrt{\frac{3}{\pi}}e^{-\frac{1}{4}\left(3r^{2}\right)}\tilde{\xi}^{\prime}_{\rho}\left(1-\frac{1}{2\tilde{\xi}_{\rho}^{2}}\right)-\frac{\tilde{\xi}^{\prime}_{\rho}\left(1-\frac{1}{2\tilde{\xi}_{\rho}^{2}}\right)^{2}e^{\frac{1}{3\tilde{\xi}_{\rho}^{2}}}\left[e^{\frac{r}{\tilde{\xi}_{\rho}}}\text{erfc}\left(\frac{1}{\sqrt{3}\tilde{\xi}_{\rho}}+\frac{\sqrt{3}r}{2}\right)-e^{-\frac{r}{\tilde{\xi}_{\rho}}}\text{erfc}\left(\frac{1}{\sqrt{3}\tilde{\xi}_{\rho}}-\frac{\sqrt{3}r}{2}\right)\right]}{2r},$
(33)

where lengths are measured in units of the radius of gyration $R_{G}$,
$\tilde{\xi}_{\rho}=1/[\sqrt{2}\left(\pi\sqrt{2}\rho+1\right)]$ and
$\tilde{\xi}^{\prime}_{\rho}=1/(2\pi\rho)$. The corresponding $S(k)$ is also
known analytically [56]. Yatsenko et al. [56] showed that for
$0.86<\rho<1.62$, Eq. (33) accurately describes the experimental radial
distributions of polyethylene and polypropylene melts at small to intermediate
$r$. However, experimental realizability of (33) outside of this density range
has not been tested. We have applied our inverse algorithm to Eq. (33) at
$\rho=0.5$ and have achieved excellent agreement between the target and
optimized pair functions by a soft-core effective potential, as shown in Fig.
10. These hypothetical pair statistics might be experimentally achievable by
polymer melts at low density, e.g. via polyethylene molecules with a smaller
number of monomers per chain than those considered in Ref. [56].

(a)

(b)

Figure 10: Target and optimized (a) $g_{2}(r)$ and (b) $S(k)$ for a
hypothetical 3D polymer motivated by the PRISM model with the total
correlation function (33) at $\rho=0.5$.

##  V Translational Order Metric and Self-Diffusion Coefficient

To demonstrate that our designed analytical pair functions enable one to
achieve a wide range of structural order and physical properties, we compute
the translational order metric $\tau$ [57] and the self-diffusion coefficient
$\mathcal{D}$ [12] for the 3D states studied here. The $\tau$ order metric is
defined by [57]

$\tau=\rho\int
h^{2}(r)d\mathbf{r}=\frac{\rho}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\tilde{h}^{2}(k)d\mathbf{k},$
(34)

where $\tilde{h}(k)$ is the Fourier transform of $h(r)$. This order metric
measures deviations of pair statistics from that of the uncorrelated (Poisson)
distribution $g_{2}(r)=1$. Since both positive and negative correlations
contribute to the integral, due to the fact that $h(r)$ is squared, $\tau$
measures the degree of translational order across all length scales. It
vanishes for the Poisson distribution and diverges for an infinite
crystal202020The metric $\tau$ can still be used to analyze crystals by
examining its growth with the system size for finite configurations [33]. and
is a positive bounded number for correlated disordered systems without long-
range order (i.e., Bragg peaks) [100].

To compute the self-diffusion coefficients, we utilize the semi-empirical
scaling laws developed by Rosenfeld [101]. The reduced self-diffusion
coefficient $\mathcal{D}\rho^{1/3}T^{1/2}$ for 3D fluids in the dilute limit
is given by [101]

$\mathcal{D}\rho^{1/3}T^{1/2}=0.37(-s_{2})^{-2/3},$ (35)

where

$s_{2}=-(\rho/2)\int_{\mathbb{R}^{d}}\\{g_{2}(r)\ln[g_{2}(r)]-g_{2}(r)+1\\}d\mathbf{r}$
(36)

is the two-body excess entropy. Equation (35) provides a good approximation
for soft-core models with $-s_{2}\leq 0.5$ [12], which is satisfied by all the
soft-core pair statistics in this study. For the ghost RSA and the
antihyperuniform packings, which have hard-core interactions, we use the
corresponding scaling law [101]

$\mathcal{D}\rho^{1/3}T^{1/2}=0.58\exp(0.78s_{2}).$ (37)

Figure 11 plots $\tau$ against the reduced self-diffusion coefficient, where
hyperuniform states are indicated with filled circles. It is evident that
$\tau$ and $\mathcal{D}$ are inversely correlated to one other. Pair functions
with smaller translational order metrics are closer to the Poisson pair
statistics. Thus, in equilibrium states with smaller $\tau$, particles
encounter less hindrance from neighboring particles, resulting in a higher
self-diffusivity. Importantly, $\tau$ has a strong positive correlation to
$-s_{2}$ [102], and in the limit that $g_{2}(r)$ is an infinitesimal
perturbation from unity (i.e., from an ideal gas), one has exactly
$\tau\sim-4s_{2}$ [102]. Consequently, $\tau$ is inversely related to
$\mathcal{D}$ via Eq. (35).

The ghost RSA and the antihyperuniform packings have hard-core interactions
and achieve higher $\tau$ and lower $\mathcal{D}$, because the hard particles
are more locally confined due to the frequent collisions with their neighbors
[103]. The soft-core hyperuniform states have intermediate $\tau$ and
$\mathcal{D}$. Figure S1 shows the forces corresponding to the effective
potentials for these states, from which it is clear that models with larger
magnitudes of small-$r$ repulsion give larger $\tau$ and smaller
$\mathcal{D}$. The OCP, which has the strongest small-$r$ repulsion among the
3D hyperuniform states, is least diffusive. This is followed by the Fourier
dual of OCP, the Gaussian pair statistics and the hyperbolic secant
$g_{2}(r)$. Interestingly, while OCP and its Fourier dual both have
$\tau=1/2$, the former has a slightly lower self-diffusivity because its
$g_{2}(r)$ scales as $r^{3}$ at small $r$, representing a stronger small-$r$
repulsion than the $r^{2}$ scaling for the Fourier dual of OCP. Our
observation that particles with soft repulsive interactions are more diffusive
than hard-core particles is consistent with the findings for macromolecules
via molecular dynamics simulations [78] and for lattice gases via Monte-Carlo
simulations [104].

Finally, the hyposurficial state, which lacks both hard-core interactions and
hyperuniformity, has the smallest $\tau$ and largest $\mathcal{D}$, as its
pair statistics are closest to those of an ideal gas. Its low translational
order is consistent with the observation in Ref. [105] that hyposurficiality
appears to be associated with spatial heterogeneities.

Figure 11: Translational order metric $\tau$ and the reduced self-diffusion
coefficient $\mathcal{D}\rho^{1/3}T^{1/2}$ for the 3D pair statistics studied
in this work and in Ref. [24]. The filled circles refer to hyperuniform states
and the unfilled ones refer to nonhyperuniform states.

##  VI Discussion

We have designed a diverse set of pair functions with varying degree of
disorder across length scales that are realizable by disordered many-particle
systems, many of which have novel structural and physical properties. We
achieved this via an efficient algorithm that determines equilibrium states
with up to pair interactions at positive temperatures. The realizability of
the prescribed pair statistics via equilibrium states with up to pair
interactions is rather strong, because it does not exclude the possibility
that the same pair statistics can be realized by equilibrium states involving
higher-body interactions or by nonequilibrium states.

We have chosen our pair functions based on their capability to generate a wide
spectrum of hyperuniform and nonhyperuniform systems with various short-range
correlations, including hyperuniform soft-core systems with different
small-$r$ repulsive forces, a nonhyperuniform 3D hyposurficial state, 1D and
2D ghost RSA packings and antihyperuniform packings. These forms of pair
functions dictate that the corresponding systems have different cluster sizes
from one another. Notably, the 3D Fourier dual of OCP possesses a larger
hyperuniformity exponent $\alpha$ than any other positive-temperature
equilibrium state known previously, which can be useful in optical
applications [81]. We have designed a diverse family of pair functions that
are self-similar under Fourier transformations, which means that their low-
density and high-density ground-state structures under the effective
potentials can be mapped to one another via duality relations [45, 46].
Furthermore, we have realized a case of pair statistics whose functional forms
are motivated by the PRISM model [56], indicating that it may be
experimentally achieved. Our designed pair functions exhibit distinctly
different translational order metric $\tau$ and self-diffusion coefficients
$\mathcal{D}$. We have shown that for our 3D models, $\tau$ and $\mathcal{D}$
are inversely correlated, the same trend of which should be true for the
corresponding one- and two-dimensional designs as well. This relationship
enables one to easily achieve desired diffusivity through designed pair
functions.

The pair potentials for many states can be experimentally achieved by polymer
or nanoparticle systems with tunable chemical components. For example,
Yatsenko et al. [56] showed that in the limit of long polymer chains, the
total distribution function of polyethylene melts with radius of gyration
$R_{G}$ can be approximated by a Gaussian-damped polynomial,

$h(\tilde{r},\xi)=-\frac{39}{16}\sqrt{\frac{3}{\pi}}\xi(1+\sqrt{2}\xi)(1-\frac{9\tilde{r}^{2}}{26})\exp(-\frac{3\tilde{r}^{2}}{4}),$
(38)

where $\tilde{r}=r/R_{G}$ and $\xi$ is a parameter inversely proportional to
the number density of polymer chains. At small to intermediate $\tilde{r}$,
Eq. (38) resembles our Gaussian pair statistics. Interestingly, it has been
suggested that hyperuniform pair statistics with soft-core interactions can be
achieved by star or bottlebrush polymers [106]. Note that the effective
potentials of many of our hyperuniform models101010These include Gaussian pair
statistics, OCP and its Fourier dual and the hyperbolic secant $g_{2}(r)$.
contain a Gaussian-damped logarithmic term for their short-ranged parts; see
the Appendix for the precise functional forms. Such logarithm-Gaussian pair
potential is commonly used to describe of intermolecular forces of star
polymers [108]. We also find that the oscillatory soft potentials in the cases
of polynomial-damped pair statistics can be achieved by “sticky”
macromolecules such as DNA-coated nanoparticles [90], and that the combination
of hard core and soft repulsion in the interactions for the ghost RSA packings
can be found in colloidal suspensions of highly charged nanoparticles [95].

The excellent agreement between the designed and simulated pair statistics in
all cases studied in this work demonstrates the power of our inverse procedure
to tackle realizability problems. It can be applied to study many prescribed
pair functions at different densities and temperatures other than those
explicitly shown in this work, provided that $g_{2}(r)$ and $S(k)$ at least
meet necessary nonnegativity conditions. For example, we have shown that
antihyperuniform critical-point systems can be designed that do not belong to
the Ising universality class. Remarkably, 1D critical-point states, which are
impossible in 1D models with Ising-like short-ranged interactions, can be
achieved as a result of the long-ranged pair potentials determined via our
inverse procedure. Knowledge of the critical point is essential in the
construction of phase diagrams fluids [40]. Previous studies on critical
points typically utilize the direct approach of statistical mechanics, i.e.,
simulations with known potentials. Our inverse approach to generate
antihyperuniform systems of different universality classes can significantly
accelerate the exploration of the phase diagram and the fabrication of systems
with fractal-like structures. One can also systematically study the
realizability of $g_{2}(r)$ with different asymptotic decay rates at large
$r$, including power-law, exponential and super-exponential rates. It is known
that the decay behavior of $g_{2}(r)$ for a fluid determines the Fisher-Widom
line that separates liquid-like and gas-like behaviors in the supercritical
region [38, 109]. Thus, $g_{2}(r)$ with prescribed decay rates may enable one
to create many-particle systems with desired phase behaviors.

Another important future project is to target disordered hyperuniform
structure factors with $\alpha\geq 3$ for their potential optical applications
[81]. Furthermore, we remark that the pair functions in Table I can be used as
basis functions to create new designs at unit density, i.e.,

$g^{*}_{2}(r)=\Sigma_{i=1}^{n}\lambda_{i}g_{2,i}(r),\quad 0\leq\lambda_{i}\leq
1,\Sigma_{i}^{n}\lambda_{i}=1.$ (39)

Equation (39) is expected to be realizable because it is realizable in all of
the extreme cases $\lambda_{i}=1$. Finally, note that if a given $g_{2}(r)$ is
realizable at $\rho=1$, then it is also realizable at all lower densities,
because the realizable density range for a given $g_{2}(r)$ is always a
continuous interval [110].

Because pair statistics based on our designed functions exhibit distinctly
different translational order and self-diffusivity, they can have potential
applications in the design of many-particle systems with tailored mass
transport properties. Importantly, while Fig. 11 shows discrete values of
$\tau$ and $\mathcal{D}$ for the 3D models in Table I, values between the
plotted points can be achieved by linear combinations of our pair functions in
the form of Eq. (39), and their corresponding potentials can be readily
determined via our inverse procedure. Thus, one can flexibly design pair
functions that yield a continuous spectrum of $\tau$ and $\mathcal{D}$ values.
In the contexts of diffusion-controlled reactions, the designed pair
statistics would guide the fabrication of nanoparticles, micelles or vesicles
loaded with reactants or catalysts to achieve desired reaction rates [111,
112]. Moreover, since our procedure can realize a wide range of short- and
long-range order, the equilibrium states can be modified into multi-phase
media with tailored long-time diffusion spreadability [113] for controlled
drug delivery purposes [114].

Our study also has important implications for the Zhang-Torquato conjecture,
which states that any realizable set of pair statistics, whether from a
nonequilibrium or equilibrium system, can be achieved by equilibrium systems
involving up to two-body interactions [32]. We have previously tested this
conjecture by accurately determining effective potentials for various
nonequilibrium pair statistics, including the unit-step function $g_{2}(r)$
from zero density to the dimension-dependent maximum realizable densities
across the first three space dimensions [24], 3D ghost RSA [24], “cloaked”
uniformly randomized lattices [33], perfect glasses and critical absorbing
states [65]. The realizability of all the prescribed pair statistics in this
study lends further support to this conjecture. In the future, it is important
to further test this conjecture with other nonequilbrium pair statistics,
including those of quantum states [45].

## Supplementary Material

The supplemental material provides a plot of the intermolecular force
$dv(r)/dr$ corresponding to the effective pair potentials for our 3D
hyperuniform models.

## Acknowledgements

This work is supported in part by the National Science Foundation CBET-2133179
and the U.S. Army Research Office 171 under Cooperative Agreement No.
W911NF-22-2-0103.

## Author Contributions

S.T. conceived of the project and research questions; H. W. and S.T. performed
the research and H.W. carried out all numerical calculations; and H. W. and
S.T. wrote the paper.

## Conflicting Interest

The authors declare no competing interests.

## Data Availability

Analytical pair functions that are known to be realizable, including the ones
in this and previous studies, are available in the PairFunctions repository,
reference number [115].

## Appendix

###  A Structure factors for antihyperuniform states

  * •

1D case:

$S(k)=1+\phi
A\sqrt{\frac{2\pi}{kD}}\left(-2\text{FresnelC}\left(\sqrt{\frac{2kD}{\pi}}\right)+1\right)-2\phi\frac{\sin(k)}{kD},$
(A1)

where we use $\phi=0.2,A=0.1$ in this work.

  * •

2D case

$S(k)=1-16A\phi\,_{1}F_{2}\left(\frac{1}{4};1,\frac{5}{4};-\frac{k^{2}}{4}\right)+\frac{4\sqrt{2}A\phi\Gamma\left(\frac{1}{4}\right)}{\sqrt{k}\Gamma\left(\frac{3}{4}\right)}-\frac{8\phi
J_{1}(k)}{k},$ (A2)

where we use $\phi=0.1,A=0.3$.

  * •

3D case:

$S(k)=1+\frac{24\sqrt{2\pi}A\phi-48\sqrt{2\pi}A\phi\text{FresnelC}\left(\sqrt{\frac{2kD}{\pi}}\right)}{\sqrt{kD}}+\frac{48A\phi\sin(kD)}{kD}-\frac{24\phi\sin(kD)}{(kD)^{3}}+\frac{24\phi\cos(kD)}{(kD)^{2}},$
(A3)

where we use $\phi=0.1,A=0.25$.

###  B Effective potentials

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 12: (a) Short-ranged part of the effective pair potentials (B1)–(B3)
for the 1D, 2D and 3D Gaussian pair statistics. (b) Short-ranged part of the
effective pair potential (B5) for the 3D OCP pair statistics. (c) Short-ranged
part of the effective potentials (B6)–(B7) for 1D and 3D Fourier dual of OCP
pair statistics. (d) Short-ranged part of the effective potential (B8) for the
1D hyperbolic secant $g_{2}(r)$. (e) Short-ranged part of the effective
potentials, Eqs. (B9) and (B10), for 2D and 3D hyperbolic secant $g_{2}(r)$,
respectively. (f) Short-ranged part of the effective potential (B11) for the
2D Gaussian-damped polynomial pair statistics. (g) Effective potentials (B12)
for the 1D hyposurficial Hermite-Gaussian $g_{2}(r)$ with $\lambda=1/20$. (h)
Effective potential (B13) for the 3D hyposurficial state. (i) Effective
potentials (B14)–(B15) for the 1D and 2D “ghost” RSA. (j) Effective potentials
(B16)–(B18) for 1D, 2D and 3D antihyperuniform states. (k) Effective
potentials for a 3D polymer model with the total correlation function (33) at
$\rho=0.5$.

  * •

1D Gaussian pair statistics:

$v(r)=-2\pi
r-\varepsilon_{1}^{(1)}\exp(-(r/\sigma_{1})^{2})(\ln(r)+\varepsilon_{1}^{(2)})+\sum_{j=2}^{4}\varepsilon_{j}\cos(r/\sigma_{j}^{(1)}+\theta_{j})\exp(-(r/\sigma_{j}^{(2)})^{2}),$
(B1)

where the optimized values of the parameters are given in Table 2.

Table 2: Parameters of the effective pair potential (B1) for 1D Gaussian pair statistics $\varepsilon_{1}^{(1)}$ | 7.368 | $\varepsilon_{2}$ | 20.15 | $\varepsilon_{3}$ | -0.2078 | $\varepsilon_{4}$ | 0.02070
---|---|---|---|---|---|---|---
$\sigma_{1}$ | 5.433 | $\sigma_{2}^{(1)}$ | 2.311 | $\sigma_{3}^{(1)}$ | 0.1564 | $\sigma_{3}^{(1)}$ | $1/(2\pi)$
$\varepsilon_{1}^{(2)}$ | 0.3824 | $\theta_{2}$ | 4.536 | $\theta_{3}$ | 2.911 | $\theta_{4}$ | 0
|  | $\sigma_{2}^{(2)}$ | 0.7983 | $\sigma_{3}^{(2)}$ | 3.305 | $\sigma_{4}^{(2)}$ | 86.07
  * •

2D Gaussian pair statistics:

$v(r)=-2\ln(r).$ (B2)

  * •

3D Gaussian pair statistics:

$v(r)=(1/r)(1-\exp(-(r/\sigma_{1})^{2}))-\varepsilon_{2}\exp(-(r/\sigma_{2})^{2})\ln(r),$
(B3)

where $\sigma_{1}=0.4672,\varepsilon_{2}^{(1)}=1.877,\sigma_{2}=0.6176$.

  * •

1D OCP pair statistics:

$v(r)=-2r-\varepsilon_{1}^{(1)}\exp(-(r/\sigma_{1})^{2})(\ln(r)+\varepsilon_{1}^{(2)})+\sum_{j=2}^{3}\varepsilon_{j}\cos(r/\sigma_{j}^{(1)}+\theta_{j})\exp(-(r/\sigma_{j}^{(2)})^{2}),$
(B4)

where the optimized values of the parameters are given in Table 3.

Table 3: Parameters of the effective pair potential (B4) for the 1D OCP pair statistics. $\varepsilon_{1}^{(1)}$ | 1.727 | $\varepsilon_{2}$ | 35.72 | $\varepsilon_{3}$ | -8.917
---|---|---|---|---|---
$\sigma_{1}$ | 14.22 | $\sigma_{2}^{(1)}$ | 6.501 | $\sigma_{3}^{(1)}$ | 1.915
$\varepsilon_{1}^{(2)}$ | 1.947 | $\theta_{2}$ | 4.438 | $\theta_{3}$ | 3.520
|  | $\sigma_{2}^{(2)}$ | 1.560 | $\sigma_{3}^{(2)}$ | 1.393
  * •

3D OCP pair statistics:

$v(r)=\frac{9}{2}\frac{3^{5/6}\Gamma\left(\frac{4}{3}\right)}{2^{2/3}\pi^{4/3}r}(1-\exp(-(r/\sigma_{2})^{2}))-\varepsilon_{2}\exp(-(r/\sigma_{1})^{2})\ln(r),$
(B5)

where $\sigma_{1}=0.4149$, $\varepsilon_{2}=3.008,\sigma_{2}=0.5624$.

  * •

Fourier dual of 1D OCP pair statistics:

$v(r)=-\ln(r)(1-\exp(-(r/\sigma_{1})^{2}))-\varepsilon_{2}\exp(-(r/\sigma_{2})^{2})\ln(r),$
(B6)

where $\sigma_{1}=0.2331,\varepsilon_{2}=1.375,\sigma_{2}=0.2331$.

  * •

Fourier dual of 3D OCP pair statistics:

$v(r)=-3\ln(r)+\ln(r/\sigma_{1}^{(1)})\exp(-(r/\sigma_{1}^{(2)})^{2}),$ (B7)

where $\sigma_{1}^{(1)}=0.4791$, $\sigma_{1}^{(2)}=0.4038$.

  * •

1D hyperbolic secant function $g_{2}(r)$:

$v(r)=-4r-\varepsilon_{1}^{(1)}\exp(-(r/\sigma_{1})^{2})(\ln(r)+\varepsilon_{1}^{(2)})+\sum_{j=2}^{3}\varepsilon_{j}\cos(r/\sigma_{j}^{(1)}+\theta_{j})\exp(-(r/\sigma_{j}^{(2)})^{2}),$
(B8)

where the optimized values of the parameters are given in Table 4.

Table 4: Parameters of the effective pair potential (B8) for the 1D hyperbolic secant function $g_{2}(r)$ $\varepsilon_{1}^{(1)}$ | 5.172 | $\varepsilon_{2}$ | 22.28 | $\varepsilon_{3}$ | -0.06283
---|---|---|---|---|---
$\sigma_{1}$ | 7.728 | $\sigma_{2}^{(1)}$ | 8.071 | $\sigma_{3}^{(1)}$ | 0.1610
$\varepsilon_{1}^{(2)}$ | 0.4883 | $\theta_{2}$ | 4.517 | $\theta_{3}$ | 3.241
|  | $\sigma_{2}^{(2)}$ | 0.6085 | $\sigma_{3}^{(2)}$ | 4.170
  * •

2D hyperbolic secant function $g_{2}(r)$:

$v(r)=-\frac{1}{G}\ln(r)(1-\exp(-(\frac{r}{\sigma_{1}})^{2}))-\varepsilon_{2}\exp(-(\frac{r}{\sigma_{2}})^{2})\ln(r),$
(B9)

where $\sigma_{1}=2.269,\varepsilon_{2}=1.367,\sigma_{2}=2.269$.

  * •

3D hyperbolic secant function $g_{2}(r)$:

$v(r)=\frac{3\sqrt[3]{\frac{2}{\pi}}}{5r}(1-\exp(-(r/\sigma_{1})^{2}))-\varepsilon_{2}\exp(-(r/\sigma_{2})^{2})\ln(r),$
(B10)

where $\sigma_{1}=0.1557,\varepsilon_{2}=1.922,\sigma_{2}=0.1174$.

  * •

2D Gaussian-damped polynomial pair statistics:

$v(r)=-\frac{6}{11}\ln(r)\left(1-e^{-\left(\frac{r}{\sigma_{1}}\right)^{2}}\right)-\varepsilon_{2}\ln(r)e^{-\left(\frac{r}{\sigma_{2}}\right)^{2}}+\varepsilon_{3}\cos\left(\frac{r}{\sigma_{3}^{(1)}}+\theta\right)e^{-\left(\frac{r}{\sigma_{3}^{(2)}}\right)^{2}},$
(B11)

where the optimized values of the parameters are given in Table 5.

Table 5: Parameters of the effective pair potential (B11) for the 2D Gaussian-damped polynomial pair statistics. $\sigma_{1}$ | 1.081 | $\sigma_{3}^{(1)}$ | 0.1843
---|---|---|---
$\varepsilon_{2}$ | 0.8532 | $\theta$ | 0.3920
$\sigma_{2}$ | 0.4884 | $\sigma_{3}^{(2)}$ | 0.9862
$\varepsilon_{3}$ | 0.7009 |  |
  * •

1D Hermite-Gaussian pair statistics with $\lambda=1/20$:

$v(r)=\sum_{j=1}^{2}\varepsilon_{j}\cos(r/\sigma_{j}^{(1)}+\theta_{j})\exp(-r/\sigma_{j}^{(2)})+\varepsilon_{3}\exp(-r/\sigma_{3}),$
(B12)

where the optimized values of the parameters are given in Table 6.

Table 6: Parameters of the effective pair potential (B12) for the 1D Hermite-Gaussian pair statistics with $\lambda=1/20$. $\varepsilon_{1}$ | 0.4143 | $\varepsilon_{2}$ | 3.174 | $\varepsilon_{3}$ | 0.2663
---|---|---|---|---|---
$\sigma_{1}^{(1)}$ | 0.4250 | $\sigma_{2}^{(1)}$ | 14.24 | $\sigma_{3}$ | 2.392
$\theta_{1}$ | 0.2047 | $\theta_{2}$ | 4.609 |  |
$\sigma_{1}^{(2)}$ | 2.829 | $\sigma_{2}^{(2)}$ | 1.544 |  |
  * •

3D hyposurficial state:

$v(r)=\sum_{j=1}^{2}\varepsilon_{j}\exp(-r/\sigma_{j})+\varepsilon_{3}\exp(-((r-\sigma_{3}^{(1)})/\sigma_{3}^{(2)})^{2}),$
(B13)

where the optimized values of the parameters are given in Table 7.

Table 7: Parameters of the effective pair potential (B13) for the 3D hyposurficial state. $\varepsilon_{1}$ | 1.646 | $\varepsilon_{3}$ | -0.1241
---|---|---|---
$\sigma_{1}$ | 0.4070 | $\sigma_{3}^{(1)}$ | 0.2020
$\varepsilon_{2}$ | 0.8742 | $\sigma_{3}^{(2)}$ | 1.312
$\sigma_{2}$ | 0.8884 |  |
  * •

1D ghost RSA:

$v(r)=\begin{cases}\infty,\quad r<0.5\\\
\varepsilon_{1}\exp(-(r/\sigma_{1}^{(1)})^{2})\cos(r/\sigma_{1}^{(2)}+\theta),\quad
r\geq 0.5\end{cases}$ (B14)

where
$\varepsilon_{1}=-2.341,\sigma_{1}^{(1)}=0.9045,\sigma_{1}^{(2)}=-0.4406,\theta=-1.686$.

  * •

2D ghost RSA:

$v(r)=\begin{cases}\infty,\quad r<1/\sqrt{\pi}\\\
\varepsilon_{1}\exp(-(r/\sigma_{1}^{(1)})^{2})\cos(r/\sigma_{1}^{(2)}+\theta),\quad
r\geq 1/\sqrt{\pi}\end{cases}$ (B15)

where
$\varepsilon_{1}=0.8107,\sigma_{1}^{(1)}=0.7860,\sigma_{1}^{(2)}=0.2867,\theta=-1.982$.

  * •

1D antihyperuniform state:

$v(r)=\begin{cases}\infty,\quad r<D\\\
-\varepsilon_{1}/\sqrt{r/D}+\varepsilon_{2}\exp(-(r/D-1)/\sigma_{2})+\varepsilon_{3}\exp(-((r/D-1)/\sigma_{3})^{2}),\quad
r>D,\end{cases}$ (B16)

where
$\varepsilon_{1}=0.3311,\varepsilon_{2}=0.0552,\sigma_{2}=70.21,\varepsilon_{3}=0.4299,\sigma_{3}=0.7560$.

  * •

2D antihyperuniform state:

$v(r)=\begin{cases}\infty,\quad r<D\\\
-\varepsilon_{1}/(r/D)^{3/2}+\varepsilon_{2}\exp(-(r/D-1)/\sigma_{2}^{(1)})\cos(r/(D\sigma_{2}^{(2)})+\theta),\quad
r>D,\end{cases}$ (B17)

where
$\varepsilon_{1}=0.1483,\varepsilon_{2}=1.647,\sigma_{2}=0.7530,\sigma_{2}^{(2)}=-4.120,\theta=4.987$.

  * •

3D antihyperuniform state:

$v(r)=\begin{cases}\infty,\quad r<D\\\
-\varepsilon_{1}/(r/D)^{5/2}+\varepsilon_{2}\exp(-(r/D-1)/\sigma_{2}^{(1)})\cos(r/(D\sigma_{2}^{(2)})+\theta),\quad
r>D,\end{cases}$ (B18)

where
$\varepsilon_{1}=0.2100,\varepsilon_{2}=0.6116,\sigma_{2}=0.7927,\sigma_{2}^{(2)}=-0.5982,\theta=7.203$.

  * •

3D hypothetical polymer motivated by the PRISM model:

$v(r)=\varepsilon_{1}\exp(-(r/\sigma_{1}^{(1)}))\cos(r/\sigma_{1}^{(2)}+\theta)+\varepsilon_{2}\exp(-r/\sigma_{2}),$
(B19)

where
$\varepsilon_{1}=3.635,\sigma_{1}^{(1)}=0.9891,\sigma_{1}^{(2)}=1.663,\theta=-0.4998,\varepsilon_{2}=-1.133,\sigma_{2}=0.3150$.

## References

  * Dyson [1962] F. J. Dyson, “Statistical theory of the energy levels of complex systems. I,” J. Math. Phys. 3, 140–156 (1962).
  * Feynman [1998] R. P. Feynman,  _Statistical Mechanics_ (Westview Press, Boulder, Colorado, 1998).
  * Van Workum and Douglas [2006] K. Van Workum and J. F. Douglas, “Symmetry, equivalence, and molecular self-assembly,” Phys. Rev. E 73, 031502 (2006).
  * Hansen and McDonald [1986] J. P. Hansen and I. R. McDonald,  _Theory of Simple Liquids_ (Academic Press, New York, 1986).
  * 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).
  * Hexner and Levine [2015] D. Hexner and D. Levine, “Hyperuniformity of critical absorbing states,” Phys. Rev. Lett. 114, 110602 (2015).
  * Torquato [2002] S. Torquato,  _Random Heterogeneous Materials: Microstructure and Macroscopic Properties_ (Springer-Verlag, New York, 2002).
  * Huang, Wang, and Holm [2021] J. P. Huang, Z. W. Wang, and C. Holm, “Swimming in circles can lead to exotic hyperuniform states of active living matter,” Proc. Natl. Acad. Sci. 118, e2100493118 (2021).
  * Nicholson  _et al._ [2021] D. M. Nicholson, C. Y. Gao, M. T. McDonnell, C. C. Sluss, and D. J. Keffer, “Entropy pair functional theory: Direct entropy evaluation spanning phase transitions,” Entropy 23 (2021), 10.3390/e23020234.
  * Bretonnet [2002] J.-L. Bretonnet, “Self-diffusion coefficient of dense fluids from the pair correlation function,” J. Chem. Phys. 117, 9370–9373 (2002).
  * Krekelberg  _et al._ [2009a] W. P. Krekelberg, M. J. Pond, G. Goel, V. K. Shen, J. R. Errington, and T. M. Truskett, “Generalized rosenfeld scalings for tracer diffusivities in not-so-simple fluids: Mixtures and soft particles,” Phys. Rev. E 80, 061205 (2009a).
  * Krekelberg  _et al._ [2009b] W. P. Krekelberg, T. Kumar, J. Mittal, J. R. Errington, and T. M. Truskett, “Anomalous structure and dynamics of the gaussian-core fluid,” Phys. Rev. E 79, 031203 (2009b).
  * Kirkwood, Buff, and Green [2004] J. G. Kirkwood, F. P. Buff, and M. S. Green, “The Statistical Mechanical Theory of Transport Processes. III. The Coefficients of Shear and Bulk Viscosity of Liquids,” J. Chem. Phys. 17, 988–994 (2004).
  * Gnan  _et al._ [2010] N. Gnan, C. Maggi, T. B. Schrøder, and J. C. Dyre, “Predicting the effective temperature of a glass,” Phys. Rev. Lett. 104, 125902 (2010).
  * Torquato [2021a] S. Torquato, “Diffusion spreadability as a probe of the microstructure of complex media across length scales,” Phys. Rev. E 104, 054102 (2021a).
  * Torquato, Scardicchio, and Zachary [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).
  * Soshnikov [2000] A. Soshnikov, “Determinantal random point fields,” Russian Mathematical Surveys 55, 923–975 (2000).
  * Hough  _et al._ [2009] J. B. Hough, M. Krishnapur, Y. Peres, and B. Virág,  _Zeros of Gaussian analytic functions and determinantal point processes_ , Vol. 51 (American Mathematical Society, Providence, RI, 2009).
  * 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).
  * Jancovici [1981] B. Jancovici, “Exact results for the two-dimensional one-component plasma,” Phys. Rev. Lett. 46, 386–388 (1981).
  * Note [15] These include $d$-dimensional ground-state free fermions [2, 16], eigenvalues for certain random matrices [1, 17, 18], 1D Lorentzian $g_{2}(r)$ [19] and 2D Ginibre ensemble [20].
  * Zernike and Prins [1927] F. Zernike and J. A. Prins, “Die Beugung von Röntgenstrahlen in Flüssigkeiten als Effekt der Molekülanordnung,” Z. Phys. 41, 184–194 (1927).
  * Torquato and Stillinger [2006a] S. Torquato and F. H. Stillinger, “Exactly solvable disordered sphere-packing model in arbitrary-dimensional Euclidean spaces,” Phys. Rev. E 73, 031106 (2006a).
  * Wang, Stillinger, and Torquato [2022] H. Wang, F. H. Stillinger, and S. Torquato, “Realizability of iso-g2 processes via effective pair interactions,” J. Chem. Phys. 157, 224106 (2022).
  * Axelrod  _et al._ [2022] S. Axelrod, D. Schwalbe-Koda, S. Mohapatra, J. Damewood, K. P. Greenman, and R. Gómez-Bombarelli, “Learning matter: Materials design with machine learning and atomistic simulations,” Acc. Mater. Res. 3, 343–357 (2022).
  * Torquato and Stillinger [2002] S. Torquato and F. H. Stillinger, “Controlling the short-range order and packing densities of many-particle systems,” J. Phys. Chem. B 106, 8354–8359 (2002), Erratum 106, 11406 (2002).
  * Torquato and Stillinger [2006b] S. Torquato and F. H. Stillinger, “New conjectural lower bounds on the optimal density of sphere packings,” Experimental Math. 15, 307–331 (2006b).
  * Yamada [1961] M. Yamada, “Geometrical study of the pair distribution function in the many-body problem,” Prog. Theor. Phys. 25, 579–594 (1961).
  * Lenard [1973] A. Lenard, “Correlation functions and the uniqueness of the state in classical statistical mechanics,” Commun. Math. Phys. 30, 35–44 (1973).
  * Kuna, Lebowitz, and Speer [2011] T. Kuna, J. L. Lebowitz, and E. R. Speer, “Necessary and sufficient conditions for realizability of point processes,” Annals Appl. Prob. 21, 1253–1281 (2011).
  * Lachieze-Rey and Molchanov [2015] R. Lachieze-Rey and I. Molchanov, “Regularity conditions in the realisability problem with applications to point processes and random closed sets,” Annals Appl. Prob. 25, 116–149 (2015).
  * 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).
  * 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 (2022).
  * Torquato and Stillinger [2003] S. Torquato and F. H. Stillinger, “Local density fluctuations, hyperuniform systems, and order metrics,” Phys. Rev. E 68, 041113 (2003).
  * Torquato [2018] S. Torquato, “Hyperuniform states of matter,” Physics Reports 745, 1–95 (2018).
  * 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).
  * Torquato [2021b] S. Torquato, “Structural characterization of many-particle systems on approach to hyperuniform states,” Phys. Rev. E 103, 052126 (2021b).
  * Widom [1965] B. Widom, “Equation of state in the neighborhood of the critical point,” J. Chem. Phys. 43, 3898–3905 (1965).
  * Kadanoff [1966] L. P. Kadanoff, “Scaling laws for Ising models near $T_{c}$,” Physics Physique Fizika 2, 263–272 (1966).
  * Fisher [1967] M. E. Fisher, “The theory of equilibrium critical phenomena,” Rep. Prog. Phys. 30, 615 (1967).
  * Wilson and Kogut [1974] K. G. Wilson and J. Kogut, “The renormalization group and the $\epsilon$ expansion,” Phys. Rep. 12, 75–199 (1974).
  * Binney  _et al._ [1992] J. J. Binney, N. J. Dowrick, A. J. Fisher, and M. E. J. Newman,  _The Theory of Critical Phenomena: An Introduction to the Renormalization Group_ (Oxford University Press, Oxford, England, 1992).
  * Mandelbrot [1982] B. B. Mandelbrot,  _The fractal geometry of nature_ (W. H. Freeman, New York, 1982).
  * Oğuz  _et al._ [2019] E. C. Oğuz, J. E. S. Socolar, P. J. Steinhardt, and S. Torquato, “Hyperuniformity and anti-hyperuniformity in one-dimensional substitution tilings,” Acta Cryst. Section A: Foundations & Advances A75, 3–13 (2019).
  * Torquato and Stillinger [2008] S. Torquato and F. H. Stillinger, “New duality relations for classical ground states,,” Phys. Rev. Lett. 100, 020602 (2008).
  * Torquato, Zachary, and Stillinger [2011] S. Torquato, C. E. Zachary, and F. H. Stillinger, “Duality relations for the classical ground states of soft-matter systems,” Soft Matter 7, 3780 (2011).
  * Torquato [2010] S. Torquato, “Reformulation of the covering and quantizer problems as ground states of interacting particles,” Phys. Rev. E 82, 056109 (2010).
  * [48] Here ${}_{p}F_{q}(a_{1},\dots,a_{p};b_{1},\dots,b_{q};z)$ is the generalized hypergeometric function and $\operatorname{ber}_{v}(z)$ is the Kelvin function.
  * [49] Here $G=0.915\dots$ is the Catalan’s constant.
  * [50] Here $\lambda$ is a parameter controlling the amplitude of oscillations in the pair statistics. The condition that $S(k)\geq 0$ for all $k$ implies that $\lambda\leq\exp(\sqrt{11/2}/2+7/4)/(8\left(\sqrt{22}+4\right)\sqrt{2\pi})=0.1067\dots$.
  * [51] Here $r^{*}=(4\pi)^{1/3}r$ and $k^{*}=(4\pi)^{-1/3}k$.
  * [52] Here Ci is the cosine integral and Si is the sine integral.
  * [53] Here $r^{*}=\sqrt{\pi}r$.
  * [54] Here $D$ is the hard-core diameter and $A+1$ is the peak height of $g_{2}(r)$ as $r\rightarrow D^{+}$.
  * Schweizer and Curro [1997] K. S. Schweizer and J. G. Curro, “Integral equation theories of the structure, thermodynamics, and phase transitions of polymer fluids,” in  _Advances in Chemical Physics_ (John Wiley & Sons, Ltd, 1997) pp. 1–142.
  * Yatsenko  _et al._ [2004] G. Yatsenko, E. J. Sambriski, M. A. Nemirovskaya, and M. Guenza, “Analytical soft-core potentials for macromolecular fluids and mixtures,” Phys. Rev. Lett. 93 (2004).
  * Torquato, Zhang, and Stillinger [2015] S. Torquato, G. Zhang, and F. H. Stillinger, “Ensemble theory for stealthy hyperuniform disordered ground states,” Phys. Rev. X 5, 021020 (2015).
  * Henderson [1974] R. L. Henderson, “A uniqueness theorem for fluid pair correlation functions,” Phys. Lett. A 49, 197–198 (1974).
  * Neidinger [2010] R. D. Neidinger, “Introduction to automatic differentiation and matlab object-oriented programming,” SIAM Rev. 52, 545–563 (2010).
  * Norgaard, Ferkinghoff-Borg, and Lindorff-Larsen [2008] A. B. Norgaard, J. Ferkinghoff-Borg, and K. Lindorff-Larsen, “Experimental parameterization of an energy function for the simulation of unfolded proteins,” Biophys. J. 94, 182–192 (2008).
  * 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 (1973).
  * Gann, Chakravarty, and Chester [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 (1979).
  * Ewald [1921] P. P. Ewald, “Die berechnung optischer und elektrostatischer gitterpotentiale,” Annalen der Physik 369, 253–287 (1921).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Math. Programming 45, 503–528 (1989).
  * Wang and Torquato [2023] H. Wang and S. Torquato, “Equilibrium states corresponding to targeted hyperuniform nonequilibrium pair statistics,” Soft Matter 19, 550–564 (2023).
  * Revels, Lubin, and Papamarkou [2016] J. Revels, M. Lubin, and T. Papamarkou, “Forward-mode automatic differentiation in Julia,” arXiv:1607.07892 [cs.MS] (2016).
  * Soper [1996] A. K. Soper, “Empirical potential Monte Carlo simulation of fluid structure,” Chem. Phys. 202, 295–306 (1996).
  * 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).
  * Flory and Bueche [1958] P. Flory and A. Bueche, “Theory of light scattering by polymer solutions,” J. Polym. Sci. 27, 219–229 (1958).
  * Koplik, Levine, and Zee [1983] J. Koplik, H. Levine, and A. Zee, “Viscosity renormalization in the Brinkman equation,” Phys. Fluids 26, 2864–2870 (1983).
  * Becker, Rosa, and Everaers [2010] N. Becker, A. Rosa, and R. Everaers, “The radial distribution function of worm-like chains,” The European Physical Journal E 32, 53–69 (2010).
  * Stillinger [1976] F. H. Stillinger, “Phase transitions in the Gaussian core system,” J. Chem. Phys. 65, 3968 (1976).
  * Louis, Bolhuis, and Hansen [2000] A. A. Louis, P. G. Bolhuis, and J. P. Hansen, “Mean-field fluid behavior of the gaussian core model,” Phys. Rev. E 62, 7961–7972 (2000).
  * Zachary, Stillinger, and Torquato [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).
  * Prestipino, Saija, and Giaquinta [2005] S. Prestipino, F. Saija, and P. V. Giaquinta, “Phase diagram of the gaussian-core model,” Phys. Rev. E 71, 050102 (2005).
  * Baus and Hansen [1980] M. Baus and J.-P. Hansen, “Statistical mechanics of simple Coulomb systems,” Phys. Rep. 59, 1–94 (1980).
  * Müller-Plathe [2002] F. Müller-Plathe, “Coarse-graining in polymer simulation: From the atomistic to the mesoscopic scale and back,” Chem. Phys. Chem. 3, 754–769 (2002).
  * Slyk, Skora, and Kondrat [2022] E. Slyk, T. Skora, and S. Kondrat, “How macromolecules softness affects diffusion under crowding,” Soft Matter 18, 5366–5370 (2022).
  * Kragset, Sudbø, and Nogueira [2004] S. Kragset, A. Sudbø, and F. S. Nogueira, “Metal-insulator transition in two- and three-dimensional logarithmic plasmas,” Phys. Rev. Lett. 92, 186403 (2004).
  * Børkje, Kragset, and Sudbø [2005] K. Børkje, S. Kragset, and A. Sudbø, “Instanton correlators and phase transitions in two- and three-dimensional logarithmic plasmas,” Phys. Rev. B 71, 085112 (2005).
  * Torquato and Kim [2021] S. Torquato and J. Kim, “Nonlocal effective electromagnetic wave characteristics of composite media: Beyond the quasistatic regime,” Phys. Rev. X 11, 021002 (2021).
  * Faliva and Zoia [2017] M. Faliva and M. G. Zoia, “A distribution family bridging the gaussian and the laplace laws, gram–charlier expansions, kurtosis behaviour, and entropy features,” Entropy 19 (2017), 10.3390/e19040149.
  * Likos  _et al._ [1998] C. N. Likos, H. Löwen, M. Watzlawek, B. Abbas, O. Jucknischke, J. Allgaier, and D. Richter, “Star polymers viewed as ultrasoft colloidal particles,” Phys. Rev. Lett. 80, 4450–4453 (1998).
  * Stellbrink  _et al._ [2000] J. Stellbrink, J. Allgaier, M. Monkenbusch, D. Richter, A. Lang, C. Likos, M. Watzlawek, H. Löwen, G. Ehlers, and P. Schleger, “Neither gaussian chains nor hard spheres-star polymers seen as ultrasoft colloids,” Trends in Colloid and Interface Science XIV , 88–92 (2000).
  * Torquato, Beasley, and Chiew [1988] S. Torquato, J. D. Beasley, and Y. C. Chiew, “Two-point cluster function for continuum percolation,” J. Chem. Phys. 88, 6540–6547 (1988).
  * Xu and Wunsch [2005] R. Xu and D. Wunsch, “Survey of clustering algorithms,” IEEE Trans. Neural Netw. 16, 645–678 (2005).
  * Donley, Rajasekaran, and Liu [1998] J. P. Donley, J. J. Rajasekaran, and A. J. Liu, “Density pair correlation functions for molecular liquids: Approximations for polymers,” J. Chem. Phys. 109, 10499–10512 (1998).
  * Yethiraj, Fynewever, and Shew [2001] A. Yethiraj, H. Fynewever, and C.-Y. Shew, “Density functional theory for pair correlation functions in polymeric liquids,” J. Chem. Phys. 114, 4323–4330 (2001).
  * Wilhelm and Frey [1996] J. Wilhelm and E. Frey, “Radial distribution function of semiflexible polymers,” Phys. Rev. Lett. 77, 2581–2584 (1996).
  * Mladek  _et al._ [2013] B. M. Mladek, J. Fornleitner, F. J. Martinez-Veracoechea, A. Dawid, and D. Frenkel, “Procedure to construct a multi-scale coarse-grained model of dna-coated colloids from experimental data,” Soft Matter 9, 7342–7355 (2013).
  * Kumar and Molinero [2017] A. Kumar and V. Molinero, “Self-assembly of mesophases from nanoparticles,” The Journal of Physical Chemistry Letters 8, 5053–5058 (2017).
  * Martelli  _et al._ [2017] F. Martelli, S. Torquato, N. Giovambattista, and R. Car, “Large-scale structure and hyperuniformity of amorphous ices,” Phys. Rev. Lett. 119, 136002 (2017).
  * Reńyi [1963] A. Reńyi, “On a one-dimensional problem concerning random space filling,” Sel. Trans. Math. Stat. Prob. 4, 203–218 (1963).
  * Widom [1966] B. Widom, “Random sequential addition of hard spheres to a volume,” J. Chem. Phys. 44, 3888–3894 (1966).
  * Liu and Luijten [2004] J. Liu and E. Luijten, “Stabilization of colloidal suspensions by means of highly charged nanoparticles,” Phys. Rev. Lett. 93, 247802 (2004).
  * Wissel [1984] C. Wissel, “A universal law of the characteristic return time near thresholds,” Oecologia 65, 101–107 (1984).
  * Dyson [1969] F. J. Dyson, “Existence of a phase-transition in a one-dimensional ising ferromagnet,” Comm. Math. Phys. 12, 91–107 (1969).
  * Cipra [1987] B. A. Cipra, “An introduction to the ising model,” Am. Math. Mon. 94, 937–959 (1987).
  * Note [20] The metric $\tau$ can still be used to analyze crystals by examining its growth with the system size for finite configurations [33].
  * Philcox and Torquato [2023] O. H. E. Philcox and S. Torquato, “Disordered heterogeneous universe: Galaxy distribution and clustering across length scales,” Phys. Rev. X 13, 011038 (2023).
  * Rosenfeld [1999] Y. Rosenfeld, “A quasi-universal scaling law for atomic transport in simple fluids,” J. Phys. Cond. Matter 11, 5415 (1999).
  * Lomba, Weis, and Torquato [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).
  * Dzugutov [2001] M. Dzugutov, “A universal scaling law for atomic diffusion in condensed matter,” Nature 381, 137–139 (2001).
  * Czarnecki  _et al._ [1996] G. Czarnecki, M. R. Dudek, A. Pekalski, and J. Cislo, “Tracer diffusion of particles with soft-core interactions studied by monte carlo simulations,” J. Phys. A 29, 3367 (1996).
  * Ma and Torquato [2017] Z. Ma and S. Torquato, “Random scalar fields and hyperuniformity,” J. Appl. Phys. 121, 244904 (2017).
  * Chremos and Douglas [2018] A. Chremos and J. F. Douglas, “Hidden hyperuniformity in soft polymeric materials,” Phys. Rev. Lett. 121, 258002 (2018).
  * Note [10] These include Gaussian pair statistics, OCP and its Fourier dual and the hyperbolic secant $g_{2}(r)$.
  * Archer, Likos, and Evans [2002] A. J. Archer, C. N. Likos, and R. Evans, “Binary star-polymer solutions: bulk and interfacial properties,” J. Phys. Cond. Matter 14, 12031 (2002).
  * Fisher and Wiodm [2003] M. E. Fisher and B. Wiodm, “Decay of Correlations in Linear Systems,” J. Chem. Phys. 50, 3756–3772 (2003).
  * Stillinger and Torquato [2005] F. H. Stillinger and S. Torquato, “Realizability issues for iso-$g^{(2)}$ processes,” Mol. Phys. 103, 2943–2949 (2005).
  * Wakao and Smith [1962] N. Wakao and J. Smith, “Diffusion in catalyst pellets,” Chem. Eng. Sci. 17, 825–834 (1962).
  * Olea and Thomas [1988] A. F. Olea and J. K. Thomas, “Rate constants for reactions in viscous media: correlation between the viscosity of the solvent and the rate constant of the diffusion-controlled reactions,” J. Am. Chem. Soc. 110, 4494–4502 (1988).
  * 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).
  * Langer and Peppas [1981] R. S. Langer and N. A. Peppas, “Present and future applications of biomaterials in controlled drug delivery systems,” Biomaterials 2, 201–214 (1981).
  * Wang [2023] H. Wang, “PairFunctions Github repository,” GitHub (2023).