Theoretical Rationale for a Thermodynamic Glass State

Article pubs.acs.org/JPCB

Theoretical Rationale for a Thermodynamic Glass State Isaac C. Sanchez* and Sean P. O’Keefe McKetta Department of Chemical Engineering, The University of Texas, Austin, Texas 78712, United States ABSTRACT: Using a chemical potential route, the square-well (SW) fluid model is solved in a quasichemical approximation (QCSW). At low temperatures, the liquid reaches a limiting density greater than the triple point density of a SW fluid but less than the equilibrium liquid-to-solid transition density for hard spheres. As this unique density is approached with decreasing temperature, the liquid entropy also approaches an asymptotic value, thus averting the “entropy catastrophe”. Mean-field models in the van der Waals (VDW) genre fail to predict this type of behavior. In VDW models, attractive force contributions to the equation of state incorrectly diverge with decreasing temperature as 1/T, whereas those for the QCSW model asymptote to a fixed value. The QCSW model posits the intuitively pleasing idea that, at high densities, attractive contributions to the configurational energy begin to saturate well before zero temperature is reached. As a consequence, the force balance between repulsive and attractive forces stabilizes the liquid density, which thereafter becomes effectively independent of temperature. This fixed density in turn fixes all other density-dependent thermodynamic properties. These lowtemperature, force-stabilized states are identified as glass states. ⎧∞ R < σ ⎪ ϕ(R ) = ⎨ −ε σ ≤ R ≤ λσ ⎪ ⎩ 0 R > λσ

I. INTRODUCTION Although ideas and models of the glass transition have been around for decades,1 it remains a subject that still engages our curiosity and attention.2,3 Interest in the nature of the glass transition was motivated by Kauzmann4 in 1948 when he noted that when the metastable liquid entropy line is extended to low temperatures it falls to the crystal entropy at what is now known as the “Kauzmann temperature”. Extension of this line reaches zero before temperature reaches zero, the so-called “entropy catastrophe”. The intervention of a liquid-to-glass transition resolves the crisis but raises a fundamental question: Is the glass transition a thermodynamic transition? An excellent review of the details of this subject exists.5 In 1958, Gibbs and DiMarzio6 developed a statistical lattice model of polymer liquids and argued that the glass transition was a thermodynamic, second-order phase transition. They proposed that the observed laboratory liquid-to-glass transition in polymeric systems is a kinetically controlled manifestation of this underlying thermodynamic transition. This dogma prevails to this day, but not without criticism.7 But how do these ideas for chain molecules translate to simple molecules such as toluene or elemental sulfur that also exhibit glass transitions? Is there a common physical principle? These questions will be addressed within the context of the quasichemical model described here. In a previous publication, an equation of state (EOS) in the van der Waals (VDW) genre was developed using scale particle theory (SPT) for the hard sphere (HS) contribution to the EOS.8 Herein, this SPT model is treated as a square-well (SW) fluid in a quasichemical (QC) approximation. A SW potential includes both attractive and repulsive interactions: © XXXX American Chemical Society

(1)

Others have developed QC-type approximations for a SW fluid. 9−11 These previous studies have concentrated on calculating the effective coordination number of a HS and use that result to calculate the free energy by isothermal integration of the energy or by simply adding the calculated energy to a HS equation of state. The present approach proceeds via a chemical potential route that differs from all previous studies and offers a more direct path to the equation of state. It also clearly defines what additional approximations must be made within the context of a QC approach.

II. QUASICHEMICAL APPROXIMATION Because of excluded volume crowding, an HS can attract at most a maximum number nmax of other spheres into its attractive domain. When λ = 1, the excluded volume around the central sphere equals 8 times the HS volume of πσ3/6. For λ > 1, the accessible shell volume (VS) around the central sphere for the placement of another sphere center is given by VS = 8(πσ 3/6)(λ 3 − 1)

(2)

The fraction of space occupied by nmax spheres within the shell volume is approximately given by Received: July 2, 2016 Revised: August 17, 2016

A

DOI: 10.1021/acs.jpcb.6b06653 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B nmax (πσ 3/6) nmax = ≡ ηmax VS 8(λ 3 − 1)

ηS =

(3)

(see Figure 1) where nmax is the maximum possible occupied volume fraction within the attraction shell. For random close

⟨n⟩(πσ 3/6) ze βε = nmax VS 1 + ze βε

(6)

