PHYSICAL REVIEW E 78, 021203 共2008兲 Simple three-state lattice model for liquid water Alina Ciach,1 Wojciech Góźdź,1 and Aurélien Perera2 1 2 Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland Laboratoire de Physique Théorique de la Matiére Condensée (UMR CNRS 7600), Université Pierre et Marie Curie, 4 Place Jussieu, F75252, Paris cedex 05, France 共Received 18 April 2008; published 19 August 2008兲 A simple three-state lattice model that incorporates two states for locally ordered and disordered forms of liquid water in addition to empty cells is introduced. The model is isomorphic to the Blume-Emery-Griffith model. The locally ordered 共O兲 and disordered 共D兲 forms of water are treated as two components, and we assume that the density of the D component is larger. The density of the sample is determined by the fraction of cells occupied by the O and D forms of water. Due to the larger density of the D state, the strength of the van der Waals 共vdW兲 interactions increases in the direction O-O ⬍ O-D ⬍ D-D. On the other hand, the H-bond interactions are assumed only for the O-O pairs. For the vdW and H-bond interaction parameters and the density ratio of the close-packed and ice forms of water compatible with experimentally known values, we find liquid-vapor and liquid-liquid transitions and the corresponding critical points in good agreement with other approaches. Water anomalies are correctly predicted within the mean-field approximation on a qualitative level. DOI: 10.1103/PhysRevE.78.021203 PACS number共s兲: 61.20.Gy I. INTRODUCTION Interactions among water molecules are characterized by the existence of the hydrogen bond, which is a strongly directional interaction. Combination of the directionality and the geometry of the water molecule gives rise to pronounced local tetrahedral ordering of the molecules in the liquid phase. The competition between the strong local order and the global disorder inherent to the liquid state is at the origin of the several anomalous properties of water, such as the maximum of density at 4C. The idea that such a competition could give rise to the coexistence of two types of water domains—a low-density ordered one and a high-density disordered one—dates back to Roentgen 关1兴 in 1892 and has known several rebirths since then, notably by Pauling 关2兴 and later by Frank and Wen 关3兴 and Robinson et al. 关4兴. However, spectroscopic experiments 关5兴 and computer simulations 关6兴 confirm only indirectly the fact that water would be formed of a coexistence of such two types of domains of water molecules, and most probably in supercooled states. Nevertheless, the two-state models for water have been very much appealing in providing a simple physical picture of its structure and in explaining its anomalies. There are other indirect indications that water can be in two different local states, such as, for example, in aqueous mixtures. There are several experimental pieces of evidence that water in the first solvation shell near ions is in a quite different state than in the bulk 关7兴, which leads to marked differences in the viscosity of these aqueous ionic solutions. Similarly, the aqueous solvation of small polar and nonpolar molecules is often discussed 关8兴 in terms of kosmotropy versus chaotropy, depending on their ability to increase or decrease the local order of water, which is also in relation to the so-called hydrophobic effect 关9兴. Water near amino acids has been found to have a different structural ordering by neutron-scattering experiments 关10兴. Temperature- and pressure-induced structural isosbestic points in neat water 关11兴 provide additional experimental support for a two-state description of liquid water. 1539-3755/2008/78共2兲/021203共13兲 The implicit existence of the coexistence of two types of liquid states in water implies possibility of the existence of a second critical point, presumably in the supercooled region 关12–14兴. The existence of such a critical point has been at the heart of much debates 关12–17兴 and has motivated the appearance of two-state models, both lattice and off-lattice, such as in Refs. 关13,16,18兴. It is worth mentioning that recent critical-point-free approaches are also able to describe the anomalies of cold water 关19,20兴. Lattice models for water have a long history, dating back to the early work of Bell in 1972 关21兴. A detailed list of such models is given in Ref. 关13兴. In these early models, water is often treated with full tetrahedral orientational dependence more or less explicitly described. The main emphasis of such models is to capture the fact that the H bonds would form a network. Recent studies 关22–24兴 of the core-softened offlattice models 关25兴 indicate that spherically symmetric interactions are equally able to capture the waterlike anomalies at low enough temperatures. These results suggest that the water anomalies should not need an explicit description of the details of the interactions even at microscopical level. The lattice- and continuum-space models introduced in Refs. 关13,14,16,18,26兴 are quite successful in predicting the thermodynamics and structure of water, but usually they are rather complex. For example, in Ref. 关16兴 orientational degrees of freedom of water molecules that occupy a bcc lattice are taken into account. Apart from the van der Waals interaction between the occupied nearest-neighbor 共NN兲 sites, the H-bond energy contribution is assumed for properly oriented NNs. The open network is promoted by an energy penalty for occupied sites around the H-bonded pair of sites. The advantage of this model is a possibility of predicting stability of fluid and solid phases. Unfortunately, on the simple meanfield 共MF兲 level the model predicts no liquid-liquid separation, which is found beyond the MF, however. Fundamental cells with 98 states 共and 98 energies兲 are considered in order to obtain the equation of state and phase diagram in Ref. 关16兴. The complexity of the model and the failure of the simple MF makes the extension of the model to aqueous 021203-1 ©2008 The American Physical Society PHYSICAL REVIEW E 78, 021203 共2008兲 CIACH, GÓŹDŹ, AND PERERA mixtures rather difficult. In Refs. 关27,28兴 extensions of the simpler, off-lattice mean-field model of water 关29兴 have been proposed to mixtures with nonpolar molecules 关27,28兴. However, we are interested in mixtures of water with polar molecules 共alcohols, amides, and other weak amphiphiles兲, because such mixtures show intriguing properties 关30兴; to some extent they resemble microemulsions formed in binary water-surfactant and ternary oil-water-surfactant mixtures. Universal properties of binary water-surfactant and ternary oil-water-surfactant mixtures are described by simple generic lattice models 关31–35兴 where water is treated as an ordinary liquid, because the effective energies involved in the selfassembly are higher than the energy involved in local ordering of water. Self-assembly of weak amphiphiles such as alcohols may compete with local ordering of water, and in the case of weak amphiphiles standard models of microemulsions are oversimplified. To understand the origin of their peculiarities it is desirable to introduce a simple model, yet capable of reproducing the special properties of pure water. To this end we need a simple model for water in the first place. Our purpose here is a construction of a generic lattice model which correctly predicts the qualitative properties of water and is as simple as possible. In general, the generic model represents the simplest system belonging to the class of systems exhibiting the same universal properties. The universal properties result from collective phenomena, where fluctuations over large length scale play the dominant role. The details of the local molecular structure are irrelevant, because these local details are averaged over large volume. The famous example of universal properties are the same values of the critical exponents for the gas-liquid critical point in simple fluids, Curie point in magnets, and demixing critical points in alloys. In the case of simple fluids such a generic model for the gas-liquid separation is the lattice-gas model, isomorphic to the Ising model for magnets. The physics of the gas-liquid critical point of water is also correctly predicted by the standard lattice-gas model. At high densities and low temperatures, however, the standard lattice-gas model is oversimplified. Here we introduce a very simple, three-state lattice model capable of predicting the special features of water. The model turns out to be of the same form as the Blume-Emery-Griffith 共BEG兲 model 关36兴. is necessary to distinguish between the locally ordered and locally disordered states, and postulate proper probability measures for each spatial distribution of these states, to be able to predict the density of the sample for given intensive parameters. The distinction between the presence and absence of the local tetrahedral order is necessary for reproducing the special properties of water. The question is if additional information about the structure is necessary for predicting the water anomalies. We shall assume that the necessary and sufficient information about the state of the sample is given when one specifies which cells are empty and selects from the remaining cells those with the tetrahedral order. Such partial information about the state of the sample is called the mesoscopic state. Our hypothesis is that for a description of collective phenomena leading to the density anomaly, no additional information about the local order is necessary, provided that the energetics of the system is properly taken into account in the probability measure ⬀ exp共−␤H兲. Here ␤ = 1 / kT, with k and T the Boltzmann constant and temperature, respectively, and H = E − ␮N is the mesoscopic Hamiltonian—a function that relates energy minus chemical potential times the number of particles to each mesoscopic state. H is derived below. The above hypothesis will be verified by comparing the predictions of the mesoscopic model with experiments. We should stress that in the mesoscopic description we lose the information about the number of microscopic states in the cells occupied by the locally ordered and locally disordered water. It is expected that the entropy associated with local ordering in the cell occupied by the locally ordered water is lower than the entropy associated with local ordering in the cell occupied by the locally disordered water. In the mesoscopic model this difference in entropy associated with local ordering is disregarded. The results obtained in the mesoscopic model should shed light on the role of the entropy of the local ordering for the water anomalies. The locally disordered and locally ordered states will be labeled 1 and 2, respectively, and the state representing the empty cell will be labeled 3. In order to introduce the Hamiltonian, we should first choose the proper size of the unit cell. For fixed volume the ratio between the number of molecules in the closely packed case 共no H bonds兲 and in the perfect tetrahedral structure is II. CONSTRUCTION OF THE GENERIC MODEL ␳HDW = 1 + 2␦ , ␳LDW In the lattice-gas models one divides space into cubic cells of a size comparable to the size of molecules. The microscopic state of the system is determined by specifying the state of each cell, instead of describing the states of all molecules. In the simplest lattice-gas model, isomorphic to the Ising model, just two states of each cell are distinguished: the cell can be either occupied or unoccupied by the center of mass of the molecule. In liquid water some molecules form H bonds with their neighbors. The H bonds are associated with formation of the local tetrahedral structure and with larger volume per molecule than in the disordered state. The density of the whole sample is thus determined by the fraction of molecules that form the local tetrahedral structure. It 共1兲 where ␳HDW and ␳LDW denote the densities in the closepacked and perfect tetrahedral structure, respectively, and ␦ is the model parameter defined by the above equation. Let us assume for a moment that only empty cells and state 2 are allowed. The natural choice for the volume v of the lattice cell is the volume per molecule in the perfect tetrahedral structure, v0. With this choice the whole density interval 共0 , ␳LDW兲 can be obtained from the model where the cells can be either unoccupied or occupied by just one molecule. The density ␳LDW corresponds to the occupancy of all cells. Let us now assume that only empty cells and state 1 are allowed. In the model where the cells can be either unoccu- 021203-2 PHYSICAL REVIEW E 78, 021203 共2008兲 SIMPLE THREE-STATE LATTICE MODEL FOR LIQUID WATER −a−h −a(1+2δ ) The mesoscopic Hamiltonian can be written in terms of the operators representing the state occupancy, which are ␳ˆ i共x兲 = 1 , 0 when the cell x is or is not in the state i, respectively. For each cell x the occupancy operators satisfy the 3 condition 兺i=1 ␳ˆ i共x兲 = 1. The above condition guarantees that each cell is in just one state, either 1, 2, or 3. The Hamiltonian is defined by 2 −a(1+2δ ) d H = − 兺 兺 兵a共1 + 2␦兲2␳ˆ 1共x兲␳ˆ 1共x + ei兲 + 共a + h兲␳ˆ 2共x兲␳ˆ 2共x x i=1 FIG. 1. Nonvanishing nearest-neighbor interactions 共upper panel兲 and two-dimensional cross sections through a part of a system in two mesoscopic states 共lower panel兲. a and h represent the van der Waals and H-bond energy parameters, respectively, and 2␦ is the relative difference between the close-packed and ice densities 共1兲. Dark-gray and light-gray squares represent cross sections through cells occupied by the closely packed and by the open tetrahedral structures, respectively, and the open squares represent empty cells. The mesoscopic state shown in the lower-left panel is typical in the 共supercritical兲 liquid, and the mesoscopic state shown in the lower-right panel is typical for vapor. Probability that a given mesoscopic state appears is proportional to the Boltzmann factor exp关−␤共E − ␮N兲兴, where the energy E is given by the sum of energies of all NN pairs 共according to the upper panel兲 and the number of molecules is N = N1共1 + 2␦兲 + N2, with N1 and N2 denoting the number of dark-gray and light-gray cells, respectively 关see 共9兲 and 共2兲兴. pied or occupied by just one molecule, the cell-volume should be v = v0␳LDW / ␳HDW ⬍ v0, so that the occupancy of all cells corresponds to the density ␳HDW. If both states and empty cells are possible, we should be able to relate the density of the liquid to the fraction of the cells in state 1 and the fraction of the cells in state 2. We choose the volume of the unit cell equal to the volume per molecule in the perfect tetrahedral structure, v0, and assume that there is one molecule per cell in state 2, whereas the cell in state 1 contains ␳HDW / ␳LDW molecules. The energy of the system depends on the number of H-bonded pairs and on the van der Waals 共vdW兲 interaction. We shall assume that only NNs interact. This assumption was verified for the vdW interaction in the context of the gas-liquid transition in the lattice-gas models in a large number of works. Clearly, an NN interaction is consistent with the short range of the H bonds. Since the H bonds are associated with the formation of tetrahedral structure, we assume that each pair of NNs that are both in state 2 contributes −h to the system energy. The pair of NNs yields no contribution to the system energy if at least one of the two cells is empty. The vdW energy of a pair of the nn cells is −a when both cells are in state 2, −a共1 + 2␦兲 if one cell is in state 1 and the other one is in state 2, and −a共1 + 2␦兲2 if both NN cells are in state 1. Here a is the vdW parameter. The above form of interaction energies directly follows from different densities in the two states, Eq. 共1兲, and the dependence of the vdW energy on the product of densities in the considered regions. Finally, the number of molecules in the system is equal to the number of cells in state 2 plus ␳HDW / ␳LDW times the number of cells in state 1. The above assumptions define the generic 共or minimal兲 model 共see Fig. 1兲. + ei兲 + a共1 + 2␦兲关␳ˆ 1共x兲␳ˆ 2共x + ei兲 + ␳ˆ 2共x兲␳ˆ 1共x + ei兲兴其 − ␮ 兺 关␳ˆ 2共x兲 + 共1 + 2␦兲␳ˆ 1共x兲兴. 共2兲 x The first sum runs over all lattice cells x, d is the space dimension, and ei denotes the unit lattice vector in direction i, where 1 艋 i 艋 d. With such a form we avoid double counting of the pairs of cells. It is convenient to introduce new variables—the local concentration and site occupancy—corresponding to two order-parameters 共OPs兲, ŝ = ␳ˆ 1 − ␳ˆ 2, ŝ = 1,− 1,0, 共3兲 ŝ2 = 1,0, 共4兲 ŝ2 = ␳ˆ 1 + ␳ˆ 2, where for the high-density water 共HDW兲, low-density water 共LDW兲, and vacuum, ŝ = 1, −1, and 0 and ŝ2 = 1, 1, and 0, respectively. The Hamiltonian in the new variables takes the familiar BEG-model form 关36兴 H = − 兺 兺 兵Jllŝ共x兲ŝ共x + ei兲 + 4Jglŝ2共x兲ŝ2共x + ei兲 x i + Q关ŝ共x兲ŝ2共x + ei兲 + ŝ2共x兲ŝ共x + ei兲兴其 − ␮␦ 兺 ŝ共x兲 − ␮共1 + ␦兲 兺 ŝ2共x兲, x 共5兲 x where h Jll = a␦2 + , 4 共6兲 h 4Jgl = a共1 + ␦兲2 + , 4 共7兲 h Q = a共1 + ␦兲␦ − . 4 共8兲 A mesoscopic state of the sample is given by the series 兵ŝ共x兲其, where ŝ共x兲 assumes a particular value 共1, −1, or 0兲 in each cell. The probability of a particular mesoscopic state is given by p关兵ŝ共x兲其兴 = ⌶−1e−␤H , 共9兲 with ⌶= 兺 e −␤H , 兵ŝ共x兲其 where the sum runs over all mesoscopic states 兵ŝ共x兲其. 021203-3 共10兲 PHYSICAL REVIEW E 78, 021203 共2008兲 CIACH, GÓŹDŹ, AND PERERA The grand-thermodynamic potential is given by the standard statistical-mechanics formula where p is pressure, V is the number of lattice cells, and v0V is volume. In the case of waterlike substance there is no direct way of measuring 具ŝ典 and 具ŝ2典 separately; the only measurable quantity is density. The density of water in this model is given by 具␳典 = ¯␳␳LDW, with 冏 冦 共11兲 ⍀ = − pv0V = − kT ln ⌶, ¯␳ = 具␳典/␳LDW = 关共1 + ␦兲具ŝ2典 + ␦具ŝ典兴 = − in the parameter space 共a , h , ␦ , ␮兲 corresponding to the stability of the three phases. From Eq. 共2兲 we obtain ⳵␤⍀ ⳵␤␮共x兲 冏 H = − 关3共a + h兲 + ␮兴V 冏 ⳵ 2␤ ⍀ + ␦ŝ共x2兲兴典 − ␳ = − ⳵␤␮共x1兲⳵␤␮共x2兲 冏 . A. Ground state Let us first focus on the low-temperature properties of the model and find the conditions corresponding to the stability of different states and to the phase equilibria. For T → 0 the probability 共9兲 becomes very narrow, and at T = 0 the probability that any state different from the ground state appears vanishes. When only the ground-state contributes to ⌶ in 共10兲 and ␮ is fixed, then ⍀ = min共H兲. Whether the H assumes a minimum for the state 1, 2, or 3 depends on the parameters. We shall calculate H for three uniform states: ŝ = 1, ŝ = −1, or ŝ = 0 in each cell. The first two states represent the close-packed and tetrahedral condensed phases, respectively, and ŝ = 0 represents the vapor. Note that at T = 0 the density of vapor tends to zero. Note also that in the latticegas models there is no transition between liquid and crystalline phases, because it is not possible to distinguish between long- and short-range order in the condensed state. We are not interested in the stable crystalline, but in the metastable liquid phases; therefore, for our purposes this deficiency of lattice models is in fact an advantage. By comparing the value of H for ŝ = 1, ŝ = −1, and ŝ = 0, we can find the regions 共14兲 共15兲 3h ␮ = − 6共1 + ␦兲 + . a 2a␦ 共16兲 Finally, the coexistence between the s = 0 and s = 1 is given by ␮ = − 3共1 + 2␦兲. a 共13兲 III. PROPERTIES OF THE MODEL 冧 At the coexistence between the s = 1 and s = −1 phases we have ␮共xi兲=␮ In the following sections we analyze three limiting cases: T → 0; next, the case of high density and finally the case of low density. In these limiting cases the model can be simplified. The simplified versions of the model can be solved or compared to previously found solutions. By analyzing the results of the simplified versions of the model, we can find for what ranges of the model parameters the waterlike behavior can be found. By requiring the consistency of the calculated quantities with the known results for water, we can find the best choice of parameters. if ŝ = 0. h ␮ =−3 1+ . a a 共12兲 ¯2 0 冉 冊 , G␳␳共x1 − x2兲 = 具关共1 + ␦兲ŝ2共x1兲 + ␦ŝ共x1兲兴关共1 + ␦兲ŝ2共x2兲 if ŝ = − 1, The coexistence between two phases is obtained by equating the values of H in these phases. In this way we find that the coexistence between the s = 0 and s = −1 phases occurs for ␮共x兲=␮ where 具¯典 denotes the average with the probability 共9兲; we replaced ␮ by the cell-dependent chemical potential ␮共x兲 in 共5兲 and used 共11兲. The 共rescaled兲 density-density correlation function G␳␳ is defined by − 共1 + 2␦兲关3共1 + 2␦兲a + ␮兴V if ŝ = 1, 共17兲 All three phases can coexist at a triple point if the hydrogenbond and vdW energies are such that h = 2␦a. The dependence of the T = 0 phase diagram on the parameter ratio h / a is shown in Fig. 2. For h ⬍ 2␦a a transition between vapor 共vacuum兲 and high-density 共s = 1兲 phases occurs at ␮ / a = −3共1 + 2␦兲. For h ⬍ 2␦a the system does not exhibit the transition between the two 共metastable兲 liquid phases. For h ⬎ 2␦a there is the transition between vacuum and the low-density phase for ␮ / a = −3共1 + h / a兲 and, next, the transition between the low- and high-density phases for ␮ / a = −6共1 + ␦兲 + 3h / 共2␦a兲. Such a coexistence between two metastable liquid phases is hypothesized for water. Thus, for parameters that satisfy the condition h ⬎ 2␦ , a 共18兲 the waterlike behavior can be expected in this model. The value of H at the coexistence between the two dense phases, ŝ = 1 and ŝ = −1, can be obtained by inserting 共16兲 into 共14兲. Because for T → 0 we have ⍀ → H, we obtain from 共11兲 the pressure at the coexistence between the two liquid phases for T → 0: pcoexv0 = 3共1 + 2␦兲 冉 冊 h −a , 2␦ T → 0. 共19兲 By equating the left-hand side 共LHS兲 of the above equation to the value expected for real water, we obtain the first relation between the model parameters. The anomalous and normal properties of water are observed for p ⬍ p0 and for p ⬎ p0, respectively, where p0 ⬇ 2000 bars. We assume that pcoex ⬇ p0 at the coexistence between the two metastable liquid phases for T → 0. 021203-4 PHYSICAL REVIEW E 78, 021203 共2008兲 SIMPLE THREE-STATE LATTICE MODEL FOR LIQUID WATER Tllc ⬇ 冉 冊 4.5Jll 4.5 h a␦2 + . = k k 4 共23兲 By assuming that the LHS of 共23兲 is equal to the temperature at the critical point of the metastable liquid-liquid separation in water, we obtain the second relation between the model parameters. The phase coexistence in the Ising model occurs for Hll1 = 0, and from 共21兲 we obtain the chemical potential at the liquid-liquid coexistence 共including the critical point兲, ␮c = −6Q / ␦. C. Model in the case of low ␳ FIG. 2. The ground-state dependence on the model parameters a, h, and ␦ and on the chemical potential ␮. Regions corresponding to stability of different phases 共close-packed, tetrahedral, and vacuum兲 at T = 0 are indicated, and the solid lines between stability regions of two phases represent the corresponding phase coexistence. Different values of the ratio between the H-bond and vdW energies, h / a, correspond to different substances. The dotted and dashed lines represent examples of the ground-state dependence on ␮ for two substances showing qualitatively different behavior. Along the dotted line the vapor and close-packed phases are stable for ␮ / a 艋 −3共1 + 2␦兲 and ␮ / a 艌 −3共1 + 2␦兲, respectively. Such “normal” behavior is found for h / a ⬍ 2␦. Along the dashed line the vapor and close-packed phases are stable for low 关␮ / a 艋 −3共1 + h / a兲兴 and high 关␮ / a 艌 −6共1 + ␦兲 + 3h / 共2␦a兲兴 values of ␮ / a, respectively, while for −3共1 + h / a兲 艋 ␮ / a 艋 −6共1 + ␦兲 + 3h / 共2␦a兲 the tetrahedral condensed phase is found. The stability of the tetrahedral condensed phase that characterizes waterlike systems is found for parameters that satisfy the inequality h / a ⬎ 2␦. In the limiting case of high densities we can reduce the model to a simplified version by eliminating the states 兵ŝ共x兲其 that occur with very low probability. In the high-density case the empty cells are very rare, and the mesoscopic states with ŝ = 0 can be neglected in calculating the average quantities. In that case ŝ = 1, −1, and by inserting ŝ2 = 1 into Eq. 共5兲 we end up with the Ising model Hll = − 兺 兺 Jllŝ共x兲ŝ共x + ei兲 − Hll1 兺 ŝ共x兲 − Hll0 兺 1, x i Hgl = − 兺 兺 4Jglŝ2共x兲ŝ2共x + ei兲 + − ␮ 兺 共1 + ␦兲ŝ2共x兲. x x i 共24兲 By changing the variables, ŝ2 = ␴ˆ + 1 , 2 共25兲 we again obtain the Ising model ˆ 共x兲 − Hgl Hgl = − 兺 兺 Jgl␴ˆ 共x兲␴ˆ 共x + ei兲 − Hgl 1兺␴ 0 兺 1, x x i x 共26兲 with B. Model in the case of high ␳ x Let us focus on the case of low densities. There is no distinction between structured and unstructured water, when the NNs of an occupied cell are typically empty, as in the gas or supercritical phase. We assume that we can average over ŝ in the Hamiltonian 共5兲. The average value of ŝ is zero. For ŝ = 0 and ŝ2 = 0 , 1 we have x gl Hgl 1 = 6J + 1+␦ ␮, 2 共27兲 gl Hgl 0 = 3J + 1+␦ ␮. 2 共28兲 The chemical potential at the critical point and along the gl phase coexistence is given by Hgl 1 = 0 and ␮c = −12J / 共1 + ␦兲. We again use the simulation result for the temperature at the critical point in the Ising model, and from 共7兲 we obtain 共20兲 Tgl c ⬇ where 冋 册 4.5Jgl 4.5 h a共1 + ␦兲2 + . = k 4k 4 共29兲 Hll1 = 6Q + ␦␮ 共21兲 By equating the LHS of the above to the temperature at the gas-liquid critical point in water, we obtain the third relation between the model parameters. Hll0 = 12Jgl + 共1 + ␦兲␮ . 共22兲 D. Model parameters for water and The critical temperature in the Ising model is known from simulations, and its value is kTc ⬇ 4.5Jll. From this result and from 共6兲 we have Note that in water the temperature of the gas-liquid critical point is much higher than the temperature of the metastable liquid-liquid critical point. We should thus require that ll in waterlike substance Tgl c Ⰷ Tc . The above condition and 共18兲 021203-5 PHYSICAL REVIEW E 78, 021203 共2008兲 CIACH, GÓŹDŹ, AND PERERA give together the range of h / a for which the model can exhibit similar properties as real water 2␦ ⬍ h 4共1 + 3␦兲共1 − ␦兲 . Ⰶ 3 a 共30兲 From 共23兲 and 共29兲 we obtain a = a0k, h = h0k, a0 = h0 = 4 ll 4Tgl c − Tc , 4.5共1 + 2␦兲 ␦兲2Tllc − 4␦2Tgl c 共1 + 4.5共1 + 2␦兲 共31兲 account when comparing the MF results with experiments. Since the MF approximation for the BEG model is well known, we describe only the main steps of derivation of the equation of state, the spinodal line, and the phase diagram. The density-density correlation function is calculated in the Gaussian approximation. In the MF approximation the grand potential 共11兲 is approximated by the minimum with respect to ␩ and s of the functional ⍀ MF关s共x兲, ␩共x兲兴 = H关s共x兲, ␩共x兲兴 − TS关s共x兲, ␩共x兲兴, 共35兲 . 共32兲 where the entropy on the lattice has the form S关s共x兲, ␩共x兲兴 = − k 兺 We assume Tgl c = 650 K, Tllc = 180 K, pcoex = 1900 bars. 共33兲 From the above and Eqs. 共31兲, 共32兲, and 共19兲 we obtain a relation between v0 and ␦. Choosing v0 = 35 Å3, we obtain ␦ = 0.12, a0 = 433.7 K 共a = 598.5⫻ 10−23 J = 3.6 kJ/ mol兲, h0 = 135.02 K 共h = 186⫻ 10−23 J = 1.2 kJ/ mol兲. These values are in semiquantitative agreement with the water parameters. The H-bond energy is small compared to the value 23 kJ/ mol, but in the mesoscopic model we do not consider the individual bond, but rather the interaction energy averaged over different orientations and distances between the molecules. Note that there are four NNs in the tetrahedral structure and six NNs on the simple cubic lattice. The H bond is strongly directional, and by averaging over all orientations we obtain 4 ⫻ 23 kJ/mol ⬇ 1.83 kJ/mol. 6 ⫻ 4␲ x + 2 ␩共x兲 − s共x兲 2 冉 ln 冉 ln ␩共x兲 + s共x兲 2 ␩共x兲 − s共x兲 2 册 冊 冊 共36兲 and H关s共x兲 , ␩共x兲兴 is given in Eq. 共5兲 with ŝ and ŝ2 replaced by s and ␩, respectively. 0 艋 ␩ 艋 1 is the fraction of occupied cells and is analogous to the number density in a mixture isomorphic to our model. Here −1 艋 s = 2x − 1 艋 1, and x is analogous to the concentration of one species 共locally disordered water兲 in the corresponding mixture. In order to simplify the notation we rewrite Eq. 共12兲 in the form 共34兲 具␳典 = 共1 + ␦兲␩ + ␦s. ␳LDW 共37兲 For uniform phases Eq. 共35兲 assumes the explicit form ⍀ MF共s, ␩兲/V = − pv0 = − 兵3Jlls2 + 12Jgl␩2 + 6Qs␩ + ␮关␦s + 共1 + ␦兲␩兴其 − TsV , where we used Eq. 共5兲 and where sV = − k IV. RESULTS IN THE MEAN-FIELD APPROXIMATION The BEG model was first solved in Ref. 关36兴 in the MF approximation. Details of our version of the MF approximation can be found for example in Refs. 关33,36兴. In the MF approximation the equation of state, the spinodal line, and the transitions between uniform phases are correctly described on a qualitative level. However, when fluctuations are included beyond the MF approximation, the stability region of the disordered phase enlarges. Thus, the critical-point temperature obtained in the MF approximation, TcMF, is overestimated. In the case of the Ising model, Tc / TcMF ⬇ 3 / 4. Also, fluctuations yield a positive contribution to the grand potential ⍀ = −pv0; therefore, the MF result for pressure, p MF, is overestimated as well. Another important source of inaccuracy of the quantitative results of the model is the underlying lattice. In particular, the symmetrical shape of the gas-liquid coexistence line, resulting from the underlying lattice, differs significantly from the shape of the corresponding line in real systems. The above features should be taken into ␩共x兲 + s共x兲 + 关1 − ␩共x兲兴ln关1 − ␩共x兲兴 , ¯␳ = Note that we could fix the model parameters based on experimental values. However, the presence of the lattice structure changes the value of the critical temperature for given energy parameters; therefore, we proceeded in the way described above. 冋 冋 ␩+s 2 ln 冉 冊 ␩+s 2 + 共1 − ␩兲ln共1 − ␩兲 + 册 ␩−s 2 ln 共38兲 冉 冊 ␩−s 2 共39兲 is the entropy per lattice cell. The explicit forms of the minimum conditions, ⳵⍀MF = 0, ⳵s 共40兲 ⳵⍀MF = 0, ⳵␩ 共41兲 are given in Appendix A. A. Spinodal line and phase diagram The function ⍀ MF共s , ␩兲 assumes a minimum for s and ␩ satisfying Eqs. 共40兲 and 共41兲 when its second derivative is positive definite. The boundary of stability is given by 021203-6 PHYSICAL REVIEW E 78, 021203 共2008兲 SIMPLE THREE-STATE LATTICE MODEL FOR LIQUID WATER 1000 1400 1200 800 1000 800 600 T(K) T 600 400 400 1 200 0.5 200 0 0 0 s 0.2 0.4 -0.5 0 0.6 0.8 η -1 1 FIG. 3. The spinodal surface T = Ts共␩ , s兲 as a function of s and ␩. The s and ␩ are not independent in equilibrium states, and in the det关⳵2␤⍀ MF兴 = 0, 册 = 2 2 冊 册 ␩ −s kT ln − 24Jgl␩ − 6Qs ␦ . 2 4共1 − ␩兲2 共43兲 The spinodal line in the 共␳, T兲 phase diagram is given as a parametric equation in Eqs. 共A3兲, 共A6兲, and 共37兲 and is shown in Fig. 4. More calculation details are given in Appendix A. From the simplified versions of the model 共Secs. III B and III C兲 we can obtain the approximate forms of the spinodal line. For low density we assume s = 0—i.e., the density is ¯␳ = 共1 + ␦兲␩—and we obtain the approximate form of the spinodal 共see Sec. III C兲: kT = 24Jgl␩共1 − ␩兲 = 24Jgl ¯␳共1 + ␦ − ¯␳兲 . 共1 + ␦兲2 6J 1.2 p(bar) 0 兵¯␳关2共1 + ␦兲 − ¯␳兴 − 1 − 2␦其. 共45兲 Equations 共44兲 and 共45兲 represent the intersection of the surface T = Ts共␩ , s兲 with the planes s = 0 and ␩ = 1, respectively. 200 (b) 400 T(K) 600 800 FIG. 4. Spinodal 共dashed兲 and binodal 共solid兲 lines in the MF approximation. The symbol shows the experimental critical point of water. Note that exact values of the critical temperatures in this model are expected to be Tc ⬇ 0.75TcMF. ␳ / ␳LDW is dimensionless; ␳LDW is the model parameter, and its order of magnitude is ␳LDW ⬃ ␳ice ⬇ 0.9 g / cm3. For ␩ → 1 the approximation 共44兲 fails, because from 共43兲 it follows that when ␩ → 1, then s → −1 at the gas-liquid branch of the spinodal 共see Fig. 3兲. Water anomalies result from this failure of the simple lattice-gas model. On the other hand, Eq. 共45兲 is a good approximation; it almost coincides with the results of the full model. The phase equilibria shown in Fig. 4 are obtained by equating the values of p 关Eq. 共38兲兴 and ␮ 关Eqs. 共40兲 and 共41兲兴 for two phases ␣ and ␤ characterized by the pair of variables 共␩␣, s␣兲 and 共␩␤, s␤兲. For more details, see Appendix B. In the simplified versions of the model 共Secs. III B and III C兲 we obtain the simple expressions kT = 12Jgl关2¯␳ − 共1 + ␦兲兴 共1 + ␦兲关ln共¯␳兲 − ln共1 + ␦ − ¯␳兲兴 共46兲 kT = 12Jll关¯␳ − 共1 + ␦兲兴 ␦关ln共¯␳ − 1兲 − ln共1 + 2␦ − ¯␳兲兴 共47兲 and ll ␦2 1 1000 共44兲 The above form is valid only for ␩ Ⰶ 1. For high density we assume ␩ = 1—hence ¯␳ = 共1 + ␦兲 + ␦s—and consider the simplified model described in Sec. III B. The spinodal for the model 共20兲 assumes the form kT = 6Jll共1 − s2兲 = 0.6 0.8 ρ/ρ LDW 共42兲 is the matrix of the second derivative of where ⳵ ␤⍀ ␤⍀MF with respect to ␩ and s. The rather lengthy explicit expression for T = Ts共␩ , s兲, following from Eq. 共42兲, is given in Appendix A 关Eq. 共A3兲兴, and the surface T = Ts共␩ , s兲 is shown in Fig. 3. In equilibrium ␩ and s are not independent, and the relation between them 关Eqs. 共40兲 and 共41兲 and Eqs. 共A1兲 and 共A2兲 in Appendix A兴 is given by ␩+s kT − 6Jlls − 6Q␩ 共1 + ␦兲 ln 2 ␩−s 0.4 2000 MF 冋 冉 冊 冋 冉 0.2 3000 MF approximation the relation between them is given in Eq. 共43兲. The cross section of the surface T = Ts共␩ , s兲 with the surface 共43兲 is shown as the thick line. This line together with Eq. 共37兲 gives the spinodal line that represents the loci of the actual instability of the model system in the 共␳, T兲 diagram. T is in kelvins; s and ␩ are dimensionless. 2 0 (a) for the gas-liquid and the liquid-liquid transition, respectively. As in the case of the spinodal, Eq. 共46兲 is reasonably good at the gas branch and fails at the liquid branch of the gas-liquid coexistence. The approximation 共47兲 for the 021203-7 PHYSICAL REVIEW E 78, 021203 共2008兲 CIACH, GÓŹDŹ, AND PERERA B. Equation of state, thermal expansion, and heat capacity The equation of state 共EOS兲 in 共␮, p, T兲 variables is given in Eq. 共38兲, with ␩ and s eliminated with the help of Eqs. 共40兲 and 共41兲. In order to obtain the EOS in the ␳, p, and T variables, we eliminate ␮ from 共38兲 by using 共40兲, and by subsequent use of Eqs. 共40兲 and 共41兲 we obtain the simple form, valid for ␩ ⬍ 1, pv0 = − 3Jlls2 − 12Jgl␩2 − 6Q␩s − kT ln共1 − ␩兲. 冉 冊 ¯␳ 12Jgl 2 ¯ , 2 ␳ = − kT ln 1 − 共1 + ␦兲 1+␦ 1000 800 600 400 200 0 0 0.2 0.4 0.6 ρ/ρLDW 0.8 1 1.2 FIG. 5. 共Color online兲 T共␳兲 isobars 共solid lines兲. From the left to the right lines p = 10, 100, 653, 1000, 1950, 2085, and 2500 bar. The dashed lines are the phase transitions. The critical pressures for the gas-liquid and liquid-liquid transitions are pgl c ⬇ 653 bars and pllc ⬇ 2085 bars, respectively. In the MF approximation the extent of the two-phase regions and pressure are both overestimated compared to exact results. In addition, the pressure is overestimated due to the lattice structure. The exact values of the critical temperatures are expected to be Tc ⬇ 0.75TcMF. T is in kelvins, and ␳ / ␳LDW is dimensionless; ␳LDW is the model parameter, and its order of magnitude is ␳LDW ⬃ ␳ice ⬇ 0.9 g / cm3. large densities pressure is significantly overestimated. This well-known deficiency of the lattice-gas models results from the hole-particle symmetry on the lattice, which in real systems is absent. For high density we assume ␩ = 1 and consequently ¯␳ = 共1 + ␦兲 + ␦s 共see Sec. III B兲. The EOS for high densities takes the form 关see 共38兲 and 共40兲兴 kTL共¯␳兲 = pv0 + R共¯␳兲, where L共¯␳兲 = 冉 共50兲 1 + 2␦ 1 ln共¯␳ − 1兲 − ln共1 + 2␦ − ¯␳兲 + ln共2␦兲 2␦ 2␦ 共49兲 which is the lattice-gas version of the vdW EOS. Note that the RHS of Eq. 共49兲 diverges for ¯␳ → 1 + ␦; therefore, for 冊 共51兲 共48兲 The EOS in the 共␳, p, T兲 variables is given in Eq. 共48兲, with s and ␩ expressed in terms of ¯␳ and T by solving Eqs. 共43兲 and 共37兲. Some calculation details are given in Appendix C. The plots of the isobars T共␳兲 for p = 10, 100, 653, 1000, 1950, 2085, and 2500 bars are shown in Fig. 5. Note that the shape of the isobars for pcoex ⬍ p ⬍ pllc results from the positive slope of the liquid-liquid coexistence line p共T兲 in our simple MF approximation. In the asymptotic cases of very low and very high density, the EOS assumes simpler forms that can be obtained from the simplified versions of the model described in Secs. III B and III C. Let us first describe the case of low density. By assuming s = 0, we obtain ¯␳ = 共1 + ␦兲␩ and pv0 + 1200 T(K) liquid-liquid branch cannot be distinguished from the exact MF result in Fig. 4. Note that the slope of the liquid-liquid coexistence line p共T兲 is very small, but positive, contrary to the common expectation that the entropy per mole in the LDW is lower than the entropy in the HDW. We verified that in the MF approximation pcoex ⬍ pllc for 0.05⬍ ␦ ⬍ 0.2, where pcoex is pressure at the coexistence at very low T 关see 共19兲兴 and pllc is pressure at the critical point of the liquid-liquid coexistence. Thus, dp共T兲 / dT ⬎ 0 is the property of the MF solution of the model for the relevant range of parameters. At the present stage we are not able to verify if beyond the MF approximation the slope of the coexistence line remains positive. As already mentioned, fluctuations usually lead to lower values of pressure and the fluctuation contribution increases with increasing temperature. We can expect that exact solution for the coexistence line in this model should correspond to a smaller slope of the line p共T兲, and since this slope is very small in the MF approximation, no definite conclusions concerning the sign of dp共T兲 / dT can be drawn before reliable simulations are performed. In any case, even if it turns out that beyond the MF approximation dp共T兲 / dT ⬍ 0 in this model, its magnitude is expected to be much smaller than in other models where orientational degrees of freedom are explicitly included 关13,16,18,29兴. We stress that our model shows that the difference in the entropy associated with the local ordering in the LDW and HDW is not necessary for obtaining the density anomaly. This is reminiscent to the presence of the density anomaly in the ramp model, where interactions are independent of orientations 关22,24兴. and R共¯␳兲 = 3Jll ␦ 2 ¯␳2 + 6Q 1+␦ ␦ 冉 冊 − 12Jgl − 3Jll 1+␦ ␦ 2 . 共52兲 The critical pressure is given by pcv0 = kTcL共¯␳c兲 − R共¯␳c兲 ⬇ 2085 bars, where Tc ⬇ 240 K and ¯␳c = 1 + ␦ are obtained from Eq. 共45兲. Note that L共¯␳0兲 = 0 for ¯␳0 ⬇ 1.111 and the pressure p0 = −R共¯␳0兲 / v0 ⬇ 2084 bars is independent of T. For p ⬎ p0, T共¯␳兲 given in Eq. 共50兲 decreases from T共¯␳0兲 = ⬁ for ¯␳ increasing from ¯␳0, whereas for p ⬍ p0, T共¯␳兲 decreases from T共¯␳0兲 = ⬁ for ¯␳ decreasing from ¯␳0. Hence, p0 is the borderline value separating the normal 关negative slope of T共␳兲兴 from the anomalous 关positive slope of T共␳兲兴 density behavior for p ⬎ p0 and p ⬍ p0, respectively. The density ¯␳0 ⬇ 1.111 is the maximum density for p ⬍ p0. Note that in the simplified model 共Sec. III B兲 no holes are present, so its validity is limited to relatively low T. At higher T the negative slope of T共␳兲 is found in the full model, as seen in Fig. 5. 021203-8 PHYSICAL REVIEW E 78, 021203 共2008兲 SIMPLE THREE-STATE LATTICE MODEL FOR LIQUID WATER 0.0004 evant range of temperatures, in agreement with the anomalous behavior of water. Recall that in calculating the entropy we disregarded the difference in the entropy of local ordering in the locally ordered and locally disordered water 共cells with ŝ = −1 and ŝ = 1, respectively兲. Evidently, the anomalous heat capacity results from the entropy associated with different spatial distributions of the locally ordered and locally disordered regions, rather than from the entropy associated with local ordering in regions on the molecular length scale. 0.0002 α 0 -0.0002 -0.0004 200 250 300 350 400 450 C. Density-density correlation function and the compressibility 500 T FIG. 6. Thermal expansivity 共53兲 for p = 100 bars 共dashed line兲 and p = 1000 bars 共solid line兲; T is in kelvins. ␣ is in inverse kelvins. Beyond the MF approximation we expect a shift of the plot towards lower T values. Since for ␩ = 1 there are no empty cells, the EOS 共50兲 shows no singularity for ¯␳ = 1 + ␦. On the other hand, 兩L共¯␳兲兩 → ⬁ for ¯␳ → 1 and ¯␳ → 1 + 2␦. Therefore we expect significantly overestimated pressures for ¯␳ → 1 and ¯␳ → 1 + 2␦. The thermal expansivity 冏 冏 ␣= − ¯ 1 ⳵␳ ¯␳ ⳵T 冏 ⳵共sVNA/¯␳兲 ⳵T with 关C兴 dewhere according to standard rules noting the matrix of the second derivative of the function ␤⍀MF共s , ␩兲. The explicit form of 关C兴 in Fourier representation is C̃␩␩共k兲 = ␩ ␩ −s C̃ss共k兲 = G̃ss共k兲 = p 5 4 300 350 400 450 500 T FIG. 7. Isobaric heat capacity 共54兲 for p = 100 bars 共dashed line兲 and p = 1000 bars 共solid line兲. T is in kelvins 共K兲; c p is in J/K/mol. Beyond the MF approximation we expect a shift of the plot towards lower T values. 1 − 4JglZ, 1−␩ 共56兲 s − QZ, ␩ − s2 共57兲 2 ␩ 2 ␩ − s2 共58兲 − JllZ, 1 + AD1Z + AD2Z2 AN0 + AN1Z , 1 + AD1Z + AD2Z2 , 共59兲 共60兲 where the coefficients ANi and ADi are given in Appendix C. Their values should be calculated for s and ␩, which satisfy the equilibrium condition 共43兲. The results of the lattice models for fluids are meaningful for distances much larger than the lattice constant, where the lattice structure is irrelevant. For k → 0 we have Z = ␤关6 − k2 3 + O共k2i k2j 兲兴, where k2 = 兺i=1 k2i , and up to O共k2兲 Eqs. 共59兲 and 共60兲 show no anisotropy resulting from the lattice structure: G̃ f f 共k兲 ⯝ 250 + ␩ − s2 − 4Jgl共␩2 − s2兲共1 − ␩兲Z G̃␳␳共k兲 = 3 2 2 3 where we introduced Z = 2␤兺i=1 cos共ki兲 and k = 共k1 , k2 , k3兲. Standard algebra yields the expressions 共54兲 , 6 200 2 C̃␩s共k兲 = − 冏 共55兲 −1 , G␣␤ = 关C兴␣␤ 共53兲 where sVNA / ¯␳ is the molar entropy, with the entropy per cell sV given in Eq. 共39兲 共ideal entropy of mixing兲, and NA is the Avogadro number. The isobaric heat capacity is shown in Fig. 7 for p = 100 and 1000 bars. Note the increase of the heat capacity upon cooling at constant pressure for the rel- Cp G␳␳ = 共1 + ␦兲2G␩␩ + ␦2Gss + 2␦共1 + ␦兲G␩s , p can be obtained by direct differentiation of the isobars shown in Fig. 5 and is shown in Fig. 6 for p = 100 and 1000 bars. As expected from the shape of the isobars, we obtain the characteristic anomalous negative expansivity for the relevant range of temperature. The isobaric heat capacity is given by cP = T The density-density correlation function 共13兲 can be written in the form Af f , + k 2兲 共␰−2 ff k → 0, 共61兲 where f denotes s or ␳. The correlation function Gss共k兲 is analogous to the correlation function for concentration fluctuations in binary mixtures. The region of excess concentration of one component in this model is analogous to the cluster of H-bonded molecules. The average spatial extent of the “cluster” is of the order of the correlation length for the concentration of the locally ordered 共or disordered兲 component. 021203-9 PHYSICAL REVIEW E 78, 021203 共2008兲 CIACH, GÓŹDŹ, AND PERERA 8 -1 χ Τ(GPa ) ξ(A) 0.25 6 4 2 0.2 0.15 240 260 280 300 T(K) 320 340 360 250 For ␩ = 1 and s = 0 the correlation functions assume the simple form kT/Jll , 共␰−2 + k2兲 k → 0, 共62兲 and in real space rGss共r兲 ⯝ kT −r/␰ e . 4␲Jll 共63兲 The correlation length ␰ = ␰ss = ␰␳␳ = 冑kT/Jll − 6 共64兲 is shown in Fig. 8 as a function of T. The above approximate forms are valid for the parameter range where the empty cells are very rare—i.e., for high densities and low temperatures. The isothermal compressibility in ␳ and T variables is given by ␹T = 350 400 450 500 T(K) FIG. 8. The correlation length for ␩ = 1 and s = 0. ␰ is in angstroms; T is in kelvins. In our MF approximation, TcMF ⬇ 240 K 共dashed line兲. G̃␳␳共k兲 = ␦2G̃ss共k兲 ⯝ ␦2 300 1 ⳵␳ G̃␳␳共0兲v0 AN0kT + 6AN1 v0 = = , ␳ ⳵p 共kT兲2 + 6AD1kT + 36AD2 ¯␳2 kT¯␳2 共65兲 where ␹T as a function of T and p can be obtained by eliminating ␩ and s from 共65兲 with the help of 共43兲 and 共48兲. The result is shown in Fig. 9 for p = 100 and 1000 bars. The anomalous behavior of the isothermal compressibility is found in this model only for high pressures 共but lower than p0 ⬇ 2000 bars兲. As discussed in the context of the EOS, pressure is significantly overestimated in the lattice models. V. SUMMARY We have introduced a very simple, three-state lattice model for waterlike fluid. The model is a special version of the BEG model 关36兴. Water is treated as a mixture of two components, whose chemical potentials are not independent. One component represents locally ordered and the other one locally disordered water. It is assumed that the density of the disordered component is larger, and the ratio between the FIG. 9. Compressibility as a function of temperature for p = 1000 bars 共solid line兲 and p = 100 bars 共dashed line兲. T is in kelvins; ␹T is in GPa−1. chemical potentials is given by the density ratio 1 + 2␦. Interaction potentials of the van der Waals type are proportional to the common parameter a and to the product of densities of the interacting species. In addition, the H-bond interaction h is assumed for the component representing the locally ordered water. The liquid-liquid critical point in this model is identified with the critical point of the demixing transition in the corresponding mixture. We fix the model parameters by requiring that the critical points for the gas-liquid and demixing transition coincide with the experimental values and the demixed components coexist at very low temperature for p = 1900 bars. There are thus three relations between four model parameters ␦, a, h, and the volume per particle in the ice structure, v0. We arbitrarily choose the volume per particle in the ice structure, v0 = 35 Å3. The chosen parameters are in semiquantitative agreement with experimental values. We solve the model in mean-field approximation. The shape of the gas-liquid coexistence line 共Fig. 4兲 is in qualitative agreement with earlier theoretical and experimental results 关16,18兴. The maximum-density temperature, however, is too high compared to real water. The shape of the phasecoexistence lines depends on the model parameters, especially on ␦, but we did not attempt to find the best choice of parameters. Improvement of the quantitative results is expected beyond the simple MF approximation. Note that in Ref. 关16兴 water anomalies are not found on the same MF level as in our case. The anomalous density increase with temperature for p ⬍ p0 ⬇ 2084 bars in this model is explained as follows. For low T and p ⬍ p0 ⬇ 2084 bars the volume fraction of the low-density component is high. When T increases, the increasing role of entropy leads first to mixing of the locally ordered and disordered components that occurs at a relatively low energy cost. Hence, density increases. Further heating leads to a still more important role of entropy. As a consequence, the two components start to mix with empty cells —this occurs at a larger energy cost than mixing of the locally ordered and locally disordered components. The presence of the empty cells leads in turn to a decrease of density, as shown schematically in Fig. 10. The EOS, thermal expansion, isobaric heat capacity, correlation length, and compressibility also show qualitative 021203-10 PHYSICAL REVIEW E 78, 021203 共2008兲 SIMPLE THREE-STATE LATTICE MODEL FOR LIQUID WATER ACKNOWLEDGMENTS This work was partially supported by the Polish Ministry of Science and Higher Education, Grant No. NN 202 006034 and partially by the POLONIUM project. p<p0 APPENDIX A: SPINODAL LINE The explicit forms of 共40兲 and 共41兲 are the following: p>p0 ␮= T1 < T2 < T3 FIG. 10. Schematic showing a sequence of typical configurations for p ⬍ p0 ⬇ 2000 bars 共top兲 and p ⬎ p0 ⬇ 2000 bars 共bottom兲 and three different temperatures. Density ratio between the dark gray and light-gray regions is 1 + 2␦ ⬎ 1. For more explanations concerning the symbolic representation of the mesoscopic states, see Fig. 1. features characteristic of water, but pressure is significantly overestimated in the simple MF approximation. Overestimated pressure is expected for the lattice models and in the MF approximation, and this deficiency will be removed 共to some extent at least兲 in the off-lattice version of the model that we plan to develop. Beyond the MF approximation pressure should be lower; the effect of fluctuations depends on T, and the slope of the coexistence lines can be significantly different in an exact solution of our model. Orientational ordering of water molecules is not taken into account directly, but only through the existence of the locally ordered water, which results from the orientational ordering. Orientational degrees of freedom certainly influence quantitative results for the grand potential 共pressure兲 and entropy; for quantitative agreement of the EOS with real water, the generic model is certainly oversimplified. On the qualitative level the only result that disagrees with common expectations is the slope of the 共metastable兲 liquid-liquid coexistence line p共T兲, which is positive 共but very small兲 in the simple MF approximation. The slope of this line beyond the MF approximation is not known; it might be negative, but large negative values are not expected. This is a direct consequence of not taking into account the difference in entropy associated with local ordering in the cells occupied by the LDW and HDW states respectively. We have shown that the model as simple as the BEG model with properly adjusted parameters is capable of predicting the distinctive anomalies of water. This result suggests that the entropy associated with different spatial distributions of the locally ordered and disordered regions is responsible for all the water anomalies, since they are found even when the entropy associated with local ordering in the locally ordered and locally disordered structures is disregarded. Finally, the relative success of this simple model in the simplest MF approximation indicates that similar approach can be used to study other associating liquids, such as alcohols, for example, as well as aqueous mixtures with polar molecules. ␮= 冋 冉 冊 册 册 ␩+s kT − 6Jlls − 6Q␩ ␦−1 , ln 2 ␩−s 冋 冉 冊 共A1兲 ␩2 − s2 kT ln − 24Jgl␩ − 6Qs 共1 + ␦兲−1 . 2 4共1 − ␩兲2 共A2兲 From Eq. 共42兲 we obtain the explicit expression for the spinodal surface in the 共␩ , s , T兲 variables: kTs共␩,s兲 = − b共␩,s兲 + 冑⌬s共␩,s兲, 共A3兲 where ⌬s共␩,s兲 = b共␩,s兲2 − 36共4JllJgl − Q2兲共␩2 − s2兲共1 − ␩兲, 共A4兲 b共␩,s兲 = − 3关Jll␩ + 4Jgl␩共1 − ␩兲 + 2Qs共1 − ␩兲 − Jlls2兴. 共A5兲 From the relation 共43兲 we can obtain temperature as a function of ␩ and s. By equating this temperature with Ts given in 共A3兲, we find that along the spinodal line the relation between ␩ and s is given by − b共␩,s兲 + 冑⌬s共␩,s兲 = 12兵关共1 + ␦兲Jll − ␦Q兴s + 关共1 + ␦兲Q − 4Jgl␦兴␩其 冋 冉 ␩ + s 2共1 − ␩兲 ln ␩−s ␩−s 冊册 2␦ . 共A6兲 APPENDIX B: PHASE EQUILIBRIA Equations 共40兲 and 共41兲 can be written in the form s = Fs共s, ␩兩, ␮,T兲, ␩ = F␩共s, ␩兩, ␮,T兲, 共B1兲 where F␩共s, ␩兩, ␮,T兲 = 2 cosh共␤Rs兲e␤R␩ , ⌶0 共B2兲 Fs共s, ␩兩, ␮,T兲 = 2 sinh共␤Rs兲e␤R␩ , ⌶0 共B3兲 ⌶0 = 2 cosh共␤Rs兲e␤R␩ + 1, 共B4兲 Rs = 6共Jlls + Q␩兲 + ␦␮ , 共B5兲 R␩ = 6共4Jgl␩ + Qs兲 + 共1 + ␦兲␮ . 共B6兲 We look for the fixed points of 共B1兲 by using iterations (si+1 = Fs共si , ␩i兲 , ␩i+1 = F␩共si , ␩i兲). The solutions give s and ␩ 021203-11 PHYSICAL REVIEW E 78, 021203 共2008兲 CIACH, GÓŹDŹ, AND PERERA in the stable or metastable phase for given T and ␮. The stable phase corresponds to a lower value of −pv0. At the phase equilibria, pv0共s␣ , ␩␣兲 = pv0共s␤ , ␩␤兲, where ␣ and ␤ denote the two coexisting phases. This method allows also for obtaining the EOS isotherms p共¯␳兲. where s1,2共␩,kT,p兲 = − 6Q␩ ⫾ 冑⌬共␩,kT,p兲 6Jll ␩ + si 2共1 − ␩兲 ␩ − si ␩ − si 冊册 2␦ − 12兵关共1 + ␦兲Jll − ␦Q兴si + 关共1 + ␦兲Q − 4Jgl␦兴␩其. 共C7兲 The latter relation follows from 共43兲. Depending on parameters, either 共C4兲 or 共C6兲 has solutions. Accordingly, the density is given either in Eq. 共C3兲 or 共C5兲, respectively. 共C1兲 APPENDIX D: CORRELATION FUNCTION and ⌬共␩,kT,p兲 = 36Q2␩2 − 12Jll关p + 12Jgl␩2 + kT ln共1 − ␩兲兴. The explicit forms of the coefficients in Eq. 共60兲 are the following: 共C2兲 AN0 = 关共1 + ␦兲2 + ␦2兴关共1 − ␩兲␩ + 2␦共1 + ␦兲兴共1 − ␩兲s The EOS is given by the parametric equations ¯␳共␩,kT,p兲 = 共1 + ␦兲␩ + ␦s1共␩,kT,p兲 冋 冉 Fi共␩,kT,p兲 = kT ln APPENDIX C: EOS Here we describe the method for obtaining the EOS isobars. From Eq. 共48兲 we obtain s = s1共␩ , kT , p兲 or s = s2共␩ , kT , p兲, where 共C6兲 F2共␩,kT,p兲 = 0, + ␦2共␩2 − s2兲, 共D1兲 共C3兲 AN1 = 关2␦共1 + ␦兲Q − 4␦2Jgl − 共1 + ␦兲2Jll兴共1 − ␩兲共␩2 − s2兲 and F1共␩,kT,p兲 = 0 = − h共1 + 2␦兲2 , 共C4兲 if Eq. 共C4兲 has a solution. When Eq. 共C4兲 has no solutions, the EOS is given by the parametric equations ¯␳共␩,kT,p兲 = 共1 + ␦兲␩ + ␦s2共␩,kT,p兲 AD1 = − 关4Jgl共1 − ␩兲␩ + Jll共␩ − s2兲 + 2Qs共1 − ␩兲兴, 共D3兲 共C5兲 AD2 = 共4JglJll − Q2兲共1 − ␩兲共␩2 − s2兲. and 关1兴 W. K. Roentgen, Ann. Phys. 45, 91 共1892兲. 关2兴 L. Pauling, in Hydrogen Bonding, edited by D. Hadzi and H. W. Thompson 共Pergamon Press, London, 1959兲. 关3兴 H. S. Frank and W.-Y. Wen, Discuss. Faraday Soc. 24, 133 共1957兲. 关4兴 G. W. Robinson, C. H. Cho, and J. Urquidi, J. Chem. Phys. 111, 698 共1999兲. 关5兴 S. Woutersen, U. Emmerichs, and H. J. Bakker, Science 278, 658 共1997兲. 关6兴 M. Sasai, Physica A 285, 315 共2000兲. 关7兴 H. J. Bakker, Chem. Rev. 共Washington, D.C.兲 108, 1456 共2008兲. 关8兴 F. Franks, Water: A Comprehensive Treatise 共Plenum Press, London, 1973兲, Vol. 3. 关9兴 B. Widom, P. Bhimalapuram, and K. Koga, Phys. Chem. Chem. Phys. 5, 3085 共2003兲. 关10兴 A. Pertsemlidis, A. M. Saxena, A. K. Soper, T. Head-Gordon, and R. M. Glaeser, Proc. Natl. Acad. Sci. U.S.A. 93, 10769 共1996兲. 关11兴 G. W. Robinson, C. H. Cho, and J. Urquidi, J. Chem. Phys. 111, 698 共1999兲. 关12兴 H. Poole, F. Sciortino, U. Essmann, and H. E. Stanley, Nature 共London兲 360, 324 共1992兲. 关13兴 S. S. Borick, P. G. Debenedetti, and S. Sastry, J. Phys. Chem. 99, 3781 共1995兲. 共D2兲 共D4兲 关14兴 H. E. Stanley, L. Cruz, S. T. Harrington, P. H. Poole, S. Satstry, F. Sciortino, F. W. Starr, and R. Zhang, Physica A 236, 19 共1997兲. 关15兴 I. Brovchenko, A. Geiger, and A. Oleinikova, J. Chem. Phys. 123, 044515 共2005兲. 关16兴 J. Roberts and P. G. Debenedetti, J. Chem. Phys. 105, 658 共1996兲. 关17兴 S. Sastry, P. G. Debenedetti, F. Sciortino, and H. E. Stanley, Phys. Rev. E 53, 6144 共1996兲. 关18兴 P. H. Poole, F. Sciortino, T. Grande, H. E. Stanley, and C. A. Angell, Phys. Rev. Lett. 73, 1632 共1994兲. 关19兴 A. Angell, Science 319, 582 共2008兲. 关20兴 For an interesting comparison of different scenarios, see K. Stoney, M. G. Mazza, H. E. Stanley, e-print arXiv:cond-mat/ 0805.3468. 关21兴 G. M. Bell, J. Phys. C 5, 889 共1972兲. 关22兴 Z. Yan, S. V. Buldyrev, N. Giovambattista, P. G. Debenedetti, and H. E. Stanley, Phys. Rev. E 73, 051204 共2006兲. 关23兴 A. Barros de Oliveira, G. Franceze, P. A. Netz, and M. C. Barbosa, J. Chem. Phys. 128, 064901 共2008兲. 关24兴 E. Lomba, N. G. Almarza, C. Martin, and C. McBride, J. Chem. Phys. 126, 244510 共2007兲. 关25兴 P. C. Hemmer and G. Stell, Phys. Rev. Lett. 24, 1284 共1970兲. 关26兴 H. Tanaka, Phys. Rev. Lett. 80, 5750 共1998兲. 关27兴 H. S. Ashbaugh, T. M. Truskett, and P. G. Debenedetti, J. 021203-12 PHYSICAL REVIEW E 78, 021203 共2008兲 SIMPLE THREE-STATE LATTICE MODEL FOR LIQUID WATER Chem. Phys. 116, 2907 共2002兲. 关28兴 S. Chatterjee, H. S. Ashbaugh, and P. G. Debenedetti, J. Chem. Phys. 123, 164503 共2005兲. 关29兴 T. M. Truskett, P. G. Debenedetti, S. Torquato, and S. Sastry, J. Chem. Phys. 111, 2647 共1999兲. 关30兴 L. Zoranic, R. Mazighi, F. Sokolic, and A. Perera, J. Phys. Chem. 111, 15586 共2007兲. 关31兴 G. Gompper and M. Schick, Self-Assembling Amphiphilic Systems, 1st ed., Vol. 16 of Phase Transitions and Critical Phe- nomena 共Academic Press, San Diego, 1994兲. A. Ciach, J. Høye, and G. Stell, J. Phys. A 21, L777 共1988兲. A. Ciach and J. Høye, J. Chem. Phys. 90, 1222 共1989兲. A. Ciach, J. Chem. Phys. 96, 1399 共1992兲. A. Ciach and W. T. Góźdź, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 97, 269 共2001兲. 关36兴 M. Blume, V. J. Emery, and R. B. Griffith, Phys. Rev. A 4, 1071 共1971兲. 关32兴 关33兴 关34兴 关35兴 021203-13