An outstanding problem in statistical mechanics is the determination of whether prescribed functional forms of the pair correlation function *g*_{2}(*r*) [or equivalently, structure factor *S*(*k*)] at some number density *ρ* can be achieved by many-body systems in *d*-dimensional Euclidean space. The Zhang–Torquato conjecture 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. To further test this conjecture, we study the realizability problem of the nonequilibrium iso-*g*_{2} process, i.e., the determination of density-dependent effective potentials that yield equilibrium states in which *g*_{2} remains invariant for a positive range of densities. Using a precise inverse algorithm that determines effective potentials that match hypothesized functional forms of *g*_{2}(*r*) for all *r* and *S*(*k*) for all *k*, we show that the unit-step function *g*_{2}, which is the zero-density limit of the hard-sphere potential, is remarkably realizable up to the packing fraction *ϕ* = 0.49 for *d* = 1. For *d* = 2 and 3, it is realizable up to the maximum “terminal” packing fraction *ϕ*_{c} = 1/2^{d}, at which the systems are hyperuniform, implying that the explicitly known necessary conditions for realizability are sufficient up through *ϕ*_{c}. For *ϕ* near but below *ϕ*_{c}, the large-*r* behaviors of the effective potentials are given exactly by the functional forms exp[ − *κ*(*ϕ*)*r*] for *d* = 1, *r*^{−1/2} exp[ − *κ*(*ϕ*)*r*] for *d* = 2, and *r*^{−1} exp[ − *κ*(*ϕ*)*r*] (Yukawa form) for *d* = 3, where *κ*^{−1}(*ϕ*) is a screening length, and for *ϕ* = *ϕ*_{c}, the potentials at large *r* are given by the pure Coulomb forms in the respective dimensions as predicted by Torquato and Stillinger [Phys. Rev. E **68**, 041113 (2003)]. We also find that the effective potential for the pair statistics of the 3D “ghost” random sequential addition at the maximum packing fraction *ϕ*_{c} = 1/8 is much shorter ranged than that for the 3D unit-step function *g*_{2} at *ϕ*_{c}; thus, it does not constrain the realizability of the unit-step function *g*_{2}. Our inverse methodology yields effective potentials for realizable targets, and, as expected, it does not reach convergence for a target that is known to be non-realizable, despite the fact that it satisfies all known explicit necessary conditions. Our findings demonstrate that exploring the iso-*g*_{2} process via our inverse methodology is an effective and robust means to tackle the realizability problem and is expected to facilitate the design of novel nanoparticle systems with density-dependent effective potentials, including exotic hyperuniform states of matter.

## I. INTRODUCTION

The relationship between interactions in many-body systems and their corresponding structural properties is fundamental in statistical physics, condensed-matter physics, chemistry, mathematics, and materials science.^{1–4} Among the commonly used structural descriptors,^{5} the pair correlation function *g*_{2}(**r**) has unique importance due to its computational simplicity, experimental accessibility through diffraction measurements, and the popularity of theoretical statistical-mechanical models with pairwise additive interactions.^{5,6} Despite the widespread application of *g*_{2}(**r**), the realizability of pair correlation functional forms by actual many-body configurations remains an outstanding problem in statistical mechanics.^{7–13} It is known that for a statistically homogeneous (i.e., translationally invariant) many-body system in *d*-dimensional Euclidean space $Rd$, knowledge of one- and two-body correlations is insufficient to determine the corresponding higher-body correlation functions *g*_{3}, *g*_{4}, ….^{8–12} Consequently, given a prescribed functional form of a pair correlation function *g*_{2}(**r**) [or equivalently, a target structure factor *S*(**k**)] and number density *ρ* ≥ 0, an uncountable number of necessary and sufficient conditions must be satisfied for *g*_{2}(**r**) to be realizable,^{9} which are practically impossible to check. The ensemble-averaged structure factor *S*(**k**) is defined as

where *h*(**r**) = *g*_{2}(**r**) − 1 is the total correlation function and $h\u0303(k)$ is the Fourier transform of *h*(**r**). While some necessary conditions can be readily checked,^{7,9,10} it remains a practically important task to probe realizability via precise numerical techniques.

A major focus of this study is the realizability of the unit-step pair correlation function for positive densities,

where *D* is the sphere diameter, hence taken to be unity. Equation (2) implies that the many-body system forms a *packing* in which no particles overlap. The *packing fraction*, i.e., the fraction of the total volume covered by the spheres is given by *ϕ* = *ρv*_{1}(1/2), where *v*_{1}(*R*) is the volume of a *d*-dimensional sphere of radius *R*,

Importantly, Eq. (2) is exactly the zero-density limit of *g*_{2}(*r*) for the equilibrium hard-sphere (HS) fluid,^{5,6} whose pair potential is given by

Hence, it is not immediately obvious that the step-function *g*_{2} given by (2) can be achieved for positive densities.

The *g*_{2}-invariant process introduced by Torquato and Stillinger^{14} aims to determine the nonvanishing density range over which a prescribed form of *g*_{2}(*r*) for a many-body system remains invariant over that range. They showed that a *g*_{2}-invariant process possesses an upper *terminal density* *ρ*_{c}, which is the highest value such that all explicitly known necessary conditions are satisfied, including *g*_{2}(*r*) ≥ 0 for all *r*, *S*(*k*) ≥ 0 for all *k*, and the Yamada condition^{7,9,10} as detailed in Sec. II D. The terminal packing fraction in the case of the step function (2) is given by^{14}

at which the structure factor associated with (2) satisfies lim_{k→0}*S*(**k**) = 0, implying that the system must be hyperuniform^{15} if realizable.^{16} Disordered hyperuniform systems anomalously suppress large-scale density fluctuations compared to typical disordered systems, such as liquids;^{15} see Sec. II B. for mathematically precise definitions. Such exotic amorphous states of matter have been receiving great attention because they are endowed with unique physical or chemical properties and connect a variety of seemingly unrelated systems that arise in physics, chemistry, materials science, mathematics, and biology.^{17} Figure 1 shows the target *S*(*k*) corresponding to the unit-step function *g*_{2} at various densities for *d* = 1, 2, 3;^{16} see Sec. II A. for their explicit functional forms.

Previous numerical studies have shown that the unit-step function *g*_{2} is realizable but only for a finite range of *r*, in one,^{8,18} two,^{8} and three^{10} dimensions up to the terminal packing fraction. However, the reverse Monte-Carlo technique^{19} used to achieve such realizations does not target small-*k* correlations in Fourier space and provides no information on interparticle interactions. For *d* = 1, Costin and Lebowitz proved that the unit-step function *g*_{2} is realizable by a special type of Markovian point processes (known as renewal processes) if and only if *ϕ* ≤ *e*^{−1} ∼ 0.370.^{9} While the results for this special *g*_{2}-invariant process do not necessarily mean that the one-dimensional (1D) unit-step function *g*_{2} is non-realizable by *any* process for *e*^{−1} < *ϕ* < *ϕ*_{c} = 1/2, they nevertheless raise the possibility that this 1D target may not be realizable at *ϕ*_{c}.