This statement is not rigorously correct because sphere centers that lie outside the attractive domain (R > λσ) can also contribute to the local occupied fraction and not all spheres within the attraction zone contribute entirely to the occupied fraction (see Figure 1). From geometry, a local occupied fraction can be rigorously defined and measured in computer simulations,13 but is not employed here. For large enough λ, the occupied volume fraction in the attractive shell must approach the global volume fraction, i.e., ηS → η = (πσ3/6)ρ where ρ is the number density. Using this as an approximation for all values of λ and taking the high-temperature limit yields a self-consistent value for the z parameter: z lim ηS = η = ηmax βε→ 0 1+z z=

η θ ≡ ηmax − η 1−θ

(7)

Substituting this result into eq 5 yields θe βε = nmax 1 + θ(e βε − 1)

where θ ≡ η/ηmax is the average occupied fraction in the attraction shell in the high temperature limit. This interpretation of θ always holds providing λ is large enough but becomes an approximation for small λ. Notice that as βε → ∞, ⟨n⟩/nmax → 1 as it should. Also, with increasing external pressure, θ → 1 and again ⟨n⟩/nmax → 1 as expected. Number fluctuations within the attractive shell are given by

Figure 1. Two-dimensional schematic illustrates the contribution of hard spheres (HS) to the occupied volume in the attraction shell around the black sphere (the gray annulus). The volume of the attraction shell is 8(λ3 − 1)v0, where v0 = (π/6)σ3 is the HS volume. Although portions of a HS can enter the excluded volume region (the white annulus), HS centers cannot. Note that the red sphere centers in the attraction shell experience an attractive energy of −ε but only partially contribute to the occupied volume in the attraction shell. In contrast, an orange sphere whose center lies outside the attraction shell can still contribute to the occupied volume in the attraction shell. As illustrated, if λ = 3/2, the volume of the attractive shell is 19v0 but shrinks to less than 8v0 for λ = 5/ 4. This latter λ probably represents a realistic lower bound for the QCSW model. As λ increases, the QCSW model morphs into the SPT mean field model. See section III.D for discussion.

⟨n2⟩ − ⟨n⟩2 =

∂⟨n⟩ θ(1 − θ )e βε = nmax ∂(βε) [1 + θ(e βε − 1)]2

⟨n2⟩ − ⟨n⟩2 ⟨n⟩

packing, ηmax ≃ 0.64, but ηmax is a parameter fixed by λ and the choice of nmax. As illustrated in Figure 1, the above approximation for ηmax becomes more exact as λ increases and breaks down as λ → 1. In the spirit of a QC approximation, the probability Pn that n sphere centers are within the attraction domain of a given sphere is given by a binomial distribution ⎛ nmax ⎞ βε n ⎟(ze ) / Q ; Pn = ⎜ ⎝ n ⎠

= e−βε /2

⎛ nmax ⎞ βε n βε n ⎟(ze ) = (1 + ze ) max ⎠ n n=0

(10)

III. QCSW THERMODYNAMIC PROPERTIES A. Configurational Energy. Let ψ = −nε be the attractive interaction energy of a HS with n other HSs in its interaction domain, then the average interaction energy/HS is ⟨ψ ⟩ = −⟨n⟩ε = −nmax ε

(4)

where as usual β = 1/kT and z is a QC parameter to be determined self-consistently. The average number of interacting spheres within the attractive shell volume VS is given by

θe βε 1 + θ(e βε − 1)

(11)

It is also straightforward to show that the constant volume heat capacity is given by CV /k = (βε)2

nmax

∂ ln Q ze βε ⟨n⟩ = ∑ nPn = = nmax ∂(βε) 1 + ze βε n=0

(1 − θ )/θ nmax

There are two interesting limits where the fluctuations vanish: first as T → 0 and then as nmax → ∞. The latter occurs at all temperatures and is characteristic of a VDW-type model (see discussion in section III.D).

nmax

∑ ⎜⎝

(9)

and fluctuations relative to the mean are

12

Q=

(8)

(5)

∂⟨n⟩ = (βε)2 [⟨n2⟩ − ⟨n⟩2 ] ∂(βε)

(12)

B. Equation of State and Chemical Potential. Now let P0(ψ) dψ be the probability that a random insertion of a HS yields interaction energy ψ and P(ψ) dψ be the probability that

The volume fraction of space occupied by the spheres in the attraction shell is given by B

DOI: 10.1021/acs.jpcb.6b06653 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B any HS already in the fluid exhibits interaction energy ψ. These two distributions are related by14 P (ψ )e−βψ P(ψ ) = 0 B

