An outstanding problem in statistical mechanics is the determination of whether prescribed functional forms of the pair correlation function g2(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-g2 process, i.e., the determination of density-dependent effective potentials that yield equilibrium states in which g2 remains invariant for a positive range of densities. Using a precise inverse algorithm that determines effective potentials that match hypothesized functional forms of g2(r) for all r and S(k) for all k, we show that the unit-step function g2, 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/2d, 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 g2 at ϕc; thus, it does not constrain the realizability of the unit-step function g2. 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-g2 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 g2(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 g2(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 , knowledge of one- and two-body correlations is insufficient to determine the corresponding higher-body correlation functions g3, g4, ….8–12 Consequently, given a prescribed functional form of a pair correlation function g2(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 g2(r) to be realizable,9 which are practically impossible to check. The ensemble-averaged structure factor S(k) is defined as
where h(r) = g2(r) − 1 is the total correlation function and 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 ϕ = ρv1(1/2), where v1(R) is the volume of a d-dimensional sphere of radius R,
Importantly, Eq. (2) is exactly the zero-density limit of g2(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 g2 given by (2) can be achieved for positive densities.
The g2-invariant process introduced by Torquato and Stillinger14 aims to determine the nonvanishing density range over which a prescribed form of g2(r) for a many-body system remains invariant over that range. They showed that a g2-invariant process possesses an upper terminal density ρc, which is the highest value such that all explicitly known necessary conditions are satisfied, including g2(r) ≥ 0 for all r, S(k) ≥ 0 for all k, and the Yamada condition7,9,10 as detailed in Sec. II D. The terminal packing fraction in the case of the step function (2) is given by14
at which the structure factor associated with (2) satisfies limk→0S(k) = 0, implying that the system must be hyperuniform15 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 g2 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 g2 is realizable but only for a finite range of r, in one,8,18 two,8 and three10 dimensions up to the terminal packing fraction. However, the reverse Monte-Carlo technique19 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 g2 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 g2-invariant process do not necessarily mean that the one-dimensional (1D) unit-step function g2 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 g2(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 process27–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 recently25 to investigate numerically the realizability problem via effective interactions. Specifically, we study the problem as an iso-g2 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 g2 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 g2(r). Density-dependent potentials are widely applied to describe coarse-grained models for polymers and macromolecules.36 The iso-g2 process is a nonequilibrium process whenever the target g2 is known to correspond to an equilibrium state: One then studies whether the same equilibrium g2 can be achieved for a range of nonvanishing densities. Since at positive densities, g2(r) for the pure equilibrium hard-sphere fluid deviates from the unit-step functional form (2) and becomes oscillatory,6 the corresponding iso-g2 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 g2(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 g2 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 methodology25 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 g2(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 potential38 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 introduced25 the following dimensionless L2-norm error:
where and DS are L2 functions given by
where g2,F(r; a) and SF(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 L2-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 methods20–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 g2 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 g2 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 g2 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) process29 is a generalization of the standard RSA process5,40,41 and is closely related to the unit-step iso-g2 process as its maximum packing fractions for all d are identical to ϕc = 1/2d. At the maximum density, g2 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 g2 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 g2. 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 g2 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 that are completely statistically characterized by the n-particle probability density functions ρn(r1, …, rn) for all n ≥ 1.6 In the case of statistically homogeneous systems, ρ1(r1) = ρ and ρ2(r1, r2) = ρ2g2(r), ρ is the number density in the thermodynamic limit, g2(r) is the pair correlation function, and r = r2 − r1. If the system is also statistically isotropic, then g2(r) is the radial function g2(r), where r = |r|.
For a single periodic configuration containing N point particles at positions r1, r2, …, rN within a fundamental cell F of a lattice Λ, the scattering intensity 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 VF 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 g2 [Eq. (2)] in dimensions 1, 2, 3 are given by16
where J1(x) is the first-order Bessel function. Note that at ϕc = 1/2d, S(k) ∼ k2 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 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 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 is hyperuniform if its variance grows in the large-R limit slower than Rd.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 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 g2(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)⟩ = ρv1(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 g2(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 fj(r/σ; aj) is the jth basis function, aj is a vector of parameters (generally consisting of multiple parameters), a = (a1, a2, …, an) is the “supervector” parameter, ɛ sets the energy scale, and σ is a characteristic length scale. Both ϵ and σ are taken to be unity. The components of aj are of four types: dimensionless energy scales ɛj, dimensionless distance scales σj, dimensionless phases θj, and dimensionless exponents pj.
The initial form of v(r; a) is informed by the small- and large-distance behaviors of the target pair statistics g2,T(r) and ST(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 cT(r) via6,53
where β = 1/(kBT) and kB is the Boltzmann constant, and the large-|r| behavior of c(r) can be extracted from the structure factor ST(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 g2,T(r) for all r or a targeted ST(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, fj(r; aj) are chosen from the possible following general forms:
- Hard core:This basis function is used for the unit-step function g2 and the RSA target as its g2,T(r) exhibits a hard core for r ≤ 1.(27)
- Exponential-damped oscillatory form:Importantly, the large-r behavior of v(r) ∼ −βc(r) for the 1D unit-step function g2 for (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.(28)
- Exponential-damped power-law form:55where κ 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 g2 in the limit for (but ϕ ≠ ϕc) is given by (29),15 with pj = 1/2 for d = 2 and pj = 1 for d = 3; see Sec. IV for further details.(29)
- Power-law damped oscillatory form:(30)
The last form is the Coulomb form (16), which corresponds to the large-r behavior of v(r) for the unit-step function g2 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 procedure56 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 and wS(k) are weight functions and g2(r; a) and S(k; a) correspond to an equilibrated N-particle system under v(r; a) at a dimensionless temperature kBT/ɛ = 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 g2(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/rd−α, 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 , wS(k) = exp(−k2/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-g2 PROCESS FOR THE UNIT-STEP FUNCTION g2
In this section, we present the realizable density ranges for the optimized effective interactions for the unit-step function g2 for d = 1, 2, 3. Taylor expansions of Eq. (11) at the point k = 0 reveal that S(k) for the unit-step function g2 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 by15,61
Furthermore, in the limit (but ϕ ≠ ϕc), the form of cT(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 g2 in one dimension exhibits exponential decay,15 i.e.,
We find that the maximum density at which the 1D unit-step function g2 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 g2(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 g2 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 L2 functions [Eqs. (7) and (8)] are , DS = 0.0010 and the L2-norm error (6) is . These errors are an order of magnitude smaller than typical errors obtained via IHNCI for realizable targets.25 Figure 2(b) also shows g2(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 g2. This demonstrates that the unit-step function-g2 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 | 0.5710 | 3.152 | |
ɛ2 | 2.809 | 0.7187 | −6.772 | −3.056 | 0.2687 | 0.3254 | |
0.7674 | 1.680 | 1.094 | 0.6015 | θ4 | 4.000 | 3.077 | |
−0.5470 | −0.3393 | 1.375 | −0.2195 | ɛ5 | 0.1934 | 0.04138 | |
θ2 | 4.059 | 0.1577 | 2.050 | 3.332 | 60.23 | 21.99 | |
ɛ3 | ⋯ | 10.74 | −1.534 | −5.844 | −425.8 | 2.073 | |
⋯ | 0.5624 | 2.211 | 1.587 | θ5 | 1.781 | 4.089 | |
⋯ | 0.7691 | 0.3275 | 1.420 | ɛ6 | 3.494 | 0.3899 | |
θ3 | ⋯ | 3.317 | 2.960 | 2.902 | 4.280 | 6.914 | |
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 | 0.5710 | 3.152 | |
ɛ2 | 2.809 | 0.7187 | −6.772 | −3.056 | 0.2687 | 0.3254 | |
0.7674 | 1.680 | 1.094 | 0.6015 | θ4 | 4.000 | 3.077 | |
−0.5470 | −0.3393 | 1.375 | −0.2195 | ɛ5 | 0.1934 | 0.04138 | |
θ2 | 4.059 | 0.1577 | 2.050 | 3.332 | 60.23 | 21.99 | |
ɛ3 | ⋯ | 10.74 | −1.534 | −5.844 | −425.8 | 2.073 | |
⋯ | 0.5624 | 2.211 | 1.587 | θ5 | 1.781 | 4.089 | |
⋯ | 0.7691 | 0.3275 | 1.420 | ɛ6 | 3.494 | 0.3899 | |
θ3 | ⋯ | 3.317 | 2.960 | 2.902 | 4.280 | 6.914 | |
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 g2(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 ST(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 methodology25 on the 1D unit-step function g2 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 ST(k) around the origin. In about 100 iterations, the trial potential became so repulsive that the system crystallized, leading to a trial g2(r) that is drastically different from g2,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 g2 in two dimensions is determined to be the exponentially screened inverse power law r−1/215 for ϕ < ϕc, i.e.,
We find that the 2D unit-step function g2 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 . 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 g2. Figure 4(c) shows that the optimized system is hyperuniform and its structure factor matches ST(k). As in the 1D case, the precision of our methodology is evident from the small values of the L2 functions [Eqs. (7) and (8)] and the L2-norm error (6), given by , DS = 0.0022, and . 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 . | . | . | θ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 . | . | . | θ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 g2, the long-ranged part of v(r; a) for the 3D unit-step function g2 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 g2 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 L2 functions [Eqs. (7) and (8)] and DS = 0.0016 as well as the small L2-norm error (6) . Figure 5(d) shows a snapshot of the system with unit-step function g2 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 g2 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 g2. Equation (19) yields immediately that τ = v1(1) for the unit-step function-g2 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 g2 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 . | . | . | θ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 . | . | . | θ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-g2 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 g2 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 process5,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 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/2d, identical to the terminal density ϕc of the iso-g2 process for the unit-step function g2.29 For finite d, the ghost RSA model is nonhyperuniform at ϕc and the corresponding g2(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 g2 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 g2,T(r) is given by29
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 by63
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 L2 functions [Eqs. (7) and (8)] , DS = 3 × 10−5, and the L2-norm error (6). Compared to the effective potential for the unit-step function g2 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 g2,T(r) does not constitute a considerable constraint to the realizability of the 3D unit-step function g2 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 g2 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 g2 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 g2(r) is given by
where σT ≥ 1 is a distance parameter in units of σ. When σT = 1, g2,T is known as the contact-δ plus step function14 and is the zero-density limit of the sticky hard-sphere (SHS) potential5 given by
which is the limit of the square well (SW) potential5,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, g2,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 g2(r) within the targeted gap 1 < r < 1.2946.
The small-k behavior of the unrealizable target is given by ST(k) ∼ 0.04689k2, 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 g2(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 g2 for d = 1, 2, 3 by determining precise density-dependent effective interactions that yield positive-temperature equilibrium states matching both target g2(r) for all r and target S(k) for all k. For d = 1, we found that the unit-step function g2 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/2d, 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 g2. Furthermore, because the unit-step function g2 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 g2 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 g2 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 g2 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 g2 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 g2, the ghost RSA model does not provide a significant constraint on the realizability of the unit-step function g2.
Importantly, we have demonstrated that our inverse methodology25 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” g2 (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” g2.11
One can also apply our methodology to study the realizability of a wide range of prescribed pair statistics, including g2 with prescribed forms of oscillations14 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 g2’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 solution67 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 nanoparticles68 as (44) consists of only short-ranged superexponential terms.
If one can find realistic chemical compositions of nanoparticles that realize the unit-step function g2 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 g2 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 applications70 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.