A powerful way to tackle the realizability problem is via inverse statistical mechanics, i.e., determining the effective interactions that attain an equilibrium system with the prescribed pair statistics.^{3,13,20–24} Recently, Zhang and Torquato conjectured that any realizable *g*_{2}(**r**) or *S*(**k**) corresponding to a translationally invariant nonequilibrium system can be attained by an equilibrium ensemble involving only (up to) effective pair interactions.^{13} They introduced a theoretical formalism that enabled them to devise an efficient algorithm to construct systematically canonical-ensemble particle configurations with such targeted *S*(**k**) at all wavenumbers whenever realizable. However, the procedure presented in Ref. 13 does not provide the explicit functional forms of the underlying one- and two-body potentials.

Very recently, Torquato and Wang tested the Zhang–Torquato conjecture for challenging target pair statistics that correspond to nonequilibrium systems of current interest,^{25,26} including a nonhyperuniform 2D random sequential addition process^{27–30} as well as exotic disordered hyperuniform states, such as a 2D perfect glass,^{26,31} a 3D “cloaked” uniformly randomized lattice,^{25,32} and a 3D critical-absorbing state.^{26,33,34} In all cases, the nonequilibrium target pair statistics was achieved accurately by equilibrium states with effective one- and two-body potentials, lending crucial support to the Zhang–Torquato conjecture. Further testing this conjecture requires the solution of inverse problems for pair statistics with other prescribed functional forms.

In this work, we use a precise inverse methodology that we introduced recently^{25} to investigate numerically the realizability problem via effective interactions. Specifically, we study the problem as an *iso-**g*_{2} *process*,^{16,18,35} which is a process corresponding to a many-body system under a density-dependent effective potential *v*(*r*; *ρ*) that yields equilibrium states in which *g*_{2} remains invariant for a positive range of densities at constant positive temperature *T*. In other words, an effective potential must negate the natural density dependence of *g*_{2}(**r**). Density-dependent potentials are widely applied to describe coarse-grained models for polymers and macromolecules.^{36} The iso-*g*_{2} process is a nonequilibrium process whenever the target *g*_{2} is known to correspond to an equilibrium state: One then studies whether the same equilibrium *g*_{2} can be achieved for a range of nonvanishing densities. Since at positive densities, *g*_{2}(*r*) for the pure equilibrium hard-sphere fluid deviates from the unit-step functional form (2) and becomes oscillatory,^{6} the corresponding iso-*g*_{2} process determines an effective potential *v*(*r*; *ϕ*) that suppresses such oscillations (both positive and negative correlations) and maintains the unit-step form (2). Approximations of the effective potential have been obtained via Percus–Yevick and hypernetted chain (HNC) approximations.^{18} Stillinger and Torquato also derived low-density approximations for the effective potential using diagrammatic expansions.^{16} However, these approximate potentials yield *g*_{2}(*r*)’s that deviate from (2), and the deviations increase for larger *ϕ*. Thus, these results do not provide a definitive conclusion on the realizability of the unit-step function *g*_{2} for *ϕ* close to *ϕ*_{c}.

Prior to the development of our inverse methodology,^{25} predictor–corrector methods,^{20–22,24} such as iterative Boltzmann inversion (IBI)^{22} and iterative HNC inversion (IHNCI),^{24,37} 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;^{25} see Sec. II B. 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. 25.