The chemical potential (μ) is given by14 βμ = βμig (T , P) + ln(ZB) = βμig (T , ρ) + ln B

where the subscript ig denotes ideal gas value. C. Critical Properties. The critical point is determined by equating ∂(ηZ)/∂η)T and ∂2(ηZ)/∂η2)T to zero and yields the following for the critical temperature

(13)

where B is the insertion factor. Equation 13 provides a useful and insightful definition for the insertion factor: ∞

B=

∫−∞ P0(ψ ) dψ ∞

∫−∞ P(ψ )e βψ dψ

0

=

∫−∞ P0(ψ ) dψ 0

∫−∞ P(ψ )e βψ dψ

=

exp[βcε] = 1 +

Pins ⟨e βψ ⟩ (14)

lim ηc =

nmax →∞

(17)

The required insertion probability can be approximated (see section VI for discussion of this approximation) using the SPT model8,15 (18)

The EOS is given by βP 1 ≡ Z = 1 − ln B + ρ η

∫0

η

Z=

ln B dη = Z HS + Zint (19)

Z int

n max

ε→0

kT * =

nmax ε ⇒ SPT model 2ηmax

This transitional behavior has also been observed in MC simulations where the phase diagram shifts from a cubic (nonclassical) to a parabolic (classical) type with increasing λ.16,17 Also note the fundamental difference between the VDW like eq 28 and the QCSW EOS given by eq 21. The attractive contribution diverges as T → 0 in eq 28, but not in eq 21.

(20)

⎧ ln[1 + θ(e βε − 1)] ⎫ 1 + η + η2 n − 1⎬ + max ⎨ 3 βε (1 − η) θ(e − 1) ⎩ ⎭

IV. MODEL PARAMETERS As illustrated in eq 1, three parameters characterize the SW potential: the well depth ε, the HS diameter σ, and the range λσ; ε establishes the temperature scale (βε), and πσ3ρ/6 defines the HS volume fraction, η. Monte Carlo (MC) studies of the SW potential10,16−19 have primarily focused on the range 5/4 ≤ λ ≤ 2 with λ = 3/2 a popular choice. For λ > 7/4, the SW model begins to exhibit mean-field (VDW) like character.16,17 The QCSW model introduces an additional parameter, nmax, the maximum number of atoms that can be accommodated in the attraction zone around a HS, also known as the coordination number. MC studies indicate that nmax → 12 as a high-density

Z = 1 + b2η + b3η2 + ... (21)

where the dimensionless virial coefficients are given by k−1 ( −1)k + 1 ⎡ e βε − 1 ⎤ 3 ⎢ ⎥ k≥2 bk = 1 + k(k − 1) + nmax ⎢⎣ ηmax ⎥⎦ 2 k

(22)

Note that b2 is the correct second virial for the SW potential: b2 = 4[1 − (λ 3 − 1)(e βε − 1)]

1 + η + η2 T* −η ; T (1 − η)3

(28)

Using eq 18 for the HS contribution yields a complete EOS Z=

(27)

n ε 1 + η + η2 − max θβ[1 + 0(βε) + 0(βε)2 + ...] 2 (1 − η)3

Z= lim →∞

where ZHS is the HS contribution to the EOS that arises from ln Pins and Zint is the contribution from attractive interactions ⎧ ln[1 + θ(e βε − 1)] ⎫ ⎬ − 1 = nmax ⎨ θ(e βε − 1) ⎩ ⎭

73 − 7 = 0.1286... 12

At the other extreme, as nmax → 1, ηc → 0 and the critical temperature Tc → 0 by eq 25; as a consequence, the liquid−vapor transition vanishes as nmax → 1. D. VDW Limit. It is interesting to note that the limiting value of ηc as nmax → ∞ corresponds to the ηc for the mean-field SPT model.8 The SPT model uses a VDW-type energy in which every HS attractively interacts with all others, or equivalently, it obtains when nmax or λ → ∞. The effect of the interaction range on the QCSW model can be seen by keeping the product nmaxε constant while letting ε → 0 and nmax → ∞. In this VDW limit, the attractive interaction range as measured by λ approaches infinity while the interaction energy per HS remains constant (geometrically, the SW area remains fixed as the range increases).

Thus,

η(14 − 13η + 5η2) 2(1 − η)3

(26)

