J. Phys. Chem. 1995, 99, 3781-3792
3781
A Lattice Model of Network-Forming Fluids with Orientation-Dependent Bonding: Equilibrium, Stability, and Implications for the Phase Behavior of Supercooled Water Steven S. Borick and Pablo G. Debenedetti” Department of Chemical Engineering, Princeton University, Princeton, New Jersey 08544-5263
Srikanth Sastry Physical Sciences Laboratory, Division of Computer Research and Technology, National Institutes of Health, Bethesda, Maryland 20892 Received: September 12, 1994@
We use a lattice model with orientation-dependent interactions to study network-forming fluids, such as water and silica. Bonds can form between pairs of molecules that have correct mutual orientation and correct separation; they can be weakened by the presence of other molecules sufficiently close to a bonded pair. These interactions give rise to competition between bonded states of low energy, density, and entropy and nonbonded states of high energy, density, and entropy. By suitable choice of parameters, the mean-field solution of the model yields the two different scenarios (stability limit conjecture and second critical point) that have been proposed to explain the anomalous behavior of supercooled water. The model’s generality suggests that similar behavior can occur in other network-forming fluids.
1. Introduction The purpose of this work is to investigate the phase behavior and limits of stability of a microscopic model of fluids that form networks held together by means of orientation-dependent bonds. The model fluid becomes denser when heated at constant pressure, like water and silica. We study in detail the relationship between such density anomalies and stability limits and its dependence upon model parameters such as bond strength and the number of distinguishable orientations that a molecule can adopt. Water and silica, the two major constituentsof the earth, form tetrahedrally coordinated networks. Water forms a network in which each oxygen atom is hydrogen-bonded to four other oxygen atoms through an 0-H-0 bridge. The oxygen atoms occupy the vertices of a tetrahedron in the resulting threedimensional arrangement. At low enough temperatures and pressures, this arrangement is regular and is called hexagonal ice (Ih). In silica, each Si atom is bonded to four other Si atoms through a Si-0-Si bridge to form a three-dimensional, tetrahedrally coordinated network. Thermal disruption causes such open networks to collapse into denser and more disordered arrangements. Macroscopically, this translates into a negative thermal expansion coefficient: the liquid contracts when heated isobarically. Such density anomalies can be observed in stable and metastable (supercooled) water and also in supercooled silica.2 Although the mechanism of bond formation is quite different in these two liquids, the competition between low-density, energetically favored and high-density, entropically favored configurations is common to both. Because the network limits the ability of molecules to move freely, its thermal or mechanical disruption causes a gain in rotational and translational entropy. Thus, an important feature that is common to these two network-forming liquids is that local regions of high density can be more disordered than open, locally tetrahedrally coordinated regions. This greater disorder (both rotational and translational) at higher
* To whom correspondence should be addressed. Abstract published in Advance ACS Abstracrs, February 15, 1995.
0022-365419512099-3781$09.0010
density is energetically unfavorable relative to network formation. Despite water’s abundance and its biological, ecological, technological, and economic significance, knowledge of its thermophysical properties is far from complete. In particular, the behavior of supercooled water is among the most intriguing open questions in the physics of liquids. Interest in the subject dates back to the early 1970s, when new experimental techniques made possible measurements on highly supercooled samples (e.g., -38 OC at atmospheric pressure). Experiments revealed surprising increases in the isobaric heat isothermal compressibility,6*’and the magnitude of the thermal expansion coefficient8-12 at low temperatures, with no comparable increases in the isochoric heat capacity. This led Speedy to postulate his stability limit ~ o n j e c t u r e . ~According ~J~ to this theory, liquid water has a continuous stability boundary (spinodal curve) bounding the superheated, supercooled, and subtriple (simultaneously superheated and supercooled) states. The superheated liquid spinodal reaches a minimum pressure (maximum tension) and retraces back toward positive pressures at low temperatures. The retracing occurs when the spinodal contacts the locus of density maxima. Speedy’s theory predicts that the isobaric heat capacity, the isothermal compressibility, and the magnitude of the thermal expansion coefficient should diverge as water is supercooled and that an absolute (thermodynamic) limit exists beyond which supercooled water cannot exist. At atmospheric pressure, for example, liquid water would become unstable below -45 “C. Debenedetti and co-workers15J6 generalized Speedy’s theory, using thermodynamic consistency arguments to conclude that the continuous stability boundary scenario should be relevant not just to water but to any liquid capable of expanding when cooled at constant pressure. Recently, a very different interpretation of the behavior of supercooled water has been proposed. In a molecular dynamics study of ST2 water,” Poole, Sciortino, Essman, and Stanley1*J9 (henceforth PSES) observed that the P-T projection of the locus of density maxima, which is negatively sloped at positive pressures, becomes infinitely sloped and then positively sloped 0 1995 American Chemical Society
Borick et al.
3702 J. Phys. Chem., Vol. 99, No. 11, 1995
with decreasing pressure and does not intersect the spinodal. Thus, upon application of sufficient tension, ST2 water again becomes a “normal” fluid that expands when heated isobarically. The simulations showed evidence of a novel metastable firstorder phase transition at low temperature, which terminates at a critical point. The binodal for this transition is negatively sloped in the P-T plane, which implies that the high-density phase has higher entropy. The cause of the increase in heat capacity and compressibility upon supercooling would then be the proximity to the metastable critical point. Because of the low temperatures at which the metastable transition occurs, it manifests itself as a transformation between two glasses.18-21 Figure 1. Representative section of the bcc lattice. Filled circles denote In the past decade, a considerable body of evidence has indeed sites on the A sublattice; open circles are sites on the B sublattice. accumulated that strongly suggests that the transition between the two known amorphous solid phases of water (low-density amorph, LDA; high-density amorph, HDA) is f i r s t - ~ r d e r . ~ ~ - ~ ~ This evidence includes measurements of the volume and heat effects that accompany the LDA HDA transition. In this work, we study two approximate solutions of a microscopic model with orientation-dependent interactions. The model fluid forms low-density bonded networks and becomes denser when heated isobarically. Both Speedy-like retracing of the spinodal and PSES-like retracing of the locus of density maxima, with a second critical point, are predicted by suitable choices of model parameters (e.g., number of distinguishable orientations of a molecule and strength of the orientationdependent interaction). The model sheds light on the microscopic basis underlying the possible behavior of supercooled water; its generality suggests that similar behavior can occur in other network-forming fluids, such as silica and germania.
-
2. The Model and Its Ground States Our model fluid is not water-like in its detailed geometry: its ground state network, for example, is 6-fold-coordinated. The goal of this work is not to describe water specifically but rather to investigate in general the stability limits of networkforming liquids and their relation to density anomalies. In this sense, strong directional bonding and network formation are essential, whereas the detailed geometry of the network is less important. The model is thus descriptive and not predictive. Several lattice models have been proposed to study the properties of ~ a t e r . ~ ~Until - ~ *r e c e n t l ~ , ~little ~ - attention ~~ has been given in such studies to the metastable states. In our lattice model, as in liquid water and silica, density anomalies appear through a competition between configurations of low energy, low density, and low entropy and configurations of higher energy, higher density, and higher entropy. To model such a competition on a lattice, we require that bonding-a condition of low energy and entropy-should only occur between molecules more than a first neighbor apart. The body-centeredcubic (bcc) lattice allows for the proper interactions in three dimensions. Figure 1 shows a section of the bcc lattice. We divide the lattice into two interpenetrating cubic sublattices A and B; sites on the A sublattice are denoted by filled circles, while empty circles are sites on the B sublattice. The possible energy interactions are shown in Figure 2. Each molecule can adopt any one of q orientations; a variable a (a = 1, 2, ..., q) describes the orientational state of each molecule. Two molecules a fist-neighbor distance apart (Le., an AB pair on the bcc lattice, A t 2 lattice units apart) contribute energy - E , regardless of their mutual orientation. Two molecules a secondneighbor distance apart (i.e., an AA or a BB pair on the bcc lattice, 1 lattice unit apart; for example, molecules i and k in Figure 2) can form a bond with energy -J if both molecules are correctly mutually oriented and the four neighbor sites that
Figure 2. Possible energy interactions between nearest and next-nearest neighbors. Sites j are the centers of the four cells that share the edge along which bond ik forms.
surround the bond are empty. These sites are labeledj in Figure 2; they are the centers of the four cells that share the edge along which the bond forms. By convention, a = 1 for a correctly oriented molecule; thus, a bond can only form in 1 out of q2 mutual orientations between two second-neighbor molecules. The entropic penalty for bonding is clear. Note that it is only the mutual orientation of a pair, and not the absolute orientation of each molecule, that determines the possibility of forming a bond. The magnitude of the bond energy between two correctly oriented second-neighbor molecules is decreased by a factor cJ/4 (0 c I1) for each nearest-neighbor site that is occupied. We attribute this nearest-neighbor energy penalty to steric hindrance. Thus, bonding requires proper mutual orientation and proper separation and is weakened by the presence of other molecules sufficiently close to a bonded pair. Energetics favor the formation of a 6-fold-coordinated network on a half-filled lattice. T h e above interactions are described by the Hamiltonian (A
B)
i
j
i
k
In eq 1, ni is an occupation (Ising) variable (ni = 1 if site i is occupied; ni = 0 if site i is empty). For brevity, we write 6il instead of 6,,1, where 6 is the Kronecker delta. Thus, 6ii is an orientation (Potts) variable (= 1 if molecule i is properly oriented for bonding, 0 if it is not). The first double sum in the Hamiltonian is the contribution of all nearest-neighbor pairs, i and j being on different sublattices. The second term in the
Phase Behavior of Supercooled Water
J. Phys. Chem., Vol. 99,No. 11, 1995 3783
Hamiltonian is the bonding contribution. In it, the first summation is over all lattice sites. For each lattice site, the second summation spans three next-nearest neighbors on the same sublattice (for example, those along the +x, +y, and +z directions); the third summation spans the four j-type sites surrounding the given ik edge (Figure 2). Thus, each occupied i j nearest-neighbor pair (Figure 2) contributes an energy -E (if site k is empty or if pair ik is not properly oriented) or cJ14 E (if site k is occupied and both i and k are properly oriented). Because there are at most three bonds originating from site i, the total energy attributable to each occupied ij nearest-neighbor pair can have one of four values: - E , cJ14 - E, cJ12 - E, or 3cJ/4 - E, depending on whether none, 1, 2, or all of these 3 bonds are active. An important feature of the model is that an open ground state can be stabilized by bonding. In the open structure, one of the two sublattices is completely filled and bonded, while the other is empty; the overall bcc lattice is thus half-occupied. Denoting the number of lattice sites by N and the number of occupied sites by M , we have, for the open and fully bonded structure, 2M=N
(2)
U = -3/2NJ = -3MJ
(3)
H = M[2PV0 - 3 4
(4)
where U and H are the configurational energy and enthalpy, P is the pressure, and vo is the volume per site. For the closepacked and fully bonded state (which is always more stable than the close-packed and nonbonded ground state when c