Our inverse methodology^{25} 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.^{25} Subsequently, a nonlinear optimization technique is utilized to minimize an objective function that incorporates *both* the target pair correlation function *g*_{2}(**r**) and structure factor *S*(**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^{38} as well as the neutralizing background one-body potential;^{25} see Sec. III for details.

To assess the accuracy of inverse methodologies to target pair statistics, we introduced^{25} the following dimensionless *L*_{2}-norm error:

where $Dg2$ and *D*_{S} are *L*_{2} functions given by

where *g*_{2,F}(**r**; **a**) and *S*_{F}(**k**; **a**) represent the final pair statistics at the end of the optimization, which depend on the supervector **a**. We have previously shown that our inverse methodology generally yields *L*_{2}-norm errors that are an order of magnitude smaller than those obtained via previous inverse methods and is able to treat challenging near-critical and hyperuniform targets,^{25} which previous methods^{20–22,24} cannot do. Importantly, for equilibrium target pair statistics, it reaches the precision required to recover the unique potential dictated by Henderson’s theorem.^{39} Thus, it is a superior method for our purpose of probing the realizability of pair statistics over all distances in direct (physical) space and all wavenumbers in Fourier space, especially as it concerns hyperuniform targets.

A major goal of this study is to ascertain the density range on which the unit-step function *g*_{2} is realizable for *d* = 1, 2, 3 in the thermodynamic limit and to determine the effective potential when it is realizable. We find that for *d* = 1, the unit-step function *g*_{2} is numerically realizable up to *ϕ* = 0.49 = 0.98*ϕ*_{c}. For *d* = 2, 3, it is realizable up to *ϕ* = *ϕ*_{c}, at which the systems are hyperuniform. The effective potentials at the maximum realizable packing fractions of the unit-step function *g*_{2} significantly suppress both positive and negative correlations that would otherwise be present in the pure equilibrium HS system, as manifested by the smaller order metrics *τ* (19) for the corresponding former systems. For *ϕ* near but below *ϕ*_{c}, the large-*r* behavior of the effective potentials are exp[ − *κ*(*ϕ*)*r*], *r*^{−1/2} exp[ − *κ*(*ϕ*)*r*], and *r*^{−1} exp[ − *κ*(*ϕ*)*r*] for *d* = 1, 2, and 3, respectively, where *κ*(*ϕ*) is the inverse screening length, and for the 2D and 3D targets at *ϕ* = *ϕ*_{c}, the potentials at large *r* are given by the pure Coulomb forms in the respective dimensions as predicted by Torquato and Stillinger.^{15} This capacity to generate hyperuniform states via interacting many-particle systems is an important advance, since it is expected to facilitate the self-assembly of tunable hyperuniform soft-matter materials and enables one to probe the thermodynamic and dynamic properties of such exotic states.^{26}

We also study the realizability of two other targets of particular theoretical importance. The ghost random sequential addition (RSA) process^{29} is a generalization of the standard RSA process^{5,40,41} and is closely related to the unit-step iso-*g*_{2} process as its maximum packing fractions for all *d* are identical to *ϕ*_{c} = 1/2^{d}. At the maximum density, *g*_{2} for the ghost RSA process contains a small positive “bump” at the contact radius *r* = 1, which may suggest that low-dimensional unit-step function *g*_{2} is unrealizable.^{29} We show that the 3D ghost RSA target at *ϕ*_{c} = 1/8 can be realized by an effective potential that is much shorter ranged than that for the unit-step function *g*_{2}. We also test the accuracy and power of our inverse methodology for a known non-realizable case that meets all of the known necessary conditions and yet is not realizable due to geometric constraints.^{11} Our method is shown to be robust as it does not reach convergence, as expected, for this known unrealizable target.

We begin by providing basic definitions and background in Sec. II. In Sec. III, we briefly review our inverse methodology to determine effective potentials for targeted pair statistics. Section IV presents the realizability range and the effective interactions for the unit-step function *g*_{2} for *d* = 1, 2, 3. Sections V and VI. present the results for the 3D ghost RSA target and the 2D non-realizable target, respectively. We provide concluding remarks in Sec. VII.

## II. DEFINITIONS AND PRELIMINARIES

### A. Pair statistics

We consider many-particle systems in $Rd$ that are completely statistically characterized by the *n*-particle probability density functions *ρ*_{n}(**r**_{1}, …, **r**_{n}) for all *n* ≥ 1.^{6} In the case of statistically homogeneous systems, *ρ*_{1}(**r**_{1}) = *ρ* and *ρ*_{2}(**r**_{1}, **r**_{2}) = *ρ*^{2}*g*_{2}(**r**), *ρ* is the number density in the thermodynamic limit, *g*_{2}(**r**) is the pair correlation function, and **r** = **r**_{2} − **r**_{1}. If the system is also statistically isotropic, then *g*_{2}(**r**) is the radial function *g*_{2}(*r*), where *r* = |**r**|.

For a single periodic configuration containing *N* point particles at positions **r**_{1}, **r**_{2}, …, **r**_{N} within a fundamental cell *F* of a lattice Λ, the *scattering intensity* $I(k)$ is defined as

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*(**k**) (1) by

where *V*_{F} is the volume of the fundamental cell and *δ* is the Dirac delta function.^{17} In simulations of many-body systems with finite *N* under periodic boundary conditions, Eq. (9) is used to compute *S*(**k**) directly by averaging over configurations. The structure factor corresponding to the unit-step function *g*_{2} [Eq. (2)] in dimensions 1, 2, 3 are given by^{16}

where *J*_{1}(*x*) is the first-order Bessel function. Note that at *ϕ*_{c} = 1/2^{d}, *S*(*k*) ∼ *k*^{2} as *k* → 0, indicating that terminal-density systems are hyperuniform if they are realizable.

### B. Hyperuniformity and nonhyperuniformity

A *hyperuniform* point configuration in *d*-dimensional Euclidean space $Rd$ is characterized by an anomalous suppression of large-scale density fluctuations relative to those in typical disordered systems, such as liquids and structural glasses.^{15,17} More precisely, a hyperuniform point pattern is one in which the structure factor $S(k)=1+\rho h\u0303(k)$ tends to zero as the wave number *k* = |**k**| tends to zero,^{15,17} i.e.,

This hyperuniformity condition implies the following direct-space sum rule:

An equivalent definition of hyperuniformity is based on the local number variance *σ*^{2}(*R*) = ⟨*N*(*R*)^{2}⟩ − ⟨*N*(*R*)⟩^{2} associated with the number *N*(*R*) of points within a *d*-dimensional spherical observation window of radius *R*, where angular brackets denote an ensemble average. A point pattern in $Rd$ is hyperuniform if its variance grows in the large-*R* limit slower than *R*^{d}.^{15}

Consider systems that are characterized by a structure factor with a radial power-law form in the vicinity of the origin, i.e.,

For hyperuniform systems, the exponent *α* is positive (*α* > 0) and its value determines three different large-*R* scaling behaviors of the number variance:^{15,17,42}

Torquato showed that for an equilibrium hyperuniform state under a pair potential whose structure factor is characterized by the power law (14), the potential *v*(*r*) must be long-ranged in the sense that its volume integral is unbounded,^{17}

Therefore, to stabilize a classical hyperuniform system at positive *T*, which is thermodynamically incompressible,^{15} one must treat it as a system of “like-charged” particles immersed in a rigid “background” of equal and opposite “charge,” i.e., the system must have overall charge neutrality.^{17} Such a rigid background contribution corresponds to a one-body potential, which is described in Sec. III.

By contrast, for any *nonhyperuniform* system, the local variance has the following large-*R* scaling behaviors:^{43}

For a “typical” nonhyperuniform system, *S*(0) is bounded. In antihyperuniform systems,^{17} *S*(0) is unbounded, i.e.,

### C. Order metric

Scalar order/disorder metrics have been profitably used to quantify the degree of order in many-particle systems, including sphere packings.^{5,50} A particularly useful measure for our purposes is the order metric *τ*,^{51} which is defined as

where *D* is a characteristic “microscopic” length scale, taken to be the sphere diameter in this work. This order metric measures deviations of two-particle statistics from that of the Poisson distribution. Since both positive and negative correlations contribute to the integral, due to the fact that *h*(*r*) or $h\u0303(k)$ is squared, *τ* measures the degree of translational order across length scales. It clearly vanishes for the uncorrelated Poisson distribution and diverges for an infinite crystal and is a positive bounded number for correlated disordered systems without long-range order (i.e., Bragg peaks).

### D. Necessary conditions for realizability

Two well-known necessary conditions must be satisfied by any pair statistics,^{8,10,16,52} viz., the nonnegativity of *g*_{2}(*r*) and its corresponding *S*(*k*),

A further necessary condition forces a lower bound on the variance associated with the number of particles within a *d*-dimensional spherical window of radius *R*,^{7}

where *θ* is the fractional (non-integer) part of the average number of particle centers contained within the window ⟨*N*(*R*)⟩ = *ρv*_{1}(*R*). The Yamada condition, in practice, is relevant in very low dimensions when (20) and (21) are satisfied and becomes less restrictive for higher dimensions.^{11} It turns out to pose no additional restriction on the realizability of all the target pair statistics considered in this study in any dimension.^{11}

## III. INVERSE METHODOLOGY

Here, we describe our methodology to determine (up to) pair interactions that yield canonical ensembles matching target pair statistics; see Ref. 25 for more details. This procedure enables one to determine equilibrium states that match both target *g*_{2}(*r*) and *S*(*k*) with unprecedented accuracy when they are realizable. Thus, it is especially suitable for probing the realizability of hypothesized functional form for the pair statistics, which heretofore has not been done. The methodology uses a parameterized family of pointwise basis functions for the dimensionless pair potential,

where *f*_{j}(*r*/*σ*; *a*_{j}) is the *j*th basis function, *a*_{j} is a vector of parameters (generally consisting of multiple parameters), **a** = (*a*_{1}, *a*_{2}, …, *a*_{n}) is the “supervector” parameter, *ɛ* sets the energy scale, and *σ* is a characteristic length scale. Both *ϵ* and *σ* are taken to be unity. The components of *a*_{j} are of four types: dimensionless energy scales *ɛ*_{j}, dimensionless distance scales *σ*_{j}, dimensionless phases *θ*_{j}, and dimensionless exponents *p*_{j}.

The initial form of *v*(*r*; **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.^{17} Specifically, under mild conditions, the exact large-*r* behavior of *v*(*r*; **a**) is enforced in (23) to be the large-*r* behavior of the targeted direct correlation function *c*_{T}(*r*) via^{6,53}

where *β* = 1/(*k*_{B}*T*) and *k*_{B} is the Boltzmann constant, and the large-|**r**| behavior of *c*(**r**) can be extracted from the structure factor *S*_{T}(**k**) using the Fourier representation of the Ornstein–Zernike integral equation,^{54}

To estimate the initial small- and intermediate-*r* behavior in (23), we simply use the hypernetted chain (HNC) approximation,^{6} i.e.,

The basis functions are chosen so that they reasonably span all potential functions that could correspond to a targeted *g*_{2,T}(*r*) for all *r* or a targeted *S*_{T}(*k*) for all *k* under the constraint that the resulting potential function *v*(*r*) satisfies the small-*r* and large-*r* behaviors dictated by statistical-mechanical theory described above. For our specific targets, *f*_{j}(*r*; *a*_{j}) are chosen from the possible following general forms:

- Hard core:This basis function is used for the unit-step function(27)$fj(r;aj)=+\u221e,r\u22641,0,r>1.$
*g*_{2}and the RSA target as its*g*_{2,T}(*r*) exhibits a hard core for*r*≤ 1. - Exponential-damped oscillatory form:Importantly, the large-(28)$fj(r;aj)=\epsilon j\u2061cosr\sigma j(1)+\theta jexp\u2212r\u2212\sigma j(2)\sigma j(3)M.$
*r*behavior of*v*(*r*) ∼ −*βc*(*r*) for the 1D unit-step function*g*_{2}for $\varphi \u2192\varphi c\u2212$ (but*ϕ*≠*ϕ*_{c}) is given by an exponential form;^{15}see Sec. IV. In this work, for simplicity and efficiency, the exponent*M*(28) is restricted to be an integer that remains fixed during the optimization process as described below. - Exponential-damped power-law form:
^{55}where(29)$fj(r;aj)=\epsilon j\u2061exp\u2212\kappa rr\u2212pj,$*κ*is the inverse screening length. This form is chosen because it has been shown that the large-*r*behavior of*v*(*r*) ∼ −*βc*(*r*) for the 2D and 3D unit-step function*g*_{2}in the limit for $\varphi \u2192\varphi c\u2212$ (but*ϕ*≠*ϕ*_{c}) is given by (29),^{15}with*p*_{j}= 1/2 for*d*= 2 and*p*_{j}= 1 for*d*= 3; see Sec. IV for further details. - Power-law damped oscillatory form:(30)$fj(r;aj)=\epsilon j\u2061cosr\sigma j(1)+\theta jr\u2212pj.$
The last form is the Coulomb form (16), which corresponds to the large-

*r*behavior of*v*(*r*) for the unit-step function*g*_{2}at the terminal packing fraction*ϕ*_{c}in all dimensions, given that the states are realizable; see Sec. II B.

Once the initial form of *v*(*r*; **a**) is chosen, a nonlinear optimization procedure^{56} is used to minimize an objective function Ψ(**a**) based on the distance between target and trial pair statistics in both direct and Fourier spaces,

where $wg2(r)$ and *w*_{S}(**k**) are weight functions and *g*_{2}(**r**; **a**) and *S*(**k**; **a**) correspond to an equilibrated *N*-particle system under *v*(*r*; **a**) at a dimensionless temperature *k*_{B}*T*/*ɛ* = 1, which can be obtained from Monte-Carlo (MC) (used here) or molecular dynamics simulations under periodic boundary conditions. The optimization procedure ends when Ψ(**a**) is smaller than some small tolerance *ϵ*. If this convergence criterion is not achieved, then a different set of basis functions is chosen and the optimization process is repeated. The additional basis functions are obtained via a fit of the difference of the potentials of mean forces between the simulated and target *g*_{2}(*r*).

For hyperuniform targets whose small-*k* behavior is given by a power law *S*(*k*) ∼ *k*^{α}, *v*(*r*; **a**) has the long-ranged asymptotic form given by *v*(*r*) ∼ 1/*r*^{d−α}, which can be regarded as a generalized Coulombic interaction of “like-charged” particles.^{17} Thus, one requires a neutralizing background one-body potential to maintain stability.^{17,57–59} 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.^{60}

In our implementation of the inverse procedure, we used $wg2(r)=exp(\u2212r2/16)$, *w*_{S}(*k*) = exp(−*k*^{2}/4), *ϵ* = 5 × 10^{−4}, and the system sizes in MC simulations were *N* = 400, 500, 1000 for *d* = 1, 2, 3, respectively. However, if a hyperuniform target is found to be numerically realizable for the aforementioned small system size, we proceed to verify its realizability for larger system sizes, i.e., *N* = 10, 000 for *d* = 1, 2 and *N* = 20, 000 for *d* = 3. If convergence is not achieved within five iterations of reselection of basis functions, the target is deemed to be potentially not realizable.

## IV. ISO-*g*_{2} PROCESS FOR THE UNIT-STEP FUNCTION *g*_{2}

In this section, we present the realizable density ranges for the optimized effective interactions for the unit-step function *g*_{2} for *d* = 1, 2, 3. Taylor expansions of Eq. (11) at the point *k* = 0 reveal that *S*(*k*) for the unit-step function *g*_{2} are analytic at the origin, implying that the power series only admits even powers of *k*.^{17} Thus, the effective potential must have exponential or superexponential decay at *ϕ* < *ϕ*_{c} and a Coulomb form (16) at *ϕ* = *ϕ*_{c}.^{17,25} Indeed, Torquato and Stillinger showed that at *ϕ* = *ϕ*_{c}, the target direct correlation function for *r* ≫ 1 is given by^{15,61}

Furthermore, in the limit $\varphi \u2192\varphi c\u2212$ (but *ϕ* ≠ *ϕ*_{c}), the form of *c*_{T}(*r*) at large *r* is given by an exponential form for *d* = 1, and by exponential-damped power-law forms for *d* ≥ 2, i.e.,^{15}

where *κ*(*ϕ*) is the inverse screening length that depends on the packing fraction,^{15}

Specifically, in the first three dimensions, one has

### A. 1D systems

The long-ranged part of effective potential for the unit-step function *g*_{2} in one dimension exhibits exponential decay,^{15} i.e.,

We find that the maximum density at which the 1D unit-step function *g*_{2} is numerically realizable is *ϕ*_{m} = 0.49 = 0.98*ϕ*_{c}. The functional form of these potentials is given by

Table I lists the values of the optimized parameters for *ϕ*/*ϕ*_{c} = 7/8, 15/16, 31/32, and 0.98, where it is shown that the larger *ϕ* is, the more oscillatory terms must be included. Figure 2(a) shows the density-dependent effective potential at the aforementioned packing fractions, which clearly shows that for *ϕ*/*ϕ*_{c} ≥ 31/32, *v*(*r*; *ϕ*) exhibits significant oscillations (or “wiggles”) to negate the oscillations in *g*_{2}(*r*). The amplitude of the sum of the oscillatory terms increases with *ϕ*. As we will show in Secs. IV B and IV C, such oscillatory behavior is insignificant in *v*(*r*, *ϕ*) for 2D and 3D unit-step function *g*_{2} due to the decorrelation principle in higher dimensions.^{17} Figures 2(b) and 2(c) show the target and optimized pair statistics at the maximum realizable packing fraction *ϕ*_{m} = 0.49 = 0.98*ϕ*_{c}. We find that the *L*_{2} functions [Eqs. (7) and (8)] are $Dg2=6\xd710\u22124$, *D*_{S} = 0.0010 and the *L*_{2}-norm error (6) is $E=0.040$. These errors are an order of magnitude smaller than typical errors obtained via IHNCI for realizable targets.^{25} Figure 2(b) also shows *g*_{2}(*r*) for the pure equilibrium hard-rod system at *ϕ*_{m} = 0.49, which exhibits significant positive and negative correlations compared to the unit-step function *g*_{2}. This demonstrates that the unit-step function-*g*_{2} system is much more disordered than the pure equilibrium hard-rod system and that the effective potential (37) significantly suppresses the correlations beyond the hard core as will be quantitatively analyzed in Sec. IV D via the order metric *τ* (19). Figure 2(d) shows a 25-particle portion of a 400-configuration at *ϕ*_{m}.

ϕ/ϕ_{c}
. | 7/8 . | 15/16 . | 31/32 . | 0.98 . | ϕ/ϕ_{c}
. | 31/32 . | 0.98 . |
---|---|---|---|---|---|---|---|

ɛ_{1} | 6.47 | 10.12 | 11.18 | 29.52 | ɛ_{4} | −4.098 | −1.828 |

κ^{−1}(ϕ) | 1.165 | 1.641 | 2.440 | 2.704 | $\sigma 4(1)$ | 0.5710 | 3.152 |

ɛ_{2} | 2.809 | 0.7187 | −6.772 | −3.056 | $\sigma 4(2)$ | 0.2687 | 0.3254 |

$\sigma 2(1)$ | 0.7674 | 1.680 | 1.094 | 0.6015 | θ_{4} | 4.000 | 3.077 |

$\sigma 2(2)$ | −0.5470 | −0.3393 | 1.375 | −0.2195 | ɛ_{5} | 0.1934 | 0.04138 |

θ_{2} | 4.059 | 0.1577 | 2.050 | 3.332 | $\sigma 5(1)$ | 60.23 | 21.99 |

ɛ_{3} | ⋯ | 10.74 | −1.534 | −5.844 | $\sigma 5(2)$ | −425.8 | 2.073 |

$\sigma 3(1)$ | ⋯ | 0.5624 | 2.211 | 1.587 | θ_{5} | 1.781 | 4.089 |

$\sigma 3(2)$ | ⋯ | 0.7691 | 0.3275 | 1.420 | ɛ_{6} | 3.494 | 0.3899 |

θ_{3} | ⋯ | 3.317 | 2.960 | 2.902 | $\sigma 6(1)$ | 4.280 | 6.914 |

$\sigma 6(2)$ | 19.36 | 1.629 | |||||

θ_{6} | 0.8264 | 4.124 |

ϕ/ϕ_{c}
. | 7/8 . | 15/16 . | 31/32 . | 0.98 . | ϕ/ϕ_{c}
. | 31/32 . | 0.98 . |
---|---|---|---|---|---|---|---|

ɛ_{1} | 6.47 | 10.12 | 11.18 | 29.52 | ɛ_{4} | −4.098 | −1.828 |

κ^{−1}(ϕ) | 1.165 | 1.641 | 2.440 | 2.704 | $\sigma 4(1)$ | 0.5710 | 3.152 |

ɛ_{2} | 2.809 | 0.7187 | −6.772 | −3.056 | $\sigma 4(2)$ | 0.2687 | 0.3254 |

$\sigma 2(1)$ | 0.7674 | 1.680 | 1.094 | 0.6015 | θ_{4} | 4.000 | 3.077 |

$\sigma 2(2)$ | −0.5470 | −0.3393 | 1.375 | −0.2195 | ɛ_{5} | 0.1934 | 0.04138 |

θ_{2} | 4.059 | 0.1577 | 2.050 | 3.332 | $\sigma 5(1)$ | 60.23 | 21.99 |

ɛ_{3} | ⋯ | 10.74 | −1.534 | −5.844 | $\sigma 5(2)$ | −425.8 | 2.073 |

$\sigma 3(1)$ | ⋯ | 0.5624 | 2.211 | 1.587 | θ_{5} | 1.781 | 4.089 |

$\sigma 3(2)$ | ⋯ | 0.7691 | 0.3275 | 1.420 | ɛ_{6} | 3.494 | 0.3899 |

θ_{3} | ⋯ | 3.317 | 2.960 | 2.902 | $\sigma 6(1)$ | 4.280 | 6.914 |

$\sigma 6(2)$ | 19.36 | 1.629 | |||||

θ_{6} | 0.8264 | 4.124 |

Figure 3(a) shows the optimized potential of the form (37) at packing fraction *ϕ* = 127/128*ϕ*_{c} > *ϕ*_{m}, obtained after five iterations of reselecting the basis functions. Figures 3(b) and 3(c) show the corresponding optimized pair statistics. It is clear that the target pair statistics are not realized by the simulated many-body system. The optimized *g*_{2}(*r*) oscillates around unity for *r* > 1 and contains a peak at *r* = 1. Thus, the system exhibits clustering of particles despite the highly repulsive nature of the potential. The optimized *S*(*k*) contains a small peak at *k* ∼ *π*, and the *S*(*k*) values at small *k* are higher than those of the target *S*_{T}(*k*).

To test whether the non-convergence of the trial potential is due to the constraints of the forms of the basis functions, we included two additional exponential-damped oscillatory basis functions (28)–(37) and repeated the optimization procedure. We observed that convergence was still not reached, which suggests that the target is potentially indeed non-realizable.

We also applied the IHNCI procedure,^{20,24} which was the most accurate inverse algorithm prior to the development of our methodology^{25} on the 1D unit-step function *g*_{2} target at *ϕ* = 127/256 = 127/128*ϕ*_{c}. We found that (1) the trial *S*(*k*) at small *k* is larger than that of the target and (2) the trial *v*(*r*) at all *r* increased in each iteration due to the fact that the algorithm attempts to match the small values of *S*_{T}(*k*) around the origin. In about 100 iterations, the trial potential became so repulsive that the system crystallized, leading to a trial *g*_{2}(*r*) that is drastically different from *g*_{2,T}(*r*). Therefore, IHNCI is definitely not suited to treat the realizability problem in this case.

### B. 2D systems

The long-ranged part of effective potential for the unit-step function *g*_{2} in two dimensions is determined to be the exponentially screened inverse power law *r*^{−1/2}^{15} for *ϕ* < *ϕ*_{c}, i.e.,

We find that the 2D unit-step function *g*_{2} is numerically realizable up to the terminal density *ϕ*_{c} = 1/4. Figure 4(a) shows the density-dependent effective potential *v*(*r*, *ϕ*), whose functional form is given by

Table II lists the optimized values of the parameters. The potential becomes increasingly long-ranged with increasing *ϕ* and the characteristic length scale *κ*^{−1}(*ϕ*) diverges to infinity as $\varphi \u2192\varphi c\u2212$. At *ϕ* = *ϕ*_{c}, the long-ranged part of the potential becomes the 2D Coulomb form given by

As shown in Fig. 4(b), the equilibrium state at *ϕ*_{c} accurately yields the unit-step function *g*_{2}. Figure 4(c) shows that the optimized system is hyperuniform and its structure factor matches *S*_{T}(*k*). As in the 1D case, the precision of our methodology is evident from the small values of the *L*_{2} functions [Eqs. (7) and (8)] and the *L*_{2}-norm error (6), given by $Dg2=6\xd710\u22124$, *D*_{S} = 0.0022, and $E=0.053$. Figure 4(d) shows a snapshot of the optimized system at *ϕ*_{c}. The particles are well-separated and do not form large clusters.

ϕ/ϕ_{c}
. | ɛ_{1}
. | κ^{−1}(ϕ)
. | ɛ_{2}
. | $\sigma 2(1)$ . | $\sigma 2(2)$ . | θ_{2}
. |
---|---|---|---|---|---|---|

7/8 | 5.254 | 0.9684 | 54.40 | 0.2182 | −0.3018 | 0.6078 |

15/16 | 5.430 | 1.455 | 595.7 | 0.1612 | −0.4543 | 0.1322 |

31/32 | 5.977 | 2.093 | 25.06 | 0.2881 | −2.857 | 3.502 |

63/64 | 6.665 | 3.425 | −20.50 | 0.3172 | 2.854 | −0.4685 |

1 | 4 | + ∞ | 6.774 | 0.4250 | 0.4834 | 1.776 |

ϕ/ϕ_{c}
. | ɛ_{1}
. | κ^{−1}(ϕ)
. | ɛ_{2}
. | $\sigma 2(1)$ . | $\sigma 2(2)$ . | θ_{2}
. |
---|---|---|---|---|---|---|

7/8 | 5.254 | 0.9684 | 54.40 | 0.2182 | −0.3018 | 0.6078 |

15/16 | 5.430 | 1.455 | 595.7 | 0.1612 | −0.4543 | 0.1322 |

31/32 | 5.977 | 2.093 | 25.06 | 0.2881 | −2.857 | 3.502 |

63/64 | 6.665 | 3.425 | −20.50 | 0.3172 | 2.854 | −0.4685 |

1 | 4 | + ∞ | 6.774 | 0.4250 | 0.4834 | 1.776 |

### C. 3D systems

As in the case of the 1D and 2D unit-step function *g*_{2}, the long-ranged part of *v*(*r*; **a**) for the 3D unit-step function *g*_{2} is found to be the Yukawa form.^{15} The target pair statistics is realizable up to *ϕ*_{c} = 1/8, and the optimized density-dependent pair potential is given by

where *κ*(*ϕ*_{c}) = 0, giving a long-ranged potential with a Coulomb tail, i.e., *v*(*r*; *ϕ*_{c}) ∼ 1/*r*.

Figure 5(a) shows the density-dependent effective potential of the unit-step function *g*_{2} for packing fractions *ϕ* in the vicinity of *ϕ*_{c}. As *ϕ* increases, *v*(*r*; *ϕ*) becomes increasingly long-ranged. Table IV lists the optimized parameters in *v*(*r*; *ϕ*), from which it is obvious that the screening length of the potential *κ*^{−1}(*ϕ*) increases with *ϕ*. Figures 5(b) and 5(c) show that the pair statistics at *ϕ*_{c} is excellently produced by the equilibrium system under the effective potential as manifested by the small *L*_{2} functions [Eqs. (7) and (8)] $Dg2=4\xd710\u22124$ and *D*_{S} = 0.0016 as well as the small *L*_{2}-norm error (6) $E=0.045$. Figure 5(d) shows a snapshot of the system with unit-step function *g*_{2} at *ϕ*_{c}.

### D. Screening lengths and order metrics across dimensions

To study the dependence of the screening length *κ*^{−1}(*ϕ*) on the packing fraction and the dimensionality, we plot in Fig. 6 the optimized *κ*^{−1}(*ϕ*) against *ϕ*/*ϕ*_{c} for *d* = 1, 2, and 3 as well as the theoretical relation (34) shown as dashed curves. In all dimensions, the theoretical and optimized *κ*^{−1}(*ϕ*) agree well, revealing the applicability of Eqs. (35) and (34) for *ϕ* near but below *ϕ*_{c}. Nevertheless, note that the discrepancy between theoretical and optimized *κ*^{−1}(*ϕ*) is larger (on the order of 10%) for *ϕ*/*ϕ*_{c} ≥ 31/32 in all dimensions due to the decreased sensitivity of pair statistics on pair potentials for higher-density states.^{62} Importantly, for the same value of *ϕ*/*ϕ*_{c}, *κ*^{−1}(*ϕ*) decreases with *d*, reflecting the decorrelation principle for higher dimensions.^{17}

To quantitatively characterize the suppression of correlations in the unit-step function *g*_{2} relative to the corresponding equilibrium system with pure HS interactions at the same *ϕ*, we compute their order metrics *τ* (19) for *d* = 1, 2, 3 at the maximum realizable packing fractions for the unit-step function *g*_{2}. Equation (19) yields immediately that *τ* = *v*_{1}(1) for the unit-step function-*g*_{2} systems. The explicit expression of *S*(*k*) for the 1D pure equilibrium HS system is given in Ref. 43, which yields *τ* = 2.43 at *ϕ*_{m} = 0.49 = 0.98*ϕ*_{c}. Note that *τ* for the unit-step function *g*_{2} at *ϕ*_{m} is only 41% of that for the equilibrium HS system (Table III), indicating that the effective potential (37) achieves significant suppression of both positive and negative correlations that would otherwise be present in the pure equilibrium HS fluid.

d
. | ϕ
. | τ_{US}
. | τ_{HS}
. | τ_{US}/τ_{HS} (%)
. |
---|---|---|---|---|

1 | 0.98ϕ_{c} = 0.49 | 1 | 2.43 | 41 |

2 | ϕ_{c} = 1/4 | π ≈ 3.14 | 3.72 | 85 |

3 | ϕ_{c} = 1/8 | 4π/3 ≈ 4.19 | 4.57 | 92 |

d
. | ϕ
. | τ_{US}
. | τ_{HS}
. | τ_{US}/τ_{HS} (%)
. |
---|---|---|---|---|

1 | 0.98ϕ_{c} = 0.49 | 1 | 2.43 | 41 |

2 | ϕ_{c} = 1/4 | π ≈ 3.14 | 3.72 | 85 |

3 | ϕ_{c} = 1/8 | 4π/3 ≈ 4.19 | 4.57 | 92 |

ϕ/ϕ_{c}
. | ɛ_{1}
. | κ^{−1}(ϕ)
. | ɛ_{2}
. | $\sigma 2(1)$ . | $\sigma 2(2)$ . | θ_{2}
. |
---|---|---|---|---|---|---|

7/8 | 3.227 | 0.926 | 13.92 | 0.2615 | −0.3368 | 5.368 |

15/16 | 3.262 | 1.352 | 15.84 | 0.2528 | −0.2850 | 0.06881 |

31/32 | 3.166 | 1.985 | −13.96 | 0.2703 | −0.3203 | 2.610 |

63/64 | 3.295 | 2.751 | −25.47 | 0.2342 | −0.3079 | 3.102 |

1 | 3.378 | +∞ | −450.3 | 0.1521 | −0.3720 | 3.502 |

ϕ/ϕ_{c}
. | ɛ_{1}
. | κ^{−1}(ϕ)
. | ɛ_{2}
. | $\sigma 2(1)$ . | $\sigma 2(2)$ . | θ_{2}
. |
---|---|---|---|---|---|---|

7/8 | 3.227 | 0.926 | 13.92 | 0.2615 | −0.3368 | 5.368 |

15/16 | 3.262 | 1.352 | 15.84 | 0.2528 | −0.2850 | 0.06881 |

31/32 | 3.166 | 1.985 | −13.96 | 0.2703 | −0.3203 | 2.610 |

63/64 | 3.295 | 2.751 | −25.47 | 0.2342 | −0.3079 | 3.102 |

1 | 3.378 | +∞ | −450.3 | 0.1521 | −0.3720 | 3.502 |

The equilibrium pair statistics for the 2D and 3D pure equilibrium HS systems are not known exactly.^{5} Thus, we numerically computed *τ* from the pair statistics obtained via MC simulations. Table III shows the *τ* order metrics (19) for the unit-step function-*g*_{2} systems and the pure equilibrium HS packings. We observe that the discrepancy in the *τ* values for the two systems decreases as *d* increases, which again reflects the decorrelation principle for higher dimensions.^{17}

## V. 3D GHOST RSA TARGET

In this section, we compare the effective potential for the unit-step function *g*_{2} with that for a closely related model known as the “ghost” RSA process.^{11,29} It is a special case of a generalization of the standard RSA process^{5,40,41} and is the only known model for which all *n*-particle correlation functions are exactly solvable for any dimension d and allowed densities.^{29} In the ghost RSA process, spherical “test” particles of unit diameter are added continually to $Rd$ during time *t* ≥ 0 according to a translationally invariant Poisson process of density *η* per unit time, i.e., *η* 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 **r** at time *t* is retained if and only if it neither overlaps any existing sphere in the packing nor any previously rejected test sphere in the time interval [0, *t*).^{29} The maximum packing fraction for the ghost RSA process is *ϕ*(*t* = +*∞*) = 1/2^{d}, identical to the terminal density *ϕ*_{c} of the iso-*g*_{2} process for the unit-step function *g*_{2}.^{29} For finite *d*, the ghost RSA model is nonhyperuniform at *ϕ*_{c} and the corresponding *g*_{2}(*r*) possesses a small “bump” at *r* = 1 and is identically unity for *r* ≥ 2. The existence of such a bump at maximum density might suggest that the unit-step function *g*_{2} is unrealizable at *ϕ*_{c} for small *d*. Here, we study the pair statistics of the 3D ghost RSA process at the maximum packing fraction (*ϕ*_{c} = 1/8), whose *g*_{2,T}(*r*) is given by^{29}

Torquato *et al.*^{63} derived the exact expression for the corresponding structure factor, which is analytic at the origin. The small-*k* expansion is given by^{63}

which implies that *v*(*r*) exhibits exponential or superexponential decay at large *r*.^{17} Thus, we used the functional forms (27) and (28) for the basis functions.

Figure 7 shows the effective pair potential as well as target and optimized pair statistics for this model. The effective potential is given by

Table V lists the optimized parameters. We achieved excellent agreement between target and optimized pair statistics, with the *L*_{2} functions [Eqs. (7) and (8)] $Dg2=5\xd710\u22124$, *D*_{S} = 3 × 10^{−5}, and the *L*_{2}-norm error (6)$E=0.023$. Compared to the effective potential for the unit-step function *g*_{2} at the same packing fraction, the potential for the ghost RSA process is much shorter ranged as the latter system is nonhyperuniform and *S*(*k*) is analytic at the origin.^{25} The fact that *v*(*r*; *ϕ*_{c}) for the 3D ghost RSA is not highly repulsive means that the bump in its *g*_{2,T}(*r*) does not constitute a considerable constraint to the realizability of the 3D unit-step function *g*_{2} at *ϕ*_{c}.

## VI. A KNOWN 2D NON-REALIZABLE TARGET

To test the accuracy and robustness of our methodology, in this section, we study a 2D target *g*_{2} whose functional form consists of a unit-step function and a delta function separated by a gap.^{11,13} For specific parameters chosen, the target *g*_{2} satisfies all known necessary conditions for realizability (Sec. II D) but is not realizable due to geometrical constraints.^{11} This model has been used to show that the densest packings in high dimensions are disordered.^{11}

The functional form of the target *g*_{2}(*r*) is given by

where *σ*_{T} ≥ 1 is a distance parameter in units of *σ*. When *σ*_{T} = 1, *g*_{2,T} is known as the contact-*δ* plus step function^{14} and is the zero-density limit of the sticky hard-sphere (SHS) potential^{5} given by

which is the limit of the square well (SW) potential^{5,35} given by

as *σ*_{0} → 1, where *ɛ*_{0} and *σ*_{0} are energy and distance parameters in units of *ɛ* and *σ*, respectively. The contact-*δ* plus step function has been shown to be realizable for a finite range of *r* for *d* = 2, 3.^{10}

On the other hand, with the parameters *σ*_{T} = 1.2946, *Z* = 4.0138, and *ϕ* = *ρπ*/4 = 0.74803, *g*_{2,T}(*r*) is not realizable in two dimensions because of the fact that *Z* > 4 implies that there are some particles that are in contact with at least five particles. However, any arrangement of the five will result in nonzero *g*_{2}(*r*) within the targeted gap 1 < *r* < 1.2946.

The small-*k* behavior of the unrealizable target is given by *S*_{T}(*k*) ∼ 0.04689*k*^{2}, implying that the large-*r* behavior of the trial *v*(*r*) must be of a Coulomb form *v*(*r*) ∼ −log(*r*). The HNC approximation suggests that the small- and intermediate-*r* behavior of the trial potential is composed of at least three oscillatory functions damped by a power law *r*^{−1} (30). The form of the trial potential is thus given by

As expected, our inverse algorithm did not find an equilibrium state that matches the target pair statistics. The inverse procedure generated an optimized potential [Fig. 8(a)] that yields Ψ = 59.5, which is much larger than the convergence criterion Ψ < 5 × 10^{−4}. The simulated *g*_{2}(*r*) [Fig. 8(b)] contains extra *δ* peaks for *r* > 1 and is nonzero in the targeted gap 1 < *r* < 1.2946. The simulated *S*(*k*) [Fig. 8(c)] contains many sharp peaks that are not present in the target. We note that the IHNCI procedure for this target yielded diverging trial potentials, i.e., |*v*(*r*)| became unbounded for all *r* values except where *v*(*r*) = 0, which again highlights the insufficiency of such methods for our purpose.

## VII. CONCLUDING REMARKS

We have numerically investigated the realizability of the unit-step function *g*_{2} for *d* = 1, 2, 3 by determining precise density-dependent effective interactions that yield positive-temperature equilibrium states matching both target *g*_{2}(*r*) for all *r* and target *S*(*k*) for all *k*. For *d* = 1, we found that the unit-step function *g*_{2} is realizable up to *ϕ*_{m} = 0.49 = 0.98*ϕ*_{c}. For *d* = 2 and 3, it is realizable up to the terminal packing fraction *ϕ*_{c} = 1/2^{d}, where the equilibrium states are hyperuniform and *v*(*r*) at large *r* behave as the Coulomb interaction in the respective dimensions. This implies that the explicitly known necessary conditions for realizability described in Sec. II D are also sufficient up though *ϕ*_{c} for 2D and 3D unit-step function *g*_{2}. Furthermore, because the unit-step function *g*_{2} is realizable up to *ϕ*_{c} for *d* = 2 and 3, it is also realizable in all higher dimensions up to their respective maximum terminal packing fractions at which the systems are hyperuniform.^{11,17} For *ϕ* near but below *ϕ*_{c}, the large-*r* behavior of the effective potentials is given by exp[−*κ*(*ϕ*)*r*], *r*^{−1/2} exp[−*κ*(*ϕ*)*r*], and *r*^{−1} exp[−*κ*(*ϕ*)*r*] for *d* = 1, 2, and 3, respectively, where *κ*(*ϕ*) is the inverse screening length, and for *ϕ* = *ϕ*_{c}, the potentials at large *r* are given by the pure Coulomb forms in the respective dimensions as predicted in Ref. 15. The effective potentials at the maximum realizable packing fractions for the unit-step function *g*_{2} significantly suppress both positive and negative correlations that would otherwise be present in the pure equilibrium HS system as manifested by the smaller order metrics *τ* (19) for the corresponding former systems. The capacity to generate hyperuniform states at positive *T* via effective potentials is expected to facilitate the self-assembly of tunable hyperuniform soft-matter materials and enables one to probe the thermodynamic and dynamic properties of such states. Due to the long-ranged nature of these potentials (16), one can use *ρ* and *T* as tuning parameters to generate equilibrium hyperuniform states or states with stronger hyperuniform forms via their *inherent structures*, i.e., local energy minima.^{26}

While we do not provide a rigorous proof that the 1D unit-step function *g*_{2} is non-realizable at *ϕ* = *ϕ*_{c} = 1/2, the accuracy of our methodology and the fact that the hyperuniform targets are realizable in higher dimensions (*d* = 2 and *d* = 3) suggests that the 1D unit-step function *g*_{2} is highly likely to be non-realizable for *ϕ* = *ϕ*_{c}, even if it is realizable for almost all packing fractions below *ϕ*_{c}, i.e., for *ϕ* ≤ 0.49. We recall that it has been shown that the realizability of a particular functional form for hypothetical correlation function corresponding to a disordered system becomes easier as the space dimension increases due to a “decorrelation” principle,^{11} and hence 1D systems are the most likely to be non-realizable in general.

We also found that the realizable 3D nonequilibrium ghost RSA target at the maximum packing fraction *ϕ*_{c} = 1/8 is realizable as an equilibrium system with an effective pair potential, which further supports the Zhang–Torquato conjecture.^{13} Interestingly, its corresponding effective pair potential is much shorter ranged than the one we determined for the 3D unit-step function *g*_{2} at *ϕ*_{c} as the former system is nonhyperuniform and *S*(*k*) is analytic at the origin. Thus, despite the fact that its maximum density coincides with *ϕ*_{c} for the unit-step function *g*_{2}, the ghost RSA model does not provide a significant constraint on the realizability of the unit-step function *g*_{2}.

Importantly, we have demonstrated that our inverse methodology^{25} is robust and generally capable of probing realizability problems as it successfully generates effective potentials in cases when the target pair statistics is realizable. Furthermore, as expected, the trial potential does not converge for the 2D “contact-delta + gap + unit-step” *g*_{2} (45), which is known to be non-realizable, despite the fact that it satisfies all known explicit necessary conditions. Note that the non-realizability of (45) in two dimensions does not imply that it is non-realizable in higher dimensions. A fascinating venue for future research is to apply our methodology to investigate the realizability of (45) for *d* ≥ 3. If (45) is shown to be realizable in relatively low dimensions, such as *d* = 3 or 4, then one has further evidence that the densest packing of identical spheres in high dimensions are disordered (rather than ordered as in low dimensions), which is based on the terminal packing fraction corresponding to the higher-dimensional versions of the “contact-delta + gap + unit-step” *g*_{2}.^{11}

One can also apply our methodology to study the realizability of a wide range of prescribed pair statistics, including *g*_{2} with prescribed forms of oscillations^{14} and hyposurficial states.^{13,64} Comparison of the realizable density range for these targets with the terminal densities derived in Ref. 14 might shed light on necessary realizability conditions that are hitherto unknown.

One could also study the thermodynamic and dynamic properties of the effective potentials that yield unit-step and ghost RSA *g*_{2}’s, such as the associated phase diagrams, inherent structures, and diffusion properties.^{6,65} Unlike ground states and jammed hyperuniform states,^{17,66} the particle diffusion rates for the positive-temperature hyperuniform fluids at *ϕ*_{c} should be positive. This could occur if the systems’ potential only allows neighboring particle pairs to exchange positions. Such pair exchanges would allow the system to explore all permutationally equivalent structures but would prevent other types of structural modifications, thus retaining the system’s vanishing isothermal compressibility while permitting a positive self-diffusion rate. One could test this hypothesis by computing particle trajectories and velocities via molecular dynamics.

Finally, we note that our work provides a challenge to experimentalists to fabricate nanoparticles that result in the effective pair potentials found in this study. Achieving potentials close to (39) and (41) for *ϕ* < *ϕ*_{c} is expected to be experimentally feasible as the exponentially damped power-law tails can be readily attained by charged nanoparticles in solution^{67} and the exponential and superexponential short-ranged terms can be achieved by polymer-grafted nanoparticle surfaces.^{68} Generating interactions close to the ghost RSA potential (44) would also be feasible via polymer-grafted nanoparticles^{68} as (44) consists of only short-ranged superexponential terms.

If one can find realistic chemical compositions of nanoparticles that realize the unit-step function *g*_{2} near or at the terminal packing fractions, then such systems could be utilized for various practical applications. For example, nanoparticles interacting via (39) in the plane near *ϕ*_{c} suppress clustering and thus can be useful substrates for reactions that are sensitive to reactant clustering.^{69} Furthermore, the packings corresponding to the 2D and 3D unit-step function *g*_{2} at *ϕ*_{c} are hyperuniform states with no short-range order beyond the hard core, unlike typical hyperuniform sphere packings. Thus, they can be useful in optical applications^{70} and controlled drug delivery.^{71} While achieving the long-ranged interactions for hyperuniform states is challenging, one could reproduce the potentials over some finite but large range of *r* to fabricate *effectively hyperuniform* states, i.e., states with very small but nonvanishing positive value of *S*(0).^{72–76} Subsequently, the deviation of such systems from perfect hyperuniformity can be characterized via the various quantitative measures described in Ref. 43.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Haina Wang**: Data curation (lead); Formal analysis (lead); Investigation (lead); Validation (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). **Frank H. Stillinger**: Conceptualization (equal); Investigation (equal); Writing – review & editing (supporting). **Salvatore Torquato**: Conceptualization (lead); Investigation (equal); Supervision (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

**75**,

**106**, 11406 (2002).

^{2}) processes

^{2}) processes in equilibrium statistical mechanics

_{c}