It should be noted that the critical volume fraction is determined by nmax, the maximum number of sphere centers that can be accommodated within the accessible shell volume. As nmax → ∞, or equivalently as λ → ∞, ηc reaches a limit of

(15)

(16)

ln Pins = ln(1 − η) −

(25)

(1 − 2ηc)3 + nmax (1 − ηc)3 [6ηc2 + 7ηc − 1] = 0

and the above average above can be calculated exactly in the QC approximation: ⎛ nmax ⎞ n n ⎜ ⎟z ∑nmax =0 ⎝ n ⎠ ⟨exp(−βεn)⟩ = = [1 + θ(e βε − 1)]−nmax Q

ln B = ln Pins + nmax ln[1 + θ(e βε − 1)]

1 ⎡ nmax (1 − ηc)4 ⎤ θc⎢ − 1⎥ ⎣ (1 + 2ηc)2 ⎦

with the critical value of the occupied volume fraction ηc given by the single positive root of the following polynomial:

where Pins is the insertion probability for a successful random insertion of a HS into a HS fluid where the average HS interaction energy is ⟨ψ⟩. The second equality above follows since P0(ψ) = 0 for ψ > 0 for a HS fluid. Thus ln B = ln Pins − ln

(24)

(23) C

DOI: 10.1021/acs.jpcb.6b06653 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B limit for λ = 3/2.20 It also agrees with a calculation using pair distribution functions for HS fluids (see the Appendix). A choice of nmax = 12 and λ = 3/2 defines the maximum packing fraction ηmax within the attraction shell: ηmax =

nmax 3

8(λ − 1)

=

12 = 0.63 19

setting Z = 0. On the liquid−vapor coexistence (COEX) curve, the saturated vapor compressibility factor Zvap → 1 with decreasing temperature, whereas for saturated liquid, Zliq → 0 as temperature decreases. The saturated vapor becomes more like an ideal gas because of decreasing vapor pressure (Psat → 0), whereas on the liquid side of the COEX curve, Zliq = Psat/ kTρliq→ 0 for the same reason. For example, for many liquids at their normal boiling points, −ln Zliq = 5.5 ± 0.2 and represents a corresponding state.8 The liquid compressibility factor rapidly becomes even smaller as the triple point is approached. An important implication of Zliq → 0 is that the negative (attractive) and positive (repulsive) contributions to the EOS must exactly cancel. Note from Figure 2 that as the temperature approaches zero the density approaches a limiting value of η = 0.478. This type of behavior is not obtained for any mean-field (VDW type) equation of state. The simple physics is that the attraction shell, which can only accommodate nmax HS centers, begins to saturate well before zero temperature is reached. For example, the average occupancy in the attraction shell already has reached 99% at βε = 4 (or T/Tc = 0.22). Once the attraction shell saturates, further lowering of temperature no longer influences the EOS. As mentioned above, long before this saturated state is approached, the vapor pressure of the liquid begins to vanish so that Z → 0 and the QCSW EOS reduces to a temperature-independent equation:

(29)

This close packing value agrees well with what is believed to be a good estimate of random close packing of spheres of the same size.12 However, it is expected that the SW fluid will crystallize well before the high-density limit is reached since MC simulations indicate that the triple point for λ = 3/2 occurs at η = 0.437.21 For comparison purposes, the normal boiling point for many simple and organic liquids is a corresponding state that occurs at η = 0.413 ± 0.005, but the triple point is not in general a corresponding state.8 Figure 1 illustrates the attraction shell for λ = 3/2.

V. EVIDENCE FOR A THERMODYNAMIC GLASS STATE Using the parameters λ = 3/2 and nmax = 12, the calculated saturated liquid density is shown in Figure 2. This requires calculation of the coexisting liquid and vapor densities by equating chemical potentials. However, for reduced temperatures of T/Tc ≤ 0.8, the QCSW saturated liquid density can be accurately calculated with an error of less than 1 part in a 1000 by

lim Z =

βε→∞

1 + η + η2 − nmax = 0 (1 − η)3

(30)

