J. Phys. Chem. 1991, 95,4540-4551
4540
by an example in Figure 6. Of course only the branch that intersects the abscissa is physically relevant. This intersection, by the way, happens at { = Sc, given by
Substitution of (17) and (18) into (16) yields f=
Jb+ k T X
or
The physically meaningful branch has two horizontal asymptotes: lim N = min(ahl,nl,, bhlflIb,chl$lc, ...) C--m
lim N = -min(ph2g2,,, qh2qn2q,r h g 2 , , ...) (-+m
This can be seen as follows. For f
-
-=, eq 19 becomes
This leads to a polynomial equation in N:
a polynomial equation in N with obvious roots ahi,ni,, bhibnib, Chigicr... The smallest of these roots is the physically relevant N(-m). Analogously, f yields the polynomial equation +
The degree of this equation equals max(a b + c + ...,p
+
+ q + r + ,,.)
+-
with roots
Solution of the equation gives the desired characteristic N( {). It often happens that the equation yields more than one real solution. This gives rise to several branches N ( f ) ,as illustrated
- p h 2 g p , -qhqnar -rh2,4,,
...
where again the value closest to zero has to be retained.
Splnodal Curve of Some Supercooled Llquids Pablo G. Debnetletti,* V. S. Raghavan, and Steven S. Borick Department of Chemical Engineering, Princeton University, Princeton, New Jersey 08544-5263 (Received: November 15, 1990)
There exist two possible limits to the extent to which a liquid can be supercooled: the K a u p n a ~temperature and the spinodal curve. The virial theorem imposes severe constraints on the type of interactionsthat can give rise to loss of mechanical stability upon supercooling and therefore to a supercooled liquid spinodal. Systemscomposed of particles interacting via pair potentials whose repulsive core has a positive curvature (such as the Lennard-Jones potential) cannot become mechanically unstable upon supercooling. Systems composed of particles interacting via potentials whose repulsive a r e is softened by a curvature change are capable of losing stability upon supercooling, and of contracting when heated isobarically. This is consistent with the idea that loss of stability upon supercooling can only occur for liquids capable of contracting when heated. Microscopically, this occurs via the formation of open structures which can be collapsed into denser arrangements through the input of thermal and mechanical energy. In the quasichemical approximation, a very simple model of a core-softened fluid, the lattice gas with attractive nearest-neighbor and repulsive next-nearest-neighbor interactions, exhibits density anomalies in one, two, and three dimensions, and a reentrant, continuous spinodal bounding the superheated, supercooled, and subtriple liquid states in three dimensions.
Introduction Liquids can be cooled below their freezing temperature without solidifying: they can be supercooled. In an ideally purified liquid devoid of any suspended impurity, small nuclei of solid phase are formed exclusively as a result of spontaneous density fluctuations. When such nuclei reach a critical size they grow, and a solid phase To whom correspondence should be addressed.
0022-365419112095-4540$02.50/0
is formed. The formation of critical nuclei (homogeneous nucleation) is a barrier which must be overcome in order for the phase transition to occur.14 In the absence of dissolved impurities, homogeneous nucleation governs the initial rate a t which the (1) Turnbull, D.; Fwher, J. C . J. Chem. Phys. 1949, 17, 71.
(2) Turnbull, D. J . Phys. Chem. 1962, 66, 609. (3) Walton, A. G. Science 1965, 148, 601. (4) Turnbull, D. Contemp. Phys. 1969, 10, 473.
0 1991 American Chemical Society
The Journal of Physical Chemistry, Vol. 95, No. 11, 1991 4541
Spinodal Curve of Some Supercooled Liquids metastable liquid relaxes toward a stable condition and solidifies. Homogeneous nucleation rates exhibit a very sharp dependence on the extent of supercooling. Consequently, when a liquid is supercooled, it evolves suddenly from a condition of apparent stability to one in which the solid phase grows spontaneously and very rapidly. As long as the sequence of events is kinetically controlled, the point at which the rapid growth of the new phase occurs can be varied by changing the design of the experiment. An example of this is the emulsification of supercooled liquid samples so that the rate of critical nuclei formation in any one droplet (nuclei per unit time) is small even though the intrinsic homogeneous nucleation rate (nuclei per unit time per unit volume) is large.s Such a limit, however, is not a true property of a particular liquid so much as it is a characteristic of the experimental technique. Supercooling limits determined by kinetic considerations, in other words, are not infinitely sharp. Nucleation is not the only kinetically controlled mechanism that can lead to solidification upon supercooling. If the rate of cooling is high enough for nucleation to be bypassed, the liquid can vitrify. Glass formation occurs when the relaxation time required for the liquid's molecules to rearrange their configurations in response to the imposed temperature change becomes comparable to the experimentally accessible Molecules then get trapped in a configuration and form an amorphous solid. The particular frozen-in configuration into which molecules get trapped is dependent on the cooling history to which the sample has been exposed. The cooling rate dependent temperature (or, more precisely, narrow temperature interval) at which the liquid solidifies is the glass transition temperature. Its lowest possible value is the temperature at which the entropy difference between the liquid and crystalline phases disappears (the Kauzmann temperature*). The Kauzmann temperature imposes a sharply defined (thermodynamic) lower limit to the possible existence of the liquid state of a given substance, since upon further supercooling the hypothetical liquid would have a lower entropy than the corresponding crystalline phase, a condition often referred to as entropy catastrophe. The Kauzmann temperature is unattainable because the slowing down of molecular motion inevitably gives rise to kinetically controlled glass transitions. There is another theoretical limit to the extent of possible liquid supercooling, namely, attainment of the spinodal locus, beyond which the liquid becomes mechanically ~ n s t a b l e . ~For a pure substance, the spinodal is defined by the simultaneous Occurrence of the conditions KT = m (1) cp = (2) where KT is the isothermal compressibility, and cp, the isobaric heat capacity. The critical point is the only stable condition along a spinodal, which is otherwise a locus of unstable limits of stability. For a stable or metastable pure substance, cp and KT are positive and finite except at the critical point. The mechanism of instability for a pure substance is the unbounded growth in density fluctuations'O (3)
In practice, a spinodal can be approached but never reached. Kinetically controlled relaxation mechanisms such as homogeneous nucleation will trigger a phase transition before a limit of stability can actually be attained. Nevertheless, it is of practical and theoretical significance to determine the location of the spinodal, (5) Rasmussen, D. H.; MacKenzie, A. P. J. Chem. Phys. 1973,59,5003. (6) Zallen, R. The Physics of Amorphous Solids: Wiley: New York,1983; Chapter I . (7) Angell, C. A. 11th IUPAC Conf. Chem. Thermodyn., Plenary Lrcr.; Como, Italy, Augusr 1990. ( 8 ) Kauzmann, W. Chem. Rev. 1948,43, 219. (9) Modell, M.;Reid, R. C. Thermodynamics and irs Applications. 2nd 4.Prentice ; Hall: Englcwood Cliffs, NJ, 1983; Chapter 9. (10) Lifshitz, E. M.; Pitaevskii, L. P. Statistical Physics, Part I, 3rd 4.; Vol. 5 of Course o Theorerical Physics; Landau, L. L., Lifshitz, E. M.; Pergamon: Oxford( 1980; Chapter 12.
since it represents a sharp limit beyond which a particular state of matter cannot exist and within which it can. The question naturally arises, therefore, as to the location of the spinodal for supercooled liquids. Molecular considerations, however, lead one to investigate critically the very existence of a supercooled liquid spinodal before attempting to predict its location. Consider the usual case of a liquid that contracts upon freezing, for which supercooling can be caused by isothermal compression. A spinodal limit would then imply a condition of infinite compressibility brought about by compression (as with supercooled vapors), but at densities such that liquid behavior is largely dominated by repulsive interactions. Clearly, we must start by exploring the conditions under which such behavior can possibly occur. In what follows, we investigate the necessary conditions for the existence of a spinodal limit to supercooling. In order to do so, we need to assume that the Kauzmann temperature is lower than the spinodal temperature, whenever the latter exists. This assumption needs to be verified for each specific substance but is required for an investigation of necessary conditions for spinodal collapse. Stability and the Virial Theorem
As a starting point for a molecular-based analysis of the existence of a supercooled liquid spinodal, we consider a fluid whose molecules interact via pairwise additive central forces. For such a fluid, we can write P = p [ k T (P/6)] (4)
+
where \k, the virial, is given by
= N(',,'f,l)
(5)
with r,, = ri - r , f,,, the force on molecule i due to j , and where angle brackets denote thermodynamic averaging. In eq 4 we have assumed N2 c N(N - 1). It follows from eq 4 that
Therefore, for a stable or metastable fluid we must have (7)
This implies the following: (i) under conditions of temperature and density such that the virial increases upon isothermal compression, loss of stability can only occur under tension ( P < 0); (ii) when loss of stability Occurs under positive pressure, the virial decreases upon isothermal compression. The former case corresponds to loss of stability in superheated liquids below a threshold temperature (TIT, = 27/32 for a van der Waals fluid) or, equivalently, above a threshold density ( p / p c = 3/2 for a van der Waals fluid). The latter corresponds to loss of stability in supercooled vapors or in superheated liquids above a threshold temperature (or below a threshold density). Consider for example the stability of a fluid whose molecules interact via a pair potential of the form
9 = e(u/r)n
(8)
P = Nne((a/r)")
(9)
for which
Since both the pressure and the virial's rate of change with respect
4542 The Journal of Physical Chemistry, Vol. 95. No. 1 1 , 1991
Debenedetti et al. 50 40
j
30
re
' 20 10
0
Figwe 1. Density dependence of the attractive (a), repulsive (r), and total virial (*/Ne) for a Lennard-Jones-typefluid below the Boyle temperature.
19
-
15 0.80
Figure 3. High-temperature density dependence of the attractive, repulsive, and total virial (*/Ne) for a Lennard-Jones fluid. kT/c = 6.5. NVT molecular dynamics of 5 0 0 atoms.
REPULSIVE ATTRACTIVE
I
I
0.82
0.84
I
0.86 DENSITY
I
I
0.88
0.90
1
TOTAL
Figure 4. A core-softened potential. DENSITY
Figure 2. Density dependence of the attractive and repulsive (a, top) and total (b, bottom) virial (*/Ne) of the Lennard-Jonesium at its triple point temperature (kT/c = 0.68) and in the vicinity of the liquid triple point density (pd = 0.86). NVT molecular dynamics of 5 0 0 atoms.
to isothermal density increase are positive, the stability inequality is never violated. There can be no supercooled liquid spinodal for a soft sphere fluid. Consider now a fluid whose molecules interact via a Lennard-Jones-type potential
4 = e[(u/r)n- (a/r)"']
(n > m )
(12)
for which the virial is given by
* = Nc[nuJyr-")- ma"(r-")] = N e [ n ( ( a / r ) " )- m ( ( u / r ) m ) ] (13) For this type of potential, Figure 1 shows the density dependence of the attractive and repulsive contributions to the fluid's virial (*/Ne) at low temperature. Above a certain temperature, the virial-density curve is monotonically increasing and has zero slope at the origin. As shown in the Appendix, this occurs when the second virial coefficient vanishes (Le., at the Boyle temperature). Returning to Figure 1, it follows from eq 7 that the only density range Ghere'the fluid can lose stability under positive pressure is below p * ; above p*, it can only lose stability under tension. Therefore, it is necessary (but not sufficient) in order for a supercooled liquid spinodal to exist at positive pressure that p' < P * , where pi is the liquid density in equilibrium with the solid at
the given temperature. This is a stringent requirement. Figure 2 shows the density dependence of the Lennard-Jones fluid's vinal at the triplepoint temperature, over a range of densities comprising the liquid's triple point (p$ = 0.86). Clearly, p' > p*, and there can be no supercooled spinodal at this temperature, since the metastable region is a d by compression and loss of stability can only occur under tension. Above the Boyle temperature, the virial behaves as shown in Figure 3. Under these conditions, there can be no supercooled liquid spinodal because the metastable region is accessed under pressure and the virial increases upon compression. Any fluid whose molecules interact via a potential whose repulsive core has a positive curvature, therefore, cannot possibly have a supercooled liquid spinodal above the Boyle temperature. Below the k y l e temperature, severe restrictions are imposed by the virial theorem on the possible divergence of a fluid's compressibility upon supercooling.
Core Softening Given the stringent limitations imposed by the virial theorem on the possible existence of a supercooled liquid spinodal, it is logical to inquire as to the type of potential, if any, that is consistent with loss of stability upon supercooling. Consider the potential shown in Figure 4, with inflection points within the repulsive core at rl and r2. For r2 > r > r , , the force between two molecules decreases as they approach each other. This behavior will henceforth be referred to as core softening, a description that also applies to the case (not shown) in which a secondary well appears between rl and r2. Mathematically, core softening means that A(rA < 0 when Ar < 0 for r2 > r > r, (i-e.,the product rfd-
The Journal of Physical Chemistry, Vol. 95, No. 11, 1991 4543
Spinodal Curve of Some Supercooled Liquids when the separation decreases). Taking into account the relationship between force and potential, this condition can be expressed as rz > I > rl r4” 4’ < 0 (14)
P
+
and also
4” > 0
r
< rl and r > r2
“h
I
(15)
where the second condition imposes that the decrease in the virial upon approach occur within the repulsive core (therefore, we require that 4’’ > 0 for some r > r2and all r < rl). Core softening has been extensively studied by Stell and co-workers.11-16 They investigated phase transitions caused by an abrupt diminution of the core diameter, which they referred to both as core softening and as core collapse. Although the emphasis of that work was on phase behavior, Stell and HemmerI2 reported density anomalies in one dimension for a linearly softened potential with long-range attractions. The contribution to the total virial due to a pair of molecules interacting via a core-softened potential does not increase monotonically as the separation decreases below r’ (potential minimum). Therefore, upon compression, the total virial does not increase monotonically with density, and the stability inequality can now be violated at high density. Consider now the thermal expansion coefficient of a coresoftened fluid. Differentiating eq 4, we obtain
Therefore, for a stable or metastable fluid, the thermal expansion coefficient will be positive so long as
(g)>-6k where the restriction to stable fluids comes from equating the signs of (aP/dT), and ( d T / d ~ )which ~ is necessarily true so long as (aP/du)Tis negative. If a given number of molecules is heated inside a rigid container, the only new contributions to the virial are those due to interpenetration of repulsive cores by pairs of energetic molecules. For a potential with positive curvature in its repulsive core (no core softening), these new contributions must necessarily lead to an increase in the virial because at the point of closest approach between two molecules during a given collision the pairwise virial is larger than for all greater separations. Consequently,eq 17 is satisfied for fluids which interact via pair potentials whose repulsive cores have only positive curvature, and such fluids cannot exhibit density anomalies. A necessary condition for a fluid to have a negative thermal expansion coefficient somewhere in its phase diagram is therefore that the isochoric rate of change of the virial with respect to temperature be negative for some condition of temperature and pressure. Core softening can lead to this condition because at the point of closest approach between two molecules during a given collision the pairwise virial is not necessarily larger than for all greater separations. Therefore, a core-softened fluid can have a negative thermal expansion coefficient and can become mechanically unstable at high density. The effective pair potentials of several liquid metals determined from experimental structure factor data exhibit core softening. Examples include AI, In, Mg;” Na, Al, Mg, K, Ca, Zn, Ga, Rb, Sr, In, Sn, Sb, Cs, Ba, T1. Pb, Bi;’* Three liquid metals (1 1) Hemmer, P.C.;Stell, G. fhys. Reu. Lerr. 1970,24, 1284. (12)Stell, 0.;Hemmer, P. C. J. Chem. fhys. 1972,56,4274. (13)Kincaid, J. M.;Stell, 0.;Hall, C. K. J . Chem. fhys. 1976,65,2161. (14)Kincaid, J. M.; Stell, G.; Goldmark, E. J. Chem. fhys. 1976,65, 2172. (15)Kincaid, J. M.; Stell, G. J. Chem. fhys. 1977,67,420. (16)Kincaid, J. M.;Stell,G. fhys. Lerr. A 3978. 65. 131. (17) Yokoyama, I.; Ono, S. J . fhys. F Mer. Phys. 1985, 15, 1215. (18) Hoshino, K.; Leung, C. H.; McLaughlin, I. L.; Rahman, S. M. M.; Young, W.H. J. Phys. F Met. fhys. 1985,17, 787.
Figure 5. Relationship between the continuous stability boundary (cgaf) within which the liquid phase can exist, and the locus of density maxima (aeg). Case in which density anomalies are confined to the metastable region. c is the critical point; b, the triple point; k, an isochore, and hb, bd, and bc, the sublimation, melting, and boiling curves. cga = spinodal; af = Kauzmann locus.
(Ga, Sn, Bi) expand upon freezing. No density anomalies have been reported for liquid metals. However, we are not aware of any volumetric measurement for supercooled liquid metals. Expansion upon freezing certainly suggests the possible existence of density anomalies. A negative thermal expansion coefficient and expansion upon freezing were both obtained by Stillinger and Weberz0 in their computer simulation study of a model coresoftened system: the Gaussian core fluid. In addition, the numerical angle-averaging of several asymmetric site-site water potentials yielded core-softened spherically symmetric potentials?’
Thermodynamic Implications The relationship between density anomalies and loss of stability upon supercooling which we have so far mentioned in connection with core softening is central to the phenomenological treatment introduced by Speedyzz*z3 to describe the stability boundaries of liquid water, and generalized by Debenedetti and DAntonioU for any fluid exhibiting density maxima Its essential features are illustrated in Figure 5 , which shows the continuous stability boundary (cgaf) of a liquid exhibiting density maxima (the locus of which is line aeg), along with the loci of solid-vapor (hb), liquid-vapor (k,c is the critical point), and solid-liquid (bd) equilibria. Figure 5 corresponds to a liquid which exhibits density maxima only in the metastable region, such as Si02.25 The picture is unchanged in its essentials in the case where aeg intersects the stable liquid region, except that the slope of line bd may then change from negative to positive at a pressure below which the liquid expands upon freezing, and above which it contracts upon freezing. Intersection with the locus of density maxima causes a change in the slope of the spinodal. At g, the liquid has attained its maximum tensile strength. Within the region bounded by aeg and the spinodal segment ag the liquid’s thermal expansion coefficient is negative. Points a and g are therefore the extremes of temperature and pressure within which contraction upon isobaric heating is possible. Line k is a liquid isochore; its slope changes sign along aeg, and it becomes tangent to the spinodal. Line af is a spinodal in both Speedy’s and Debenedetti and D’Antonio’s treatments.=-% However, because the liquid has a positive thermal expansion coefficient upon approaching this limit, it is more appropriate to view line af as a Kauzmann locus. (19)March, N.H. Can. J . fhys. 1987,65,219. (20)Stillinger, F. H.;Weber, T. A. J . Chem. fhys. 1978,68,3837. (21)D’Antonio, M.C.Ph.D. Thesis, Princeton University, 1989. (22) Speedy, R. J. J. fhys. Chem. 1982,86,982. (23) Speedy, R. J. J. fhys. Chem. 1982,86,3002. (24)Debcnedetti, P. G.;D’Antonio, M . C. AIChE J . 1988,34, 347. (25)Angell, C.A.; Kanno, H. Science 1976,193. 1121.
Debenedetti et al.
4544 The Journal of Physical Chemistry, Vol. 95, No. 11, 1991
U
-E
Figure 6. A pair potential with nearest-neighbor attraction and nextnearest-neighbor repulsion.
Experimental evidence in support of spinodal collapse in supercooled liquids includes measurements of the isobaric heat capacity of supercooled water and heavy water26*27 and of the isothermal compressibility of supercooled water.% Both properties showed anomalous increases; for water, the compressibilitywas fitted with a power law of the form KT = At7 [ t = ( T -T,)/T,] and yielded a spinodal temperature of 228 K at atmospheric pressure. An anomalous increase in the isobaric heat capacity of supercooled tellurium (Te) has recently been reported.29 The slope of the pressure-temperature projection of the homogeneous nucleation temperatures of supercooled water and heavy water exhibits high-pressure d i s c o n t i n u i t i e ~similar ~ ~ * ~ ~to that occurring at point a in Figure 5. In the case of heavy water, furthermore, the relation between the locus of density maxima and the homogeneous nucleation temperature locus appears to be similar to that shown in Figure 5 between the locus of density maxima and the spinodal.25 Both observations are consistent with the existence of strong pretransitional fluctuations due to the proximity to a spinodal having the same underlying shape. Angell and co-workers” have recently used the general picture described in Figure 5 (continuous stability boundary, relationship between loss of stability and density anomalies) to explain and interpret their measurements of water and aqueous solutions at highly negative pressures. Figure 5 , which we do not rederive here, follows only from thermodynamic consistency arguments and the assumption of an analytic Helmholtz energy along the spinodal. It shows that density anomalies provide a mechanism for a continuous stability boundary spanning the supercooled, superheated, and simultaneously superheated and supercooled, or subtriple*I conditions, and hence also for a supercooled liquid spinodal. From a continuum, macroscopic perspective, there appears a fundamental connection between loss of stability upon supercooling and the ability of a liquid to contract when heated isobarically. This agrees with essential aspects of the behavior associated with core softening. This connection, however, needs to be established more rigorously, starting from a microscopic model. We address this question in the next section via a deliberately simple model which we use not to predict the stability boundaries of a specific substance but to demonstrate that the behavior shown in Figure 5 follows also from microscopic considerations.
Microscopic Model and Implications We choose to formulate the problem in a lattice. The simplest case of core softening is a lattice fluid with single occupancy of (26) Rasmussen, D. H.;MacKenzie, A. P.; Angell, C. A.; Tucker, J. C. Science 1973, 181. 342. (27) Angell, C. A.; Shuppert, J.; Tucker, J. C. J . fhys. Chem. 1973, 77. 3092. (28) Speedy, R. J.; Angell, C. A. J . Chem. fhys. 1976,65. 851. (29) Angell, C. A. NATO Advanced Summer Institute on Correlation and
Connectivity in Biophysics, Corsica, France, July 1990. de Neufvillc, J. Private communication to C. A. Angell. (30) Kanno, H.;Speedy, R. J.; Aogell, C. A. Science 1975, 89, 880. (31) Green, J. L.; Durben, D. J.; Wolf, G. H.; Angell, C. A. Science 1990, 249, 649.
sites (hard core repulsion), nearest-neighbor attraction (-e; e > 0), and next-nearest-neighbor repulsion (At, A > 0; see Figure 6 ) . Stell and 00-workers studied phase transitions due to core softening by repulsive shoulders with and without long-range attractions.l1-l6 We are interested in density anomalies and stability limits in addition to phase behavior. In magnetic terminology, the present problem is equivalent to solving the king model with ferromagnetic nearest-neighbor and antiferromagnetic next-nearest-neighbor interactions in nonzero field. Previous work on aspects of the Ising problem with next-nearest-neighbor interactions includes that of Domb and Potts?* Fan and Wu,” Dalton and Wood,w Stephenson and Bett~,’~ L a n d a ~ , ~ ’ .Campbell ~* and S~hick,’~ Nauenberg and Nienhuis,40Meijer and Cunningham,” Taka~e,’z’~ Selke,44Binder and and Binder.& Domb and P ~ t t used s ~ ~series expansions to investigate the dependence of the critical temperature, energy, and entropy upon the sign and magnitude of X for a square, two-dimensional lattice. Fan and Wu3’ investigated the square Ising lattice with nextnearest-neighbor interactions in zero field and obtained expressions for the free energy, energy, specific heat, and critical temperature. Dalton and Wood34used series expansions to investigate the dependence of the critical temperature, critical exponents for the magnetization and susceptibility, critical energy, and critical entropy upon the ratio of next-nearest-to-nearest-neighbor interactions (both ferromagnetic) in two and three dimensions. Gibberd35investigated the two-dimensional problem on a square lattice and calculated the critical temperature and magnetization critical exponent for ferromagnetic nearest- and next-nearest interactions. Stephenson and B e t t ~calculated ~~ the value of X for which the phase transition disappears, as well as the X-dependence of the critical temperature for a variety of lattices. L a n d a used ~ ~ Monte ~ ~ ~ Carlo ~ simulations to calculate the critical temperature, as well as the behavior of the energy, specific heat, susceptibility, and magnetization in the vicinity of the critical point for the two-dimensional, zero-field next-nearest-neighbor problem. Campbell and Schick39studied the plane triangular lattice with nearest-neighbor repulsion and next-nearest-neighborattraction in the Bethe-Peierls approximation. Nauenberg and Nienhuis“ used the renormalization group method to obtain the Xdependence of the critical temperature for a square Ising lattice in the range -1 C X C 1. Meijer and Cunningham4I compared the BethePeierls results of Campbell and Schick with their own Monte Carlo calculations. T a k a ~ used e ~ Monte ~ ~ ~ Carlo simulations and series expansions to explore the X-dependenceof the critical temperature (as well as of the energy and the magnetization in the vicinity of the critical point) for a square Ising lattice with repulsive next-nearest-neighbor interactions. SelkeU performed Monte Carlo calculations to study the threedimensional Ising model with next-nearest-neighbor repulsion along one direction only. Binder and Landa~,4~ using Monte Carlo calculations, obtained the phase diagrams and critical behavior of two-dimensional antiferromagnets with next-nearest-neighbor interactions in nonzero field. BindeP used Monte Carlo methods to investigate the energy, free energy, entropy, and magnetization of the facecentered cubic Ising antiferromagnet with next-nearest-neighbor interactions in nonzero field. The present model is also to be viewed in the context of several lattice-based studies of water and water-like behavior. These (32) (33) (34) (35) (36) (37) (38) (39) (40) (41) (42) (43) (44) (45) (46)
Domb, C.; Potts, R. B. froc. R. Soc. (London) 1951, A210, 125. Fan, C.; Wu,F. Y . fhys. Rev. 1969, 179, 560. Dalton, N. W.;Wood,D. W. J. Murh. fhys. 1969, IO, 1271. Gibberd, R. W. J . Murh. fhys. 1969, 10, 1026. Stephenson, J.; Betts, D. D. fhys. Rev. B 1970. 2, 2702. Landau, D. P. J . Appl. fhys. 1971, 42, 1284. Landau, D. P. fhys. Reu. B 1980, 21, 1285. Campbell, C. E.; Schick, M. fhys. Rev. A 1972, 5, 1919. Naunberg, M.; Nienhuis, B. fhys. Rev. Let?. 1974, 33, 944. Meijer, P. E.;Cunningham, G. W. fhys. Rev. B 1977, 15, 3436. Takase, S. J . fhys. Soc. Jpn. 1976.40, 1240. Takase. S. J. fhvs. Soc. J m . 1977. 42, 1819. Selke, W. 2.Phis. B 197$, 29, 133. Binder, K.; Landau, D. P. fhys. Rev. B 1980, 21, 1941. Binder, K. 2.fhys. B 1981, 45, 61.
Spinodal Curve of Some Supercooled Liquids
The Jcwnal of Physical Chemistry, Vol. 95, No. 11, 1991 4545
include the work or' Bell, Lavis, and co-worker~,)~-~~ Meijer and co-worker~,~)-~~ and Fleming and G i b b ~ . ~Bell , ~and ~ La~is~~" introduced two plane (two-dimensional) models, the i n t e r ~ t i t i a l ~ ~ and the orientable moleculea models on a triangular lattice, which they solved through a first-order (quasichemical) approximation.58 The Bell-Lavis models are both characterized by a competition between a low-energy open structure and a high-energy closepacked structure brought about by the possibility of forming bonding and nonbonding interactions. In the interstitial model, there is a division of the original lattice into two sublattices (honeycomb and interstitial sublattices), with nearest-neighbor molecules on the honeycomb lattice forming bonded interactions (lower energy), while nearest-neighbor pairs belonging to different sublattices are unbonded. The open structure is then the completely filled honeycomb sublattice with an empty interstitial sublattice, both sublattices being filled in the close-packed structure. In the orientable molecule model, each molecule has three possible bonding directions and two possible orientations. If nearest-neighbor molecules have bonds pointing toward each other, bonding results (lower energy). A competition results between an open structure characterized by a honeycomb arrangement with 1/3 of the sites vacant, and the close-packed structure in which the lattice is completely filled. In both the interstitial and orientable molecule models, the open structure is stable at low temperatures and presdures, giving rise to water-like behavior (density maxima). Bell49extended these ideas to three dimensions. He considered a water model consisting of a bodycentered lattice in which molecules can form tetrahedrally oriented bonds with their nearest-neighbors, when appropriately oriented. In addition, he introduced a (repulsive) energy penalty against close packing. Bell solved the problem in the quasichemical approximation. Density anomalies arise as a consequence of the competition between an open tetrahedrally bonded structure on a half-filled lattice and a close-packed network of interpenetrating tetrahedrally bonded structures, which Bell compared to ice I(c) and ice VII, respectively. Young and Lavism and Southern and Lavis5' used real-space renormalization to investigate the critical behavior of the two-dimensional Bell-Lavis model. The threedimensional model of Bell49 was solved in the zeroth-order (Bragg-Williams) approximation by Bell and Saltes2This caused the disappearance of density anomalies in the liquid phase, but enabled the incorporation of two long-range ordered (icelike) solid phases, of which the open one melted into a denser liquid. Meijer and worker^^"^ used the cluster variation method to solve Bell's three-dimensional water model and studied its phase behavior and stability limits. They pointed out explicitly the importance of next-nearest-neighbor repulsions. Fleming and GibbsS6vs7used a perturbation scheme to calculate the free energy and several thermodynamic properties of a model similar to Bell's three-dimensional model. They obtained density maxima, as well as a water-like vapor-liquid coexistence curve, with a maximum in the saturated liquid density. The present treatment is not aimed at modeling water. Rather, we seek to establish a general microscopically based connection between density anomalies and loss of stability upon supercooling. The phenomenologicalarguments on which Figure 5 was originally b a ~ e d ~ *will - ~ ' thus be given a molecular basis. Core softening provides the simplest and most general microscopic model consistent with the complex behavior described in Figure 5 . The lattice gas with attractive nearest-neighbor and repulsive next(47) Bell, G. M.; Lavis, D. A. J . Phys. A: Gen. Phys. 1970, 3, 427. (48) Bell, G. M.; Lavis, D. A. J. Phys. A: Gen. Phys. 1970, 3, 568. (49) Bell, G. M. J . Phys. c: Solidsrure Phys. 1972, 5, 889. (50) Young, A. P.; Lavis, D. A. J . Phys. A: Murh. Gen. 1979, 12, 229. (51) Southern, B. W.; Lavis. D.A. J. Phys. A: Murh. Gen. 1980,13,251. (52) Bell, G. M.; Salt, D. W. J . Chem. Soc., Furuduy Truns. 2 1972,72, 76.
(53) Meija, P. H. E.; Kikuchi. R.; van Royen, E. Physicu 1982,115A, 124. (54) van Royen. E.; Meijer, P. H. E. Physicu 1984, 127A, 87. (55) van Royen, E.; Meijer, P. H. E. Physicu 1986, 139A, 412. (56) Fleming, P. D.; Gibbs, J. H. J. Sruf. Phys. 1974, IO, 157. (57) Fleming, P. D.; Gibbs, J. H. J . Slur. Phys. 1974. 10. 351. (58) Guggenheim, E. A.; McGlashan, M. C. Proc. R. Soc. 1951,206,335.
BB (b)
(0)
Figure 7. Formation of open (a) and close-packed (b) structures in a two-dimensional lattice fluid in the low-temperature limit.
Figure 8. Pressure dependence of the low-temperature limits of the enthalpy and energy changes associated with the transition from opento-close-packed structure in two dimensions. M is the number of molecules in the lattice. Ae and Ah are the energy and enthalpy changes per molecule.
nearest-neighbor interactions (Figure 6) is the simplest case of core softening. It will now be shown that the competition between nearest-neighbor attraction and next-nearest-neighbor repulsion is enough to give rise to density anomalies and hence to the possibility of loss of stability upon supercooling. The underlying mechanism, as in the case of water, will be shown to be the competition between open structures which can 'melt" into denser, high-energy, close-packed configurations through the input of thermal or mechanical energy. To see this, consider the perfectly ordered two-dimensional arrangements shown in Figure 7. In the open structure, there are no next-nearest-neighbor (diagonal) interactions, and half of the sites are occupied. Then, with N denoting the number of sites in the lattice, and M the number of molecules,
H = -Mt
2M=N
(18)
E = -Me
(19)
[ 2y
+ 2PvoM = M t -- 1 1
(20)
where E and H a r e the configurational energy and enthalpy, P is the (two-dimensional) pressure, and uo, the two-dimensional volume per site. For the compact structure, on the other hand,
M=N E = 2M4X - 1)
(21) (22)
1
where h is the height of the repulsive barrier in units of the attractive well (see Figure 6). Therefore, at absolute zero, the close-packed structure will be stable if Po0 e > 2(x -
1)
4546
Debenedetti et al.
The Journal of Physical Chemistry, Vol. 95, No. 11, 1991
TABLE I: Configurations, Energies, and Degeneracies for the One-Dimensional Model configurn no.
energy
degeneracy
1, 0-0-0
(A - 2)t
2,o-0-0
-E
1 2
3,o-0-0
At
4,o-0-0 5,o-0-0
0 0 0
6,O-0-0
1 2 1 1
$i"
P3 P2(1 - P ) P2(1 - P ) P(1 - P)2 P(1 - d2 ( 1 - PI3
TABLE 11: Configurations, Energies, and Degeneracies for the Twoand Three-Dimensional Models configurn no. energy degeneracy 9; 1,
0
-
1
(1
4
P(1
-d3
Z
4
Y
1f.X < 1/2, the denser structure will always be favored at equilibrium when T = 0. Furthermore, the close-packed structure is energetically unfavorable only if X > 1/2. These considerations are illustrated in Figure 8, where the vertical axis is the enthalpy (or energy) change associated with the transition open-to-closepacked as a function of pressure, at T = 0. The open structure X is stable at pressures lower than that at which the enthalpy line Figure 9. Tessellations for the quasichemical solution of the one-, two-, intersects the horizontal axis. Under these conditions, isobaric and three-dimensional problems. The tessellations preserve nearestheating leads to "melting" into denser configurations. This neighbor contacts. possibility of forming open structures stabilized by next-nearest-neighbor repulsion is obviously not limited to two dimensions. In a one-dimensional lattice, the ordered structure with no nearest-neighbor interactions occurs with pairs of neighboring sites alternatively vacant and full, and in three dimensions, by stacking planes such as the one shown in the left-hand side of Figure 7 separated by empty planes. In both cases, the dense structure is the completely filled lattice. The minimum A-value below which the denser structure is always favored at T = 0 is also 1/2 in one and three dimensions; the minimum pressures (in units of c/oo) Figure 10. Bond-moving step required to preserve next-nearest-neighbor above which the denser structure is favored at T = 0 are (A - 1/2) interactions upon tessellation. and 4(X - 1 /2)/3, respectively. The preceding analysis is identical with the one used by Bell and Lavis47*48 and by Bell49to explore of cells is determined by requiring that the total number of the low-temperature stability of open and close-packed structures nearest-neighbor contacts be conserved upon tessellating the lattice in their two- and three-dimensional models. with elementary cells The equation of state and chemical potential for the lattice gas Ncells = W 2 r l (26) with attractive nearest-neighbor and repulsive next-nearestneighbor interactions were obtained through the application of where z is the coordination number (2,4,6 in one, two, and three the quasichemical approximati~n;~~ the technical details of the calculations follow closely the approach of Bell and c o - ~ o r k e r s . " ~ ~ ~dimensions, respectively), and 7, the number of nearest-neighbor contacts in a cell (2, 4, 4 in one, two, and three dimensions, In essence, the method involves choosing an elementary cell and respectively). Tessellations for the one- two- and three-dimensional approximating the number of configurations accessible to M cases are shown in Figure 9. Note that cells share sites, but not molecules on N lattice sites ( p = M / N ) as interactions (contacts). The number of next-nearest-neighbor contacts is N, 2N, and 6N in the linear, square, and simple cubic lattices. The tessellations shown in Figure 9, however, conserve only nearest-neighbor contacts, giving rise to N/2, N, and 3N/2 next-nearest-neighbor contacts in one, two, and three dimensions, respectively. Therefore, where Ncellsis the number of cells into which the lattice has been the tessellation implies a bond-moving step (shown in Figure 10 divided, + i ( p , T ) is the probability of occurrence of the ith confor the two-dimensional case), such that X/X,,,,, = 2, 2 , 4 in one, figuration of the elementary cell, wi is the degeneracy associated two, and three dimensions, respectively (e.g., one-dimensional with the ith configuration, m is the number of distinguishable results for X = 1 would represent the quasichemical approximation configurations which can be assigned to an elementary cell, and to the properties of a "true" system with X = 1/2). Since our goal Qo is a normalization constant. Configurations, degeneracies, is not to model a particular true fluid, the bond-moving step is energies, and +i" [ = \ t i ( p , ~ ) ]values for the one- two- and threeonly formal. Two-dimensional phase envelopes for several values dimensional problems (the latter two on square and simple cubic of A, were obtained via grand canonical Monte Carlo simulations lattices, respectively) are given in Tables I and 11. The number
The Journal of Physical Chemistry, Vol. 95, No. 11, 1991 4547
Spinodal Curve of Some Supercooled Liquids
on a 40 X 40 lattice with periodic boundary conditions; they were in good qualitative agreement with the corresponding quasichemical phase envelopes for X = 2X,,, except in the critical region. The normalization constant in eq 25 is calculated by requiring that Q(p,T)be exact in the high-temperature limit. In this case, molecules are distributed randomly in the lattice, and the number of configurations is N!/((Np)![N(l- p ) ] ! ]
which gives N-' In 51, = y [ p In p + ( 1 - p ) In ( 1 - p ) ] (28) where y is 112, 1, and 2 in one, two, and three dimensions. The probabilities $i must satisfy the material balance constraints l m P = -Cwi$ini (29) ai-1
where
P* = P u o / c P*
= r/e
TI = k T / c x = exp(-c/kT) (44) In the above equations, r is a parameter related to dimensionless
density p (fractional coverage) by the relations iIX-(2+A) y3[4x-(.'+A) 21 r where a is the number of sites per cell and n,, the number of occupied sites in the ith configuration. Equations 29 and 30 imply the normalization condition
P =
+
+
+ + + + 3r + 9 x A
r-lx-(2+A) !1~[4x-(I+~) 21
(45)
1 + #/2[2x-A/2 + x('+W)]+ 3r
+ x(V2-1)9/2 r-1/2x(V2-l)+ 4 + 2rI/2[2x-V2 + x(I+A/2)]+ 4r + x(A/2-1)9/2
(46) where eq 45 is valid in one dimension and eq 46 in two and three dimensions. The quantity r is defined in the Appendix (eqs A.6, A.11). The structure of the equations is therefore P* = P*(p,r,T) (47)
The energy is given by
and the canonical partition function, Q,by A(N,p,T) = -kT In Q(N.p,T) =
P = p(r,T)
-kT In (o(p,T) exp[
(33)
where A is the Helmholtz energy, and we have replaced the true partition function by its (still undetermined) maximum term. By use of eqs 25 and 32, the Helmholtz energy becomes ( a = N-'A I pa'; a' I M-'A)
4PlT) = m
- r W p In P + ( 1 - P ) In ( 1 - P ) I + GCwi$i(a, + kT In $J 1-1 (34) where 6 is 112, 112, and 314 in one, two, and three dimensions, respectively, and is a (still undetermined) function of density and temperature. We minimize the Helmholtz energy by finding the probability distribution ($i) that maximizes (In 51 - E / k T ) . Then, as shown in the Appendix, we obtain, for the one-dimensional fluid's pressure and chemical potential
P * =~ - *Tl 1 n- P (
P =
~ )
(48) = cL*(p,r,T) (49) The pressure and chemical potential are given as functions of density and temperature by two pairs of parametricequations ((47) and (48) for the pressure, (49) and (48) for the chemical potential). r is obtained numerically from eq 48 for any given ( p , T ) . Figure 1 1 shows isobars for the one-dimensional fluid. Note that some curves have positive slope at low temperature, indicating that the fluid contracts when heated at constant pressure, and therefore has a negative thermal expansion coefficient. The region where density anomalies occur is bounded from above and from below in tcmperature, pressure, and density. Qualitatively very similar behavior was obtained by Bells9 via a one-dimensional model in which next-nearest-neighbor molecules can form bonded (low-energy) interactions when the intermediate site is empty. For pressures lower than 0.2, the first-order approximation predicts that as T tends to 0 the fraction of nearest-neighbor contacts approaches 5%, and the fraction of next-nearest-neighbor contacts vanishes (Figure 12). In the quasichemical approximation, this open structure is realized by tessellating the lattice with cells having configuration 2 (Table I) but alternating in direction (e',f,f',f,e',f,f',f,e', ...; where e denotes empty, f, full, and primes denote the beginning and end of a cell). For pressures above 0.067 (Figure 13a), this structure gives way, upon isobaric heating, to arrangements with a progressively larger fraction of completely filled cells (configuration 1; Table I), and hence to a density increase. For pressures below 0.067 (Figure 13b), the open structure gives way, upon isobaric heating, to arrangements with a progressively larger fraction of empty cells (configuration 6; Table I), and hence to a density decrease. Note the invsrsion in the relative stability of configurations 1 and 6 between parts a and b of Figure 13, which results in the disappearance of density P*
(59) Bell, G. M.J. Morh. Phys. 1969, 10, 1753.
Debenedetti et al.
4548 The Journal of Physical Chemistry, Vol. 95, No. 11, 1991 1. o
1.Oh
I
1-d
Q
u
0
T' Figure 11. Density-temperature relationship for the one-dimensional fluid. Labels on each curve are dimensionless pressures. .7
1
h=0.7 P*=0.05
h 0
1 -d
a
2a
)\=0.7
.41&
3
0
n .5
0
nn,0.15
0
0
.2
OO
Figure 12. Fraction of active nearest-neighbor (nn) and next-nearestneighbor (nnn) contacts for the one-dimensional fluid. Labels on the curves are dimensionless pressures.
.BO I
anomalies at low pressure. This is illustrated in Figure 14, which shows loci of density maxima for different values of A. The fractions of active nearest- and next-nearest-neighbor contacts (fnn, fnnn) were calculated as follows
.76
m
I- 1
(50)
m
fnnn = 2Co,$,nnnI/s 1=1
.2
rr*
.3
.5
.4
Fig& 13. Fraction of cells in configurations 1, 2, and 6 (see Table I) for the one-dimensional fluid. P* = 0.15 (a, top); 0.05 (b, bottom).
T*
fnn = CwI$,nni/q
.1
(51)
where 7 has already been defined (see eq 26) and nn, and nnnI represent the number of active nearest- and next-nearest-neighbor contacts in the ith configuration. From Figures 11-14 it can be seen that indeed core softening gives rise to density anomalies in one dimension and that this behavior results from the possibility of forming an open structure which can be disrupted into closer packed configurations by input of thermal or mechanical energy. The normalization constant in eq 25, however, is exact in the high-temperature limit. This means that the quasichemical approximation becomes progressively less accurate a t low temperatures. Important aspects of the behavior shown in Figures 11-14 as T tends to 0 need to be critically examined. In the first place, the perfectly regular open structure described above has a true fractional occupancy of 3/4, whereas the quasichemical approximation predicts a value of 2/3 (100% of the cells are in configuration 2, whose fractional occupancy is 2/3). Secondly, there are no next-nearest-neighbor contacts in configuration 2, whereas 1/2 of the next-nearestneighbor contacts in the tessellation e',f,f',f,e', ..., are active. These discrepancies arise from applying a cell-based configuration counting to a perfectly ordered lattice. Such inconsistencies disappear at higher temperatures; a t low temperatures, they become progressively less severe in higher dimensional lattices, as
I
1
I
\I
i;/,$d!i .66
.06
,lo
.14
,I 8
k1.2
.22
k1.5
.26
30
T'
Figure 14. Loci of density maxima for the one-dimensional fluid for different values of X. will be shown below in connection with the two-dimensional fluid. Figure 15 shows the coexistence region and spinodal curves for the two-dimensional fluid for different values of A. There is no phase transition for X 2 1. The critical temperature for X = 0 is 0.696 (in units of elk); the exact value for the two-dimensional Ising model is 0.567 27, and the prediction of the quasichemical approximation using a two-site linear cell, 0.721 347. Stable density maxima occur only for X 2 1. In the low-temperature limit, the minimum and maximum pressures between which density maxima are possible are 0.125z(X - 1) and 0.25z(X - l ) , respectively. These bounds are valid in two and three dimensions. Figure 16a shows the P* = 0.23 isobar for the two-dimensional fluid with X = 1.2. Since this pressure exceeds 0.2, there are no density maxima. In the low-temperature limit the lattice is completely filled, and the fraction of active nearest- and next-
The Journal of Physical Chemistry, Vol. 95, No. 11, 1991 4549
Spinodal Curve of Some Supercooled Liquids
1.
.70
.65
A = 1.2 P*=O. 2-d 18
h
.60
C
'
(d
.55
a
I
¶
0
.45
:q; I
.10
.
0
0
.1
~~,
,
.2
.3
,
I
, , ,
,
.5 P*
.6
.4
,
I
,
.7
I
.8
-11'
T'
.9 1.0
1 .o
Figure 15. Binodal and spinodal curves for the two-dimensional fluid for different values of A.
2-d A = 1.2 P*=O. 1
t
e h 0
5
h=1.2
.8
V
.6
0 &
20
.4
.I
4
V
2
.2
w
0 '
.2
.4
.6
.8
1.0
T* Figure 17. Density, fractional occupancies of the various typcs of cells (see Table 11), and fraction of active nearest-neighbor (nn) and nextnearest-neighbor (nnn) contacts for the twedimensional fluid. P = 0.18 (a, top); 0.10 (b, bottom). All nn, but only active nnn contacts shown in squares.
T*
0 '
.2
.4
.6
.8
1.0
T* Figure 16. Density, fractional occupancies of the various types of cells (see Table 11). and fraction of active nearest-neighbor (nn) and nextneerest-neighbor (nnn) contacts for the twedimensional fluid. P = 0.23 (a, top); 0.20 (b, bottom). All nn, but only active nnn contacts shown
in squares.
nearest-neighbor contacts approaches 100% (configuration 6). Upon increasing the temperature isobarically, the lattice at first becomes progressively populated with cells having configuration 3. The fraction of these cells peaks slightly below T* = 0.2and decreases thereafter. Figure 16b shows the P* = 0.2 isobar. In two dimensions and for A = 1.2,this is the upper limit of pressure above which density maxima no longer occur. Note that in the low-temperature limit 75% of the cells are in configuration 3, and 25% in configuration 6, giving an overall fractional coverage of 62.5%. Below P* = 0.2 (Figure 17a). the lattice is populated completely by cells having configuration 3, giving rise to a fractional coverage of 1/2 and fractions of active nearest- and next-nearest-neighborcontacts of 50% and 25% respectively. This
is exactly the open structure shown in Figure 7. It can be generated by tessellating the lattice with cells in configuration 3 with orientations which alternate from row to row. Note that there is perfect agreement between the cell-based predicted fractions of occupied sites and of active nearest- and next-nearest-neighbor contacts (1/2,1/4, 0) and the true values, in contrast to the one-dimensionalcase. Density maxima occur as a consequence of the gradual population of the lattice by completely filled cells (configuration 6) in agreement with the discussion of Figure 7. Below P* = 0.1 density anomalies disappear due to the inversion in the relative stability of configurations 1 and 6 (Figure 17b). Although in the low-temperature limit the structure of the lattice is identical for P* = 0.18 and 0.1 (Figure 17, a and b), isobaric heating at or below P = 0.1 now leads to a preferential population with empty cells (configuration l), and hence density anomalies disappear. Figure 18 shows the coexistence region for the three-dimensional fluid, for different values of A. The critical temperature for A = 0 is 1.225;the best accepted numerical value for the three-dimensional Ising model is 1.12428,and the prediction of the quasichemical approximation using a two-site linear cell, 1.2333. As A is increased above 0.9, a second phase transition occurs. This is an example of the type of phase transition due to core softening investigated by Stell and co-workers.11-16In the present case, both transitions have the same critical temperature and are symmetric about p = 1 /2,the triple-point liquid density. The phase diagram is shown in pressure-temperature coordinates in Figure 19. It can be seen that the effect of increasing A is to decrease the critical temperature, decrease the critical pressure of the low-density transition, and increase the critical pressure of the high-density transition. The appearance of a second phase transition is significant in the context of core softening, because it allows the liquid to become supercooled with respect to a denser phase. The nature
Debenedetti et al.
4550 The Journal of Physical Chemistry, Vol. 95, No. 1I, 1991 12.5
1,'
0'
.1
.2
.3
.4
.5 P'
.6
.7
.8
.9 1.0
Figure 18. Phase envelopes for the three-dimensional fluid for different values of A.
,121
i
.10
1
-12.5
0
d
I
10.0
!
I
.05
.10
.15 T'
.20
.25
.30
.35
Figure 20. Continuous stability boundary (cgaf) and locus of density maxima (aeg) for the three-dimensional liquid; A = 1. b and c are triple and critical points. cn and m are portions of the supercooled vapor and superheated high-density fluid spinodals.
associated with the second fluid-fluid transition. It is a consequence of core softening but bears little relation to a limit of stability associated with a solid-liquid transition.
3-d
_-01k%j---zm
'
'
..A5
T* Figure 19. Pressurttemperature Drojection of the phase diagram of the three-dimensionalfluid for diffei. I t values of A.
of the latter, however, plays no role in the mechanical stability (or lack thereof) of the supercooled phase. Hence, it will be possible to test the stability of a supercooled core-softened fluid with considerably more generality than would at first appear possible from the fluid-fluid character of the high-density transition. Figure 20 shows the relationship between the phase boundaries, stability limits, and density anomalies for the three-dimensional core-softened fluid with X = 1. The lettering corresponds to that of Figure 5. The superheated liquid spinodal of the low-density transition (cg) exhibits a change in slope upon encountering the (entirely metastable) locus of density maxima (aeg). Curve af is the spinodal along which the supercooled core-softened fluid becomes unstable. Portions of the superheated high-density fluid spinodal (m, corresponding to transition bd) and the supercooled vapor spinodal (cn, comsponding to transition bc) are also shown. The similarity between Figures 5 and 20 is obvious. The behavior shown in Figure 20 occurs for 0.99 IX I1.001. Since the maximum pressure at which density maxima occur in the low-temperature limit is (z/4)(X - 1) in this model,@ the high-pressure intersection of density maxima and spinodal loci (point a, Figure 5 ) occurs at pressures very close to zero. Except for this limitation, all the other features shown in Figure 5 are reproduced by the lattice model. The complex predicted behavior (density maxima in conjunction with a tensile strength maximum or reentrant spinodal, loss of stability upon supercooling) is therefore in qualitative agreement with phenomcnological thermodynamic arguments and experiments. The most significant aspect of Figure 20 is the reentrant spinodal (cag) and its relationship to density anomalies. Curve af is simply the spinodal (60)Raghavan, V. S. MS Thesis, Princeton University, 1990.
Conclusion There are two possible absolute limits to the extent of supercooling to which liquids can be subjected: the Kauzmaiin temperature and the spinodal condition. In this paper we have investigated the necessary conditions for spinodal collapse. A consistent picture of this phenomenon is gradually emerging. In this view, spinodal collapse is possible only for liquids capable of contracting when heated isobarically. Microscopically, this occurs via the formation of open structures which are stabilized by repulsions and which can be collapsed into denser arrangements through the input of thermal or mechanical energy. Within the context of a pairwise additive and spherically symmetric description, both negative thermal expansion coefficients and loss of stability at high density can be explained by core softening. Effective pair potentials for several liquid metals exhibit core softening. To the best of our knowledge, no stable density maxima for liquid metals have been reported. Nevertheless, Ga, Bi, and Sn expand on freezing. Density measurements in the supercooled region as well as heat capacity and isothermal compressibility data for these elements should shed considerable light on the general problem of the stability boundaries of supercooled liquids which we have addressed in this paper. The possibility of spinodal collapse in supercooled Ge and Te has recently been suggested, evidence of heat capacity anomalies in supercooled Te having recently been reported.29 A lattice model with nearest-neighbor attraction and nextnearest-neighbor repulsion exhibits density ancmalies in one, two, and three dimensions, and a reentrant spinodal in three dimensions. This model is deliberately simple. Consequently, the phase with respect to which the suptrcooled liquid becomes metastable is not solid in the extended first-order approximation. However, the purpose of the model is not to mimic the behavior of a particular real fluid, much less to predict the location of stability boundaries of, say, water. Rather, it shows how a system capable of forming stable open structures at low temperatures and pressures can exhibit density anomalies and become mechanically unstable when supercooled, leading to a reentrant spinodal. It is hoped that a combination of more sophisticated models, together with the availability of data on the behavior of supercooled liquid metals such as Ga, Sn, Bi, Te, and Ge, will contribute toward an improved understanding of the important yet currently poorly understood problem of determining the absolute boundaries within which the liquid state of matter can exist. Acknowledgment. We gratefully acknowledge the financial support of the National Science Foundation (Presidential Young
The Journal of Physical Chemistry, Vol. 95, No. 11, 1991 4551
Spinodal Curve of Some Supercooled Liquids Investigator Award CBT-8657010 to P.G.D.), the Camille and Henry Dreyfus Foundation (1989 Teacher-Scholar Award to P.G.D.), and the US.Department of Education for a Fellowship for Graduate Assistance in Chemical Engineering to SSB. We are indebted to Ariel Chialvo for the simulation results of Figures 2 and 3. Appendix
We first prove that, for a potential such as that defined by eq 12, the virial is a monotonically increasing function of density above the temperature at which the second virial coefficient vanishes (the Boyle temperature). To this end, we write the virial equation of state
P = pkT(1
In one dimension, we choose I) to obtain
and $4 as independent (see Table
where n,,x , and el have already been defined. Then, using eqs 29 and 30 P = $30
1 - p = $3@
(A.7)
and, therefore,
+ pE + p2C + ...) = pkT + pkT(pE + p2C + ...) (A. 1)
Comparison with eq 4 yields $ = 6kT(pE
+ p2C + ...)
(A.2)
and, therefore,
I(
") aP
T
= EkT + 2pkTC + ...7EkT
(A.3)
Equation A.3 shows that the sign of the slope of an isotherm in the virial-density plane is that of the second virial coefficient in the zero density (ideal gas) limit. Therefore, the initial slope of the virial-density curve decreases with density below the Boyle temperature and increases with density at higher temperatures. Since isotherms necessarily have positive curvature in the virial-density plane for a potential such as that defined by eq 12 (see Figures 1,2, and 3), it follows that above the Boyle temperature the virial is a monotonically increasing function of density at a given temperature. We now give the derivation of the equilibrium values for $,(p,T), the equation of state, and the chemical potential. The goal is to find the probability distribution of cells that maximizes the Helmholtz energy, In fl - E/kT. To this end, we write
(A.14) The equation of state for M molecules in N lattice sites is obtained from the identities
and the chemical potential, from b=-G =-=
where gk and $, are the two independent probabilities. The existence of two independent probabilities follows from the two material balance conditions (eqs 29 and 30). The partial derivatives of the independent probabilities with respect to $i are obtained from solving eqs 29 and 30 for $k and In addition, from eq 34 we obtain m
= bCw,(c, p,T
'-I
+ kT(l + In $,)I
(A.5)
A+PV
-(a 1
+ Pvo)
(g)
(A.16)
M M P T where uo is the volume of a unit cell, and p = M/N.Note that a is related to the Helmholtz energy per molecule (a') by a = pa' (A.17) The partial derivative in eq A.16 is evaluated as
where ($1 denotes constancy of all $,. Equations 35-40 then follow after straightforward, if tedious, algebra. Equations A.8-A. 10 yield eq 45; eqs A.12-A.14 yield eq 46.
