Ensemble Theory for Stealthy Hyperuniform Disordered Ground States

It has been shown numerically that systems of particles interacting with"stealthy"bounded, long-ranged pair potentials (similar to Friedel oscillations) have classical ground states that are, counterintuitively, disordered, hyperuniform and highly degenerate. Disordered hyperuniform systems have been receiving recent attention because they are distinguishable exotic states of matter poised between a crystal and liquid with novel properties. The task of formulating an ensemble theory that yields analytical predictions for the structural characteristics and other properties of stealthy degenerate ground states in $d$-dimensional Euclidean space is highly nontrivial because the dimensionality of the configuration space depends on the number density $\rho$ and there is a multitude of ways of sampling the ground-state manifold, each with its own probability measure. The purpose of this paper is to take some initial steps in this direction. Specifically, we derive general exact relations for thermodynamic properties that apply to any ground-state ensemble as a function of $\rho$ in any $d$, and show how disordered degenerate ground states arise as part of the ground-state manifold. We then specialize our results to the canonical ensemble by exploiting an ansatz that stealthy states behave remarkably like"pseudo"equilibrium hard-sphere systems in Fourier space. Our theoretical predictions for the structure and thermodynamic properties of the stealthy disordered ground states and associated excited states are in excellent agreement with computer simulations across the first three space dimensions. The development of this theory provides provide new insights regarding our fundamental understanding of the nature and formation of low-temperature states of amorphous matter. Our work also offers challenges to experimentalists to synthesize stealthy ground states at the molecular level.