The solution of this equation for nmax = 12 is η = 0.478. This limiting density is the global density; the local density in the attraction shell around a HS at saturation is higher, 12/8(λ3 − 1) = 12/19 = 0.63. However, if the central atom is included, the density of a 13 HS cluster becomes 13/(8λ3) = 13/ 27 = 0.481, which to 2 significant figures is identical to the global density of 0.478. This reflects well on the self-consistency of the QCSW model for the choices λ = 3/2 and nmax = 12. For any completely repulsive system, such as a HS fluid, density is controlled by βP where P is the pressure. In a system with attractive interactions of strength ε, both βP and βε control the density with density increasing as either variable increases. However, for a liquid in equilibrium with its vapor and below its normal boiling point, only βε controls the density with repulsive forces acting to keep the liquid from completely collapsing. But what happens when the liquid becomes so dense that its density no longer responds to an increase in βε? For example, if attractive interactions were added to a metastable HS fluid (η > 0.494), would its density increase? The intuitive expectation is that it would change very little as long as the system remained disordered (it would begin to densify if crystallization intervenes). So from a qualitative viewpoint, the QCSW model is consistent with the intuitive idea that attractive interactions become ineffective in controlling liquid density at high densities and a hypothetical limiting density, defined by the force balance, might be reached prior to absolute zero. As noted above, the triple-point density for a SW fluid for λ = 3/2 is 0.437.21 Thus, the observed limiting density of 0.478 is in the metastable liquid regime. It represents a density where the liquid gets “stuck”. Equation 30 defines what “stuck” means thermodynamically. With decreasing temperature, liquid density is determined by a balance of attractive and repulsive (packing) forces. Once attractive forces become ineffective the liquid

Figure 2. Variation of density (occupied volume fraction) with reduced temperature for the QCSW model that asymptotes to 0.478 (nmax = 12 and λ = 3/2). All densities are saturated liquid densities. Glass states are operationally defined as those that have reached at least 99.9% of the saturated configurational energy. MC simulations indicate that the triple point for the SW model occurs at η = 0.437.21 States with greater densities are metastable/glassy. Interestingly, extension of the dashed straight line to zero temperature intercepts the density axis at η = 0.545, which is the estimated Kauzmann density for the putative fluid-to-glass transition for hard spheres.13 D

DOI: 10.1021/acs.jpcb.6b06653 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

changing the density or temperature,8 or equivalently, it can be viewed as the entropy loss in transferring a molecule from an ideal gas to the liquid with both gas and liquid at the same temperature and density.23 For comparison, this entropy change for a HS fluid is given by

approaches a density that defines this balance. If the more accurate Carnahan−Starling EOS21 is used, this density becomes 0.488 or about 2% greater.22 Some insight can also be gained by calculating the liquid entropy −ln B = β[μ − μig (T , ρ)] = β[HE − TS E]

(31)

S E/k = ln Pins +

where the superscript E denotes an excess property relative to an ideal gas at the same temperature and density. It is straightforward to show that βHE = β⟨ψ ⟩+Z − 1

= −4.7 for η = 0.478

(33)

At low temperatures where Z → 0, it is easy to show that SE becomes a function of density only:

VI. APPROXIMATE NATURE OF THE QCSW MODEL Although the QCSW model provides a new physical insight, it remains an imperfect model. As an example, using eqs 25 and 26, the QCSW critical properties for λ = 3/2 and nmax = 12 are ηc = 0.1043 and kTc/ε = 1.136 as compared to the exact MC values of ηc = 0.16 and kTc/ε = 1.218.19 Equation 15 is exact for the SW potential. The QC approximation used to calculate ⟨exp (−nβε)⟩ in eq 15 takes into account short-range fluctuations in energy but fails to account for the long wavelength density fluctuations that occur near the critical point. From eq 14 it is seen that the insertion probability for a SW fluid will have a finite number of contributions:

lim S E/k = ln Pins + nmax ln θ − 1

βε→∞

= −20.0 for η = 0.478

(35)

This reveals that adding attractive interactions to a HS fluid causes an additional and significant lowering of the entropy as expected. In the jargon of statistical mechanics, attractive interactions act to restrict the available configuration space. This loss of configuration space corresponds to an entropy loss, which is compensated by the lowering of the configurational energy (energy−entropy compensation).

(32)

Thus, the excess entropy of the fluid relative to an ideal gas at the same density and temperature is given by S E/k = ln B + β⟨ψ ⟩+Z − 1

1 + η + η2 −1 (1 − η)3

(34)

Figure 3 illustrates the entropy behavior as a function of temperature. Equation 34 indicates that the configurational

Pins = P0(0) + P0( −ε) + P0( −2ε) + ...

(36)

