Localization, Disorder, and Entropy in a Coarse-Grained Model of the Amorphous Solid

We study the role of disorder in producing the metastable states in which the extent of mass localization is intermediate between that of a liquid and a crystal with long-range order. We estimate the corresponding entropy with the coarse-grained description of a many-particle system used in the classical density functional model. We demonstrate that intermediate localization of the particles results in a change of the entropy from what is obtained from a microscopic approach using for sharply localized vibrational modes following a Debye distribution. An additional contribution is included in the density of vibrational states g(ω) to account for this excess entropy. A corresponding peak in g(ω)/ω2 vs. frequency ω matches the characteristic boson peak seen in amorphous solids. In the present work, we also compare the shear modulus for the inhomogeneous solid having localized density profiles with the corresponding elastic response for the uniform liquid in the limit of high frequencies.


Introduction
In isotropic fluids, the constituent particles move randomly, with the corresponding time-averaged density being constant at all points. This random movement in space is absent at lower temperatures or higher densities, and the crystalline state forms. The latter has a characteristic inhomogeneous density ρ(x) having the symmetry of the corresponding crystal. The particles vibrate around the respective lattice points. Even in computer simulations of a system of hard spheres, it has been seen [1] that a fluid to crystal transition occurs. The hard-sphere fluid transforms into an inhomogeneous state with the particles vibrating around the points of an fcc lattice. This transformation is a result of the competition between the energetic and entropic contributions to the free energy and occurs at the packing ρ 0 σ 3 = 0.502, which is much lower than the close pack fcc structure with ϕ = 0 .740. The classical density functional theory (DFT) presents a theoretical model [2] for understanding the freezing transition by treating the coarse-grained one-particle density ρ(x) as an order parameter [3]. A typical parametrization of the inhomogeneous density function for the crystalline state involves superposition of Gaussian density profiles centered on a lattice with the long-range order of the crystalline state [4].
where the {R i } denotes the underlying lattice. The function φ is taken as the isotropic Gaussian φ(r) = ( α π ) 3 2 e −αr 2 . The parameter α characterizes the inverse width or spread of the Gaussian density profile and represents the degree of mass localization in the system. The homogeneous liquid state density is described by Gaussian profiles of very large widths in the limit α→0. Metastable states with aperiodic structures [5][6][7] have also been studied using DFT. They have free energies intermediate between a crystal and a homogeneous liquid state. These metastable states consist of localized density profiles around the respective points {R i } of an amorphous lattice. Such local minima of the free energy functional exist for hard-sphere interactions and hard-core Hertzian potentials. The thermodynamic properties of the inhomogeneous amorphous solid are computed, assuming the system to be in a single phase. The amorphous nature of the underlying structures plays a crucial role in producing the metastable free-energy minima. Breaking of the isotropic symmetry of the liquid state leads to the development of the transverse sound modes in the crystal which are of the Goldstone modes [8]. Transverse sound modes occur in the amorphous solid state as well. However, the vibrational states of glassy materials are beyond the simple plane-wave phonon picture for crystals.
The crystalline and the amorphous (termed as glass) states of matter have a characteristic high degree of mass localization compared to that of the liquid state. In the present context of DFT, this refers to free-energy minima for the system corresponding nonzero values of the parameter α. The localized particles in the disordered glassy state are on a metastable amorphous lattice structure without any long-range order. The individual particles vibrate around the sites, which remain localized on a disordered lattice over long times. For the amorphous solid, density of vibrational modes is modified from the predictions of the Debye distribution. The excess in the density of states (DOS) g(ω) for a disordered system is at the THz frequency (ω) range and appears as a peak in the reduced DOS representation as in g(ω)/ω 2 vs. ω plot. This peak is analogous to the so-called boson peak, whose height decreases and location shifts towards the higher frequencies with the increase of pressure or density. Inelastic scattering studies of light and neutron from amorphous solids [9,10] has established the presence of the boson peak.
Several authors have studied various models for understanding the boson peak formation. These include analysis based on (a) potential energy landscape description to describe a transition from a saddle-point-dominated phase without phonons to a minima-dominated phase with low-energy phonons [11], (b) localized phonons in cohesive clusters [12,13], (c) an-harmonic interaction potentials [14] for a cluster of particles, (d) domain walls between configurations forming a mosaic of disordered states [15], (e) descriptions in terms of crystalline lattices of springs with randomly assigned stiffness [16], (f) randomly fluctuating elastic constants [17,18], and (g) the disordered network solids with fixed bond connectivity [19]. The model for a disordered solid is studied [16] in terms of a geometrically perfect crystal with random interactions between nearest neighbors or in terms of a crystal having particles with different masses. A common aspect of the boson peak models is that they are based on the existence of localized modes in the amorphous solid [20] and are a manifestation of the disorder. The system's free energy reaches a minimum in the equilibrium state with contributions from the vibrational modes in the solid.
In the present work, we demonstrate the role of disorder in producing the metastable states with a degree of mass localization, intermediate between a liquid and a crystal with long-range order. Using coarse-grained description in terms of the density field, we obtain that localization of the particle over intermediate length scales makes a change in the system's entropy. This modified entropy is obtained by including the excess density of states in the form of a boson peak. For the inhomogeneous states using the standard formulation of DFT, thermodynamic properties like pressure and elastic constants like shear modulus are also obtained. For this, we compute the first and second derivatives of the free energy with respect to the average density. In the present work, we also consider how the solid-like properties, like the shear modulus for the inhomogeneous disordered states (with mass localization), change compared to the similar property for the uniform liquid in the high-frequency limit. We organize the paper as follows. In Section 2, we discuss the calculation of entropy of the hard-sphere system using a continuum model in terms of the coarse-grained density. We consider a model for the solid-like system in terms of vibrational modes and a choice of the appropriate density of states g(ω). In the next section, we discuss the calculation of the total free energy taking into account the role of the interaction and identify the amorphous metastable states with inhomogeneous density ρ(x). We analyze the appropriate density of states for the solid-like state with a low degree of mass localization (compared to a crystal) and identify the characteristic boson peak. In Section 4, we calculate the elastic constants for the amorphous solid-like state by extending the DFT to include the weighted density approximation. The paper ends with a discussion of the key points.