I. INTRODUCTION
The equilibrium structure and phase behavior of soft matter systems span from the relatively simple, as found in strongly repulsive colloidal particles, to the highly complex, as seen in microemulsions and polymers [1][2][3][4][5][6][7][8][9][10]. Soft matter has been fruitfully microscopically modeled as classical many-particle systems in which the particles (or metaparticles) interact with effective pair potentials. Bounded (soft) effective interactions have been particularly useful in modeling polymer systems, and display zero-temperature ground states with a rich variety of crystalline structures, depending on the composition of the constituents and interaction parameters [1,4,[7][8][9][10]. * torquato@electron.princeton.edu We have previously used a "collective-coordinate" approach to generate numerically exotic classical ground states of many particles interacting with certain bounded isotropic long-ranged pair potentials in one-, two-and three-dimensional Euclidean space dimensions [11][12][13][14][15][16][17]. It was shown that the constructed ground states across dimensions are the expected crystal structures in a lowdensity regime [11,12,14,18], but above some critical density there is a phase transition to ground states that are, counterintuitively, disordered (structurally isotropic with no long-range order), hyperuniform, and highly degenerate [19]. These unusual amorphous states of matter have been shown to be endowed with novel thermodynamic and physical properties [14][15][16][20][21][22] and belong to the more general class of disordered hyperuniform systems, which have been attracting recent attention, as detailed below.
The disordered ground states are highly degenerate with a configurational dimensionality that depends on the density, and there are an infinite number of distinct ways to sample this complex ground-state manifold, each with its own probability measure. For these reasons, it is theoretically very challenging to devise ensemble theories that are capable of predicting structural attributes and other properties of the ground-state configurations.
A new type of statistical-mechanical theory must be invented to characterize these exotic states of matter. The purpose of this paper is to take some initial steps in this direction. However, to motivate the theoretical formalism, it is instructive to first briefly review the collectivecoordinate numerical procedure that we have used to achieve disordered ground states.
In the simplest setting, we previously examined pairwise additive potentials v(r) that are bounded and integrable such that their Fourier transformsṽ(k) exist. If N identical point particles reside in a fundamental region F of volume v F in d-dimensional Euclidean space R d at positions r N ≡ r 1 , . . . , r N under periodic boundary conditions, the total potential energy can be expressed in terms ofṽ(k) as follows: whereñ(k) = N j=1 exp(−ik · r j ) is the complex collective density variable, which can be viewed as a nonlinear transformation from the finite set of particle coordinates r 1 , . . . , r N to the complex functionsñ(k) that depend on the infinite set of wave vectors k in reciprocal space appropriate to the fundamental cell F . The crucial idea is that ifṽ(k) is defined to be bounded and positive with support in the radial interval 0 ≤ |k| ≤ K and if the particles are arranged so that |ñ(k)| 2 , a quantity proportional to the structure factor S(k), is driven to its minimum value of zero for all wave vectors whereṽ(k) has support (except k = 0), then it is clear from relation (1) that the system must be at its ground state or global energy minimum. We have referred to these ground-state configurations as "stealthy" [14] because the structure factor S(k) (scattering pattern) is zero for 0 < |k| ≤ K, meaning that they completely suppress single scattering of incident radiation for these wave vectors, and thus, are transparent at the corresponding wavelengths [23]. Various optimization techniques were employed to find the globally energy-minimizing configurations within an exceeding small numerical tolerance [11][12][13][14][15][16][17]. Generally, a numerically obtained ground-state configuration depends on the number of particles N within the fundamental cell, initial particle configuration, shape of the fundamental cell, and particular optimization technique employed.
As the number of k vectors for which |ñ(k)| is constrained to be zero increases, i.e., as K increases, the dimensionality of the ground-state configuration manifold D C decreases. Because |ñ(k)| is inversion-symmetric, the number of wave vectors contained in a sphere of radius K centered at the origin must be an odd integer, say 2M (K) + 1, and thus M (K) is the number of independently constrained wave vectors [24]. A parameter inversely proportional to density provides a measure of the relative fraction of constrained degrees of freedom compared to the total number of degrees of freedom d(N − 1) (subtracting out the system translational degrees of freedom). We will show in Sec. IV B that the dimensionality of the configuration space per particle is given by d(1 − 2χ) in the thermodynamic limit.

Compression Dilation
Direct Space Reciprocal Space S(k) FIG. 1. Schematic illustrating the inverse relationship between the direct-space number density ρ and relative fraction of constrained degrees of freedom χ for a fixed reciprocalspace exclusion-sphere radius K (where dark blue k-points signify zero intensity with green, yellow and red points indicating increasingly larger intensities) for a stealthy ground state. A compression of a disordered ground-state configuration with a fixed number of particles N in direct space leads to a dilation of the lattice spacing in reciprocal space. This means that during the compression process, the k-points for which |ñ(k)| is zero associated with the initial uncompressed system move out of the exclusion zone, i.e., the value of M (K) [cf. (2)] decreases. Since there are fewer constrained degrees of freedom (dimensionality of the ground-state configuration manifold increases), the disordered direct-space configuration becomes less spatially correlated. For a fixed N in the limit ρ → ∞ (i.e. system volume vF → 0), every k-point (except the origin) is expelled from the exclusion zone, and the system tends to an ideal-gas configuration [25], even if not an ideal gas thermodynamically, as will be shown in Sec. IV B. χ = 0, it is a non-interacting classical ideal gas. While it is not usual to think of classical ideal-gas configurations as ground states, at T = 0, they indeed are global energyminimizing states that are highly degenerate and typically disordered for N large enough. While the groundstate manifold contains periodic configurations (e.g., Bravais lattices and lattices with a basis), these are sets of zero measure in the thermodynamic limit. Clearly, if χ is made to be positive but very small, the ground states remain disordered and highly degenerate, even if the dimensionality of the configuration space D C has now been suddenly reduced due to the imposed constrained degrees of freedom, the number of which is determined by the radius K. From relation (2), we see that if K is fixed, configurations with ideal-gas-like pair correlation functions correspond to the limit χ → 0 or equivalently, to the limit ρ → ∞ [25]. The latter situation runs counter to traditional understanding that ideal-gas configurations correspond to the opposite zero-density limit of classical systems of particles. The reason for this inversion of limits is due to the fact that a compression of the system in direct space leads to a dilation of the lattice spacing in reciprocal space, as schematically shown in Fig. 1. While it is not surprising that the configuration space is fully connected for sufficiently small χ, quantifying its topology as a function of χ for all allowable χ is an outstanding problem, which is discussed further in the Conclusions.

FIG. 2. A disordered nonhyperuniform configuration (top)
and a disordered hyperuniform configuration (bottom). The latter is arrived at by very small collective displacements of the particles above. (Each particle on average moves a rootmean-square distance that is about an order of magnitude smaller than the mean-nearest-neighbor distance as measured by the configuration proximity metric [26].) These two examples show that it can be very difficult to detect hyperuniformity by eye, and yet their large-scale structural properties are dramatically different ("hidden order").
It is noteworthy that stealthy point patterns (disordered or not) constitute a special class of so-called hyperuniform states of matter. Hyperuniform systems are characterized by vanishing (normalized) density fluctuations at large length scales, i.e., the structure factor S(k) tends to zero in the limit |k| → 0 [27]; see Sec. II A for details. The hyperuniformity concept provides a means of categorizing crystals, quasicrystals and special disordered systems according to the degree to which large-scale density fluctuations are suppressed [27,28]. Disordered hyperuniform patterns, of which disordered stealthy systems are special cases, behave more like crystals in the manner they suppress large-scale density fluctuations, and yet also resemble typical statistically isotropic liquids and glasses with no Bragg peaks. In this sense, they have a "hidden order" on large length scales that is not apparent at small length scales, even if short-range order is present; see Fig. 2 for a vivid illustration. During the last decade, a variety of disordered hyperuniform states have been identified that exist as both equilibrium and nonequilibrium phases, including maximally random jammed particle packings [29][30][31][32], jammed athermal granular media [33], jammed thermal colloidal packings [34,35], cold atoms [36], transitions in nonequilibrium systems [37,38], surface-enhanced Raman spectroscopy [39], terahertz quantum cascade laser [40], wave dynamics in disordered potentials based on supersymmetry [41], avian photoreceptor patterns [42], and certain Coulombic systems [43]. Moreover, disordered hyperuniform materials possess novel physical properties potentially important for applications in photonics [20-22, 44, 45] and electronics [46][47][48].
The well-known compressibility relation from statistical mechanics [49] provides some insights about the relationship between temperature T and hyperuniformity for equilibrium systems at number density ρ: We see that any ground state (T = 0) in which the isothermal compressibility κ T is bounded and positive must be hyperuniform because the structure factor S(k = 0) must be zero. This includes crystals as well as disordered ground states such as stealthy ones. However, in order to have a hyperuniform system at positive T , the isothermal compressibility must be zero, i.e., the system must be incompressible [17]. Subsequently, we will use relation (3) to draw some conclusions about the excited states associated with stealthy ground states. Our general objective is the formulation of a predictive ensemble theory for the thermodynamic and structural properties of stealthy degenerate disordered ground states in arbitrary space dimension d that complements previous numerical work on this topic [11][12][13][14][15][16]. After providing basic definitions and describing a family of stealthy potentials (Secs. II and III), we derive general exact relations for thermodynamic properties (energy, pressure and isothermal compressibility) that apply to any well-defined ground-state ensemble as a function of the number density or, equivalently, χ in any space dimension d (Sec. IV) . We subsequently derive some exact integral conditions that both the pair correlation function g 2 (r) and structure factor S(k) must obey (Sec. V). The existence of periodic stealthy ground states enables us to show how disordered degenerate ground states arise as part of the ground-state manifold for sufficiently small χ (Sec. VI). Subsequently, we derive analytical formulas for the pair statistics [g 2 (r) and S(k)] for sufficiently small χ in the canonical ensemble in the limit that the temperature T tends to zero (Sec. VII) by exploiting an ansatz that stealthy states behave like "pseudo" equilibrium hard-sphere systems in Fourier space. Our theoretical predictions for g 2 (r) and S(k) are in excellent agreement with computer simulations across the first three space dimensions. These results are then used to predict, with high accuracy, other structural characteristics of stealthy ground states across dimensions, such as order metrics, local number variance and nearest-neighbor functions (Secs. VIII and IX). Subsequently, we derive analytical formulas for the structure factor and thermal expansion coefficient for the associated excited states for sufficiently small temperatures (Sec. X). Finally, we provide concluding remarks in Sec. XI.

II. DEFINITIONS AND PRELIMINARIES
Roughly speaking, a point process in d-dimensional Euclidean space R d is a distribution of an an infinite number of points in R d with configuration r 1 , r 2 , . . . at a well-defined number density (number of points per unit volume). For a statistically homogeneous point process in R d at number density ρ [50], the quantity ρ n g n (r n ) is the probability density associated with simultaneously finding n points at locations r n ≡ r 1 , r 2 , . . . , r n in R d [49]. With this convention, each n-particle correlation function g n approaches unity when all of the points become widely separated from one another. Statistical homogeneity implies that g n is translationally invariant and hence only depends on the relative displacements of the positions with respect to any chosen system origin, e.g., g n = g n (r 12 , r 13 , . . . , r 1n ), where r ij = r j − r i .
The pair correlation function g 2 (r) is a particularly important quantity. If the point process is also rotationally invariant (statistically isotropic), then g 2 depends on the radial distance r ≡ |r| only, i.e., g 2 (r) = g 2 (r). Thus, it follows that the expected number of points Z(R) found in a sphere of radius R around a randomly chosen point of the point process, called the cumulative coordination function, is given by where s 1 (r) = 2π d/2 r d−1 /Γ(d/2) is the surface area of a d-dimensional sphere of radius r. The total correlation function h(r) is trivially related to g 2 (r) as follows: When there are no long-range correlations in the system, h(r) → 0 or, equivalently, g 2 (r) → 1 as |r| → ∞. The structure factor S(k), which plays a prominent role in this paper, is related to the Fourier transform of h(r), denoted byh(k), via the expression A lattice Λ in R d is a subgroup consisting of integer linear combinations of vectors that constitute a basis for R d and thus represents a special subset of point processes.
Here, the space can be geometrically divided into identical regions F called fundamental cells, each of which contains just one point specified by the lattice vector p = n 1 a 1 + n 2 a 2 + · · · + n d−1 a d−1 + n d a d , where a i are the basis vectors for a fundamental cell and n i spans all the integers for i = 1, 2, · · · d. We denote by v F the volume of F . A lattice is called a Bravais lattice in the physical sciences. Unless otherwise stated, we will use the term lattice. Every lattice has a dual (or reciprocal) lattice Λ * in which the lattice sites are specified by the dual (reciprocal) lattice vector q · p = 2πm for all p , where m = 0, ±1, ±2, ±3 · · · . The dual fundamental cell F * has volume v F * = (2π) d /v F . This implies that the number density ρ Λ of Λ is related to the number density ρ Λ * of the dual lattice Λ * via the expression Some common d-dimnesional lattices are mathematically defined in Appendix A.
A periodic point process (crystal) is a more general notion than a lattice because it is obtained by placing a fixed configuration of N points (where N ≥ 1) within a fundamental cell F of a lattice Λ, which is then periodically replicated. Thus, the point process is still periodic under translations by Λ, but the N points can occur anywhere in F ; see Fig. 3.

A. Hyperuniform Point Processes
Consider uniformly sampling the number of points that are contained within a spherical window of radius R of a point process in R d . A hyperuniform point process has the property that the local number variance σ 2 (R) grows more slowly than R d [27]. Because σ 2 (R) is exactly related to a d-dimensional volume integral of the structure factor S(k) [27] (see Sec. VIII), this implies that infinite-wavelength density fluctuations (appropriately normalized) vanish, i.e., S(k) obeys the condition lim |k|→0 S(k) → 0.
For a Poisson (spatially uncorrelated) point process and many disordered point patterns, including typical liquids and structural glasses, the number variance grows like the volume of the window, i.e., σ 2 (R) ∼ R d , implying that S(k) is positive at k = 0. All perfect crystals and quasicrystals are hyperuniform such that σ 2 (R) ∼ R d−1 ; in other words, the variance grows like window surface area. By contrast, it is much more unusual to find disordered systems that are also hyperuniform. In recent years, evidence has been emerging that disordered hyperuniform many-particle systems can be regarded to be new distinguishable states of disordered matter (see examples given in the Introduction). Whenever the structure factor goes to zero with the power-law form S(k) ∼ |k| α , the number variance has the following large-R asymptotic scaling that depends on the value of the exponent α [17,28]: Since disordered as well as ordered stealthy states can be viewed as systems in which α tends to infinity, we see from (10) that they have the asymptotic scaling σ 2 (R) ∼ R d−1 . We give theoretical predictions for the variance of disordered stealthy ground states in Sec. VIII.

III. FAMILIES OF STEALTHY PAIR POTENTIALS
As we will see in the next section, the specific form of a stealthy potential does not affect the ground-state energy manifold, but it can affect other thermodynamic properties, such as the pressure. This has consequences in simulations of such properties, especially with respect to convergence issues. Hence, it is instructive to remark on some mathematical aspects of the long-range nature of the direct-space stealthy potentials, which are very similar to the weakly decaying Friedel oscillations of the electron density in a variety of systems, including molten metals as well as graphene [51,52]. As we will see, in some cases, stealthy potentials may mimic effective interactions that arise in certain polymer systems [10].
Here we will limit ourselves to pair potentials v(r) that are radial functions in R d , where r = |r|, and therefore their Fourier transformsṽ(k) are also radial functions in R d , where k ≡ |k| is a wavenumber. The d-dimensional Fourier transform of any integrable radial function f (r) in R d is given by [27] and the inverse transform off (k) is given by where J ν (x) is the Bessel function of order ν.
Consider the class of "stealthy" radial potential functionsṽ(k) in R d that are bounded and positive with compact support in the radial interval 0 where for simplicity V (k) is infinitely differentiable in the open interval [0, K) and is the Heaviside step function. The corresponding directspace radial pair potential v(r) is necessarily a delocalized, long-ranged function that is integrable in R d . Moreover, without any loss of generality, it will be assumed that V (k) ≤ v 0 . For concreteness and purposes of illustration, we will examine properties of two specific families of potentials that fall within the aforementioned wide class of stealthy interactions: "power-law" and "overlap" potentials.

A. Power-Law Potentials
The power-law potentials are defined in Fourier space as follows:ṽ where the exponent m can be any of the whole numbers. The corresponding direct-space potential v(r) will depend on d for any given m and is exactly given by is a special case of the generalized hypergeometric function p F q (a 1 , ..., a p ; b 1 , ..., b q ; x) [53]. Because the potential (16) is derived from the Fourier power-law potential (15), we will refer to (16) as the direct-space power-law potential. In the instance when m = 0 in (15) (i.e., simple step function), this expression for v(r) simplifies as follows: for which the large-r asymptotic behavior is given by (18) For any fixed d and m ≥ d, the direct-space power-law potential has the asymptotic form where s(r; m, d) is a bounded function (a sinusoidal function or constant of order one). For any fixed d and 1 ≤ m < d, the long-range oscillations of v(r) are controlled by an envelope that decays like 1/r β , where In Fig. 4, we plot the Fourier power-law potential for selected values of m (which applies in any dimension) and the corresponding direct-space potentials for d = 3. In all cases, we set v 0 = K = 1. It is to be noted that the amplitudes of the oscillations in v(r) decrease as m increases for a fixed dimension.

B. Overlap Potentials
Let α(r; R) represent the intersection volume of two identical d-dimensional spheres of radius R (scaled by the volume of sphere) whose centers are separated by a distance r. This quantity is known analytically in any space dimension, and has a variety of representations [54], including the following: where c(d) is the d-dimensional constant given by For d = 1, 2, 3 and 4, we respectively have Consider the class of "overlap" potentials, which for any dimension is given bỹ Note that for d = 1, the overlap potential is identical to power-law potential when d = 1 and m = 1. The thermodynamics of the ground-state manifold of this potential in the special case d = 2 was numerically investigated in Refs. [15] and [16]. It follows from (26) that the corresponding direct-space overlap potential is given by which is clearly nonnegative for all r. Its large-r asymptotic behavior is given by revealing that the long-ranged decay of the direct-space overlap potential has an envelope controlled by the inverse power law 1/r d+1 . Figure 5 depicts the overlap potentialṽ(k) for the first three space dimensions and the corresponding directspace overlap potentials v(r), the latter of which vividly shows the increasing decay rate of v(r) with increasing dimension. The direct-space overlap potential v(r) is similar in functional form to effective positive pair interactions that arise in mulitlayered ionic microgels [10].

IV. ENSEMBLE THEORY FOR STEALTHY DISORDERED GROUND STATES: EXACT RESULTS FOR THERMODYNAMIC PROPERTIES
Our general objective is the formulation of an ensemble theory for the thermodynamic and structural properties of stealthy degenerate disordered ground states that we previously investigated numerically [11][12][13][14][15][16]. In this section, we derive general exact relations for thermodynamic properties that apply to any well-defined ensemble as generated by a particular way to sample the stealthy disordered ground-state manifold as a function of number density ρ. In the subsequent section, we will derive some exact results for the pair statistics for general ensembles.

A. Preliminaries
To begin, consider a configuration of N identical particles with positions r N ≡ r 1 , . . . , r N in a large region of volume V in d-dimensional Euclidean space R d . For particles interacting via a pair potential v(r), the total potential energy Φ(r N ) is given by where r ij = r j − r i . Of particular interest are classical ground states, i.e., those configurations that minimize the energy per particle Φ(r N )/N . The ensemble average of the energy (29) per particle u in the thermodynamic limit can be written in terms of the pair correlation function g 2 (r): where angular brackets denote an ensemble average and ρ is the number density in the thermodynamic limit.
Because the collective-coordinate approach relies on the Fourier representation of the energy, we recast (30) in terms of the structure factor S(k) by applying Parseval's theorem to the second line of Eq. (30): (31) whereṽ(k) andh(k) are the Fourier transforms of v(r) and h(r), respectively, both of which are assumed to exist, k is a wave vector, and S(k) is the structure factor defined in relation (6). Note that the structure factor is a nonnegative, inversion-symmetric function, i.e.,

B. Ground-State Energy and Dimensionality of Its Configuration Space
Consider a stealthy potential functionṽ(k) with support in 0 ≤ k ≤ K of the class specified by (13). In light of (31), it clear that whenever particle configurations in R d exist such that S(k) is constrained to achieve its minimum value of zero for 0 ≤ k ≤ K, the system must be at its ground state or global energy minimum. This follows because the integrandṽ(k)S(k) in the nontrivial term on the right side of (31) is identically zero due to the conflicting demands of the step functions. When such configurations exist , the average ground-state energy per particle in any well-defined ensemble is given exactly by which is a structure-independent constant, depending on the density ρ, v 0 ≡ṽ(k = 0), v(r = 0) and d, as explicitly shown below. Importantly, because u is a constant independent of the structure, its value has no effect on the ground-state manifold, which is generally degenerate, and hence this manifold is invariant to the specific choice of the stealthy functionṽ(k) at fixed ρ and d.
We would like to express the energy (33) in terms of the parameter χ, defined by relation (2), which measures the relative number of independently constrained degrees of freedom for a finite system under periodic boundary conditions. Note that, in the thermodynamic limit, M (K) in (2) is simply half of the volume of a sphere of radius K [due to the inversion symmetry of S(k)] multiplied by the density ρ Λ * of the dual lattice [cf. Eq. (8)], i.e., where we have used the fact that ρ → N ρ Λ in this distinguished limit. Hence, from (2), we obtain the following expression for χ in the thermodynamic limit: where is the volume of a d-dimensional sphere (hypersphere) of radius R. We see that for fixed K and d, which fixes the potential, χ is inversely proportional ρ, which is the situation that we usually consider in this paper [55]. Hence, as χ tends to zero, ρ tends to infinity, which configurationally corresponds counterintuitively to the uncorrelated ideal-gas limit (Poisson distribution), as discussed in the Introduction. As χ increases from zero, the density ρ decreases and the dimensionality of the groundstate configuration manifold D C decreases. The configurational dimensionality per particle in the thermodynamic limit, d C , can easily obtained from the relation D C = dN − 2M (K) for a finite system [56]; specifically, where d C = lim DC →∞,N →∞ D C /N . Equations (33) and (35) yields the average groundstate energy per particle to be where ρ * min is defined in Sec. VI, and is a constant whose value depends on the specific form of the stealthy-potential classṽ(k) defined by (13), and hence must lie in the interval (0, 1], where γ = 1 corresponds to the step-function choiceṽ(k) = v 0 Θ(K − k). While the system in the limit χ → 0 (ρ → ∞) corresponds configurationally to an ideal gas in so far as the pair correlation function is concerned, as we will explain in detail in Sec. VII A, thermodynamically, it is nonideal; see Eq. (38) for u and Eq. (41) for the pressure.

C. Energy Route to Pressure and Isothermal Compressibility
The pressure in the thermodynamic limit at T = 0 can be obtained from the energy per particle via the relation Therefore, for stealthy potentials, we see from (33) that the ground-state pressure, for all possible values of ρ or χ, is given by the following simple expression: where Hence, the isothermal compressibility κ T ≡ ρ −1 ∂ρ ∂p T of such a ground-state is We see that as ρ tends to infinity, the compressibility tends to zero. Two important remarks are in order. First, estimates of the pressure obtained from simulations that we previously performed for d = 2 [15,16] as well as those carried out in the present study across the first three space dimensions (Appendix B) are in very good agreement with the exact result (41) across a wide range of densities, thus validating the accuracy of the simulations. Second, the fact that the pressure (41) is a continuous function of density implies that any phase transition that may take place could be a continuous one, the implications of which are discussed in the Conclusions.

D. Virial Route to Pressure and Isothermal Compressibility
An alternative route to the pressure for a radial pair potential function v(r) is through the "virial" equation, which at T = 0 in the thermodynamic limit is given by Although the pressure obtained via the virial route is generally expected to be equivalent to that obtained from the energy route (as described in the previous section), we will show that, for a certain class of stealthy potentials, the pressure obtained from (43) is either ill-defined or divergent. This has practical implications for what types of stealthy potentials can be used in constant-pressure simulations.
It is convenient to re-write the virial relation (43) in the following form: (44) whereF (k) is the Fourier transform of F (r) ≡ rdv/dr, when it exists, and we have used Parseval's theorem and definition (5) for the total correlation function h(r).
To continue with this analysis, we will make use of the following lemma.
Consider a bounded radial functionz(k) with compact support on the radial interval [0, K] in R d that is an infinitely differentiable in the open interval [0, K). Therefore, its Fourier transform z(r) exists.
Lemma 1: The Fourier transform of the radial function w(r) = rdz/dr in R d is given bỹ Proof: Differentiation ofz(k) [defined via (11)] with respect to k leads to the following identity: The Fourier transform of w(r) is given bỹ Integrating relation (47) by parts and use of (46) proves the lemma. Corollary: It immediately follows from Lemma 1 that w(k) has the same support asz(k) and meaning that the volume integral of rdz/dr over all space is proportional to the corresponding volume integral of z(r).
Note that by the Corollary of Lemma 1,F (0) = −d · v(0) and hence we can re-write the virial relation (44) as The second term in the second line of (49) follows because ρh(k) = −Θ(K − k) inside the exclusion sphere of radius K [see also Eq. (52) below] and has support in this exclusion zone by the Corollary of Lemma 1. But this second term must vanish in light of the trivial identity We see that the virial ground-state pressure (49) agrees with that of the pressure obtained via the energy per particle [c.f. (41)] provided thatF (k) exists. Since the latter is a stronger condition than the existence ofṽ(k), it is possible to devise a stealthy functionṽ(k) for which F (k) does not exist and hence a virial pressure that either diverges or is nonconvergent. For example, this problem occurs for the power-law potential (15) with m = 0 (step function) for any dimension d. By contrast, the virial pressure is always well-defined for the overlap potential (26) in any dimension [57]. This example serves to illustrate the mathematical subtleties that can arise due to the long-ranged nature of stealthy potentials in direct space.

V. ENSEMBLE THEORY FOR STEALTHY DISORDERED GROUND STATES: EXACT INTEGRAL CONDITIONS ON THE PAIR STATISTICS
Here we derive some exact integral conditions that must be obeyed by both the pair correlation function g 2 (r) and the structure factor S(k) for stealthy ground states that apply to general ensembles. These analytical relations can be profitably employed to test corresponding computer simulation results.

A. General Properties
In any stealthy ground state, the structure factor attains its minimum value S(k) = 0 for 0 < k ≤ K, and hence has the form where Θ(x) is the Heaviside step function defined by (14) andQ(k) = S(k)− 1 is a function that obeys the inequal-ityQ(k) ≥ −1. Therefore, from (6), we have that the Fourier transform of the total correlation function h(r) has the form wheref It is noteworthy that the functionf (k) is identical to the Mayer-f function for an equilibrium hard-sphere system in direct space. Taking the inverse Fourier transform of (52) yields the direct-space total correlation function to be given by where f (r) = − K 2πr and where the lower bound on P (r) indicated in (57) follows from the fact that h(r) ≥ −1 for all r for any point pattern. It trivially follows that since ρh(k = 0) = ρ R d h(r)dr = −1, the volume integral of P (r) must be zero, i.e., Less trivially, because the productṽ(k)P (k) is zero for all k, then by Parseval's theorem, we have the integral condition Thus, the functions v(r) and P (r) are orthogonal to one another. The exact integral conditions (58) and (59) can be used to test the accuracy of numerical methods that yield estimates of the pair correlation function.

B. Behavior of the Pair Correlation Function Near the Origin
It is instructive to determine the behavior of the pair correlation function g 2 (r) for small r. Substitution of the general form (11) forh(k) into the definition of the total correlation function h(r) as obtained from (12), and expanding h(r) in a Taylor series around r = 0 through second-order in the radial distance r yields h(r) = h(r = 0) + 1 2 where and the corresponding curvature is Therefore, from (61), we see that the pair correlation function at the origin is given by Since g 2 (r) must be nonnegative for all r, we have the following integral condition onQ(k): Hence, this integral must be positive for We also conclude from (62) that for h(r) or g 2 (r) to have positive curvature at the origin,Q(k) must obey the additional integral condition: Finally, we note that when g 2 (r = 0) = 0, the results above yield the equality and, because the curvature must be positive in this instance, the inequality (66) must generally be obeyed. The inequality (64), conditional inequality (66) and conditional equality (67) provide integral conditions to test the accuracy of numerical methods that yield estimates of the structure factor.

VI. EXISTENCE OF STEALTHY DISORDERED DEGENERATE GROUND STATES
It is noteworthy that any periodic crystal with a finite basis is a stealthy ground state for all positive χ up to its corresponding maximum value χ max (or minimum value of the number density ρ min ) determined by its first positive Bragg peak k Bragg [minimal positive wave vector for which S(k) is positive]. Tables I-IV list the pair χ max ,ρ min for some common periodic patterns in one, two, three and four dimensions, respectively, all of which are part of the ground-state manifold; see Appendix A for mathematical definitions. (The crystals denoted by Dia d and Kag d are d-dimensional generalizations of the diamond and kagomé crystals, respectively, for d ≥ 2 [58].) While the mere existence of such periodic ground states does not provide any clues about their occurrence probability in some ensemble, we will use these results here to show how disordered degenerate ground states arise as part of the ground-state manifold for sufficiently small χ.       Tables  I-IV, which we call ρ * min , corresponds to the dual of the densest Bravais lattice in direct space, and represents the critical density value below which a stealthy ground state does not exist for all k ≤ |k * Bragg |. The fact that ρ * min corresponds to the body-centered-cubic (bcc) lattice for d = 3 was initially shown analytically in Ref. [18] and subsequently numerically in Ref. [14]. We note that the values of ρ min for the simple hexagonal lattice and hexagonal close-packed crystal for d = 3 reported in Ref. [18] are incorrect because those calculations were based on the erroneous assumption that the structure factors at the corresponding shortest reciprocal lattice vectors have non-vanishing values.
Observe that in the case d = 1, there is no non-Bravais lattice (periodic structure with a basis n ≥ 2) for which χ max is greater than 1/2, implying that the ground-state manifold is nondegenerate (uniquely the integer lattice) for 1/2 < χ ≤ 1. This is to be contrasted with the cases d ≥ 2 where the ground-state manifold must be degenerate [60] for 1/2 < χ < χ * max and nondegenerate only at the point χ = χ * max , as implied by Tables II-IV. Here χ * max is the largest possible value of χ max in some fixed dimension.
Lemma: At fixed K, a configuration comprised of the union (superposition) of m different stealthy groundstate configurations in R d with χ 1 , χ 2 , . . . , χ m , respectively, is itself stealthy with a value χ value given by which is the harmonic mean of the χ i divided by m. Proof: Formula (68) is a direct consequence of the fact that χ is inversely proportional to the number density ρ = m i=1 ρ i of the union of the configurations in R d , where ρ i is the number density associated with ith configuration, which is inversely proportional to χ i . This Lemma together with the fact that any periodic crystal with a finite basis is a stealthy ground state can be used to demonstrate rigorously how complex aperiodic patterns can be ground states, entropically favored or not. A sketch of such a proof would involve the consideration of the union of m different periodic structures in R d with densities ρ 1 , ρ 2 , . . . , ρ m , respectively, each of which are randomly translated and oriented with respect to some coordinate system such that m is very large but bounded and ρ i = ρ j for all i and j. It is clear that the resulting configuration will be a highly complex aperiodic structure in R d that tends toward a disordered stealthy pattern with a value of χ that is very small but positive according to relation (68).

VII. PAIR STATISTICS IN THE CANONICAL ENSEMBLE: "PSEUDO" HARD SPHERES IN FOURIER SPACE
The task of formulating an ensemble theory that yields analytical expressions for the pair statistics of stealthy degenerate ground states is highly nontrivial because the dimensionality of the configuration space depends on the density (or χ) and there is a multitude of ways of sampling the ground-state manifold, each with its own probability measure for finding a particular ground-state configuration. Therefore, it is desirable to specialize to equilibrium ensembles with Gibbs measures because the characterization of both the ground states (as well as the corresponding excited states) would be most tractable theoretically. In particular, our objective is to derive analytical formulas for the pair statistics of stealthy disordered ground states for sufficiently small χ in the canonical ensemble as temperature T tends to zero, i.e., the probability of observing a configuration is proportional to exp[−Φ(r N )/(k B T )] in the limit T → 0. We show here that under such circumstances the pair statistics in the thermodynamic limit can be derived under the ansatz that stealthy ground states behave remarkably like "pseudo" equilibrium hard-sphere systems in Fourier space. This ansatz enables us to exploit well-known accurate expressions for the pair statistics in direct space. As will be shown, agreement with computer simulations is excellent for sufficiently small χ.

A. "Pseudo" Hard-Sphere Ansatz
We have already noted that the step-function contribution to ρh(k) for stealthy ground states, denoted byf (k) in relation (52), is identical to the Mayer-f function for an equilibrium hard-sphere system in direct space. This implies that the corresponding contribution to S(k) is a simple hard-core step function Θ(k − K), which can be viewed as an equilibrium hard-sphere system in Fourier space with "spheres" of diameter K in the limit that χ tends to zero. Why? Because such a step-function is exactly the same as the pair correlation g 2 (r = k) of an equilibrium hard-sphere system in direct space in the limit that ρ tends to zero. That the structure factor must have the behavior S(k) → Θ(k − K) in the limit χ → 0 is perfectly reasonable, since a perturbation about the ideal-gas limit [where S(k) = 1 for all k] in which an infinitesimal fraction of the degrees of freedom are constrained should only introduce an infinitesimal change in S(k) of zero inside the exclusion zone (constrained region). We call this the weakly-constrained limit where a step-function S(k) is expected on maximum entropy grounds; it corresponds to the most disordered (decorrelated) form of S(k) subject to the impenetrability condition in Fourier space. We refer to this phenomenon as equilibrated "pseudo" hard spheres in Fourier space because there are actually no points in that space that have a hard-core repulsion like true hard spheres do in direct space.
It is expected on the same maximum entropy grounds that a perturbation expansion about the weaklyconstrained limit χ = 0 will lead to a perturbation expansion in χ for S(k) that can be mapped to the low-density expansion of g 2 (r) for equilibrium hard spheres. More generally, we make the ansatz that in the canonical en-semble as T → 0 this hard-sphere analogy continues to hold as χ is increased from zero to positive values, provided that χ is small enough, implying that the collective coordinate variablesñ(k) (defined in the Introduction) are weakly correlated. Though the pseudo hard-sphere picture must breakdown in some intermediate range of χ, for d = 1 and d = 2, this hard-sphere mapping is again exact when χ = χ max , which corresponds to the maximal value of the packing fraction η in these dimensions (see Tables I and II). This exact correspondence with the maximal value of η when χ = χ max does not hold for d = 3 or d = 4, however. Thus, one should only expect that χ and η are proportional to one another, even at small χ values.
Under the pseudo-hard-sphere ansatz, the direct-space pair correlation function g HS 2 (r; η) of a disordered hardsphere system at a packing fraction η for sufficiently small η can be mapped into the structure factor S(k; χ) for a disordered stealthy ground state derived from the canonical ensemble at fixed χ for sufficiently small χ as follows: As alluded to above, the parameter χ can be viewed as an effective packing fraction for pseudo hard spheres of diameter K in reciprocal space that is proportional to η, i.e., where b(d) is a d-dependent parameter that is to be determined. Let h HS (r) be the total correlation function of a disordered equilibrium hard-sphere system in direct space and let us define for stealthy ground states The ansatz is also defined by the alternative mapping This mapping then enables us to exploit the well-known statistical-mechanical theory of equilibrium hard-sphere systems. In particular, we can employ a generalized Ornstein-Zernike convolution relation that defines the appropriate direct correlation functionC(k), namely, where the symbol ⊗ denotes the convolution operation in R d . Therefore, in direct space, H(r) is given by the relation of the following form For example, for d = 1, Inverting this function yields C(r) = −r sin(r) + (r[sin(r) + cos(r)] − 1)η For d = 2 and d = 3, one can use the Percus-Yevick closure of the Ornstein-Zernike integral equation [49], which is highly accurate for low to intermediate densities along the liquid branch, or when mapped to the stealthy problem, for low to intermediate values of χ. It is noteworthy that the exact low-density expansion of h HS (r) for practical purposes is sufficient to produce accurate estimates ofH(k) = ρh(k) and its counterpart ρh(r) for low to intermediate values of χ or η. In particular, using the mapping (72), we obtain for any dimension d, the following low-χ expansion of ρh(k) : where b(d) is the proportionality constant in (70) and α(k; K) is the scaled intersection volume of two identical d-dimensional spheres of diameter K whose centers are separated by a distance k [cf. (20)] [61]. This formula indicates that S(k) develops a peak value at k = K (over and above the value of unity due to the step function in the limit χ → 0) and then monotonically decreases until k = 2K where its achieves its long-range value of unity for all k > 2K, which we will see is verified by computer simulations. Fourier inversion of (77), division by ρ, and use of (35) yields a corresponding low-χ expansion of the total correlation function h(r) through second-order in χ and hence has an error term of order χ 3 .
To get an idea of the large-r asymptotic behavior of the pair correlations, consider the limit χ → 0 for any d. In this limit, the total correlation function h(r) for any r obtained from (77) is given by which for large r is given asymptotically by (79) Thus, the longed-ranged oscillations of h(r) are controlled by the power law −1/r (d+1)/2 . Equation (78) indicates that in the limits χ → 0 and ρ → ∞, h(r) → 0, and therefore the pair correlation function tends to the ideal gas even though the structure factor [Eq. (69)] cannot tend to the ideal-gas form because of its stealthy property. This is in contrast to the situation considered in Fig. 1, where we take the ρ → ∞ limit by fixing N and letting v F → 0. In that case, both g 2 (r) and S(k) tend to the associated ideal-gas forms, i.e., g 2 (r) = 1 for all r and S(k) = 1 for all k.

B. Comparison of Theoretical Predictions to Simulations
In order to test our theoretical results for the pair statistics of stealthy ground states in the canonical ensemble, we have carried out computer simulations to generate and sample such configurations, the details of which are described in Appendix B. In all cases, we take K = 1, which sets the length scale. Our simulation results reveal that the functional trends for S(k) and g 2 (r) predicted by the ansatz of pseudo hard spheres in Fourier space with an effective packing fraction χ are remarkably accurate for a moderate range of χ about χ = 0. Because it is theoretically highly challenging to ascertain the proportionality constant b(d) in (70) that arises in (77), we must rely on the simulations to guide us in its determination. First, we observe that for d = 1, the mapping between χ and η is one-to-one, i.e., b(1) = 1. Second, the simulation data suggest that to an excellent approximation, b(d) for d ≥ 2 is given by assuming that the peak value of S(k) or ρh(k), achieved at k = K for sufficiently small χ, is invariant with respect to this peak value in the one-dimensional case, and consequently b(d) = [α(K; K)2 d ] −1 . Here K = 1. Figure 6 shows that the structure factor S(k), as obtained from (71) and (77), is in excellent agreement with the corresponding simulated quantities for χ = 0.05, 0.1 and 0.143 for d = 3. In Fig. 7, we compare our theoretical results for the pair correlation function g 2 (r), as obtained by Fourier inversion of (77), to corresponding simulation results across the first three space dimensions. Again, we see excellent agreement between theory and simulations, which validates the pseudo hard-sphere Fourier-space ansatz. Figure 8 depicts our theoretical predictions for g 2 (r) for χ = 0.15 across the first four space dimensions. It is seen that increasing dimensionality increases short-range correlations.  Here K = 1.

C. Translational Order/Disorder Metric
We have seen that both short-and long-scale correlations increases as χ increases. A useful scalar positive order metric that captures the degree to which translational order increases with χ is given by where we have used Parseval's theorem and D is some characteristic length scale [62]. Note that for an ideal gas (spatially uncorrelated Poisson point process), τ = 0 because h(r) = 0 for all r. Thus, a deviation of τ from zero measures translational order with respect to the fully uncorrelated case. Because τ diverges for any perfect crystal, it is a quantity that is better suited to distinguishing the degree of pair correlations in amorphous systems.
In the case of stealthy ground-state configurations, τ is given explicitly by the relation whereH(k) is given by (71). Substitution of the leading order term in the χ-expansion (77) into (81) yields where we have taken D = K −1 . Thus, for stealthy ground states, the order metric τ grows quadratically with χ for small χ. Since the error is of order χ 3 , we expect that this quadratic form will be a very good approximation of τ up to moderately large values of χ. Indeed, this is confirmed by our simulations up to χ = 0.25.

VIII. LOCAL NUMBER VARIANCE FOR STEALTHY HYPERUNIFORM DISORDERED GROUND STATES
Here we investigate theoretically the local number variance for stealthy disordered ground states as a function of χ and then use these results to extract an order metric [27] that describes the extent to which large-scale density fluctuations are suppressed as χ increases in these hyperuniform systems (see Sec. II A). The local number variance σ 2 (R) associated with a general statistically homogeneous and isotropic point process in R d at number density ρ for a spherical window of radius R is determined entirely by pair correlations [27]: where v 1 (R) is the d-dimensional volume of a spherical window [cf. (36)], h(r) is the total correlation function [cf. (5)], α(r; R) is the scaled intersection volume of two spherical windows of radius R, as given by (20), and α(k; R) is the Fourier transform of α(r; R), which is explicitly given by [27] α(k; R) We have already noted that the stealthy ground states considered in the present paper are hyperuniform, i.e., S(k) → 0 as k → 0 (see Sec. II A). This means that such systems obey the sum rule ρ R d h(r)dr = −1 and, because of the rapid manner in which S(k) vanishes in the limit k → 0, the number variance has the following large-R asymptotic behavior [27]: where Λ(R) is a bounded function that oscillates around an average value The scaling (85) occurs for a broader class of hyperuniform systems, as specified by relation (10). The parameter Λ is an order metric that quantifies the extent to which large-scale density fluctuations are suppressed in such hyperuniform systems [27]. To compare different hyperuniform systems, Torquato and Stillinger used the following rescaled order metric: which is independent of the density, where φ = ρv 1 (1/2). Among all hyperuniform point patterns having the scaling (85), B is minimized (greatest suppression of largescale density fluctuations) for the integer, triangular, bcc and D 4 lattices for d = 1, 2, 3 and 4, respectively [27,28]. Using the analytical results for g 2 (r) or S(k) described in the previous section, we have computed relation (83) for σ 2 (R) versus R for selected values of χ across the first three dimensions and compared them to our corresponding simulation results; see Fig. 9. The analytical and numerical results are in excellent agreement with one another. Table V lists the order metric B for various values of χ across the first four dimensions for disordered stealthy ground states, as obtained from the analytical estimates of σ 2 (R) and Eq. (87). These results are also compared to the corresponding optimal values. As expected, B decreases as χ increases for fixed d.

IX. NEAREST-NEIGHBOR FUNCTIONS
Here we obtain theoretical predictions for the nearestneighbor functions of stealthy disordered ground states.
Nearest-neighbor functions describe the probability of finding the nearest point of a point process in R d at some given distance from a reference point in the space. Such statistical quantities are called "void" or "particle" nearest-neighbor functions if the the reference point is an arbitrary point of the space or an actual point of the point process, respectively [63]. Our focus here is on the particle nearest-neighbor functions.
The particle nearest-neighbor probability density function H P (r) is defined such that H P (r)dr gives the probability that the nearest point to arbitrarily chosen point lies at a distance between r and r + dr from this chosen point of the point process. The probability that a sphere of radius r centered at a point is empty of another points, called the exclusion probability E P (r), is the associated complementary cumulative distribution functions and so E P (r) = 1 − r 0 H P (x) dx and and hence H P (r) = −∂E P /∂r.
The nearest-neighbor functions can be expressed as infinite series whose terms are integrals over n-body correlation function defined in Sec. II [63,64]. In general, an exact evaluation of this infinite series is not possible because the g n are not known accurately for n ≥ 3, except for simple cases, such as the Poisson point process. Theoretically, one must either devise approximations or rigorous bounds to estimate nearest neighbor quantities for general models [63,65].
Torquato has given rigorous upper and lower bounds on the so-called canonical n-point correlation function H n for point processes in R d [64,66]. Since nearestneighbor functions are just special cases of H n , then we also have strict bounds on them for such models [63,64]. Here we will employ upper and lower bounds on E P (r) that relies on knowledge of the pair correlation function: where Z(r) is the cumulative coordination number [cf. (4)]. The upper bound (89) was presented in Ref. [43].
These bounds are evaluated for stealthy ground-state configurations using the analytical expression for the pair correlation function given in Sec. VII. Figures 10 and 11 compare these bounds to our numerical results for both χ = 0.05 and χ = 0.1 for d = 2 and d = 3, respectively. We see that the bounds on E P (r) provide the correct qualitative trends as a function of r; the upper bound being the sharper of the two bounds for these cases.  The mean-nearest neighbor distance λ is defined as the first moment of H P (r) or, equivalently, zeroth moment of E P (r), i.e., For an ideal gas (Poisson point process) at number density ρ, the mean nearest neighbor can be explicitly given in any dimension [64]: Using the upper bound (89) and relation (90), we plot in Fig. 12 upper bounds on the mean nearest-neighbor distance λ, scaled by the corresponding ideal-gas quantity obtained from (91), as a function of χ for the first four space dimensions. For fixed χ, the upper bounds on λ/λ Ideal decrease as the space dimension increases, as expected, and tends to unity in the large-d limit consistent with the so-called "decorrelation" principle [43,54].

X. EXCITED STATES: STRUCTURE FACTOR AND THERMAL EXPANSION COEFFICIENT
Here we derive accurate analytical formulas for the structure factor and thermal expansion coefficient for the excited states associated with stealthy ground states at sufficiently small temperatures. We see from the compressibility relation (3) that if the isothermal compressibility κ T is bounded, then S(0) must be zero for any ground state, stealthy or not. Recall that for stealthy ground states κ T is bounded according to (42). Now consider excited states infinitesimally close to the stealthy ground states, i.e., when temperature T is positive and infinitesimally small. Under the highly plausible assumption that the structure of such excited states will be infinitesimally near the ground-state configurations for sufficiently small χ and T , then to an excellent approximation, the pressure is given by where first term is the ideal-gas contribution and the second term is the configurational contribution, which under the stated conditions, is effectively the same as the ground-state expression (41), where we have set k B = v 0 = K = 1. Thus, relation (92) yields the isothermal compressibility to be κ T = [ρ(ρ + T )] −1 , which when substituted into (3) for large ρ (small χ) and small T yields that S(0) varies linearly with T for such excited states: Figure 13 shows that the prediction of relation (93) is in excellent agreement with our MD simulation results (Appendix B) in the case d = 2. It is expected that this positive value of S(0) will be the uniform value of S(k) for 0 ≤ k ≤ K for the special case of the step-function power-law potentialṽ(K) [the case m = 0 in (15)] for small χ. This behavior of S(k) has indeed been verified by our simulation results in various dimensions. For other stealthy-potential function choices, S(k) will no longer be a constant for 0 ≤ k ≤ K. An interesting conclusion to be drawn from this analysis is that S(0) for a system can be made to be arbitrarily close to zero at positive temperature, even if T is itself arbitrarily small. This means that, for all practical purposes, such systems at positive T are effectively hyperuniform. Perfect hyperuniformity is not necessarily required in order to achieve novel physical properties in technological applications.
Using the cyclic identity (∂ρ/∂T ) p (∂T /∂p) ρ (∂p/∂ρ) T = −1 and approximation (92), it immediately follows that the thermal expansion coefficient α ≡ −ρ −1 (∂ρ/∂T ) p , for sufficiently small χ and T is given by We see that the thermal expansion is positive under such conditions, which is to be contrasted with the anomalous negative thermal expansion behavior for sufficiently large χ over a low temperature range demonstrated in our earlier numerical work [15,16].

XI. CONCLUSIONS AND DISCUSSION
Stealthy hyperuniform disordered ground states in R d are infinitely degenerate and arise from a class of bounded long-ranged pair potentials with compactly-supported Fourier transforms. Such exotic many-particle states of matter were previously studied only numerically. Because the configurational dimensionality depends on the density (or χ), a highly unusual situation, and there are an infinite number of distinct ways to sample the groundstate manifold, each with its own probability measure, it has been theoretically very challenging to devise predictive ensemble theories. A new type of statisticalmechanical theory needed to be be invented. This paper has initiated such a theoretical program.
Specifically, we have derived general exact relations for the ground-state energy, pressure and isothermal compressibility that apply to any ensemble as a function of the number density ρ in any dimension d. It was demonstrated how disordered degenerate ground states can arise as part of the ground-state manifold. We also obtained exact integral conditions that both the pair correlation function g 2 (r) and structure factor S(k) must satisfy in any ensemble. Then we specialized our results to the canonical ensemble in the zero-temperature limit by exploiting an ansatz that stealthy states behave like "pseudo" equilibrium hard-sphere systems in Fourier space. The resulting theoretical predictions for g 2 (r) and S(k) were shown to be in excellent agreement with computer simulations across the first three space dimensions for sufficiently small χ. These results were used to theoretically obtain order metrics, local number variance and nearest-neighbor functions across dimensions. We also derived accurate analytical formulas for the structure factor and thermal expansion coefficient for the excited states associated with stealthy ground states at sufficiently small temperatures. Our analyses provide new insights on our fundamental understanding of the nature and formation of low-temperature states of amorphous matter. Our work also offers challenges to experimentalists to synthesize stealthy ground states at the molecular level, perhaps with polymers, as suggested in Sec. III.
There are many remaining open theoretical problems. While the pseudo hard-sphere system picture for the canonical ensemble is almost surely exact in the limit χ → 0, a future challenge would be to provide rigorous justification for this picture for positive but small χ. One possible avenue that could be pursued is the formulation of an exact perturbation theory for the pair statistics about the weakly-constrained limit (χ → 0). As noted in the Introduction, while the configuration space is fully connected for sufficiently small χ, quantifying its topology as a function of χ up to χ * max is an outstanding open problem. At some intermediate range of χ, the topology of ground-state manifold undergoes a sequence of one or more disconnection events, but this process is poorly understood and demands future study. In the limit χ → χ * max , the disconnection becomes complete at the unique crystal ground state [67]. The simple constraint/degrees-of-freedom counting arguments described earlier lead to definite predictions for the entropically favored stealthy ground states derived form the canonical ensemble in the limit T → 0. For χ between 0 and 1/2, the ground states are disordered and possess a configurational dimension per particle d C of d(1 − 2χ); see Eq. (37). At χ = 1/2, the configurational dimensionality per particle collapses to zero and there is a concomitant phase transition to a crystal phase. The fact that the pressure is a continuous function of density [cf. (41)] for all χ up to χ * max (see Tables I-IV) implies that any phase transition could be continuous. While this eliminates a first-order phase transition in which the phase densities are unequal, it does not prohibit such a phase transition in which the two distinct phases possess the same number density. At χ = 0.5, our numerical evidence indicates that the structurally distinct fluid and crystal phases have equal free energies. This does not necessarily imply that the two phases could coexist side by side within the system separated by an interface. This phase diagram is depicted in Fig. 14, which applies to the first four space dimensions. For d = 1, the only crystal phase allowed is the integer lattice (see Sec. VI), and hence there can be no phase coexistence. For d = 2, our simulations [14][15][16]68] indicate that the crystal phase is the triangular lattice for 1/2 ≤ χ ≤ χ * max . However, for d = 3, it is possible that there may be more than one crystal phase. For example, while for sufficiently high χ up to χ * max , we expect the stable crystal to be the bcc lattice, our simulations cannot eliminate the possibility that the fcc lattice is a stable phase for some χ in the range 1/2 ≤ χ ≤ 0.9068 . . .; for χ > 0.9068 . . ., the fcc lattice cannot be ground state (see Table III). Since four dimensions is more similar to two dimensions in that the lattice corresponding to χ * max is equivalent to its dual, we would expect that the D 4 lattice is the stable crystal for 1/2 ≤ χ ≤ χ * max , but this remains to be confirmed. All of our previous and current simulations for the first three space dimensions [12][13][14][15][16] strongly suggest that all of the energy minima attained were global ones for χ < 0.5, but when χ > 0.5, the topography of the energy landscape suddenly exhibits local minima above the groundstate energies. The possible configurations that can arise as part of the ground-state manifold for χ > 1/2, regardless of their probability of occurrence, not only include periodic crystals for d ≥ 2, as discussed in Sec. VI, but aperiodic structures, reflecting the complex nature of the energy landscape. For example, for d = 2, the mani-fold includes aperiodic "wavy crystals," which have been shown to arise via numerical energy minimizations from random initial conditions with high probability in a range of χ where the triangular lattice is entropically favored [12,14]. A deeper understanding of such aspects of the ground-state manifold would undoubtedly shed light on the topography of the energy landscape.