P0(−nε) is the probability that a random insertion of a HS into the SW fluid system of density η will result in n attractive interactions. These probabilities will depend on the average interaction energy, ⟨ψ⟩, at the given fluid density. At high densities, which are the most interesting, Pins is approximated as the HS fluid insertion probability given by eq 18. This should be a good approximation at high densities since the void distribution in a SW fluid at high densities should not differ much from that of a HS fluid. However, at low densities, this approximation will underestimate the insertion probability because clustering effects favor the formation of larger cavities. On the positive side for the QCSW model, the critical compressibility factor, Zc = 0.338, which is an improvement over the mean-field SPT model (0.360) and the classical VDW model (3/8). In a previous publication, three new critical properties have been identified as corresponding states: fugacity coefficient, thermal pressure coefficient, and internal pressure.24 The strongest of these corresponding states is the critical fugacity coefficient, which has a near-universal value of 2/3 ± 1% for 65 classical fluids, which includes water and other hydrogen bonding fluids as well as three quantum fluids. The QCSW calculated value of the fugacity coefficient is 0.671. The calculated values of the other two critical properties are also in good agreement with experiment.

Figure 3. Plot of excess configurational entropy calculated from eq 33 against reduced temperature. Calculated saturated liquid densities, illustrated in Figure 2, were used at every temperature to calculate the entropy.

liquid entropy asymptotes to a finite value at low temperatures and avoids the “entropy catastrophe.” At reduced temperatures of less than about 0.1, the liquid becomes a glass. The excess entropy discussed above is also known as the selfsolvation entropy. It can be viewed as the entropy loss that results from turning on repulsive and attractive interactions without

VII. SUMMARY Experience, especially with computer simulations, suggests that as a liquid becomes dense repulsive forces dominate its structure and density with attractive forces playing a minor role. Meanfield models in the VDW genre do not yield this important E

DOI: 10.1021/acs.jpcb.6b06653 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B qualitative result. Instead, attractive force contributions to the EOS increase with decreasing temperature and incorrectly diverge as T → 0 as 1/T. This negative divergence requires the repulsive contribution to the EOS to diverge positively and usually at a nonphysical value of density. In contrast, the QCSW model indicates that this tug-of-war between attractive and repulsive forces results in a unique and physically realistic limiting density at low temperatures with no divergences in the EOS. For the QCSW model considered herein, this limiting density is η = 0.478, which is greater than the triple point density determined from MC simulations20 (0.437) but less than the equilibrium fluid-to-solid transition density of a HS system24 (η = 0.494). Using the MC value of the critical density (ηc = 0.16) yields a physically realistic reduced QCSW density of 3.0.25 Many simple liquids such as the inert elements, methane, and LennardJones fluids26 begin to freeze at reduced densities of about 2.6. A reduced density of 3 would place the QCSW limiting density in the metastable liquid range for simple fluids. This intuitive and desirable saturation behavior obtains because the attractive range in the QCSW model is restricted. For the specific SW model considered herein, the maximum number of possible attractive interactions of one HS with its nearest neighbors within a range of λ = 3/2 is nmax = 12. The attraction shell, which can only accommodate 12 HS centers, begins to saturate well before zero temperature is reached. Once the attraction shell begins to saturate, further lowering of temperature has little influence on the EOS. However, if the range (λ) is allowed to increase, the QCSW model morphs into a VDW-type model with its undesirable qualitative features (see eq 28 and associated discussion). Keeping attractive forces short ranged yields a low temperature range where attractive forces become ineffective in controlling density. At all temperatures, repulsive (packing) forces balance attractive forces. At low temperatures, eq 30 defines the unique density where this balance is achieved. As mentioned above, the limiting density corresponds to a reasonable reduced density of 3.0. Other checks on the physical reasonableness of the calculated densities can also be made. The limiting density of 0.478 is the global density; the local density in the attraction shell around a HS at saturation is higher, 12/8(λ3 − 1) = 12/19 = 0.63. However, if the central atom is included, the density of a 13 HS cluster becomes 13/(8λ3) = 13/ 27 = 0.481, which is very close to the global density of 0.478. This reflects well on the self-consistency of the QCSW model for the choices λ = 3/2 and nmax = 12. It has also been observed that the saturated liquid density in the normal liquid range (NLR) is linear for many liquids with only water and helium as known exceptions.27 The five stable liquid points in Figure 2 more than span the NLR and form the dotted straight line shown in Figure 2. If this line is extended to zero temperature, it intersects at η = 0.545 or at a reduced density of 0.545/.16 = 3.4. Extrapolation of experimental density data for Ar, Kr, Xe, CH4, S, O2, and N2 to zero temperature yields an average intercept of 3.6227 in reasonable agreement with the QCSW value of 3.4. As an additional observation of interest, the density η = 0.545 also corresponds to the estimated Kauzmann density for the putative fluid-to-glass transition for hard spheres.13 What separates liquid from glass? The QCSW model does not exhibit a second-order phase transition that might define a clear boundary between glass and liquid states. However, with decreasing temperature both the thermal expansion coefficient and heat capacity decrease dramatically and asymptote to zero at zero temperature. The heat capacity is illustrated in Figure 4.