The Coarse Grained Model
Let us consider a system of N identical particles of mass m and the position and momentum coordinates {r α , p α }, for α = 1, ..N. The whole set of phase space variables are to be denoted as {r N , p N }. At the microscopic level, the one particle densityρ(x) is formally defined as a function of the phase space variables (r N , p N ), The microscopic Hamiltonian for the system of N particles is given by where H 0 N is the intrinsic part of the Hamiltonian N particle system, and U(r 1 , ...., r N ) is the interaction energy of the N particle system. U ext is the external field contribution to the Hamiltonian and is obtained as an interaction term with a local field arising from a one-body potential φ, The coarse-grained density function corresponding to the distribution function f is then obtained as ρ(x) = Trρ(x) f , where the operation "Tr" refers to statistical mechanical average. For the equilibrium density, we choose the standard grand canonical ensemble. The free energy of the liquid is obtained as a functional of the coarse-grained density function ρ(x) as a sum of two parts-the ideal gas term and the interaction term: The ideal gas part F id is obtained exactly by setting the interaction potential U = 0 in the Hamiltonian. In this case, Hamiltonian H o N of the N particle system only has the kinetic energy term K = ∑ α p 2 α /2m and thus explicitly integrating out the 3N momentum variables, the grand canonical partition function Ξ is obtained as where u(x) = µ − φ(x), in terms of the external field φ and Λ 0 = h/ √ 2πmk B T, is the thermal de Broglie wavelength for the liquid particles. Using the expression (6), the coarse grained density is obtained for noninteracting case in the exact form, where βu(x) = ln[Λ 3 0 ρ(x)] . The free energy functional F[ρ(x)] for the noninteracting system is obtained in terms of the coarse grained density by generalization of the corresponding thermodynamic relation, We have used the ideal gas equation of state in writing the second term on the right-hand side. Therefore, using the expression (7) we obtain For the noninteracting system of N particle, the partition function is and thus the ideal gas free energy is obtained from the logarithm of Z N for the equilibrium state. Therefore, βF id = Vρ 0 (ln(ρ 0 Λ 3 0 ) − 1) for the uniform density fluid ρ 0 , which is the limit ρ(x)→ρ 0 in Equation (9). From Equation (9), it is easy to see that the entropy drops upon localization of the particles due to the restriction of available phase space. The ideal gas contribution for the N particles system changes by an amount ∆F id when ρ 0 →ρ(x): whereρ(r) = ρ(r)/ρ 0 . As by definition the quantityρ(r) is always positive, using the Gibbs inequality {x ln x − x + 1} ≥ 0 for positive x, it follows that ∆F id ≥ 0. The relation (7) between the chemical potential and the density is inverted here, giving an exact expression for the free energy functional. In this regard, it is essential to note that the corresponding free energy for the system, expressed as a function of the density ρ(x), is exact. It is, however, possible only for the so-called ideal-gas contribution for the noninteracting system. The role of interactions are accounted for by introducing the direct correlation functions c (i) 's for i = 1,2,3, which are defined as successive order functional derivative of an excess contribution F ex to the free energy due to interactions. We will be discussing the interacting system in the next section. Two length scales are inherent in the coarse-grained picture for the solid presented above. The interaction potential between the constituent particles has the characteristic length scale, which is the hard-sphere diameter σ in the present case. The other length is = 1/ √ α used in parametrization of the inhomogeneous density function (1). In the coarse-grained description, the scale signifies the degree of localization of mass in the system, and the limit ασ 2 → 0 or >> σ represents density profiles of uniform liquid with an average ρ 0 = N/V. In this case, the sum in the RHS of Equation (1) has matching contributions from a large number of terms, each of which corresponds to a Gaussian profile centered on a corresponding site ∈ {R i }. On the other hand, the limit ασ 2 >> 1 or << σ corresponds to sharply localized density profiles like that in a crystal. The ratio /σ is generally referred to as the Lindeman parameter for the crystal. For ασ 2 ≥ 50, the density ρ(r) is constructed in terms of the contributions from non-overlapping density profiles and Equation (9) is well approximated by its asymptotic value for large α, For lower values of α where the overlapping of the Gaussian profiles from different sites are relevant in the sum in Equation (1), f id is evaluated numerically from the integral The free energy in the right-hand side of (12) is averaged over a random set of Gaussian centers {R i } for the density profiles, and is obtained in terms of the site-site correlation function w(R). However, it is helpful to note here that even if the Gaussian centers are on a regular crystalline lattice, the result for the non-interacting system is not changed. It follows directly by expressing the inhomogeneous density function ρ(r) in terms of reciprocal lattice vector (RLV) expansion: where {G i } denotes the RLV [21], and for the fcc structure the solid density ρ s = 4/a 3 0 in terms of the lattice constant a 0 . The RLV expansion for the density reduces to the form (1) when we identify For smaller values of α, the sum in the expansion for density has to include contributions from increasingly the larger number of RLVs. Figure 1 shows the variation of the total entropy S (in units of k B ) obtained by evaluating the ideal gas part of the free energy f id from the expression (9) using two different approximations for the inhomogeneous density ρ(x). First, in Equation (1), ρ is described with in terms of Gaussian density profiles centered on a random lattice {R i } in real space. Second, in Equation (13), density is expressed in terms of an expansion in reciprocal wave vectors {G} corresponding to a regular lattice with long range order. We use the fcc structure going up to Gσ = 14. The ρ(x) in the respective cases are used in expression (9) to obtain the f id . These results shown in the main panel correspond to a hard-sphere system with packing-fraction ϕ = πρ 0 σ 3 /6 = 0.576. The thermal wavelength is kept a constant at Λ = 0.025σ. For a crossover value of α≤α 0 , the entropic contribution is evaluated numerically from the integral given in Equation (12). The pair function w(R) is approximated here with pair function for the Bernal's random structure [22] and is generated through Bennett's algorithm [23]. We have obtained w(R) using the following relation: where ϕ denotes the average packing fraction and ϕ 0 is a scaling parameter for the structure [24] such that at ϕ = ϕ 0 , Bernal's structure is obtained. The important thing to note here that for α < α 0 , the free energy curve (shown as the solid line) deviates completely from the asymptotic result. This line extends to the correct α→0 limit (ln{ρ 0 Λ 3 } − 1) corresponding to the uniform liquid state. In the DFT model with coarse-grained density ρ(x) expressed in terms of Gaussian profiles of inverse width α, a corresponding α 0 is identified such that for all α<α 0 the asymptotic formula (11) deviates from the entropy curve (shown in Figure 1) of the coarse-grained model. The inset of Figure 1 shows the weak dependence of α 0 on density ρ 0 . Note, here, that the modified entropic contribution discussed above, corresponding to small values of α or more spread out density distributions, be it around a random structure or a regular set of lattice points, it is a manifestation of the delocalization of the particles and in the α→0 it goes to the case of an isotropic liquid. . Arrow marks location of α 0 , such that for α < α 0 the asymptotic formula deviates from the exact result. Inset shows α 0 ( in units of σ −2 ) vs. density ρ 0 σ 3 for the hard sphere system.

Microscopic Description
As an alternative to the coarse-grained model, evaluating the corresponding partition function (for a microscopic model) obtains the free energy for the system. The microscopic description is introduced here to focus on the density of vibrational modes in the amorphous state with inhomogeneous mass distribution. The description involving the vibrational modes is needed to study the corresponding modifications of the density of states in the metastable liquid. Though no long-range symmetry breaking occurs, the glassy state still supports transverse sound waves. These transverse modes are described in terms of vibrational excitations over an amorphous structure constituting the phonon modes. The respective free energies obtained from the continuum and microscopic approaches closely agree in the form of parameters of the (density) field-theoretic model that are identified with the density of states and characteristic cut-off frequency of the vibrational modes. The set of vibrational modes for the density of states g(ω)dω between frequencies ω and ω + dω contributes to this. For an average ρ 0 number of particles per unit volume, we have a constraint on g maintaining the total number of vibrational modes in three dimensions, The cut-off frequency ω m is a characteristic property of the material and sets the shortest time scale for the vibrational modes. For fixed ρ 0 , the upper cut-off ω m [g] in the integral in Equation (15) is treated as a functional of the density of states g(ω). In the constant NVT ensemble, the partition function for 3N harmonic oscillators is where ω i is the frequency of the i-th vibrational mode in the system. The free energy of the non-interacting system is obtained in units of β −1 aŝ whereĒ denotes the average kinetic energy of the N particles and scaled energy ε m = βhω m . We write the free energy asf with a hat to indicate that it is obtained from microscopic considerations. The density of states g(ω) is written in Equation (17) in a scaled formḡ(x), as a function of the variable . In this case,ḡ(x) is only nonzero for 0 ≤ x ≤ 1. The function κ(y) in the right hand side of Equation (17) is obtained as Please see Appendix A for details of obtaining the above result for the function κ(y). The function κ(ε m x) determines the functional dependence off [g] on the scaled density of statesḡ as given in Equation (17). A link between the density of states g(ω) and the mass localization parameter α, is obtained by equating the results (9) and (17) for the ideal gas part of the free energy following respectively from the coarse-grained model and the microscopic calculation.