Figure 4. Variation of the QCSW constant volume heat capacity with temperature. See eq 12. Liquid densities used in the calculation are the same as in Figure 2.

This behavior might be termed a pseudo-second-order phase transition and tends to mimic what is often observed experimentally. But what clearly emerges from the QCSW model is that on cooling the liquid gets “stuck” at a density where the configurational energy has saturated. This causes other thermodynamic properties, such as the configurational entropy, to also saturate. So an operational definition is adopted for the glass transition: the liquid-to-glass transition occurs when the conf igurational energy reaches 99.9% of its saturated value. This occurs near βε = 10 (or about T/Tc = 0.1). This view of the glass transition differs from the prevailing idea that the ideal thermodynamic liquid-to-glass transition occurs when the configurational entropy reaches zero (or a constant) via a second-order phase transition. Similarity exists in the sense that the configurational entropy of the QCSW model reaches a limiting value, but the cause and effect relationship is different. The limiting entropy property of the QCSW model is a direct consequence of configurational energy saturation. For shortrange attractive interactions, excluded volume limits the number of other spheres within the attractive domain of a given sphere. At high densities, this limiting number is quickly approached. In real molecular systems with short-range interactions, whether monomeric or polymeric, it is expected that the same limited attractive interaction between atoms or group of atoms will prevail. Extending the concept to chain molecules follows because if short-ranged monomer−monomer attractive interactions in polymer models are treated in a QC approximation, F

DOI: 10.1021/acs.jpcb.6b06653 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(7) Stillinger, F. H. Supercooled Liquids, Glass Transitions, and the Kauzmann Paradox. J. Chem. Phys. 1988, 88, 7818−7825. (8) Sanchez, I. C. Dimensionless Thermodynamics: A New Paradigm for Liquid State Properties. J. Phys. Chem. B 2014, 118, 9386−9397. (9) Guo, M.; Wang, W.; Lu, H. Equations of State for Pure and Mixture Square-Well Fluids. II. Equations of State. Fluid Phase Equilib. 1990, 60, 221−237. (10) Heyes, D. M. Coordination Number and Equation of State of Square-Well and Square-Shoulder Fluids: Simulation and QuasiChemical Model. J. Chem. Soc., Faraday Trans. 1991, 87, 3373−3377. (11) Haghtalab, A.; Mazloumi, S. H. A New Coordination Number Model for Development of a Square-Well Equation of State. Fluid Phase Equilib. 2009, 280, 1−8. (12) Speedy, R. J. On the Reproducibility of Glasses. J. Chem. Phys. 1994, 100, 6684−6691. (13) Sanchez, I. C.; Lee, J. S. On the Asymptotic Properties of a Hard Sphere Fluid. J. Phys. Chem. B 2009, 113, 15572−15580. (14) Widom, B. Potential-Distribution Theory and the Statistical Mechanics of Fluids. J. Phys. Chem. 1982, 86, 869−872. (15) Reiss, H.; Frisch, H. L.; Lebowitz, J. L. Statistical Mechanics of Rigid Spheres. J. Chem. Phys. 1959, 31, 369−380. (16) Vega, L.; deMiguel, E.; Rull, L. F.; Jackson, G.; McLure, I. A. Phase Equilibria and Critical Behavior of Square Well Fluids with Variable Width by Gibbs Ensemble Monte Carlo Simulation. J. Chem. Phys. 1992, 96, 2296−2305. (17) Orea, P.; Varga, S.; Odriozola, G. A Heuristic Rule for Classification of Classical Fluids: Master Curves for Mie, Yukawa, and Square-Well Potentials. Chem. Phys. Lett. 2015, 631−632, 26−29. (18) Orkoulas, G.; Panagiotopoulos, A. Z. Phase Behavior of the Restricted Primitive Model and Square-Well Fluids from Monte Carlo Simulations in the Grand Canonical Ensemble. J. Chem. Phys. 1999, 110, 1581−1590. (19) Orkoulas, G.; Fisher, M. E.; Panagiotopoulos, A. Z. Precise Simulation of Criticality in Asymmetric Fluids. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 63, 1−14. (20) Khanna, K. N. A New Coordination Number Model and Heat Capacity of Square-Well Fluids of Variable Width. Fluid Phase Equilib. 2006, 243, 101−106. (21) Liu, H.; Garde, S.; Kumar, S. Direct Determination of Phase Behavior of Square-Well Fluids. J. Chem. Phys. 2005, 123, 174505. (22) Carnahan, N. F.; Starling, K. E. J. Chem. Phys. 1969, 51, 635−636. (23) Ben-Naim, A. Solvation Thermodynamics; Plenum Press: New York, 1987. (24) Sanchez, I. C.; Boening, K. L. Universal Thermodynamics at the Liquid-Vapor Critical Point. J. Phys. Chem. B 2014, 118, 13704−13710. (25) Hoover, W. G.; Ree, F. H. Melting Transition and Communal Entropy for Hard Spheres. J. Chem. Phys. 1968, 49, 3609−3617. (26) Mastny, E. A.; de Pablo, J. J. Melting Line of the Lennard-Jones System, Infinite size, and Full Potential. J. Chem. Phys. 2007, 127, 104504. (27) Sanchez, I. C.; O’Keefe, S. P.; Xu, J. F. New Zeno-Like Liquid States. J. Phys. Chem. B 2016, 120, 3705−3712. (28) Barker, J. A.; Henderson, D. Theories of Liquids. Annu. Rev. Phys. Chem. 1972, 23, 439−484.