Debye Distribution
If the density of states is the Debye distribution, we obtain where the upper cut-off ω m is the Debye frequency ω D = (6πρ 0 ) 1/3 c, in terms of the speed of sound c. The corresponding scaled distributionḡ D is obtained as a function of the scaled variableḡ D (x) = 3x 2 where x = ω/ω D . In terms of the reduced variables, we obtain ε m ≡ε D =βhω D . In the limit, ε D << 1, when the density of states g(ω) is the Debye distribution g≡g D in Equation (17), we have to the leading orders in ε D : The result (20) for the free energy per particle calculated from the partition function of the microscopic model is identical to the asymptotic formula (11) where we obtain the (reduced) Debye frequency as By comparing numerical evaluation of the integral on the right hand side of Equation (17) (for the choiceḡ≡ḡ D ), with the asymptotic result (20) for different values of ε D , a maximum value ε max D of ε D is identified. It is ε D up to which the two results agree. Figure 2 shows this comparison. The point at which the two curves (approximate and exact) as given by Equations (20) and (17), respectively, deviate is marked with an arrow at ε max D = 2.5 in the main panel. In the inset of the same figure, we show the corresponding agreement between the two formulas for entropy (calculated from DFT) shown in Figure 1. The arrow shown in this inset corresponds to ε D = 0.069, i.e., using the relation (21), it follows that (for Λ/σ = 0.025) the corresponding α = 24 ≈ α 0 . From this, it is clear that over the entire range of α values in the DFT model, the Formula (11) is in agreement with the results of the corresponding microscopic model with a Debye density of states. However, as the actual entropy, S, deviates from the asymptotic Formula (11) for α<α 0 , the corresponding density of states g(ω) in the microscopic model gets modified from the Debye form.  (17) and the asymptotic form (20). The arrow indicates the value ε D = 2.5 at which the asymptotic result separates from the exact value. Inset focus on the part of the curve for the α range given in Figure 1 and arrow indicates the location of α 0 in terms of ε D . The width parameter α is related to ε D through Equation (21).

Excess Density of States
To summarize the previous section, the free energy curve using DFT for α ≥ α 0 follows the asymptotic form (11), and this is identical with the form (20) obtained in the microscopic framework using the Debye distribution g D (ω) as density of states g(ω). However, for α < α 0 , the departure of the entropy curve from the asymptotic form (11) as shown in Figure 1 requires modifying in the microscopic model the corresponding density of states for g(ω) from the Debye form g D (ω). The normalization condition (15) for the modified g(ω) is maintained with a corresponding upper cut-off in frequency ω m or equivalently ε m (=βhω m ) which is different from that of the Debye distribution. The upper cut-off ω m (and thus ε m ) will depend on the form of function g(ω) as well as the value of width parameter α(<α 0 ), as the corresponding entropy, calculated in the coarse grained model has to match with the microscopic results (17)-(18) with the modified g(ω). We require that the function ε m (α) → ε D (α) (of Equation (21)) as α → α 0 . The correction part in this functional relation is denoted by C(α), such that ε m (α) = ε D (α 0 ) + C(α) with C(α)→0 as α→α 0 . For the scaled distribution functionḡ(x) we make a modification over the corresponding Debye formḡ D (x) = 3x 2 in terms of the function ∆, such thatḡ(x, x m (α)) = 3x 2 (1 + ∆(x, α)). As α→α 0 , we must have g →ḡ D for all frequencies x, and we express ∆ with the separation of variables: ∆(x, α) = B(α)∆(x), where B(α)→0 as α→α 0 . In order to maintain the positivity of the density of states, we also require that the function∆(x) is in the range ±1. The x dependence of the∆(x) is chosen to have an intermediate peak over the whole frequency range, and the position of the peak of∆(x) is kept fixed at x = 0.22. In Figure 3, we show the results obtained for the appropriate density of statesḡ(x) which reproduce the f id for α < α 0 . We study a hard sphere system at a fixed density ϕ = 0.576. The reduced density of states g(ω)/ω 2 vs. ω/ω D is shown in Figure 3 for five different values of the localization parameter = 1/ √ α. The variation of the corresponding upper cutoff ω m with width parameter α is shown in the main panel of Figure 4. The peak frequency ω m is smaller than the Debye frequency ω D for α < α 0 and approaches ω D as α→α 0 .