one expects the same energetic saturation behavior as observed in the monomeric QCSW model.

VIII. CONCLUSION The QCSW model endorses the intuitively pleasing concept that at high densities, the short-ranged attractive contribution to the configurational energy begins to saturate well before zero temperature is reached. As a consequence, the force balance between repulsive and attractive forces defines a unique liquid density, which thereafter effectively becomes temperature independent. This fixed density in turn establishes the value of the configurational entropy and all other density-dependent thermodynamic properties. These low temperature, forcestabilized states are identified as glass states. This new physical insight offered by the QCSW model appears to be widely applicable to liquids with short-range attractive forces, which includes both monomeric and polymeric organics. It might even be applicable to ionic liquids with highly screened columbic interactions. Molten metals might be excluded because of the long-range attractive effects of valence electron delocalization.

■

APPENDIX In the high density limit, the QCSW model behaves structurally like a HS fluid. In this limit, ⟨n⟩ can be approximated as the number of HS centers between σ and λσ of a given central sphere. Using the pair distribution g(r;ρ), we have ⟨n⟩ = 4πρσ 3

∫1

λ

(r /σ )2 g[r /σ ; ρσ 3] d(r /σ )

Tabulated values of g(r;ρ) are available in the literature28 as a function of r/σ at selected values of ρσ3. The above integral can be evaluated numerically for λ = 3/2 with the following results: η 0.340 0.393 0.445 0.484 ⟨n⟩ 9.0 10.5 11.8 12.6 Therefore, at the limiting density of 0.48, the calculated coordination number of 12.6 is in good agreement with the chosen nmax = 12.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS Financial support from the William J. Murray, Jr., Endowed Chair in Engineering is gratefully acknowledged. REFERENCES

(1) Berthier, L.; Biroli, G. Theoretical Perspective on the Glass Transition and Amorphous Materials. Rev. Mod. Phys. 2011, 83, 587− 645. (2) Langer, J. The Mysterious Glass Transition. Phys. Today 2007, 60, 8−9. (3) Berthier, L.; Ediger, M. Facets of Glass Physics. Phys. Today 2016, 69, 40−46. (4) Kauzmann, W. The Nature of the Glassy State and the Behavior of Liquids at Low Temperatures. Chem. Rev. 1948, 43, 219−256. (5) Debenedetti, P. G.; Stillinger, F. H. Supercooled Liquids and the Glass Transition. Nature 2001, 410, 259−267. (6) Gibbs, J. H.; DiMarzio, E. A. Nature of the Glass Transition and the Glassy State. J. Chem. Phys. 1958, 28, 373−383. G

DOI: 10.1021/acs.jpcb.6b06653 J. Phys. Chem. B XXXX, XXX, XXX−XXX