The Free Energy Landscape
In the density functional theory, for identification of the equilibrium state of the fluid, a proper thermodynamic potential or free energy is minimized with respect to the density ρ(x). Free energy is expressed as a sum of two parts which are, respectively, the free energy of the non-interacting system F id and the contribution F ex due to interactions between the particles. The part F id we have already discussed in details above. The interaction part F ex of the free energy is obtained in a functional Taylor series expansion around the corresponding free energy for the uniform density state in terms of the density fluctuation δρ(x) = ρ(x) − ρ 0 . For the N particle system we obtain, the Ramakrishnan-Yussouff (RY) functional, At low densities, the state with spatially uniform mass distribution has the lowest free energy. However, at higher densities, the crystalline state with highly localized density profiles centered around the points of a lattice with long-range order has lower free energy and represents the equilibrium state. Metastable minima corresponding to the above free energy functional, intermediate between the isotropic liquid and the crystalline state, have been identified in several extensions of traditional DFT methods [5,6]. Those metastable states have characteristic ρ(r) represented in terms of Gaussian profiles centered on the points of an amorphous lattice. However, the inverse width parameter α for these metastable states is generally lower than those for the sharply peaked density profiles depicting the crystalline state. A qualitatively different free energy minimum, with the optimum density function characterized by much smaller α values, has also been obtained subsequently [7]. The low α minimum appears at α ≈ 18 for ρ 0 = 1.12. We have chosen here to describe the amorphous lattice {R i } in terms of the Bernel structure. However, metastable minima with the less localized density profiles also occur for other choices of the random structure {R i } [25,26]. The expansion (22) for the free energy of the liquid works as a better approximation when the coarse-grained density corresponds to a low degree of mass localization compared to that for the crystalline state with sharp density profiles.
The role of the interaction term enters here through the direct correlation function c(r). The local minima of the free energy, signifying a low degree of mass localization, have been obtained for various interaction potentials. These include hard sphere interaction [7,27], Lennard-Jones [28], soft-sphere interaction [29], and Hertzian potentials [30]. A typical result obtained for the hard-sphere system in which the solution of the Percus-Yevick equation [31] has been used to compute the direct correlation function c(r), is displayed in Figure 5. In this case, minimization of the total free energy shows that disorder plays a vital role in producing (in the free energy landscape) the local minima, which signify the metastable states [32] in the partially delocalized region (α ≤ α 0 ). The modified entropy is matched with changed density of vibrational states, and is apparent from the fact that the interaction part of the free energy F ex is very different in the respective cases of ordered and disordered structures. If the Gaussian profiles for ρ(r) have their respective centers on a lattice with long-range order similar to the crystalline state, no metastable minimum for the total free energy is obtained in the small α (< α 0 ) region. Figure 6 shows this case with no minima appearing in the low α region. These low α or partially localized state exists only if F ex is computed with the lattice points {R i } on a random structure. Thus, disorder is essential for the metastable minimum of the free energy in the delocalized region (α ≤ α 0 ). For disordered systems with weaker mass localization than the crystal, the corresponding density of states g(ω) differs from the Debye distribution. This modified density (over Debye) produces the entropy contribution appropriate for the delocalized amorphous state. This density of states is consistent with a boson peak seen in amorphous solids.
For a specific interaction potential, density (ϕ) dependence of properties like peakheight (B h ) and peak-position (ω p ) of g(ω) is determined in terms of the optimum α min . In Figure 7, we show for the hard sphere potential, the dependence of α min on the packing fraction ϕ. Using this dependence, we obtain the corresponding boson peak curve (see Figure 3, for example). In each case, the properties B h and ω p for the distribution curve at the α min corresponding to packing fraction ϕ are obtained. Linking of α min to ϕ is done using Figure 7. The density (ϕ) dependence of height and position of the boson peak are displayed, respectively, in Figures 8 and 9. Insets of Figures 8 and 9, respectively, display the corresponding results from experiments. The trends seen from the experimental data are the same as the theoretical model shown in the respective main panels. With the increase of density, the boson peak height decreases, and the peak shifts to higher frequencies. Thus accounting for the entropy for delocalized states by modifying (in the microscopic model) the corresponding density of states in the form of boson peak agree with the experimental data for amorphous metastable systems. Further studies with various interaction potentials are needed to better understand this link between microscopic and coarse-grained models.   . Results for properties of the modified density of states g(ω), e.g., shown in Figure 3 for the (amorphous) metastable hard sphere system with a low degree of mass localization. Main panel: boson peak frequency ω p in units of c/σ (c is the sound speed) vs. packing fraction ϕ. Inset: experimental data [33] boson peak height ω p (in units of 3.2 meV) vs. density ρ (in units of 3.66 gm/cm 3 ).

Elasticity of the Localized State
The elastic constants for the amorphous metastable state is an important property to understand its solid-like nature. In the DFT, equilibrium free energy of the inhomogeneous state is obtained as a function of the average density ρ 0 . The bulk modulus K of the isotropic solid is obtained in terms of the second derivatives of the free energy [24] with respect to ρ 0 .
The pressure in the solid is also obtained from the first derivative of the free energy as We use the modified weighted density functional approximation (MWDA) [34] for calculating the free energy of the inhomogeneous liquid in the metastable state. This is an effective medium approach in which the nonuniform solid is mapped to an equivalent homogeneous liquid of lower density. In calculating the excess part per particle f ex = F ex /N using MWDA [35][36][37][38] in the canonical ensemble, a self-consistent integral equation [34] is obtained for the density of the effective liquidρ. The corresponding packing fraction ϕ = πρσ 3 /6 in terms of the suitably chosen free energy function f ex (φ). (25) where the integral I is defined as I is evaluated using the density in the parametric form (1). The latter involves the set of points {R i } at which the Gaussian density profiles are respectively centered. Averaging over the different choices of this amorphous lattice, we express the final result in terms of a site-site pair correlation function w(r) for the lattice points.
The single and double primes over f ex (x), in the above Equation (25), respectively, denote the first and second derivatives of the function with respect to the argument x. For solving Equation (25), the free energy f ex (ϕ) is taken from the standard expression of excess free energy of a hard-sphere system [39], in the form The Percus-Yevick solution for the hard sphere system [31,40] is used for the direct correlation function c(r) in the integral Equation (25). For a fixed value of ρ 0 , the total free energy (sum of the respective ideal and excess parts of the free energy) is obtained over a range of the width parameter α values by solving the MWDA equation in each specific case. Metastable amorphous states, distinct from the uniform liquid state, are identified by locating the intermediate minima of the corresponding free energy with respect to the mass localization parameter α at α = α min for different values of the packing fraction ϕ. The quantity = 1/ √ α min is the localization parameter scaled with respect to σ. With increasing ρ 0 or ϕ, the particles are more localized, and thus the amplitudes of vibration of the particles around their respective mean position fall.
A class of minima corresponding to heterogeneous structures characterized by weak mass localization for low values of α is detected [27]. Close to freezing, these delocalized structures are more stable than the highly localized "hard-sphere glass". However, at high densities, or packing fraction ϕ > 0.500, the highly localized states corresponding to large values of α (which signifying strong localization) become more stable. In Figure 10, we show the free energy minima for different ϕ at 0.617, 0.581, and 0.554. As the packing fraction increases, the curvatures of the free energy plots with respect to the width parameter α keeps changing. The solid-like behavior of these amorphous states with inhomogeneous density distributions is manifested in the corresponding elastic constants. The elastic constants are calculated by analyzing the nature of the local free-energy minima. We use the formulas (23) and (24) to compute the pressure and the bulk modulus by computing first and second derivatives of the free energy at the two respective minima shown in Figure 10. Using the K and P, in the Cauchy relation, For the isotropic solid, we obtain the corresponding shear modulus G. These DFT results for K, P, and G are, respectively, shown in the main panel as well as in the two insets of Figure 11. We show the results for the two type minima depicted in Figure 9. For low packing fractions (ϕ < 0.580), the less localized (low α) state has higher values for the elastic constants. For high packing fractions (ϕ > 0.580), the corresponding elastic constants are higher for the sharply localized (high α) state. The corresponding pressure, calculated from the same DFT results, however, does not show any cross over in trend like the elastic coefficients and is always lower for the less localized state, i.e., the pressure for the less localized state is always more than that for the sharply localized state. This is shown in the inset of Figure 10. For the elastic constants, qualitative changes in relative behavior between high and low alpha minima occur with increasing packing fraction due to subtle changes in the curvatures of the corresponding free energy curves (shown in Figure 9). Even a uniform liquid behaves like a solid over very short time scales, or equivalently in the high-frequency limit [41]. We studied this short-time elastic response of the uniform hard-sphere liquid in terms of the high-frequency elastic constants. The high-frequency elastic constants for a many-particle system with pairwise interaction are expressed in the well-known Mountain-Zwangig formulas [42,43]. However, these formulas do not go to a finite limit (at a fixed density) for the pure hard-sphere interaction. Therefore, to calculate the high-frequency elastic constant, we use for the elastic moduli of the hard-sphere system, formulas obtained through an analysis of the stress tensor [44]. For local equilibrium, the stress tensor is expanded in terms of strains [45] assuming only instantaneous binary collision between the hard spheres. Three particle and higher-order collisions are being ignored here. From the long-wavelength expansion of the stress tensor, the corresponding shear modulus, G ∞ , and bulk modulus, K ∞ , are obtained (in units of ρ 0 k B T) in the form [46], where P/(ρ 0 k B T) = Z(ϕ). In reaching the above relations we have used for thermodynamics pressure P for the hard-sphere system the Carnahan-Starling approximation [39] as, We assume that the three quantities K, G, and P are related through the Cauchy linear relation (29) as was the case for the DFT model. Therefore, the bulk modulus K is obtained using the same linear relation. Results for these short time or high-frequency quantities are shown with dashed lines in Figure 12. We express these three thermodynamic quantities in units of ρ 0 k B T where T is the temperature. Main panel: high-frequency limit of the shear modulus G ∞ (solid) and bulk modulus K ∞ (dashed) for the uniform density hard sphere system vs. packing fraction ϕ. Inset: pressure P for the uniform hard sphere liquid vs. packing fraction ϕ. Here G ∞ , K ∞ , and P (all expressed in units of ρ 0 k B T) are respectively obtained [44,46] using formulas (30)- (32) in the text.

Discussion
In the present work, we obtained the entropy of the metastable hard-sphere liquid with inhomogeneous (coarse-grained) density, a (functional) Taylor series expansion of the free energy functional around the homogeneous state used in classical DFT models. The result from coarse-grained density functional model is matched with that from the microscopic model by using a modified distribution for the vibrational modes in the amorphous solid-like state. By choosing the modified density of states in the form of a boson peak, seen in amorphous solid-like states, we estimated the characteristic properties, like the height of the peak or position of the peak frequency. The theoretical predictions obtained from the model agree with the corresponding trends seen in experiments. With increasing density, the height B h of the peak decreases, while the position ω p of the peak shifts to higher frequencies.
An essential characteristic of the boson peak is that it gets weaker with increasing fragility, which is related to the long-time relaxation behavior of the glass-forming material. This observed behavior follows naturally in the present model. The stronger the liquid is, the more it displays an increasing tendency to form network structures. The density profiles are more localized for the stronger glasses. Therefore the increase of fragility is synonymous with decrease of α, i.e., in the coarse-grained DFT models discussed here. Due to the fact that, as discussed above, the boson peak height decreases with decreasing α, it also implies a weaker boson peak for more fragile systems [47,48].
The coarse-grained model of DFT considered here can be linked to more traditional models of disordered solids built in terms of springs. The phenomena of boson peak in the disordered system are also studied in terms of a geometrically perfect crystal having random interactions between the neighbors [16]. For large values of α in the DFT model, the sharply localized density profiles are non-overlapping and are interpreted as individual harmonic oscillators having spring constant κ proportional to α [49,50]. Including fluctuations of α at different sites on the amorphous structure will be a natural extension of the present model to make it more appropriate to describe the heterogeneous glassy state.
For calculating free energy, we have used the MWDA of DFT to account for the inhomogeneous density distributions of the amorphous solid with purely hard-sphere interaction. The weighted density approximation [51], keeping up to second-order density fluctuations, accurately maps a purely repulsive hardcore system in terms of an equivalent low-density fluid. In this regard, it is useful to note that the hard-sphere solid is somewhat anomalous in the usual descriptions of lattice dynamics. In a microscopic level description, no expansion in terms of displacements from equilibrium sites exists for the Hamiltonian with purely hard-sphere interactions. Collisions entirely control the system, and between the collisions, motions of the hard spheres lose coherence very rapidly. This ballistic motion of the freely moving hard spheres in the crystal between collisions is quite analogous to the corresponding motion of the particles in the low-density fluid. In MWDA, the thermodynamic properties of the hard-sphere solid is successfully computed in terms of an equivalent liquid of much lower density.
The equivalent densityρ of the crystalline state in MWDA is generally much smaller than ρ 0 , and thus for the low-density system, the PY approximation (28) for the free energy is appropriate to use. However, when MWDA is applied to describe the amorphous solid state, two qualitatively different types of minima occur, as discussed in the previous section. First, the low α minimum, theρ comes out to be close to ρ 0 , and is not small in the metastable region beyond the freezing point. Thus, the approximation (28) is not well applicable in this case. On the other hand, for the highly localized state (for large α), theρ is much smaller than ρ 0 , and, in this case, the approximation (28) is more appropriate. For the low α minimum in the MWDA results for free energy, the boson peak is identified by accounting for the new entropy in terms of a modified density of states. This was described in the previous section. For the other minimum (see Figure 9) at α = α min on the higher side, the corresponding entropyŜ calculated from the partition function of the microscopic model with a Debye density of vibrational states, will agree (as shown in Figure 2) with the corresponding DFT result for the entropy obtained using the asymptotic Formula (11). This matching of entropy from the coarse-grained and microscopic models indicates that any correction for the density of vibrational modes over the Debye distribution g D (ω) will imply a zero correction to entropy S. How this will affect the height and position of the boson peak for amorphous states with a high degree of mass localization (α min >> α 0 ), will be studied elsewhere.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In a constant NVT ensemble comprising of 3N harmonic oscillators, the logarithm of the partition function is obtained as In writing the last equality, we change from discrete to continuum limit and express the sum over possible states in terms of the density of states g(ω). The latter is assumed to be nonzero up to an upper cutoff of frequency ω m . Using the standard thermodynamic relations, we obtain the entropy from the logarithm of Z N as, Derivative of Equation (A1) with respect to −β obtains for average energy as, Taking the average kinetic and potential energies to be the same (Ē) and equal to half of the total energy E, we obtain Since the total number of vibrational modes in the system is 3N, we have the following constraint in terms of the scaled frequency x = ω/ω m , where the scaled density of states has been defined asḡ(x) = (3ρ 0 ) −1 ω m g(ω). Now, Equation (A5) gives the free energy density for the noninteracting system in terms of the scaled frequency as, where we have defined the function κ(y) as, κ(y) = 1 4 y + 4 ln 1 − e −y − 2y e y − 1 , and m = βhω m .