ABSTRACT Title of dissertation: TWO-STATE THERMODYNAMICS OF SUPERCOOLED WATER John W. Biddle, Doctor of Philosophy, 2016 Dissertation directed by: Professor Mikhail A. Anisimov Water has been called the “most studied and least understood” of all liquids, and upon supercooling its behavior becomes even more anomalous. One particularly fruitful hypothesis posits a liquid-liquid critical point terminating a line of liquid- liquid phase transitions that lies just beyond the reach of experiment. Underlying this hypothesis is the conjecture that there is a competition between two distinct hydrogen-bonding structures of liquid water, one associated with high density and entropy and the other with low density and entropy. The competition between these structures is hypothesized to lead at very low temperatures to a phase transition between a phase rich in the high-density structure and one rich in the low-density structure. Equations of state based on this conjecture have given an excellent ac- count of the thermodynamic properties of supercooled water. In this thesis, I extend that line of research. I treat supercooled aqueous solutions and anomalous behavior of the thermal conductivity of supercooled water. I also address supercooled water at negative pressures, leading to a framework for a coherent understanding of the thermodynamics of water at low temperatures. I supplement analysis of experimen- tal results with data from the TIP4P/2005 model of water, and include an extensive analysis of the thermodynamics of this model. TWO-STATE THERMODYNAMICS OF SUPERCOOLED WATER by John W. Biddle Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2016 Advisory Committee: Professor Mikhail A. Anisimov, Chair/Advisor Professor Jan V. Sengers Professor Theodore L. Einstein Professor Theodore R. Kirpatrick Professor John D. Weeks c© Copyright by John W. Biddle 2016 Dedication To Katie, my wife. ii Acknowledgments I wish to thank first of all my advisor, Mikhail Anisimov, who has provided me with genuine mentorship, encouragement, patience, a steady stream of creative ideas, and journeys across the ocean during my four years as his student. Another special thank-you is due to Jan Sengers, who served as an additional mentor to me in my research. The research group of Anisimov and Sengers has been very supportive to me, and I wish to acknowledge my group-mates: Deepa Subramanian, Daphne Fuentevilla, Jan Leys, and in particular Vincent Holten, with whom I collaborated closely and who showed me countless small ways to be more thorough and efficient in my research. I thank Fre´de´ric Caupin for hosting me for several months in Lyon, France and his research group for making me feel so welcome. I also thank my co- authors, Fernando Bresme, Rakesh Singh, and Pablo Debenedetti, who made much of this research possible. This research was supported by Grant No. 52666-ND6 of the American Chem- ical Society Petroleum Research Fund, by a Chateaubriand Fellowship from the French Republic, and by an Ann G. Wylie Dissertation Fellowship from the Gradu- ate School at the University of Maryland. At other times the research was supported by my work as a Teaching Assistant, work that was provided by the Physics Ed- ucation Research Group at the University of Maryland. I thank them for many semesters of interesting and innovative courses and professional development. More personally, my wife Katie has been incredibly supportive throughout the process of preparing and writing this dissertation. In addition, she has tolerated iii both long working hours and a period of our engagement during which we were on opposite sides of the Atlantic Ocean. I cannot thank her enough. The friends that I have made in graduate school have made my time in Mary- land much better, and I will fondly remember the bike rides, fishing trips, and flag football championships. My parents have always encouraged and supported my studies, and I thank them for that, and for showing me how an academic career can be part of a happy family life. iv Table of Contents List of Figures vii List of Abbreviations ix 1 Introduction 1 2 The Two-Structure Equation of State 9 2.1 Phase Transition and Criticality . . . . . . . . . . . . . . . . . . . . . 12 2.2 Implementation of the TSEOS . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Relationship to Scaling Theory of the Critical Point . . . . . . . . . . 17 2.3.1 The Widom Line . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3 Supercooled Aqueous Solutions 23 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 Overview of the Experimental Situation in Supercooled Aqueous So- lutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3 Thermodynamics of Criticality in Dilute Solutions . . . . . . . . . . . 26 3.3.1 Isomorphism . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.2 Implications of Constant Composition . . . . . . . . . . . . . 29 3.4 Aqueous Solutions of Sodium Chloride . . . . . . . . . . . . . . . . . 35 3.5 Liquid-Liquid Transition in Water-Glycerol . . . . . . . . . . . . . . . 42 3.6 Appendix to Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.6.1 Heat Capacity at Constant Composition . . . . . . . . . . . . 44 3.6.2 Width of the Two-Phase Region at Constant Composition . . 48 4 The LLCP in TIP4P/2005 50 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2 Computational Model and Simulation Details . . . . . . . . . . . . . 54 4.3 Description of thermodynamic properties of TIP4P/2005 water . . . . 55 4.4 Water-Like Models versus Real Water . . . . . . . . . . . . . . . . . . 62 4.5 Discussion: Does a Metastable LLPT Exist in TIP4P/2005? . . . . . 64 v 5 Negative Pressures 69 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.2 Experimental Situation . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.3 Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.5 Appendix: Fitting Parameters . . . . . . . . . . . . . . . . . . . . . . 88 6 Thermal Conductivity and its Relationship to Thermodynamics 90 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.2 Thermal Diffusivity and Thermal Conductivity . . . . . . . . . . . . . 91 6.3 Predictions of Mode-Coupling Theory . . . . . . . . . . . . . . . . . . 97 6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.5 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.7 Appendix: Details of the Mode Coupling Calculations . . . . . . . . . 113 Bibliography 117 vi List of Figures 1.1 Experimentally measured response functions in supercooled water. . . 2 1.2 A proposed phase diagram for cold and supercooled water. . . . . . . 3 3.1 An example of phase boundaries at constant composition in a super- cooled aqueous solution . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2 Schematic T -x diagram of a supercooled aqueous solution . . . . . . . 34 3.3 A comparison of our equation of state with the simulation results . . 38 3.4 Temperature of maximum density in aqueous solutions of supercooled water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.5 Suppression of the anomaly of the heat capacity in aqueous solutions of sodium chloride. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.1 Isochores in the P -T plane TIP4P/2005 . . . . . . . . . . . . . . . . 57 4.2 Densities along isobars in TIP4P/2005 . . . . . . . . . . . . . . . . . 58 4.3 Isothermal compressibility in TIP4P/2005 . . . . . . . . . . . . . . . 59 4.4 Isochoric heat capacity of TIP4P/2005 . . . . . . . . . . . . . . . . . 61 4.5 Density isobars in ST2 and mW . . . . . . . . . . . . . . . . . . . . . 63 4.6 Density in real water . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.1 Isochores in TIP4P/2005 over broad range . . . . . . . . . . . . . . . 78 5.2 Density isobars in TIP4P/2005 from negative to positive pressure . . 79 5.3 Density along isotherms in TIP4P/2005 including negative pressures . 80 5.4 Isothermal compressibility along isobars in TIP4P/2005 including negative pressures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.5 Isothermal compressibility along isotherms in TIP4P/2005 including negative pressures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.6 Isobaric heat capacity along isotherms in TIP4P/2005 including neg- ative pressures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.7 Isobaric heat capacity along isobars in TIP4P/2005 including negative pressures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.8 Lines of thermodynamic extrema in TIP4P/2005 compared to the extended TSEOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.9 Lines of thermodynamic extrema in the ST2 model . . . . . . . . . . 86 vii 5.10 A comparison between the overall behavior of the ST2 and TIP4P/2005 models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.1 Thermal-diffusivity data for supercooled water . . . . . . . . . . . . . 92 6.2 Heat-capacity data for supercooled water . . . . . . . . . . . . . . . . 93 6.3 Thermal conductivity data for supercoled water . . . . . . . . . . . . 95 6.4 Viscosity data for supercooled water . . . . . . . . . . . . . . . . . . 99 6.5 A comparison of the background thermal conductivity of supercooled water and the critical enhancement predicted by mode-coupling theory102 6.6 Density, compressibility, and speed of sound of TIP4P/2005 at 0.1, 70, and 120 MPa as a function of temperature . . . . . . . . . . . . . 104 6.7 Normalized heat-flux correlation functions in TIP4P/2005 . . . . . . 107 6.8 Thermal conductivity of TIP4P/2005 at 0.1, 70, and 120 MPa . . . . 108 6.9 Comparison of the simulated temperatures of minimum thermal con- dutivity in TIP4P/2005 with theoretical predictions . . . . . . . . . . 109 6.10 Translated thermal conductivity of TIP4P/2005 . . . . . . . . . . . . 110 6.11 Comparison of the Bridgman formula in TIP4P/2005 and in real water111 viii List of Abbreviations cxc Coexistence Curve GROMACS Groningen Machine for Chemical Simulations HDA High-Density Amorphous HDL beta High-Density Liquid IAPWS International Association for the Properties of Water and Steam LDA Low-Density Amorphous LDL Low-Density Liquid LINCS Linear Constraint Solver LLCP Liquid-Liquid Critical Point LLPT Liquid-Liquid Phase Transition LSI Local Structure Index LVS Liquid-Vapor Spinodal mW Monatomic Water Model ST2 Revised Pair Potential TIP4P Four-Point Transferable Intermolecular Potential TIP5P Five-Point Transferable Intermolecular Potential TMD Temperature of Maximum Density TSEOS Two-Structure Equation of State UCM Universidad Complutense de Madrid ix Chapter 1: Introduction Water has been called the “most studied and least understood” of all liquids [1]. Cold water is particularly noted for its anomalous behavior, and probably the best- known of water’s thermodynamic anomalies is the maximum in water’s density. This maximum occurs just above the freezing point (4 ◦C at atmospheric pressure) [2] and has proven crucial for the development of life on earth. The thermodynamic anomalies of water become more pronounced—and new ones manifest themselves— when water is supercooled, that is, cooled below its melting temperature while remaining liquid. This intriguing metastable state of water occurs in a variety of natural settings, including clouds, where liquid water has been observed at temperatures approaching −40◦C [3], and the bodies of insects that use supercooling to remain active in the winter [4]. Its properties are therefore of interest to a wide range of scientists and engineers. Supercooled water can also be created in the laboratory, and for physicists and chemists, perhaps the greatest promise lying in the study of supercooled water is the light that it can shed on the nature of liquid water generally. The most striking thermodynamic anomalies in supercooled water are those of the thermodynamic response functions. The isobaric specific heat capacity cp 1 Figure 1.1: Experimentally measured response functions in supercooled water. The for- mulation of the International Association of the Properties of Water and Steam (IAPWS) is shown in the stable region by way of comparison. [5,6] and isothermal compressibility κT [7] increase dramatically upon supercooling, as can be seen in Fig. 1.1. The correlation length characterizing fluctuations of density increases markedly as well [8]. In 1992, Poole et al. proposed a coherent and particularly fruitful account of the thermodynamic anomalies of water: that below the line of homogeneous nucleation is a first-order phase transition between two liquid states distinguished by their different densities and entropies, called the high-density liquid (HDL) and low-density liquid (LDL). This transition line would terminate at a critical point, and proximity to this critical point could explain the apparent divergence of the response functions [9]. This liquid-liquid phase transition (LLPT) and liquid-liquid critical point (LLCP) are hypothesized to lie in the region of the phase diagram that is inaccessible to experiment due to homogeneous ice nucleation, as shown in Fig. 1.2. It should be noted that while kinetic factors make it impossible to carry out experiments on the liquid state in this region of the phase diagram, the liquid state is not thermodynamically unstable there. The hypothesized LLPT is a phase transition between two liquid states, i. e., it does not entail any long-range symmetry breaking. Rather, like the transition 2 Figure 1.2: A proposed phase diagram for cold and supercooled water. The melting line and the line of homogeneous nucleation are experimentally known; the LLPT and LLCP are hypothesized, with the positions proposed in Ref. [10] used in this figure. Figure prepared by Vincent Holten. 3 between the liquid and vapor phases of a single-component fluid, the LLPT is con- ceptualized as a line of first-order transitions characterized by discontinuous changes in the density and entropy of water, terminating at a (LLCP) beyond which all thermodynamic properties vary smoothly. In order to account for this, researchers have proposed the existence of two distinct hydrogen-bonding structures in water: a high-density, high-energy, high-entropy structure (“structure A”) that is prevalent at room temperature and above, and a low-density, low-energy, low-entropy structure (“structure B”) that is favored at lower temperatures (and pressures).1 Depending on the nature of the interactions between molecules in the different structures— in thermodynamic terms, the non-ideality of the mixture—there could be a phase transition at which the liquid changes discontinuously from a liquid rich in state A (high-density liquid or HDL) to a liquid rich in state B (low-density liquid or LDL). On the other hand, the fraction of water in each state might vary continuously throughout the phase diagram, in which case no phase transition or critical point would be observed. It is now broadly agreed, based on evidence derived primarily from scattering experiments, that two distinct local structures exist in cold and supercooled wa- ter [13–16]. Simulation studies on a variety of simulation models yield similar find- ings. The precise nature of the low-density structure remains a subject of debate, 1In fact, the idea of two distinct structures in liquid water is an old one: in the late 19th century both Harold Whiting [11] and Wilhelm Roentgen [12] suggested that the density maximum in water could be explained as a kind of “pre-freezing” effect associated with the (also anomalous) fact that water expands upon crystallization. They each suggested that in cold water, there might be small pockets of ice, too small to destroy the liquid nature of cold-stable water as a whole (we might anachronistically say that no long-range symmetry-breaking has yet occurred at this stage), but numerous enough to effect a decrease in the density upon cooling. This is not plausible if taken literally, but it has remained useful and intriguing as a heuristic. 4 but an increase in the number of molecules with four nearest neighbors, and indeed upon local tetrahedral arrangements (without any long-range symmetry-breaking) is seen upon cooling in the TIP4P/2005 [17], ST2 [18], and mW [19] models of water. The existence of two competing local structures does not by itself guarantee the existence of a phase transition, however. The ST2 [20], TIP4P/2005 [17], and mW [21] models are all described very well by two-structure thermodynamics, and while the LLPT and LLCP have been confirmed in the ST2 model [22, 23], they are absent from the mW model [21]. Liquid-liquid separation in the TIP4P/2005 model continues to be debated in the literature, with this Dissertation presenting further evidence for the existence of the LLPT and LLCP in that model in Chapter 4. Experimental evidence for the LLPT and LLCP in real water has so far been indirect. Mishima’s finding of a phase transition in glassy water between two amor- phous states of different density has lent further support to this thesis [24], as have Mishima and Stanley’s intriguing experiments on the melting lines of metastable phases of ice [25, 26]. In short, these latter experiments investigate the melting lines of metastable phases of ice, and find kinks in those lines just below the ho- mogeneous nucleation line in the liquid. Morevoer, equations of state based on the two-structure model and including a liquid-liquid critical point have provided in- creasingly accurate accounts of the thermodynamic properties of water over the last several years [27,28], and the recent model of Holten and Anisimov shows excellent agreement with the thermodynamic data with fewer adjustable parameters than any model so far [10]. In 2014, Holten et al. published an extended form of this model that fits pressures up to 400 MPa, and this version is now the basis of the most 5 recent formulation for supercooled water employed by the International Association for the Properties of Water and Steam (IAPWS) [29]. Other hypotheses than the LLCP have been proposed to account for the anomalies of water at low temperatures. In the singularity-free hypothesis, CP and κT pass through finite maxima at temperatures too low to be accessed by ex- periment, but do not diverge, and there is no first-order phase transition. Although it lacks the sensational punch of the LLPT and LLCP, this hypothesis is in fact con- sistent with the two-structure hypothesis, and two-structure thermodynamics might still be invoked to explain water’s anomalies (as in the case of the mW model of water [21]). A very different suggestion was made by Speedy in 1982: that the absolute limit of stability of the liquid state retraces to positive pressures, eventually acquiring a steep negative slope and lying just below the homogeneous nucleation limit at atmospheric pressure [30]. He proposed that the locus of temperatures of maximum density (TMD locus) and the apparent divergences in the response functions κT and CP were signs of an approach to this spinodal, at which κT and CP are expected to diverge. Debenedetti and Stanley have since argued convincingly that the liquid- vapor spinodal (LVS) does not reach a pressure above the saturation pressure at low temperatures because the intersection of the LVS and the binodal must be a critical point, and a second liquid-vapor critical point, occurring at a temperature far below the triple-point temperature, is fantastically improbable [31]. However, the related hypothesis of a line of absolute instability of the liquid state other than the LVS (Speedy’s original work [30] did not specify the phase toward which the liquid 6 might become unstable) remains a possibility as an explanation for the anomalies of supercooled water. It is also possible that the LVS goes through a minimum and returns to higher but still negative pressures, a prospect discussed further in Chapter 5. This Dissertation takes the ideas and basic structure of Holten and Anisimov’s TSEOS to investigate a new set of anomalies in supercooled water, both dynamic and thermodynamic. In this way, I hope that it will contribute to a coherent, unified explanation of anomalies in liquid water at low temperatures in terms of the conjecture of two distinct hydrogen-bonding structures in liquid water. The Dissertation also seeks to shed light on the question, “is there a liquid-liquid phase transition and a second critical point in supercooled water?” That question will not be answered conclusively in this Dissertation, but the proof of the pudding is in the eating, and the hypothesis of the liquid-liquid critical point continues its success in explaining the anomalies of supercooled water, as shown below. After introducing the equation of state (Chapter 2) the Dissertation presents an analysis of aqueous solutions, showing that the TSEOS can describe even counter-intuitive behavior in aqueous solutions with only slight modifications (Chapter 3). Chapter 4 examines the thermodynamics of the TIP4P/2005 model and presents evidence that the LLCP and LLPT exist in this model. Chapter 5 addresses questions surrounding water at negative pressure. By incorporating a liquid-vapor spinodal into the TSEOS, it presents an equation of state that matches new, negative-pressure data on the TIP4P/2005 model and sheds light on the behavior of real water. Chapter 6 presents calculations arguing that the anomaly in the thermal conductivity of supercooled 7 water is of thermodynamic origin, and presents simulation studies to support this claim. 8 Chapter 2: The Two-Structure Equation of State Water is a single-component fluid in the most intuitive and important sense: it comprises molecules of only one type, H2O. However, the two-structure hypothesis conjectures that it is also in a certain sense a mixture. More specifically, H2O can form two distinct supramolecular structures, and clusters of molecules in each of the supramolecular structures are able to mix with each other. Consequently, considera- tions on a length-scale slightly longer than the molecular but still microscopic suggest that the thermodynamics of water will be similar to that of a binary mixture. The research presented in this Dissertation comprises several phenomenological studies of supercooled water. In order to model actual data, whether from experiment or simulation, each of the studies presented here makes use of an equation of state based on the two-structure conjecture. The basic structure of these equations of state is the same, and that structure is described in this Chapter. Small variations in structure are made for each study, and the values of several parameters are fitted to best match the data at hand in each case (the values of the parameters will be very different in a study of simulation data than in one of real water, for example). This family of equations of state, as well as each individual version, is called the two-structure equation of state (TSEOS). 9 Because experiments on water are typically carried out at a given temperature and pressure, it is most convenient to implement the TSEOS as an expression for the Gibbs energy of water. The Gibbs energy is also advantageous because a field- dependent potential lends itself easily to an analysis in terms of scaling fields, as will be discussed in a later section. The Gibbs energy of a binary mixture can be described (as described in, e. g., Ref. [32]) as the sum of the Gibbs eneries of each of the two species as a function of temperature and pressure, weighted by the molar fraction of each species, and a Gibbs energy of mixing. Using φ to represent the molar fraction of state B, one can write (1− φ)GA(T, P ) + φGB(T, P ) +Gmixing. (2.1) Here, GA and GB are the Gibbs energy of pure state A and pure state B at (T, P ) (if they could somehow be isolated and studied alone), and not the chemical potential of the respective species in actual water at (T, P ). In practice, for reasons that will become clear, it is easier to work with GA(T, P ) + φGBA(T, P ) +Gmixing, (2.2) where GBA is the difference in Gibbs energies GB − GA. In an ideal mixture, the entropy of mixing takes the Lewis-Randall form [32] SLR = φlnφ+ (1− φ)ln(1− φ), (2.3) 10 giving rise to a term RT [φlnφ+ (1− φ)ln(1− φ)] (2.4) in the Gibbs energy, a term which will sometimes be denoted by I. But the mixture is not, in general, ideal. The non-ideal part of the mixture is called W , and it will be a function W (T, P, φ). The φ-dependence of the W is modeled with a simple, symmetric form: W = ω(T, P )φ(1− φ). In general, then, the TSEOS has the form: G = GA(T, P )+φGBA(T, P )+RT [φlnφ+ (1− φ)ln(1− φ)]+ω(T, P )φ(1−φ). (2.5) The Gibbs energy has so far been constructed using the thermodynamic theory of binary mixtures as G(T, P, φ). But in order to implement the TSEOS, one must also consider the interconvertibility of the two species. It is obvious that distinct forms of hydrogen bonding within the same liquid are, in principle, interconvertible. The outstanding question is the timescale over which this interconversion takes place. Experimental evidence is clear that any relaxation timescale in supercooled water is many orders of magnitude shorter than the experimental timescale for the measure- ment of thermodynamic properties or transport properties [33–36]. Consequently, φ is not a thermodynamic variable on a level with T and P . Rather, for a given (T, P ), the possible values of φ are sampled ergodically, and the thermodynamic value of φ at a given (T, P ) is the value that minimizes the Gibbs energy. When applying, for example, the Gibbs phase rule to supercooled water as described by the TSEOS, the number of intensive degrees of freedom is two, not three. 11 The derivatives of the Gibbs energy with respect to temperature and pressure provide expressions for the thermodynamic properties (let us call a generic ther- modynamic property X) as functions X(T, P ;φ). In light of the considerations in the previous paragraph, they must be evaluated in two steps. First, one solves the equation ( ∂G ∂φ ) T,P = 0 (2.6) to find the equilibrium value of φ, φe. Next one evaluates the expression X(T, P ;φe). In the event, Eq. (2.6) is a transcendental equation in φ and does not have a closed- form solution, so numerical methods must be used in order to implement the TSEOS. 2.1 Phase Transition and Criticality The first derivatives of the Gibbs energy evaluated at constant φ, i. e., ( ∂G(T,P ;φ) ∂T ) P,φ and ( ∂G(T,P ;φ) ∂P ) T,φ , are continuous everywhere. There can, however, be discontinu- ities in the function φe(T, P ), and a locus of such discontinuities constitutes a line of first-order transitions, possibly terminating in a critical point. In particular, a sufficiently strong non-ideality in the mixing term of the Gibbs energy may lead to a function G(φ) with more than one local minimum. If this is the case, the global minimum of course defines the value of φe. When the two minima are of equal depth, however, two-phase coexistence is possible, which defines the LLPT. A first-order phase transition occurs upon crossing such a point of two-phase coexistence such that one one side, one of the minima is the global minimum and on the other side, the other is the global minimum, leading to a discontinuity in φe(T, P ). 12 The possibility of a phase transition and a critical point can be determined from the derivatives of the Gibbs energy as a function of φ: ( ∂G ∂φ ) T,P =GBA(T, P ) +RT · ln ( φ 1− φ ) + ω(T, P )(1− 2φ), (2.7)( ∂2G ∂φ2 ) T,P =RT ( 1 φ(1− φ) ) − 2ω(T, P ), (2.8)( ∂3G ∂φ3 ) T,P =RT ( 1− 2φ φ2(1− φ)2 ) . (2.9) The conditions of criticality, then, are ω = 2RT, (2.10) GBA = 0, (2.11) and under these conditions, φc = 0.5. (2.12) As discussed above, two-phase coexistence and a LLPT are possible only if the the function G(φ) has two minima, which requires ω > 2RT. (2.13) Because I and W are symmetric in φ, two distinct values of φ can only yield the same value of the Gibbs energy if condition 2.11 is also met, so GBA = 0 also locates the LLPT. Where ω < 2RT there is no phase transition. In this region GBA = 0, 13 the extension of the LLPT into the one-phase region, is the locus of points at which φ = φc and fluctuations of φ are at a maximum, as can be seen by minimizing the RHS of Eq. (2.8). This locus is known as the Widom line, and is discussed further in the Section on scaling theory. One could construct a version of the TSEOS in which the non-ideality is never strong enough to produce criticality or a first-order phase transition, as was done in Ref. [21] for the mW water model. But the TSEOS was developed to explore the idea of the LLCP and LLPT, and both real water and the TIP4P/2005 model are best described by an equation that includes both a LLPT and a LLCP. Thus, since the values of certain parts of the TSEOS at criticality are determined by thermo- dynamics (to wit: Eqs. 2.10 and 2.11), it is convenient to work in dimensionless variables that are given with respect to the critical parameters. In particular Tˆ = T/Tc, (2.14) Pˆ = P/ρcRTc, (2.15) ∆Tˆ = (T − Tc)/Tc, (2.16) ∆Pˆ = (P − Pc)/ρcRTc. (2.17) In general, properties made dimensionless by the appropriate combination of the critical parameters and the universal gas constant are represented by a circumflex, so Gˆ = G/RTc, etc. From now on, ω will be presented in units of RTc. 14 2.2 Implementation of the TSEOS The Gibbs energy of pure structure A, GA, is treated phenomenologically: GˆA = ∑ m,n cmn∆Tˆ m∆Pˆ n, (2.18) with {cmn} as adjustable coefficients. c10 is associated with the absolute value of the entropy at the critical point and has no physical significiance, so it is implicitly set to zero. c01 is related to the volume at the critical point, so c01 and ρc are not independent parameters. As discussed in the previous section, The LLPT, the LLCP, and the Widom line are all characterized by condition 2.11, so it is convenient also to implement GˆBA as a power series in ∆Tˆ and ∆Pˆ . In the most general form that is used for the TSEOS, GˆBA = λ(∆Tˆ + a∆Pˆ + b∆Tˆ∆Pˆ + d∆Pˆ 2 + f∆Tˆ 2), (2.19) or, in some versions, GˆBA = Tˆ λ(∆Tˆ + a∆Pˆ + b∆Tˆ∆Pˆ + d∆Pˆ 2 + f∆Tˆ 2). (2.20) For convenience, the rest of this Dissertation will refer to Eq. (2.19) as the normative form of GBA. In this case, λ gives the difference in molar entropy Sˆ between the two struc- 15 tures at the critical point. a gives the slope −(dPˆ /dTˆ ) of the of the LLPT and Widom line at the critical point, and the product λa gives the difference in mo- lar volume Vˆ between the two structures at the critical point, and b, d, and f are proportional to the difference in give the difference in thermal expansivity αˆP , isothermal compressibility κˆT , and isobaric heat capacity CP at the critical point. 1 In addition, b, d, and f determine the curvature of the LLPT and Widom line in the (T, P ) plane. The parameters λ, a, b, d, and f are adjustable parameters whose values are chosen based on the best available fit to the data. In practice, not all of the terms may be necessary; real water and indeed aqueous solutions at posi- tive pressures can be fitted with only λ, a, and b, as can TIP4P/2005 at positive pressures [17]. However, extension of the equation of state to negative pressures, especially if low-temperature data are to be considered, requires that the differences in compressibility and heat capacity between the two states be accounted for. As discussed above, the non-ideality of mixing between the two structures is assigned the simiple, symmetric form ω(T, P )φ(1− φ). (2.21) According to Eq. (2.10), we must have ω = 2 at the critical point. In putting the TSEOS in roughly the form that I have used in this research, Holten and Anisi- mov conjectured that the non-ideality might take the form2 ω = Tˆ (2 + ω0∆Pˆ ). 1It should be clarified that these parameters give the differences in properties of each of the pure structures A and B at the critical point. All differences in properties between HDL and LDL vanish at the critical point. 2This is expressed according to the definition ω relative to the Gibbs energy that is used in this 16 Later work found that the ST2 model of water was better matched if ω took the form 2 + ω0∆Pˆ . In the first case, the non-ideal Gibbs energy is proprotional to the non-ideality in the entropy of mixing and there is no non-ideal enthalpy of mixing. This has been dubbed an “entropy-driven” LLPT. In the second case, the non-ideal Gibbs energy is equal to the non-ideal enthalpy of mixing, and there is no non-ideal entropy of mixing. This has been dubbed an “energy-driven” LLPT. In both cases, however, non-ideal energy and volume of mixing are both present. In a particular implementation of the TSEOS, either form for ω, or a linear combination of the two, can be used. In this case, ω = [ Tˆ (1− δ) + δ ] (2 + ω0∆Pˆ ). (2.22) 2.3 Relationship to Scaling Theory of the Critical Point For a lucid and complete treatment of critical pheonomena, in particular scal- ing and universality, I recommend in particular Michael Fisher’s lectures, collected in Ref. [37]. Landau’s treatment in Ref. [38] is also worthwhile. This section of the dissertation presumes some knowledge of critical phenomena on the part of the reader and outlines the main features of the TSEOS in terms of the scaling theory of phase transitions. In the theory of critical phenomena, the thermodynamic potential can be sepa- rated into a regular background part and a critical part, the latter of which contains dissertation. Ref [10] uses a definition of ω by a factor of RT . 17 all of the singular behavior associated with the critical point. The critical part of the potential is associated with the dependent scaling field h3, which can be expressed in terms of two independent scaling fields: the ordering field h1 and the second, “thermal” field h2. h3 = Φ1dh1 + Φ2dh2, (2.23) where Φ1 is the order parameter [38] and Φ2 is the second (weakly fluctuating) scaling density. Because we use the molar Gibbs energy G(T, P ) as the thermodynamic po- tential, we start with h3 = Gˆ− GˆA. (2.24) Once the scaling densities have been identified, regular terms can be subtracted from h3 so that the scaling densities are set to zero at the critical point. In a “simple-scaling” model, each thermodynamic field is associated with one of the scaling fields [39]. Such a model is indeed too simple for the present case. In principle, all physical fields must be considered for a complete treatment of a fluid critical point (“complete scaling” [40]). However, if the analysis is confined to a mean-field approximation, then a mixture of two thermodynamic fields suffices for a simple fluid. For the reasons of convenience mentioned in the introduction, we choose to express the independent scaling fields in terms of T and P . This is consistent with an intermediate approach between simple and complete scaling known as “revised 18 scaling” [10,39,41]. Thus we have: h1 = a1∆Pˆ + a2∆Tˆ , (2.25) h2 = b1∆Tˆ + b2∆Pˆ . (2.26) The order parameter Φ1 in this model is associated with the fraction of molecules in the low-density state, such that Φ1 = φ− φc. (2.27) Condition 2.10 locates the LLPT and LLCP, and in the one-phase region gives the locus of maximum fluctuations of the order parameter, from which we find h1 = G BA, (2.28) which is consistent with the condition Φ1 = ( ∂h3 ∂h1 ) . (2.29) In practice, a linear approximation of GˆBA may be used. The inclusion of a dependence upon pressure in the expression for ω—this indicates non-ideality in the volume of mixing between the two species, and is a necessary part of the TSEOS—means that the thermal field h2 includes a term proportional to the pressure. In the athermal case, h2 = −∆Pˆ , while for the regular- 19 solution or a mixed case h2 is a linear combination of pressure and temperature. In general, h2 = −∆Pˆ + 2δ ω0 ∆Tˆ . (2.30) This is a difference between the thermodynamic environments of the liquid-liquid and the vapor-liquid critical points. In the latter case, simple scaling models have h2 = ∆Tˆ , while revised or complete scaling models include contributions from the other fields. In particular, the lattice-gas model of the vapor-liquid critical point has h2 = ∆Tˆ [39]. From these definitions, expressions can be derived for the scaling susceptibili- ties, χ1, χ2, and χ12, and the physical response functions can be expressed as linear combinations of these susceptibilities, expressions that are used in the analysis of Chapter 3. Of particular interest will be χ1 = ( ∂Φ1 ∂h1 ) , the strongly divergent scal- ing susceptibility. χ1 diverges strongly at the critical point, and κT , and CP each contain terms proportional to χ1, which accounts for the strong divergence of these properties. 2.3.1 The Widom Line The definition of the Widom line deserves special consideration here. Accord- ing to scaling theory, the Widom line corresponds to h1 = 0 and Φ1 = 0 in the one-phase region, and the definitions of the Widom line and h1 are thus consis- tent. In the literature, however, various definitions of the Widom line are used. The reason for this is that in the one-phase region, κT and CP go through finite 20 maxima but do not diverge, and asymptotically, the loci of these maxima approach the line of maximum fluctuations in the order parameter. Because of the relative ease of measuring these maxima, one or the other is often referred to as the Widom line. In the vicinity of the liquid-vapor critical point, where the order parameter is ρ− ρc to an excellent approximation, the locus of κT maxima is a suitable proxy for the Widom line. The thermodynamics of the liquid-liquid critical point are more complicated, however. Moreover, the TSEOS aims to describe water over a broad range of temperatures and pressures. Far from the critical point the loci of CP and κT maxima may have very different shapes, and both may differ from the locus of maximum fluctuations in the order parameter, as is discussed in detail Chapter 5. Consquently, neither is a suitable proxy for the Widom line. In this Dissertation, the Widom line refers to the locus of maximum fluctuations of the order parameter, which in the TSEOS corresponds to the condition 2.10 and P < Pc. 2.4 Conclusion This Chapter has laid out the structure of the TSEOS and the thermodynamic considerations that motivate it. The Chapters that follow document applications of the TSEOS to particular systems and problems, namely: the thermal conductivity of supercooled water, supercooled aqueous solutions, the TIP4P/2005 model in the critical region, and the TIP4P/2005 model at negative pressures. Each project has required at least slight modifications or extensions of the TSEOS, as documented in the respective Chapters, most notably the addition of a liquid-vapor spinodal to 21 study negative pressures in Chapter 5. The core of the TSEOS has remained the same throughout, however. 22 Chapter 3: Supercooled Aqueous Solutions Author’s Note: this chapter incorporates large selections of text from Ref. [42], which I published in 2014 along with my colleague Vincent Holten and my thesis advisor Mikhail Anisimov. “We” refers to myself and these co-authors. –JWB 3.1 Introduction With the hypothesized liquid-liquid critical point lying just out of reach of experiment in pure water, the addition of solutes to water shows obvious promise as a way to shed light on the question. The addition of a solute to water depresses the freezing point, and this disruption of the freezing process also lowers the tem- perature of homogeneous nucleation. As a result, temperature and pressures that are inaccessible in the pure fluid due to homogeneous ice nucleation become acces- sible in aqueous solutions. In the best-case scenario, addition of a solute would bring the LLPT and LLCP themselves into the experimentally accessible region, where they could be directly observed. Now, even were a LLPT to be observed directly and incontrovertibly, a thermodynamic study would still be needed in or- der to determine whether it emanated from a liquid-liquid transition in water, or whether it came about due to non-ideality in the mixing between the water and the 23 solute. Moreover, claims of direct observation of a liquid-liquid transition remain controversial. With that said, studies of solutions in which the LLPT and LLCP cannot be observed can still shed light on the nature of water. The thermodynamics of dilute aqueous solutions near the vapor-liquid critical point is a well-developed science, and is easily applicable to aqueous solutions near the LLCP by the principle of critical- point universality. In this Chapter, I will present an account of the experimental data on NaCl(aq), showing that the response of water anomalies to the addition of NaCl is best accounted for by two-structure thermodynamics, in particular by the hypothesis of a LLPT and LLCP. These concepts yield a convincing qualitative picture of the behavior of a variety of aqueous solutions. Moreover, the TSEOS can give a striking quantitative account of the experimental density and heat-capacity data in supercooled NaCl(aq) with only simple changes to the background terms. This chapter also analyzes experimental claims regarding the direct observation of LLPT in supercooled water-glycerol solutions. 3.2 Overview of the Experimental Situation in Supercooled Aqueous Solutions It has been known since the mid-19th century that the addition of sugars and salts to water depresses the temperature of maximum density. Moreover, the rela- tionship between temperature of maximum density and concentration of solute is nearly linear, a relationship known as “Despretz’s Law”. Upon addition of NaCl, 24 homogeneous nucleation moves to lower temperatures and pressures, while the ho- mogeneous nucleation curve in the P -T plane keeps roughly the same shape. The experimental information on the thermodynamic properties of supercooled aqueous solutions of salts, in particular of NaCl, is very limited. The available data are the isobaric heat capacity measurements of Archer and Carter [6], and the density measurements of Mironenko et al. [43], both at atmospheric pressure. The density data confirm the depression of the TMD and show an increase in density upon addition of a solute [43]. Archer and Carter found that upon addition of NaCl, the anomalous increase in the heat capacity moves to lower temperatures and decreases in magnitude, practically disappearing for concentrations greater than 2 mol/kg. In the initial publication announcing these results, Archer and Carter interpreted them as evidence against the existence of the liquid-liquid critical point. However, our study reaches the opposite conclusion: in light of what can be inferred about the movement of the locus of liquid-liquid phase transitions upon addition of NaCl, this behavior of the heat capacity is precisely what thermodynamics predicts if the anomaly in the heat capacity is indeed associated with a liquid-liquid critical point. In addition, Murata and Tanaka have reported direct visual observation of a liquid-liquid transition in supercooled aqueous solutions of glycerol [44]. They have argued that the formation of a more stable liquid phase in this solution may occur by two alternative types of kinetics: nucleation and spinodal decomposition. They have also claimed that the transition is mainly driven by the local structuring of water rather than of glycerol, suggesting a link to the hypothesized liquid-liquid transition in pure water. However, they did not observe two-phase coexistence, 25 leading them to claim that the transition is “isocompositional” and the nucleation and spinodal decomposition occurs “without macroscopic phase separation.” Murata and Tanaka’s results have been questioned by other researchers [45], and controversy concerning the liquid-liquid transition in supercooled water-glycerol persists. Still, it is worthwhile to examine the observations and analysis of Murata and Tanaka in light of two-structure thermodynamics. 3.3 Thermodynamics of Criticality in Dilute Solutions 3.3.1 Isomorphism There is a well-developed approach to treating the thermodynamics of mixtures near their critical points, closely related to the principle of critical-point universality, known as “isomorphism” [46]. It has been postulated that upon the addition of solute, the form of the equation of state remains unchanged under the condition of constant thermodynamic fields, including chemical potentials [46–52]. The molar Gibbs energy of a binary system is expressed through the chemical potentials of the two components, solvent and solute, µ1 and µ2 as G = (1− x)µ1 + xµ2 = µ1 + x(µ2 − µ1), (3.1) where x is the mole fraction of solute and δ = (∂G/∂x)T,P = µ2 − µ1 is the ther- 26 modynamic field conjugate to x, and dG = V dP − SdT + δdx. (3.2) In the theory of isomorphism, the chemical potential of the solvent in solution, µ1 = G − xδ, which is the same in the binary-fluid coexisting phases, replaces the concentration-dependent Gibbs energy as the relevant thermodynamic potential such that dµ1 = V dP − SdT + xdδ. (3.3) There are two alternative cases of fluid-fluid separation in a binary solution. One is caused by non-ideality of mixing between the two species. The other is the offspring of a transition in the pure solvent. The former case is typical for liquid- liquid separation in weakly compressible binary solutions, while the latter case is observed as fluid-fluid transitions stemming from the vapor-liquid transition in the pure solvent. For the second case, the mixing of the two species in the solution does not need to be non-ideal, as the phase separation in the solution is a continuation of the phase-separation in the pure solvent [53]. This work models liquid-liquid transitions in supercooled aqueous solutions as instances of this latter case. The stability criterion in fluid mixtures can be written in a form convenient for the latter case: ( ∂P ∂V ) T,δ = ( ∂P ∂V ) T,x + ( ∂P ∂x )2 T,V ( ∂x ∂δ ) T,V ≤ 0. (3.4) 27 We assume that the form of the isomorphic thermodynamic potential, which is the chemical potential of water in solution µ1, is the same as that of the Gibbs energy of pure solvent (water) given by the TSEOS: µ1 = µA + φµBA = µ1A + φµ1BA+ RT [φ lnφ+ (1− φ) ln(1− φ) + ωφ(1− φ)] (3.5) where µ1A is the chemical potential of pure state A, and µ1BA = µ1B − φµ1A, is the difference in the chemical potentials between the pure states. There are two key differences difference from the thermodynamics of the pure solvent. One is that is that the critical parameters, Tc and Pc, are functions of the chemical potential difference δ = µ2−µ1. The other is that the dependence of the Gibbs energy of pure structure A (the “background”) upon concentration is accounted for by making the {cmn} of Eq. 2.18 dependent on the concentration x of NaCl, so that GA = ∑ m,n cmn(1 + βmnx+ γmnx 2)∆TˆmPˆ n. (3.6) The {cmn} and all other parameters for pure water used in this research are un- changed from the mean-field version of the TSEOS that is presented in Ref. [10]. Note that the chemical potentials of the two components, solvent and solute, are not the chemical potentials µA and µB of two alternative structures in water; µ1 and µ2 are controlled by the concentration of solution and interactions between the 28 solvent and solute molecules. 3.3.2 Implications of Constant Composition Properties at constant composition are found from the derivatives of the Gibbs energy. In order to obtain them, we must perform a Legendre transformation G = µ1 +µx. In addition, we use an approximation, called the critical-line condition [52], that requires along the critical line in solution x = xc = e δ/RT . (3.7) When a critical line emanates from the critical point of the pure solvent, the prin- cipal thermodynamic property that controls the behavior of solutions at constant composition is the so-called Krichevskii parameter defined in the dilute-solution limit as [54] K = lim x→0 ( dP dx ) c,cxc = dTc dx [ dPc dTc − ( dP dT ) c,cxc ] . (3.8) where (dP/dT )c,cxc is the slope of the line of liquid-liquid coexistence at the critical point of pure solvent. The derivatives dTc/dx and dPc/dTc determine the initial slopes of the critical line in the (T, P, x) space. Thus, since in the limit of the solvent critical point (dP/dx)c,cxc = (∂P/∂x)T,V , the absolute stability limit (spinodal) can be formulated through the Krichevskii parameter as ( ∂P ∂V ) T,δ = ( ∂P ∂V ) T,x +K2 ( ∂x ∂µ ) T,V = 0. (3.9) 29 Along the critical line the inverse compressibility at constant composition vanishes in the pure solvent limit as ( ∂P ∂V ) T,x = −xK2 → 0. (3.10) Note that if the fluid phase separation in solutions stems from the transition in pure solvent, the stability criterion of a pure fluid smoothly transforms into the stability criterion of a solution. If the solute dissolves more favorably in the higher-pressure phase, then Le Chatelier’s principle indicates that the phase transition at a given temperature will move to lower pressures, and vice versa. The direction in which the phase transition pressure moves at constant temperature, positive or negative, indicates the sign of the Krichevskii parameter. This does not necessarily mean that the sign of dPc/dx determines the sign of the Krichevskii parameter, as this derivative is influenced both by the movement of the transition line in the P -T plane and the movement of the critical point along that line. In binary fluids, the critical point becomes a critical line and the phase- transition line becomes a surface of two-phase coexistence in the “theoretical” (T, P, δ) space. However, in the “experimental” (T, P, x) space the behavior of thermody- namic properties evaluated at constant composition will, in general, be different from that of the corresponding properties in the pure solvent and from that of the isomorphic properties in solutions (evaluated at constant chemical-potential differ- ence δ). Remarkably, the nature and magnitude of this difference depend primarily 30 on the value of the Krichevskii parameter [51, 52]. In particular, the concentration gap in the (T, x) plane at constant temperature can be found from Eq. (3.3) as ( dP dδ ) T,cxc = ∆x ∆V , (3.11) where ∆V is the difference in volume of the coexisting phases. In the dilute-solution approximation, dx/dδ = x/RT , so the concentration gap at the first-order transition and constant pressure can be evaluated through the Krichevskii parameter, ∆V , and the slope of the transition line as ∆x ' −xK∆V RT , (3.12) where x = xc in accordance with the critical-line condition (3.7). Correspondingly, the temperature gap at constant composition can be evaluated as (see the Appendix to this chapter) ∆T ' xK2 ∆V RT ( dT dP ) c,cxc . (3.13) In the solvent-critical-point limit, ∆V ∝ x and the phase diagram develops a so- called “bird’s beak” where the concentration gap vanishes to first order in x and the two branches of the transition merge with the same tangent [54,55]. The above-described thermodynamics explains possible phase behavior of a supercooled aqueous solution with a critical line emanating from the pure solvent (water) critical point, as shown in Figs. 3.1 and 3.2. Only in a special case, when 31 p re s s u re temperature Pc' Tc Pc Tc' C C' 1 2 4 3 Figure 3.1: An example of phase boundaries at constant composition in a supercooled aqueous solution exhibiting a liquid-liquid transition between HDL and LDL. The black curve is the liquid-liquid transition in pure water, terminated at the critical point C. The critical line is shown by solid red with the critical point of the solution labeled C′. The blue curve shows the appearance of the first droplet of of LDL. The green curve shows the disappearance of the last droplet of HDL. The blue and green dashed curves are the thermodynamic stability limits of HDL and LDL, respectively. The shaded region shows where HDL forms by nucleation. The dotted lines labeled 1,2,3, and 4 show different thermodynamic paths as explained in the text. 32 the the critical line and the liquid-liquid transition line merge with the same slope in the (P, T ) plane, the Krichevskii parameter is zero, and the liquid-liquid transition in solution will be isocompositional. That case corresponds to the so-called critical azeotrope [51, 52]. The case demonstrated in Figs. 3.1 and 3.2 corresponds to a negative value of the Krichevskii parameter. The sign of the Krichevskii parameter determines the partition of the solute between the coexisting phases. The negative sign of the Krichevksii parameter means that HDL has a higher concentration of the solute. The existence of phase separation in two-structure thermodynamics, caused by coupling of the order parameter with density and entropy, raises an interesting question concerning the path dependence of the character of spinodal decomposition in such systems. Conventionally, spinodal decomposition in fluids is observed along the critical isochore which, asymptotically close to the critical point, merges with the Widom line. For this path, the final equilibrium state will be the two-phase coexistence between liquid and vapor. However, if a fluid, initially (for example) in the gaseous state, is quenched at constant pressure to the liquid state, the forma- tion of the new equilibrium state may occur by two alternative mechanisms, either nucleation or spinodal decomposition, both without macroscopic phase separation. The same will be true for the liquid-liquid transition in water. This is illustrated in Fig. 3.1. The conventional spinodal decomposition toward macroscopic phase sepa- ration will be observed upon quenching along paths 1 (pure water) and 3 (solution). However, if the final equilibrium state is located in the shaded region between the spinodal (the absolute stability limit of the high-temperature liquid) and the phase 33 c o n c e n tr a ti o n o f s o lu te temperature 0 C C' P2 P1 P3 P1 = Pc' < water Pc P2 = water Pc ("Bird's Beak") P3 > water Pc Figure 3.2: Schematic T -x diagram of a supercooled aqueous solution exhibiting a liquid- liquid transition between HDL and LDL. The red line is the critical line with the critical point of pure water labeled C and the critical point of the solution at a certain concen- tration C’. Solid blue and green lines show the coexistence between two phases, HDL and LDL, respectively. Blue and green dashed lines show the thermodynamic stability limits of HDL and LDL, respectively. 34 transition line, the new state will be formed by nucleation without macroscopic phase separation. If the final state is reached beyond the spinodal, the process will be similar to spinodal decomposition, but without macroscopic phase separation. These events are illustrated in Fig. 3.1 by thermodynamic paths 2 and 4. 3.4 Aqueous Solutions of Sodium Chloride The experimental information on the thermodynamic properties of supercooled aqueous solutions of salts, in particular of NaCl, is very limited. The available data are the isobaric heat capacity measurements of Archer and Carter [6], and the density measurements of Mironenko et al. [43], both at atmospheric pressure. Archer and Carter observed that for small NaCl concentrations, upon lowering the temperature, the heat capacity increases in the supercooled region. As the salt concentration is increased, this anomalous rise in heat capacity moves to lower temperatures and decreases in magnitude, virtually disappearing for salt concentrations greater than 2 mol/kg. Mironenko et al. found that as the concentration of NaCl was increased, the density of the solution increased while the density maximum moved to lower temperatures [43]. About forty years ago, Angell observed the suppression of the heat capacity anomaly in supercooled water upon addition of lithium chloride [56], qualitatively similar to the effects reported by Archer and Carter for sodium chloride [6]. In light of what can be inferred about the movement of the locus of liquid- liquid phase transitions upon addition of NaCl, this behavior of the heat capacity is precisely what thermodynamics predicts if the anomaly in pure supercooled water 35 is indeed associated with a liquid-liquid critical point. Homogeneous ice nucleation in solutions of NaCl is shifted to lower tempera- tures as the concentration of salt increases [57, 58], with the lines of homogeneous nucleation keeping nearly the same shape in the P -T plane as in pure water [58]. The kinks in the melting lines of metastable phases of ice in aqueous solutions of LiCl, observed by Mishima, suggest that the liquid-liquid transition also moves to lower temperatures and pressures as the salt is added, remaining just below the temperature of homogeneous nucleation for any given concentration of solute [59]. Mishima has also observed that the transition in amorphous water between the HDA and LDA phases moves to lower pressures upon addition of LiCl [60]. The location of the critical point in pure water along the liquid-liquid transition is uncertain; according to the analysis of Ref. [10], one can currently only say that the critical pressure is smaller than 30 MPa, and could even be negative. Above the lines of homogeneous ice formation, negative pressures are experimentally accessible and correspond to doubly metastable liquid water, with respect to both the solid and vapor states [61]. Liquid-liquid transitions at these pressures are an intriguing possibility [62–65]. Simulations on the TIP4P [66] and mW [67] models of water suggest that hy- drophilic solutes dissolve more easily in HDL than in LDL, the tetrahedral structure of which they tend to disrupt. This further corroborates the hypothesis that the liquid-liquid transition and Widom line will move to lower pressures (at constant temperature) and to lower temperatures (at constant pressure) as the concentra- tion of salt increases. Corradini and Gallo examined the slope of the liquid-liquid 36 phase transition line and the position of the liquid-liquid critical point in TIP4P water at several concentrations of NaCl [66]. From these results we can estimate the derivatives in Eq. (3.8) as follows: dTc/dx = 770 K, dPc/dTc = −13.1 MPa/K, and (dP/dT )cxc = −3.1 MPa/K, yielding a value of -7700 MPa for the Krichevskii parameter in this water model. Such a large magnitude of the Krichevskii param- eter indicates that the critical anomalies will be greatly suppressed even for small concentrations of NaCl, and its sign indicates that NaCl dissolves better in HDL than in LDL. As can be seen in Fig. 3.3, our equation of state is in qualitative agreement with simulation studies on the TIP4P model of water. Both our equation of state and the simulations of Corradini and Gallo [66] display a large, negative value of the Krischevksii parameter, driven primarily by the movement of the critical point to lower pressures. Corradini and Gallo also find a slight increase in the critical temperature as NaCl is added, and they find a smaller slope for the LLPT. A re- scaling of the transition line obtained for the TIP4P model to match the slope of the transition line in our equation of state suggests an almost vertical critical line in real NaCl solutions (Fig. 3.3). A vertical critical line is adopted in our equation of state and, as shown below, is also supported by further analysis of the heat capacity data. Simulations, experiments on the metastable ices in aqueous LiCl, and experi- ments on the homogeneous nucleation in NaCl are thus in agreement that the locus of liquid-liquid transitions moves rapidly to lower temperatures and pressures upon addition of NaCl, yielding a negative Krichevskii parameter on the order of 103 MPa. 37 0T/Tc T'/Tc P / c R T c P '/  c R T c 0 Figure 3.3: A comparison of our equation of state with the simulation results of Corradini and Gallo for the TIP4P model [66]. Because the systems have very different critical pres- sures, the features are presented in terms of difference from the critical point in variables reduced by the critical parameters of the system, as indicated. The blue solid line and blue dashed line show linear approximations of the LLPT and Widom line, respectively, for our equation of state. The green solid line and green dashed line show a linear approx- imation of the LLPT and Widom line, respectively as reported in Ref. [66] for the TIP4P model. Green circles show the critical points calculated at different mole fractions of NaCl in simulation with the thin green line as a guide to the eye, whle the red line shows the critical line for our equation of state. The critical point of H2O is shown as a red circle. Top: our data and that of Corradini et al. Bottom: data of Corradini et al. re-scaled so that the the LLPT has the same slope as in our equation of state. In order to form a more precise estimate for our model, we take note of Mishima’s evidence that in solutions of LiCl, the liquid-liquid transition remains just below the line of homogeneous ice nucleation as both move to lower temperatures and pressures [59]. With a linear approximation for the curve comprising the locus of liquid-liquid phase transitions and the Widom line, the Krichevskii parameter can be calculated based on the movement of this curve, regardless how the critical point 38 might move along it. Thus, taking the behavior of the line of homogeneous nucle- ation as a proxy for that of the line of liquid-liquid phase transitions, Ref. [58] gives a Krichevskii parameter of K = −2230 MPa and Ref. [57] gives K = −2860 MPa. Within that range, the value that we adopt for the Krichevskii parameter makes only a small difference in the fit of the model to the data, and the slope of the critical line makes little difference provided that the dominant contribution to the Krichevskii parameter comes from the movement of the critical point to lower pres- sures, as suggested by [66]. We find that a vertical critical line and a Krichevskii parameter of K = −2860 MPa provides the most accurate calculations of the heat capacity, and accordingly adopt these parameters. For dilute solutions of simple electrolytes such as NaCl, there is a linear rela- tionship between the concentration of the solute and the depression of the tempera- ture of maximum density of water, a relationship known as Despretz’s law [68–70]. For those salinities at which data exist, our equation of state reproduces this phe- nomenon adequately and matches the experimental data, as shown in Fig. 3.4. 0.00 0.01 0.02 265 270 275 Equation of State Wright (1919) Cawley et al. (2006) T M D ( K ) Mole Fraction of NaCl Figure 3.4: Temperature of maximum density in aqueous solutions of supercooled water. The black line shows the equation of state used in this work, while the red squares and blue circles show the measurements of Refs. [68] and [71] respectively. 39 To calculate the isobaric heat capacity at constant composition we use the thermodynamic relation between this experimentally available property and the “theoretical” (isomorphic) heat capacity CP,δ = T (∂S/∂T )P,δ : CP,x = CP,δ − T (∂x/∂T )2P,δ (∂x/∂δ)P,T , (3.14) yielding CP,x R = T R ( ∂S ∂T ) P,x = aˆ2 χ1 1 + x(φ1Lˆ+ Kˆ)2χ1 +B, (3.15) where the background heat capacity B is approximated as a polynomial function of T and x, and aˆ = (ρcR) −1(dP/dT )c,cxc, Kˆ = K/ρcRTc, χ1 is a strongly divergent susceptibility, and Lˆ = (dPc/dx) /ρcRTc The details of this calculation are provided in the first section of the Appendix to this Chapter. Eq. (3.15) describes the crossover of the heat capacity between two limits. In the limit x → 0 one recovers the expression for the heat capacity of pure water, diverging at the critical point as Cp R = aˆ2χ1. (3.16) As the solution critical point is approached, χ1 →∞, φ→ 0, and the heat capacity approaches a finite value, growing with decreasing concentration: CP R → aˆ 2 xKˆ2 +B. (3.17) 40 The large negative value of the Krichevskii parameter for this system, Kˆ ' −30 is mainly responsible for the significant suppression of the heat capacity anomaly even in dilute solutions of NaCl. The results of fitting Eq. (3.15) to the experimental data of Archer and Carter are shown in Fig. 3.5. The agreement between the theory and experiment is remarkable. 200 220 240 260 280 300 3 4 5 6 C P , x ( k J k g - 1 K - 1 ) T (K) NaCl Mol Fraction 0 .00091 .00178 .00298 .00448 .00657 .00892 .01768 .03474 .05123 .07495 .09750 T M T Figure 3.5: Suppression of the anomaly of the heat capacity in aqueous solutions of sodium chloride. Symbols: experimental data of Archer and Carter [6]. Solid curves: predictions based on two-state thermodynamics. Dashed curve shows the positions of the melting temperatures as given by the IAPWS formulation for saltwater [72]. Dashed-dotted curve shows the temperatures of homogeneous ice formation [58]. 41 3.5 Liquid-Liquid Transition in Water-Glycerol Addition of glycerol lowers the temperature of homogeneous nucleation. Glyc- erol stabilizes the liquid state, because hydrogen bonding between water and glycerol increases the nucleation barrier for ice formation [44]. At mole fractions of glycerol x ≥ 0.135, Murata and Tanaka have reported phase transitions between two liq- uid states in the solution. They observed two alternative types of kinetics in the formation of the low-temperature liquid state: nucleation and spinodal decomposi- tion. They also found that the transition is mainly driven by the local structuring of water rather than of glycerol, suggesting a link to the hypothesized transition between LDL and HDL in pure water. However, Murata and Tanaka have also claimed that the transition between two liquids in supercooled water-glycerol solu- tions is “isocompositional”, i. e., at the transition point, LDL and HDL have the same concentration of glycerol. They also argue that the transition occurs without macroscopic phase separation. Furthermore, they relate these putative features of the phase transition to the non-conserved nature of the order parameter. We suggest an alternative interpretation of the experiments of Murata and Tanaka. As we have shown above, the HDL-LDL transition in aqueous solutions stemming from the transition in pure water cannot be isocompositional, except for the case of a special behavior of the critical line, yielding the Krichevskii parameter to be zero. Moreover, it cannot take place without macroscopic phase separation if there exists a coupling between the order parameter and density and concentration. As suggested by Murata and Tanaka, in HDL glycerol molecules destabilize 42 hydrogen bonding as pressure does in pure water, whereas in LDL cooperative inter- water hydrogen bonding and the resulting enhancement of tetrahedral order promote clustering of glycerol molecules. This suggests that the critical line emanating from the critical point of pure water continues down to negative pressures, while the critical temperature decreases. Therefore, the difference in interaction of glycerol molecules with the two alternative liquid structures practically rules out the possi- bility that the critical point moves tangent to the phase transition line (dP/dT )c,cxc. Thus the coexisting phases will not have the same composition and the Krichevskii parameter will not be zero. Adopting the extrapolation of Murata and Tanaka for atmospheric pressure, the critical point of the solution is found at x ' 0.05 and T ' 225 K, and locating the critical point of pure water at 13 MPa and 227 K as suggested in Ref. [10], we obtain from Eq. (3.8) the Krichevskii parameter to be K ' −600 MPa. The temperature gap for the transition at constant concentration, e. g. x = 0.165 and atmospheric pressure can be evaluated from Eq. (3.13). The difference in the molar volumes ∆V/Vc of the coexisting phases can be estimated as about 0.05 based on the distance between the transition at atmospheric pressure and the critical point at the same concentration of glycerol. Then we find ∆T ' 5 K. In light of this, it is unsurprising that Murata and Tanaka observed the for- mation of LDL alternatively by spinodal decomposition and by nucleation without observing macroscopic phase separation at x = 0.165. As illustrated in Fig. 3.1, the transition should occur through spinodal decomposition if it takes place below the absolute stability limit of the HDL phase, and by nucleation and growth if it takes place between the point where the last drop of HDL vanishes in the metastable state 43 and the absolute stability limit. The slow kinetics in supercooled water-glycerol and the narrow width ∆T of the two-phase region make this scenario worthy of consid- eration. However, available experimental data on the phase behavior of supercooled glycerol-water solutions are still inconclusive. Other interpretations of the results reported by Murata and Tanaka [44], in particular involving partial crystallization, have been proposed by researchers in recent years [45, 73, 74] and the discussion is ongoing in the literature. 3.6 Appendix to Chapter 3 3.6.1 Heat Capacity at Constant Composition Eq. (3.14) in the main text of this Chapter is an expression for the isobaric heat capacity at constant composition. Eq. (3.15) expresses the same property in terms of scaling susceptibilities, so that the reader can more easily see the suppression of the heat-capacity anomaly upon the addition of solute. This Appendix shows the derivation of Eq. (3.15) from Eq. (3.14) and the principles of the scaling theory of the critical point. As discussed in Chapter 2, in this model the independent scaling fields can be expressed in linear approximation as combinations of the temperature T and pressure P , expressed as [10,41] h1 = a1∆P + a2∆T, (3.18) h2 = b1∆T + b2∆P. (3.19) 44 For pure water in this version of the TSEOS, we find: a1 = 1 ρcRTc a2 = − 1 ρcRTc ( dP dT ) c,cxc , (3.20) b1 = 0 b2 = 1 ρcRTc , (3.21) where (dP/dT )c,cxc is the slope of the phase transition line at the critical point. The condition b1 = 0 corresponds to an entropy-driven phase separation [10]. In a two-component mixture there is an additional thermodynamic degree of freedom to consider, and the scaling fields should be generalized to [51,52] h1 = a1∆P + a2∆T + a3∆δ, (3.22) h2 = b1∆T + b2∆P + b3∆δ. (3.23) According to the principle of critical-point universality, the dependent scaling field h3 must depend on the independent scaling fields h1 and h2 in the same way for a mixture as for a pure fluid. Our approximation that the isomorphic Gibbs energy µ1 = G− xδ retains the same form in mixtures as in pure water entails that the coefficients in the scaling fields, a1, a2, b1, and b2 remain unchanged. With respect to an arbitrary point on the critical line, the scaling fields can 45 be expressed to linear order as h1 = a1∆P + a2∆T − ( a1 dPc dδ ∆δ + a2 dTc dδ ∆δ ) , (3.24) h2 = b2∆P − ( b2 dPc dδ ∆δ ) . (3.25) So we can approximate a3 = − ( a1 dPc dδ + a2 dTc dδ ) , (3.26) b3 = −b2dPc dδ . (3.27) The critical-line condition [52] implies that (∂δ/∂x)T,P = RTc/x, therefore a3 = − x ρc(RTc)2 [ dPc dx − ( dP dT ) c,cxc dTc dx ] , (3.28) b3 = − x ρc(RTc)2 ( dPc dx ) . (3.29) Thus a3 is associated with the Krichevskii parameter in accordance with Eq. (3.8); b3 is associated with the parameter K2 = dPc/dx, which plays a secondary role in the behavior of response functions at constant composition. We now evaluate the response functions entering Eq. (3.14). The critical parts of these response functions can be expressed in terms of the scaling susceptibilities, 46 which are defined as follows in the mean-field approximation: χ1 = ( ∂2h3 ∂h21 ) h2 , (3.30) χ2 = ( ∂2h3 ∂h22 ) h1 = Φ21χ1, (3.31) χ12 = ( ∂2h3 ∂h1∂h2 ) = Φ1χ1. (3.32) With b1 = 0, the critical parts of the response functions, denoted with a superscript c, read: 1 RTc ( ∂S ∂T )c P,δ = a22χ1, (3.33) 1 ρcRTc ( ∂x ∂T )c P,δ = a2a3χ1 + a2b3χ12, (3.34) 1 ρcRTc ( ∂x ∂δ )c P,T = a23χ1 + 2a3b3χ12 + b 2 3χ2. (3.35) We approximate the the regular parts of the response functions, denoted by a su- perscript r, as ( ∂x ∂T )r P,δ = 0, (3.36)( ∂x ∂δ )r P,T = x RTc , (3.37)( ∂S ∂T )r P,δ = ( ∂S ∂T )r P,x . (3.38) Then, from Eq. (3.14) we have 47 CP,x R = CrP,x R + aˆ2 χ1 1 + x(Φ1Lˆ+ Kˆ)2χ1 . (3.39) 3.6.2 Width of the Two-Phase Region at Constant Composition This section of the Appendix gives the derivation of Eq. (3.13) in the main text. This is a first-order expression for the difference in temperatures between the appearance of the first drop of LDL and the disappearance of the last drop of HDL as an aqueous solution is cooled through the LLT. A linear approximation of the coexistence surface, P − ( P 0c + dPc dδ ∆δ ) = ( dP dT ) c,cxc ( T − ( T 0c + dTc dδ ∆δ )) , (3.40) gives ( dP dδ ) T,cxc = dPc dδ − ( dP dT ) c,cxc dTc dδ (3.41) ' x RTc [ dPc dx − ( dP dT ) c,cxc dTc dx ] (3.42) ' x RTc K, (3.43) where here and below x = xc, as follows from the critical-line condition (3.7). Thus from Eq. (3.11), ∆x ' −xK∆V RT . (3.44) 48 To a first approximation, the width of the two-phase region can be estimated as ∆T ' ∆x ( dT dx ) P,x , (3.45) where (dT/dx)P,x is the slope of the line of symmetry (along the average concentra- tion x) of the two-phase region in a (T, x) plane. The approximate slope of this line is ( dT dx ) P,x = x RTc ( dT dδ ) P,cxc . (3.46) From Eq. (3.40), ( dT dx ) P,x ' K ( dT dP ) c,cxc . (3.47) Therefore, ∆T ' xK2 ∆V RT ( dT dP ) c,cxc . (3.48) 49 Chapter 4: The LLCP in TIP4P/2005 Author’s Note: This Chapter incorporates large selections of text from Ref. [17]. The simulations were carried out by Rakesh Singh at Princeton University. I am primarily responsible for the thermodynamic analysis presented here. “We” refers to myself and Mikhail Ansimov (my thesis advisor), and Rakesh Singh and Pablo Debenedetti of Princeton University (the four authors of Ref. [17].) 4.1 Introduction Several computer simulation studies directly or indirectly suggest the existence of a metastable LLPT for some molecular models of water [9,20,22,23,75–80] and for tetrahedral network-forming models [64, 81–84]. Recent state-of-the-art free-energy computations convincingly confirm the LLPT for the ST2 model [23, 80, 85] and some coarse-grained water-like network-forming models [83]. Of special relevance is the work of Smallenburg and Sciortino [85]: these authors showed that the LLPT in the ST2 model persists upon making the crystal phase metastable with respect to the liquid by tuning the hydrogen bond angular flexibility. This disproves interpre- tations according to which the LLPT is generically a misinterpreted crystallization transition [86–89]. On the other hand, metastable liquid-liquid separation is not 50 observed in the coarse-grained mW model [21, 86, 90, 91], a re-parametrized version of the Stillinger-Weber (SW) model [92]. Placing phenomenological (equation of state) calculations on firmer theoretical [20] and computational ground [80], under- standing the molecular basis underlying the existence [80] or absence [91] of a LLPT in specific models, and elucidating the model-dependent time and length scales over which a metastable LLPT can be observed [93–95] are current objects of activity and robust discussion. The main focus of this Chapter is the TIP4P/2005 water model [96]. This model reproduces satisfactorily the thermodynamics of liquid water and the com- plex, experimentally observed phase diagram of water in its numerous crystalline phases, and is considered to be one of the most accurate classical molecular models of liquid water. The existence of a LLPT in the TIP4P/2005 model is still a subject of debate. Abascal and Vega [77] reported the existence of a LLPT with critical param- eters T = 193 K, ρ = 1012 kg/m3 and P = 135 MPa from molecular dynamics (MD) simulations in the NPT ensemble. The more recent study of Sumi and Sekino [97] in the NPT ensemble also suggests a LLPT for the TIP4P/2005 model, but the critical parameters (T ≈ 182 K, ρ ≈ 1020 kg/m3 and P = 1580−1620 bar) were found to be significantly different than those proposed by Abascal and Vega [77]. An equation of state based on the concept of the presence of two different local structures, proposed by Russo and Tanaka [18], also suggests the existence of a metastable LLPT for this model. Recently, Yagasaki et al. [79] have carried out MD simulations in the NV T ensemble and have observed a spontaneous low- and high-density liquid-liquid phase separation in three models: ST2, TIP5P, and TIP4P/2005. The critical tempera- 51 ture and density of TIP4P/2005 reported by Yagasaki et al. [79] are in agreement with the predictions of Sumi and Sekino [97]. These authors also observed a clear separation of time scales between LLPT and crystallization. However, studies of Limmer and Chandler [87, 88] and Overduin and Patey [94, 95] found no evidence for two metastable liquid phases around the temperature-pressure range suggested by Abascal and Vega [77]. Overduin and Patey [94] also argued that the simulations of Abascal and Vega [77] are too short to obtain converged results. In recent studies, Limmer and Chandler [88] and Overduin and Patey [95] have also challenged the results of Yagasaki et al. Specifically, for TIP4P/2005 [96] and TIP5P [98], Overduin and Patey [95] show that the spontaneous liquid-liquid phase separation reported by Yagasaki et al. [79] exhibits a strong system-size dependence. For a system size of 4000 molecules, both studies, Ref. [79] and Ref. [95], observe regions of different densi- ties separated by well-defined planar interfaces. However, Overduin and Patey [95] also observed that the density difference between these regions was sharply reduced with increasing system size, and disappeared for a system size of 32000 molecules. These authors further argue that, as the appearance of regions of low density is always accompanied by an excess of local ice-like molecules, the regions of different densities observed by Yagasaki et al. [79] are likely associated with appearance and coarsening of local ice-like structures, rather than with liquid-liquid phase separa- tion. This argument supports the conclusion of Limmer and Chandler [88], who also argued that the density differences observed by Yagasaki et al. [79] are due to ice coarsening. 52 There is another aspect of the physics of supercooled water that is closely related to the discussion of the possibility of a metastable LLPT. This has to do with the physical nature of the thermodynamic anomalies, in particular, the trend toward diverging response functions. Currently, there is broad consensus based on the experimental [13–16, 99, 100] and simulation [16, 18, 101–103] studies, that in supercooled water two competing local structures indeed exist. Could this competi- tion, which is assumed to be responsible for the thermodynamic anomalies, be sharp enough to trigger a metastable LLPT? This is the central question. Definitely, this possibility is strongly model-dependent and could also depend on specific experi- mental/simulation conditions. In this work, we have carried out extensive computer simulations in order to explore the nature of the thermodynamic anomalies and, consequently, the possibil- ity of a metastable LLPT in TIP4P/2005. To describe the computed properties, we have applied two-structure thermodynamics, viewing water as a non-ideal mixture of two interconvertible local structures. The thermodynamic behavior of the model in the one-phase region is fully consistent with the existence of an energy-driven LLPT in this model (at least for the simulated length and time scales). We have compared the behavior of TIP4P/2005 [96] with the mW [90] and ST2 [104] models, and with real water. 53 4.2 Computational Model and Simulation Details We performed molecular dynamics (MD) simulations of 216 water molecules interacting via the TIP4P/2005 pair potential [96] in a cubic box at constant tem- perature and volume (NV T ensemble). We computed the properties of liquid water at approximately 200 state points at densities ranging from 1120 − 960 kg/m3 in steps of 20 kg/m3 and temperatures ranging from 300 K down to 185− 180 K (de- pending on the density of the system) in steps of 5 K. This choice of ensemble was partly motivated by the possibility of observing van der Waals loops in the two- phase region (below LLCP), in case they exist. It turned out that we were not able to relax the system in the region where one would expect to observe van der Waals loops. However, using many state points in the NV T ensemble enables us to follow different isochores throughout the one-phase region, and to extrapolate them into the region where the slow relaxation of the system impedes reliable computation. We may thus distinguish between a system with a LLPT, in which the isochores are projected to cross, and a system with competition between two structures but without a LLPT, in which the isochores do not cross. We have also performed MD simulations in the NPT ensemble with 512 water molecules at 0.1 MPa and temperatures ranging from 300 K to 200 K. All simula- tions were performed with use of GROMACS 4.6.5 molecular dynamics simulation package [105]. In all cases, periodic boundary conditions were applied, and a time step of 2 fs was used. The short-range interactions were truncated at 8.5 A˚ for 216 water molecule system and 9.5 A˚ for 512 water molecule system. Long-range 54 electrostatic terms were computed by particle mesh Ewald with a grid spacing 1.2 A˚. Long-range corrections were applied to the short-range Lennard-Jones interaction for both energy and pressure. Bond constraints were maintained using the LINCS algorithm [106]. To maintain constant temperature we used a Nose´-Hoover thermo- stat [107, 108] with 0.2 ps relaxation time. Constant pressure was maintained by a Parrinello-Rahman barostat [109] with 2 ps relaxation time. Molecular models of water are notorious for extremely slow structural relax- ation in the superooled state. This slow structural relaxation often leads to con- troversy over thermodynamic behavior of supercooled water observed in computer simulation studies [22, 86, 87, 94]. In this work, in order to ensure the relaxation of the system at each state point, we computed and carefully monitored the decay of the self part of the intermediate scattering function (Fs(k, t)) [110] with time t. To ensure the relaxation of the system at each thermodynamic condition investigated in this work, MD trajectories are at least 400 times as long as the structural relaxation time (defined as the time at which Fs(k ∗, t) = 1/e, k∗ is the wavenumber corre- sponding to the first peak of structure factor). Depending on the thermodynamic condition, MD trajectory lengths vary between 20 ns and 15 µs. 4.3 Description of thermodynamic properties of TIP4P/2005 water In Fig. 4.1, we present the results of our NV T simulations (open circles) along with isochores predicted by the TSEOS (solid lines). The densities range from 960 to 1120 kg/m3 in steps of 20 kg/m3 and the temperatures range from 300 K 55 down to 180− 190 K (depending on the density) in steps of 5 K. The error bars of the simulation data points are approximately equal to the size of the circles. The shape of the isochores in the supercooled region strongly suggests the existence of a LLCP for this model as predicted by the TSEOS. In Fig. 4.2, we have compared TSEOS predictions for the densities along isobars with the previously reported data by Sumi and Sekino [97] as well as by Abascal and Vega [77] obtained by NPT simulations. We observe quantitatively good agreement between predictions of the TSEOS and simulation data for T > 200 K. However, for very low temperatures, the densities predicted by the TSEOS deviate significantly from previously reported data [97]. This discrepancy most likely arises due to the approximations used in the current form of the TSEOS. In any case, both the new simulation data and the TSEOS strongly imply the existence of a liquid-liquid critical point near 182 K and 170 MPa, consistent with the recent simulation studies by Yagasaki et al. [79] and by Sumi and Sekino [97]. In Fig. 4.1, and Fig. 4.2 we show the TSEOS prediction for the LLPT and LLCP in the P -T and ρ-T planes, respectively. The TSEOS parameters are reported in Table 4.1. In order to gain deeper insight into the thermodynamic behavior of the TIP4P/2005 model in the supercooled state, in Fig. 4.3(a) and Fig. 4.3(b) we demonstrate the behavior of the isothermal compressibility(κT ) and the corresponding predictions of the TSEOS. The isothermal compressibility, defined as, κT = −(1/V )(∂V/∂P )T = 〈(δV )2〉/kBTV (kB is Boltzmann’s constant, V is the volume), is a measure of the mean-square volume fluctuations 〈(δV )2〉 at constant temperature. The compress- 56 160 180 200 220 240 260 280 300 T (K) −100 0 100 200 300 400 P (M Pa ) 1000 kg m−3 980 kg m−3 960 kg m−3 1120 kg m−3 1100 kg m−3 1080 kg m−3 1060 kg m−3 1040 kg m−3 1020 kg m−3 Figure 4.1: Isochores in the P -T plane for the TIP4P/2005 model. The open circles indicate simulation data, while the solid lines show the same iscohores according to the TSEOS. The LLPT line, the LLCP, and the Widom line are shown as the solid black line, large red circle, and black dashed line, respectively. The thin dotted line is the melting line of TIP4P/2005 as reported in Ref. [77]. 57 160 180 200 220 240 260 280 300 T (K) 950 1000 1050 1100 1150 ρ (k g m −3 ) Figure 4.2: Densities along isobars computed by Sumi and Sekino [97] (open squares), Abascal and Vega [77] (open diamonds), in this work (open circles; 0.1 MPa), and fits by the TSEOS (solid lines). The black dashed line bounds the two phase region as predicted by the TSEOS, and the red circle shows the predicted location of the critical point. Isobars shown, from top to bottom, are 300, 200, 175, 150, 125, 120, 100, 70, 40, and 0.1 MPa. 58 Table 4.1: Parameters for the two-structure equation of state Parameter Value Parameter Value Tc 182 K c20 −5.3481 Pc 170 MPa c12 0.000493 ρc 1017 kg/m 3 c21 0.1094 λ 1.407 c30 1.3293 a 0.171 c22 −0.02129 b −0.100 c31 −0.02446 ω0 0.0717 c40 −0.13173 c01 0.8617 c23 0.003687 c02 −0.003412 c32 0.01229 c11 0.01351 c33 −0.003513 ibility as a function of density along isotherms, shown in Fig. 4.3(a), is computed from the simulation data. In Fig. 4.3(b)) we show the compressibility data along the isobars reported by Abascal and Vega [77] to compare with the predictions of the TSEOS. The results presented in these figures show that the two-structure thermo- dynamics successfully describes the observed anomalous behavior of thermodynamic response functions in the supercooled region. 960 980 1000 1020 1040 1060 1080 1100 ρ (kg m−3) 0.002 0.004 0.006 0.008 0.010 κ T (M Pa −1 ) 185 K 190 K 195 K 200 K 210 K B 180 200 220 240 260 280 300 T (K) 0.0000 0.0005 0.0010 0.0015 0.0020 κ T (M Pa −1 ) 120 MPa 100 MPa 70 MPa 40 MPa 0.1 MPa C Figure 4.3: (a) Isothermal compressibility along isotherms. Symbols are simulation data and the curves are predictions of the TSEOS (this work). (b) Isothermal compressibility along isobars. Symbols are simulation data by Abascal and Vega [77] (open circles) along with our work at 0.1 MPa (open squares). The curves are the predictions by the TSEOS. 59 In Fig. 4.4, we demonstrate the temperature dependence of the heat capac- ity at constant volume, CV , along different isochores. The isochoric heat capacity CV , defined as (∂E/∂T )V = 〈(δE)2〉/kBT 2, is a measure of the total energy (E) fluctuations of the system. We observe excellent agreement between two-state ther- modynamics and computer simulation predictions for higher densities (greater than 1040 kg/m3) and reasonably good agreement (considering the larger uncertainties in- volved in computing energy fluctuations) at lower densities and temperatures. From the figure, it is evident that unlike κT , CV does not show any significant anomaly on supercooling down to 185 K. This is not surprising, as the anomaly of the isochoric heat capacity near the critical point is very weak. It originates solely from fluctu- ation effects associated with the divergence of the correlation length and does not exist in the mean-field approximation. According to scaling theory [37], the weak divergence of CV should only be noticeable in the close vicinity of the critical point (practically, within 1-2 degrees, i. e. at (T − Tc)/Tc < 10−2 [50]). Moreover, in a finite-size system, the correlation length cannot exceed the size of the box, and so the critical anomalies are rounded. Our system contains only 216 molecules, which is too small for weak (fluctuation-induced) anomalies to be observed. A crossover TSEOS that incorporates fluctuation effects upon approaching the critical point has recently been applied for the description of the ST2 model [20], and it was shown that within the accuracy of simulation data for that model, fluctuation effects are negligible. 60 4000 4500 5000 5500 6000 C V (J kg −1 K −1 ) ρ = 960 kg m−3 ρ = 980 kg m−3 ρ = 1000 kg m−3 4000 4500 5000 5500 6000 C V (J kg −1 K −1 ) ρ = 1020 kg m−3 ρ = 1040 kg m−3 ρ = 1060 kg m−3 200 220 240 260 280 300 T (K) 4000 4500 5000 5500 6000 C V (J kg −1 K −1 ) ρ = 1080 kg m−3 200 220 240 260 280 300 T (K) ρ = 1100 kg m−3 200 220 240 260 280 300 T (K) ρ = 1120 kg m−3 Figure 4.4: Temperature dependence of the specific heat capacity at constant volume (CV ) along different isochores. Symbols are simlutation data computed from total energy fluctuations and solid lines show the predictions of the TSEOS. 61 4.4 Water-Like Models versus Real Water Our study, together with three previously published simulation results [77,79, 97], shows that the TIP4P/2005 model in the range of pronounced thermodynamic anomalies behaves similarly to the ST2 model. This is clearly seen from the equally sharp behavior of isobars in the vicinity of the projected critical point as demon- strated in Figs. 4.5(a) and 4.2. Even without computational data obtained for the two-phase region, such van der Waals-like behavior of the isobars suggests the proximity of the critical point. Contrary to the ST2 and TIP4P/2005 models, in the mW model of water the isobars, shown in Fig. 4.5(b), only weakly change with changing pressure and never become steep enough to suggest criticality. Indeed, the presence of a LLPT is model-dependent. While in the mW model the non-ideality in mixing of the two structures never becomes strong enough to cause a metastable LLPT, the thermodynamics of the ST2 and TIP4P/2005 models strongly implies the existence of a metastable LLPT. The thermodynamics of real supercooled water is more ambiguous. Properties of bulk supercooled water in the experimentally accessible region are well described by two-structure thermodynamics (for example, density data and theoretical pre- dictions along isobars are presented in Fig. 4.6(a)). However, the projected phase separation is located so far below the homogeneous ice nucleation limit that the location of a LLPT and even its very existence becomes uncertain [10]. This prob- lem with real water is clearly illustrated by comparison of Fig. 4.1, showing the convergence of isochores in TIP4P/2005 at a point that is interpreted as the LLCP, 62 B7 . E Figure 4.5: Simulated data from the ST2 and mW model (symbols) are compared with the predictions (curves) of two-structure thermodynamics, as adapted for the respective models. (a) Temperature-dependent density (ρ) along different isobars computed for (a) the ST2(II) model. The thick black curve indicates two-phase coexistence (dashed: mean field equation, solid: crossover equation) and black dots represent the critical point. The isobar pressures vary from 100 MPa to 200 MPa in steps of 10 MPa. Figure adapted with permission from Ref. [20], c© 2014, American Institute of Physics. (b) Temperature- dependent density (ρ) along different isobars (from the top: 1013, 810, 507, 228, 159, 57, 0.1 MPa) computed for the mW model. Figure adapted with permission from Ref. [21], c© 2013, American Institute of Physics. 180 200 220 240 260 280 300 T (K) 960 980 1000 1020 1040 1060 1080 1100 ρ (k g m −3 ) 0.1 MPa 20 40 60 80 100 120 140 160 180 TH B    7 .   3 03 D      70' 7+            ȡ NJP  E Figure 4.6: (a) Density of cold and supercooled water as a function of temperature along different isobars (black lines are the predictions of an extended version of the TSEOS [29]). Symbols are experimental data reported in Refs. [111] (crosses), [112] (open red circles), and [113] (filled blue diamonds). TH indicates the homogeneous nucleation line. The data from Ref. [112] have been adjusted by at most 0.3% to correct for small systematic errors, as explained in Ref. [41]. (b) Isochores of cold and supercooled water computed with an extended version of the TSEOS [29]. The dashed curve is the homogeneous nucleation line and the blue curve is the TMD locus. 63 and Fig. 4.6(b) for real water in which such convergence is in principle allowed but far from certain. Obviously, the real-water dilemma cannot be resolved with the experimental data that are currently available. Future studies will need to ei- ther penetrate into “no-man’s land” or bring the critical point into experimentally accessible conditions by adding a solute [42,74,114]. 4.5 Discussion: Does a Metastable LLPT Exist in TIP4P/2005? We have investigated the thermodynamic behavior of the TIP4P/2005 water model in the supercooled region. The convergence of the isochores around a density of about 1020 kg/m3 at about 180 − 185 K suggests the presence of a metastable LLPT in the TIP4P/2005 model. Our results are supported by the data of Sumi and Sekino [97] and consistent with the conclusions of Yagasaki et al. [79]. The substantiation of this viewpoint will require free-energy calculations such as those that have yielded unambiguous evidence [80,85] of a liquid-liquid transition in the ST2 model of water. Because the phenomenon under scrutiny is metastable, the question of how sampling time and system size constrain the possibility of observing a liquid-liquid transition arises in addition to the question of its existence in a free-energy or equation-of-state calculation. However, the most recent extensive study of the TIP4P/2005 and TIP5P mod- els by Overduin and Patey [95], which reported simulations in the projected two- phase region for systems ranging in size from 4000 to 32000, found density differences 64 between the regions of low and high densities to decrease with increasing system size. The difference finally disappeared for a system composed of 32000 molecules. Over- duin and Patey further argued that, as the appearance of regions of low density is always accompanied by small ice-like crystallites, the regions of different densities observed by Yagasaki et al. [79] might be associated with the appearance and coars- ening of local ice-like structures, rather than with liquid-liquid phase separation. This argument is similar to that of Limmer and Chandler [88], who also argued that the density differences observed by Yagasaki et al. [79] are due to ice coarsening, rather than to spontaneous liquid-liquid phase separation. This argument deserves serious consideration. However, we must note that separated liquid states observed in NV T simulations are always metastable with respect to ice formation. Consequently, as Overduin and Patey note [95], the mere presence of ice-like crystallites (6 − 8% for TIP4P/2005 model at the lowest tem- perature studied by Overduin and Patey [95]) having finite lifetime in the system does not provide unambiguous proof for the ice-coarsening hypothesis proposed by Limmer and Chandler [86–88]. Also, the computed fraction of ice-like particles or crystallites is very sensitive to the definition adopted for classifying a water molecule as ice-like. On the contrary, the observed excess local density of ice-like crystallites and strong correlations among them in low-density regions can also be understood without invoking the ice-coarsening hypothesis. Liquid-liquid phase separation leads to spatial heterogeneity in water, and it is to be expected that the ice-like fluctu- ations or crystallites will be more stable in the low-density regions due to a lower surface free-energy cost. In this context, the recent simulations of Smallenburg and 65 Sciortino [85] have clearly demonstrated that the liquid-liquid transition in the ST2 model is not a misinterpreted crystallization transition, as had been claimed [86,87]. Moreover, the fact that phenomena observed in the metastable region depend on the system size and on the duration of observation time is not surprising. This is, in fact, an essential characteristic of metastability. A metastable phase separation cannot exist at all in the thermodynamic limit (infinite size and infinite time). If we denote by τrelax the internal reaxation time in the metastable state, and τout the time it takes for the system to exit the metastable state and form the stable phase (i. e. a characteristic crystallization time in our case), then the metastable state is well defined if τrelax << τout. When this condition is met, thermodynamics can be applied to a metastable state. As Overduin and Patey note, there are several reasons why the metastable LLPT might not be manifested in large enough systems [95]. In particular, as has been emphasized by Binder [115], the divergence of the correlation length at the critical point causes the relaxation time to diverge, an effect known as critical slowing-down. Increasing the system size, on the other hand, decreases the lifetime of metastability, and thus at certain conditions prevents the manifestation of metastable phase separation. The formation of two liquid phases can also be impeded by the unfavorable interfacial energy between them. Consequently, the extent of phase separation not only depends on the choice of initial density of the system but also on the aspect ratio of the simulation box. Due to the large surface energy cost for the formation of well-defined stable interfaces, phase separation is not observed in cubic boxes, even in systems far below the LLCP. In order to observe phase separation one always 66 simulates rectangular boxes (1 : 1 : 4 in case of Yagasaki et al. [79] as well as Overduin and Patey [95] for 4000 molecules) to minimize interfacial free energy cost for formation of the LDL-HDL interface. It is thus very plausible that the observation of two different metastable liquid densities in water-like models, such as TIP4P/2005 and TIP5P, would involve length and time scale constraints that would also influence the pathway to homogeneous ice nucleation. As explained above, attempts to directly observe metastable liquid-liquid sep- aration in NV T simulations are subject to non-trivial limitations. We have used an alternative approach to evaluate the hypothesis of the metastable LLPT in super- cooled water. We have studied a relatively small system of hundreds of molecules and performed a series of simulations (about 200) to obtain reliable information on the thermodynamic surface. Our study does not support one of the scenarios discussed by Overduin and Patey [95] in which “liquid-liquid coexistence is simply not a possibility” for the TIP4P/2005 water model. On the contrary, the clear con- vergence of the isochores around 1020 kg/m3 and the behavior of thermodynamic properties demonstrate the tendency to criticality. Furthermore, the equation of state that is built on the assumption of the existence of LLPT fits the simulation data very well. An alternative hypothesis to the competition between two liquid structures would be to attribute supercooled water anomalies (the sharp increases of the re- sponse functions) to pre-crystallization effects [88,95]. Indeed, the theory of so-called “weak crystallization”, which accounts for translational-order fluctuations, describes the properties of the supercooled mW model as well as two-structure thermodynam- 67 ics does [21]. However, pre-crystallization effects cannot explain the convergence of isochores that is clearly observed in the ST2 and TIP4P/2005 models. There is another puzzling result of Overduin and Patey [95] that requires fur- ther studies. The correlation length characterizing fluctuations of density increases sharply upon supercooling in real water [8, 116]. In Ref. [95], Overduin and Patey examine this correlation length in both TIP5P and TIP4P/2005 and claim that it apparently diverges along the critical isochore in TIP5P, but does not exhibit such an anomaly in TIP4P/2005. We note that in our simulations the isothermal com- pressibility increases by an order of magnitude along the critical isochore, which is a strong effect, especially in view of a relatively small size of the system (about 2 nm). The correlation length of density fluctuations is approximately proportional to the square root of the compressibility. Accordingly, the correlation length should increase by about three times, the effect indeed observed for TIP5P [95]. In conclusion, the results of our study strongly support the presence of a liquid-liquid critical point in the TIP4P/2005 model, and are consistent with the possiblity of a liquid-liquid phase transition for this model. Our study does not answer the questions regarding conditions under which the metastable LLPT can or cannot be observed in the region below the projected critical point. Systematic studies at various simulation conditions are required to further our understanding of this deep and important problem. As far as the one-phase metastable liquid region is concerned, investigation of finite-size effects on the shape of the thermodynamic anomalies would be highly desirable. 68 Chapter 5: Negative Pressures 5.1 Introduction The strong intermolecular forces in water enable it to exist as a liquid not only at temperatures below the freezing point, but also at pressures far below the saturation pressure, and indeed at negative pressures. Water that is both stretched and supercooled is thus doubly metastable—metastable against both the vapor and the crystal. Recent experimental breakthroughs have greatly expanded the exper- imentally accessible region of the phase diagram by probing supercooled water at greatly negative pressures [61]. The failure of extrapolations of current equations of state to account for these new data (discussed in Refs. [61] and [117]) highlights the need for a new equation of state based on theoretical considerations relevant to low pressures. In particular, current equations of state do not include a liquid-vapor spinodal (LVS), and so are ipso facto inaccurate at very low pressures. Moreover, the behavior of stretched wa- ter could shed light on questions surrounding the anomalies in supercooled water at ambient and high pressures, including the relevance of two-structure thermodynam- ics and the hypothesis of the LLCP. Of particular interest are the loci of minima and maxima in various thermodynamic properties. The TMD locus is perhaps the 69 signature anomaly of cold water and is a key to understanding the thermodynamics behind water’s odd behavior. The shape of this line as it continues to lower pres- sures can help distinguish among the various scenarios that have been proposed to explain the anomalies of supercooled water. It is also hoped that other maxima and minima that are predicted by two-structure thermodynamics will emerge into the experimentally accessible region at negative pressures, confirming these predictions and allowing, among other things, a more precise location of the Widom line. Al- ready, a minimum in the speed of sound—predicted by the TSEOS to occur just out of the reach of experiment at ambient pressures—has been observed along n isochore at negative pressure [61]. Furthermore, a series of papers by D’Antonio and Debenedetti [118–120], building on [30], has laid out necessary thermodynamic relationships among the lines of extrema in various thermodynamic properties and between them and the LVS. These relationships clarify the important role that the behavior of LVS plays in the thermodynamics of supercooled water and highlight that it must be included in a complete understanding of the subject. As yet, however, the experimental data at large negative pressures are too sparse to form the basis for a reliable equation of state. Consequently, I have examined a wide range of simulation data for the TIP4P/2005 model, with the goal of constructing an extended version of the TSEOS. With the inclusion of a liquid- vapor spinodal, this new equation of state matches the data in TIP4P/2005 over a very broad range of temperatures and pressures without sacrificing the quality achieved in the fit to the critical region by the Equation of State in Chapter 4. Moreover, patterns characterizing the lines of minima and maxima that one 70 finds in TIP4P/2005—and that this extension of the TSEOS reproduces so well—also appear in other water models, both with and without a LLPT or LLCP. Thus it seems likely that the often anomalous behavior of cold and supercooled water can be explained by two key concepts: the interconversion of alternative hydrogen-bonding structures, and the liquid-vapor spinodal. It is my hope that as more experimental data become available in real water, this most recent extension of the TSEOS can be used to give a comprehensive picture of the thermodynamic behavior of water at low temperatures. 5.2 Experimental Situation There is fairly extensive literature on cavitation and the cavitation limit in stretched water, which is reviewed well in Ref. [117]. Experimental data on the thermodynamic properties, however, is more sparse: Henderson and Speedy have measured the TMD at pressures down to -20.3 MPa (where it is found to occur at 281.45 K) [121, 122], Davitt et al. have measured the density and speed of sound down to -26 MPa at 296.45 K [123]. In recent years, much greater tension has been achieved by cooling inclusions of water in quartz at constant density. Azouzi et al. measured the TMD at a pressure near -120 MPa [124]. Doubly metastable water—both supercooled and stretched—was first observed practically as a nov- elty (only down to -0.02 MPa, and reported on comically) by Hayward [125]. The moderate-tension measurements of Ref. [122] include supercooling to 255.15 K. The experimental phase diagram was greatly expanded by measurements of Pallares et 71 al. along two isochores: they measured density and speed of sound in water down to 258 K and pressures of around -100 MPa. Their density data supported the conclusion that the TMD locus remains monotonic down to -100 MPa. The speed of sound results were more surprising: their measured values of cs were nearly dou- ble the extrapolations of the IAPWS equation of state, and they observed a clear minimum in the speed of sound along the isochore ρ = 933kg/m3. The handful of isochores that have yielded usable data on the thermodynamic properties in the doubly metastable region are both intriguing and informative, but they are far from adequate to the task of constructing an equation of state. For this reason, I have returned to the TIP4P/2005 model. As part of a collaboration between the University of Maryland, Universite´ Claude Bernard Lyon-1 in France, and Universidad Complutense de Madrid, the research group of J. L. F. Abascal and C. Valeriani carried out extensive simulations on this model, including data at deep supercooling and negative pressures. The experimental data on this topic are intriguing to say the least, but as yet they are too sparse to permit the construction of reliable equation of state. In order to proceed with this project while experimental data remain unavailable, J. L. F. Abascal and C. Valeriani’s group at the Universidad Complutense de Madrid carried out simulations of the TIP4P/2005 model over a wide range of pressures. M. A. Gonza´lez, at that time a student at UCM, carried out additional simulations along isotherms in an attempt to locate the LVS of the model. At the time of this submission, these authors, together with F. Caupin, have published data as Ref. [126]. A journal article by them is pending review. I have combined this data 72 with the data that were analyzed in the previous chapter and published in Ref. [17] so that a broad picture of the EOS can be seen. 5.3 Equation of State In order to fit this broad range of data, I found it necessary to extend the TSEOS in two ways. First, by the inclusion of terms proportional to ∆Tˆ 2 and ∆Pˆ 2 in the expression for GBA, so that it reads GBA = λ ( ∆Tˆ + a∆Pˆ + b∆Tˆ∆Pˆ + f∆Tˆ 2 + d∆Pˆ 2 ) . (5.1) Physically, these terms account for the difference in heat capacity and compress- ibility between the two structures, respectively. One would speculate that both the heat capacity and the compressibility of the low-density structure would be smaller than those of the high-density structure, and the results of the fitting support this speculation: all good fits to the data had positive values for both d and f . The second addition to the equation of state, and one which proved crucial for the fitting of data at negative pressures, was the explicit addition of a liquid-vapor spinodal to the equation of state. In order to implement the LVS, I took my cue from Ref. [30]. In that work, Speedy argues that because (∂P/∂ρ)T must vanish at the LVS, P in the vicinity of the LVS can be expanded as a function of ρ in a Taylor series whose first non-vanishing, non-constant term will be second-order in ρ, such that P = Ps +B(ρ− ρs)2 + ..., (5.2) 73 where Ps and ρs are the pressure and density of the liquid at the LVS, and B is identified as B = 1 2 ( ∂2P ∂ρ2 ) . (5.3) The original application of this expansion in Ref. [30] was to locate the LVS in real water. It does not seem to be up to that task: the predicted location of the spinodal depends heavily on the choice of the choice of ρ vs. V as an expansion variable and on the range of data that is used for the fit. However, if the position of the spinodal can be located by other means, such an expansion can be used to implement its effect on the thermodynamic properties, as it is in this work. This work uses V rather than ρ as an expansion variable, so that lim P→Ps ( ∂P ∂V ) T = 0, (5.4) and, asymptotically, P − Ps ∼ (V − Vs)2 (5.5) This work meets those conditions in its implementation of the spinodal by adding term of the form Gσ(T, P ) = A(T )(P − Ps(T ))3/2 (5.6) in the expression for the Gibbs energy of structure A. It must be noted here that the expression for GB, i. e., the Gibbs energy of pure structure B, does not appear 74 explicitly in the equation of state, and must be calculated from the expression GB = GBA +GA. Therefore, since there is no compensating term in the expression for GBA, this form of the spinodal term entails that the effects of the LVS appear equally in the Gibbs energy of each species–or to put it another way, that the effects of the spinodal occur irrespective of the fraction φ. The presence of the spinodal does not directly affect φ nor, as a result, does it affect the position or shape of the LLPT or Widom line. Note also that the contribution of the spinodal to a thermodynamic property X is denoted by Xσ, while the value of that thermodynamic property on the spinodal is denoted by Xs. From this expression we find the contributions to the thermodynamic proper- ties as follows. V σ = 3 2 A(P − Ps)1/2 (5.7)( ∂V ∂P )σ = 3 4 A(P − Ps)−1/2 (5.8) − ( ∂S ∂T )σ = 3 4 A ( dPs dT )2 (P − Ps)−1/2 + 3 ( dA dT )( dPs dT ) (P − Ps)1/2 (5.9) + ( d2A dT 2 ) (P − Ps)3/2( ∂V ∂T )σ = 3 4 A ( dPs dT ) (P − Ps)−1/2 + 3 2 ( dA dT ) (P − Ps)1/2 (5.10) From these expression it can be seen that κT will diverge as (P −Ps)−1/2, and we can identify A = 2 √ 2 3 ( ∂2P ∂V 2 )−1/2 T ;P=Ps , (5.11) and CP and αP will diverge with the same exponent. The exception is an extremum 75 in the spinodal pressure as a function of temperature, where (dPs/dT ) = 0, and therefore CP and αP do not diverge. In fact, as Ref. [30] has shown, αP = 0 at such a point, and therefore if the spinodal pressure goes through a minimum, the LMD terminates at that minimum. For this work, A(T ) takes the form A(T ) = A0 + A1∆Tˆ , (5.12) where A0 and A1 are optimized to fit the data on thermodynamic properties. It is likely that the variation of the thermodynamic properties as the spinodal is approached—in other words, the “range” and “strength” of the pre-spinodal effects—differs between structure A and structure B. However, I did not find that the data were fit any better when such an explicit dependence was built into A, so I continued to use a linear temperature dependence as a phenomenological approach that implicitly incorporates this effect inter alia. As explained in [30], if the LVS retraces to higher pressures, then there must exist a TMD locus terminates terminate at the minimum of the LVS. Moreover, αp must have the same sign as the slope of the LVS in the (P, T ) plane in the neighborhood of the LVS itself. For TIP4P/2005 water, this means that the points of minimum density observed along isobars by Abascal and Valeriani’s group is inconsistent with a retracing spinodal except in the improbable event that there is a completely separate TMD locus of which no sign has yet been seen. Consequently, a monotonic LVS is used to model TIP4P/2005. 76 The position of the LVS, Ps(T ), is given a quadratic dependence on tempera- ture: Ps(T ) = S0 + S1∆T + S2∆T 2. (5.13) The parameters {Sn} are derived from a least-squares fit to the cavitation data of M. A. Gonza´lez. S1 and S2 have exactly the values derived from the least-squares fit. In practice, though, one always observes cavitation in simulations before the LVS can be reached, so S0 is adjusted down by 25 MPa. Thus the spinodal has the same shape in the (T, P ) plane as the cavitation line, but lies at somewhat lower pressures. This shift yields a superior fit to the thermodynamic properties. 5.4 Results and Discussion Fig. 5.1 shows PV T data along isochores along with the corresponding iso- chores plotted from the TSEOS. These are the same data that were used for the previous chapter. However, that fit, which did not include a LVS, was unable to account either for the two lowest-density isochores of 920 and 940 kg/m3 or the two highest of 1140 and 1160 kg/m3. These isochores are now matched very well with- out any sacrifice of quality in the critical region. Density data are also shown along isobars (Fig. 5.2) and isotherms (Fig. 5.3). Data from Ref. [97] were not included in the fitting process, but they are shown in Fig. 5.2 nonetheless for the purpose of comparison. The pre-spinodal effects are most clearly visible in the behavior of the higher-temperature isotherms, and the extended TSEOS accounts well for these isotherms. However, the improvement of the extended TSEOS over previous models 77 160 180 200 220 240 260 280 300 T (K) −200 −100 0 100 200 300 400 500 P (M P a) 1000 kg/m3 980 kg/m3 960 kg/m3 940 kg/m3 920 kg/m3 1160 kg/m3 1140 kg/m3 1120 kg/m3 1100 kg/m3 1080 kg/m3 1060 kg/m3 1040 kg/m3 1020 kg/m3 Figure 5.1: Symbols represent simulation data of Abascal’s group on the minima and maxima of various thermodynamic properties as shown in the legend. Solid lines of the same color represent the predictions of the extended TSEOS. The solid black line, open red circle, and dashed black line are the LLPT, LLCP, and Widom line, respectively. is especially noticeable in the low-pressure, low-temperature region. This region is further from the LVS and its behavior is less directly affected by pre-spinodal effects, but the inclusion of an explicit LVS in the model is necessary in order to fit it. This is probably because previous attempts to model the behavior at higher temperatures and very low pressures ignored the LVS and relied on polynomial “background” 78 180 200 220 240 260 280 300 T (K) 800 850 900 950 1000 1050 1100 1150 1200 ρ (k g/ m 3 ) ) 300 MPa 200 MPa 175 MPa 150 MPa 125 MPa 100 MPa 1 MPa -40 MPa -80 MPa -125 MPa -170 MPa Figure 5.2: Density along isobars, from Ref. [97] (squares), this work (diamonds), and fits by the extended TSEOS (solid lines). terms, which led to over-fitting in that region and poor predictions elsewhere. A more theoretically grounded approach to the higher-temperature region, incorporat- ing an LVS, solves this problem. The effects of the two key features that the extended TSEOS aims to capture– the LVS and the LLCP–can both be seen clearly in the κT data along isotherms: the compressibility goes through a maximum in the vicinity of the Widom line, 79 −300 −200 −100 0 100 200 300800 850 900 950 1000 1050 1100 1150 1200 ρ (k g/ m 3 ) T = 200 K −300 −200 −100 0 100 200 300800 850 900 950 1000 1050 1100 1150 1200 T = 210 K −300 −200 −100 0 100 200 300800 850 900 950 1000 1050 1100 1150 1200 T = 230 K −300 −200 −100 0 100 200 300800 850 900 950 1000 1050 1100 1150 1200 T = 240 K −300 −200 −100 0 100 200 300800 850 900 950 1000 1050 1100 1150 1200 ρ (k g/ m 3 ) T = 247 K −300 −200 −100 0 100 200 300800 850 900 950 1000 1050 1100 1150 1200 T = 260 K −300 −200 −100 0 100 200 300800 850 900 950 1000 1050 1100 1150 1200 T = 270 K −300 −200 −100 0 100 200 300800 850 900 950 1000 1050 1100 1150 1200 T = 280 K −300 −200 −100 0 100 200 300 P (MPa) 800 850 900 950 1000 1050 1100 1150 1200 ρ (k g/ m 3 ) T = 290 K −300 −200 −100 0 100 200 300 P (MPa) 800 850 900 950 1000 1050 1100 1150 1200 T = 300 K −300 −200 −100 0 100 200 300 P (MPa) 800 850 900 950 1000 1050 1100 1150 1200 T = 310 K −300 −200 −100 0 100 200 300 P (MPa) 800 850 900 950 1000 1050 1100 1150 1200 T = 320 K Figure 5.3: Density along isotherms, simulated for this work (symbols) compared with fits by the extended TSEOS (solid lines). decreases and then begins to increase once again as the LVS is approached. This effect is captured beautifully by the extended TSEOS, as is shown in Fig. 5.5. Fig. 5.4 shows κT along isobars. The isobaric heat capacity CP is also strongly affected by the presence of the spinodal at low temperatures, and is matched well by the TSEOS, with the exception of a few data points at very low T and P, as shown in Figs. 5.6 and 5.7. The immediate goal of the extended TSEOS has been to give a complete pic- ture of the anomalous thermodynamics of supercooled TIP4P/2005 water. With that in mind, the most concise demonstration of its success is its striking agreement with the simulation data concerning the lines of thermodynamic maxima and min- 80 180 200 220 240 260 280 300 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 κ T (J /(k g K )) P = -170 MPa 180 200 220 240 260 280 300 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 P = -125 MPa 180 200 220 240 260 280 300 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 P = -80 MPa 180 200 220 240 260 280 300 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 κ T (J /(k g K )) P = -40 MPa 180 200 220 240 260 280 300 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 P = 1 MPa 180 200 220 240 260 280 300 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 P = 100 MPa 180 200 220 240 260 280 300 T (K) 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 κ T (J /(k g K )) P = 150 MPa 180 200 220 240 260 280 300 T (K) 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 P = 200 MPa 180 200 220 240 260 280 300 T (K) 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 P = 300 MPa Figure 5.4: Isothermal compressibility along isobars, simulated for this work (symbols) compared with fits by the extended TSEOS (solid lines). ima, as shown in Fig. 5.8. The implications of these lines of extrema for the shape of the LVS and more generally the underlying causes of water’s anomalies deserve further comment. The stability-limit hypothesis [30] explains the thermodynamic anomalies of water as the result of a “re-entrant spinodal”. In the experimentally accessible region of the phase diagram, the LVS in water goes monotonically to lower pressures as the temperature is decreased. But the negative slope of the TMD 81 −300 −200 −100 0 100 200 3000.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 κ T (J /(k g K )) T = 200 K −300 −200 −100 0 100 200 3000.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 210 K −300 −200 −100 0 100 200 3000.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 230 K −300 −200 −100 0 100 200 3000.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 240 K −300 −200 −100 0 100 200 3000.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 κ T (J /(k g K )) T = 247 K −300 −200 −100 0 100 200 3000.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 260 K −300 −200 −100 0 100 200 3000.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 270 K −300 −200 −100 0 100 200 3000.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 280 K −300 −200 −100 0 100 200 300 P (MPa) 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 κ T (J /(k g K )) T = 290 K −300 −200 −100 0 100 200 300 P (MPa) 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 300 K −300 −200 −100 0 100 200 300 P (MPa) 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 310 K −300 −200 −100 0 100 200 300 P (MPa) 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 T = 320 K Figure 5.5: Isothermal compressibility along isotherms, simulated for this work (symbols) compared with fits by the extended TSEOS (solid lines). −300 −200 −100 0 100 200 300 P (MPa) 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 C p (J /(k g K )) T = 210 K −300 −200 −100 0 100 200 300 P (MPa) 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 T = 230 K −300 −200 −100 0 100 200 300 P (MPa) 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 T = 240 K −300 −200 −100 0 100 200 300 P (MPa) 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 T = 247 K Figure 5.6: Isobaric heat capacity along isotherms, simulated for this work (symbols) compared with fits by the extended TSEOS (solid lines). locus and positive slope of the LVS put them on course to intersect, and Speedy demonstrated [30] that such an intersection could only occur where the pressure of the stability limit of the liquid phase goes through a minimum. Approach to such a stability limit upon deep supercooling would explain the increase in κT and CP , (κT must diverge at the spinodal and CP contains a term proportional to κT ) without 82 180 200 220 240 260 280 300 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 C p (J /(k g K )) P = -170 MPa 180 200 220 240 260 280 300 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 P = -125 MPa 180 200 220 240 260 280 300 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 P = -80 MPa 180 200 220 240 260 280 300 T (K) 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 C p (J /(k g K )) P = -40 MPa 180 200 220 240 260 280 300 T (K) 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 P = 100 MPa 180 200 220 240 260 280 300 T (K) 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 P = 150 MPa Figure 5.7: Isobaric heat capacity along isobars, simulated for this work (symbols) com- pared with fits by the extended TSEOS (solid lines). the necessity of postulating a LLCP or LLPT. Such a hypothesis is not consistent with the data in TIP4P/2005 water, however. If the stability limit of the liquid does retrace to higher pressures upon deep supercooling, then (1) a TMD locus must emanate from the point of its minimum pressure, and (2), in the neighbor- hood of a negatively sloped stability limit, αP must be negative. The TMD locus in TIP4P/2005 water, however, does not terminate at the LVS or any other stability limit. Rather, it bends back to lower temperatures, avoiding the LVS, and joins a locus of minimum density at a point where each locus has zero slope, as can be seen in Fig 5.8. Furthermore, κT is expected to diverge at the stability limit. However, 83 160 180 200 220 240 260 280 300 320 T (K) −300 −200 −100 0 100 200 300 P (M P a) ρmax ρmin Cmaxp Cminp κmaxT κminT Figure 5.8: Symbols represent simulation data of Abascal’s group on the extrema of various thermodynamic properties as shown in the legend. Solid lines of the same color represent the predictions of the extended TSEOS. The solid black line, open red circle, and dashed black line are the LLPT, LLCP, and Widom line, respectively. The red dashed line is the LVS, and the red stars show the points at which cavitation was observed in the simulations of M. Angel. in TIP4P/2005, κT decreases monotonically upon cooling for pressures lower than -102 MPa, and for pressures above this, a finite maximum is observed for pressures up to 120 MPa [77]. A retracing spinodal, therefore, would also entail the existence of an additional TMD locus and an additional locus of minima in the isothermal 84 compressibility, occurring at extremely low temperatures. Therefore, while techni- cally possible, the hypothesis of a retracing spinodal is extraneous to the analysis of the observed anomalies of TIP4P/2005, and cannot explain the anomalies of any liq- uid whose thermodynamic properties show a qualitatively similar pattern of minima and maxima. The success of a version of the TSEOS that incorporates a LVS in matching the lines of extrema in the thermodynamic functions is particularly significant in light of the current state of research into other models. Poole et al. have investigated the lines of thermodynamic minima and maxima in the ST2 model at low temperatures for a wide range of pressures [75]. A figure from their publication is reproduced here as Fig. 5.9. The similarity of these results to those found in TIP4P/2005 is suggestive to say the least, and Fig. b5.10 shows a direct comparison between the two models. In that figure, the temperatures and pressures are shifted according to the critical parameters of the model, and then re-scaled empirically. Upon such a re-scaling, the patterns of loci of minima and maxima are remarkably similar. Furthermore, private communications with the Author show similar patterns in other water models, both with and without a LLPT and LLCP. It seems likely, therefore, that the extended TSEOS can be fit to these models as well, with similar success. A project implementing this is currently being planned. Moreover, if such behavior is indeed ubiquitous in water models, one may spec- ulate that it is also a feature of real water. In this case, it will be possible to create a comprehensive picture of the anomalous thermodynamics of low-temperature wa- ter based on two essential concepts: the interconversion of two alternative liquid 85 structures and the liquid-vapor spinodal. The extended TSEOS can then hopefully be adapted to give a comprehensive equation of state for real water at low temper- atures, especially in the metastable region. However, this project must await the availability of further experimental data at negative pressures. Figure 5.9: Lines of thermodynamic minima and maxima in the ST2 model. This figure is reproduced from Ref. [75], with that reference’s notation intact. Temperatures of max- imum and minimum density are thick and thin black lines, respectively. Temperatures of maximum and minimum isothermal compressibility are thick and thin blue lines, re- spectively. Pressures of maximum and minimum isobaric heat capacity are thick and thin green lines, respectively. Red diamonds indicate the estimated position of the LVS. The open blue circle is the LLCP, and the up and down blue triangles are the spinodals of the LDL and HDL phases, respectively. 86 0.0 0.5 1.0 1.5 T (Rescaled) −1.2 −1.0 −0.8 −0.6 −0.4 −0.2 0.0 0.2 P (R es ca le d) TIP4P/2005 ρmax ρmin Cmaxp Cminp κmaxT κminT Figure 5.10: A comparison of the lines of thermodynamic extrema in the TIP4P/2005 (large symbols) and ST2 (small symbols) models. Temperatures and pressures are shifted and displayed relative to the LLCP in each model, and then rescaled so that the loci of extrema coincide. 87 5.5 Appendix: Fitting Parameters The values of the parameters used in the TSEOS are reported in Tables 5.1 and 5.2. The critical parameters are the same as those used in Chapter 4: Tc = 182 K, Pc = 170 MPa, ρc = 1017 kg/m 3. Table 5.1: Background Parameters c02 -0.00261876 c30 2.18819 c11 0.257249 c13 -0.000994166 c20 -6.30589 c22 -0.00840543 c03 0.000605678 c31 0.0719058 c12 0.0248091 c40 -0.256674 c21 -0.0400033 Table 5.2: Special Parameters λ 1.55607 S0 -5.40845 a 0.154014 S1 0.0305542 b 0.125093 S2 -7.61× 10−5 d 0.00854418 A0 -0.0547873 f 1.14576 A1 -0.0822462 ω0 0.03 The value of c01 is not provided because it is not a separate parameter, but is specified by specifying ρc. It should be noted, however, that the formula used for c01 in previous versions of the TSEOS is no longer valid. Rather, the spinodal term must be taken into account when evaluating c01. This is done as follows. Because properties are reduced by critical temperature and critical volume (or density), is is obvious that Vˆc = 1. (5.14) We can also evaluate the general formula for volume at the critical point, where 88 ∆Tˆ = 0, ∆Pˆ = 0, and xe = 0.5 and find that Vˆc = c01 + 0.5λa+ 0.25ω0 + 1.5A0 √ Pc − S0 ρcRTc , (5.15) so c01 = 1− 0.5λa− 0.25ω0 − 1.5A0 √ Pc − S0 ρcRTc . (5.16) For the values of the parameters used in this work, c01 = 1.09621. 89 Chapter 6: Thermal Conductivity and its Relationship to Thermo- dynamics Author’s Note: This section quotes extensively from two published papers, Refs. [127] and [128]. In the section concerning simulations of the TIP4P/2005 model, the simulations were performed by Fernando Bresme at Imperial College, London (the first author of Ref. [128]). I am primarily responsible for the theoretical analysis. –JWB 6.1 Introduction While most of the phenomenology surrounding two-structure models generally and the liquid-liquid critical point in particular has focused on thermodynamics, there are important implications for dynamics as well. To take one example, the viscosity of water decreases upon compression, which is anomalous. A suggested explanation for this anomaly is that compression forces a greater fraction of water into the HDL state, which has greater mobility than the tetrahedrally ordered LDL state [31]. Furthermore, the dispersion of sound at high frequencies seems likely to reflect viscoelastic behavior associated with a structural relaxation in supercooled water [129]. 90 This Chapter examines the thermal conductivity of supercooled water in light of the two-structure conjecture, and examines the possible effects that critical fluc- tuations might have on thermal conductivity. According to our calculations, any effect of critical fluctuations associated with the virtual liquid-liquid critical point on the thermal conductivity would be too small to be measured at experimentally accessible temperatures. Remarkably, the behavior of thermal conductivity can be fully explained by the anomalies of the thermodynamic properties. The difference between the behavior of the thermal conductivity in the vapor-liquid and in the hypothesized liquid-liquid critical regions is the result of differences in both the dynamic and thermodynamic environments in the respective regions. 6.2 Thermal Diffusivity and Thermal Conductivity The thermal diffusivity a of water has been measured by Taschin et al. down to 256 K [36] and by Benchikh et al. down to 250 K [130], and it decreases steadily with decreasing temperature, as shown in Fig. 6.1. Thermal conductivity λ can be calculated from these data by means of the formula λ = ρCPa, since isobaric specific heat capacity CP and mass density ρ are both known experimentally in the relevant temperature range. The heat capacity changes little in the range for which thermal- diffusivity data are available (Fig. 6.2), so in that temperature range the thermal diffusivity and thermal conductivity are nearly proportional. At atmospheric pres- sure, thermal conductivity decreases with temperature, from the boiling point to the lowest temperatures at which thermal conductivity has been measured [131]. 91 Water’s behavior in this regard is unique among non-metallic liquids of low molec- ular weight, as all other such liquids show an increase in thermal conductivity upon cooling [132]. 220 240 260 280 300 0.00 0.05 0.10 0.15 0.20 0.25 Taschin et al. Benchikh et al. Bridgman Formula with TSEOS Eyring Formula with TSEOS T h e r m a l D i f f u s i v i t y a ( m m 2 s - 1 ) Temperature (K) 220 240 260 280 300 0.1 0.2 0.3 0.4 0.5 Figure 6.1: Thermal-diffusivity data of Benchikh et al. [130] and Taschin et al. [36], com- pared with the thermal diffusivity calculated from Bridgman’s and Eyring’s formulae. We use the TSEOS to evaluate the thermodynamic properties in these formulae [10]. The inset shows simulation results of Kumar and Stanley for the TIP5P model [133]. Benchikh et al. have noted a strong correlation between the thermal con- ductivity of both supercooled and stable liquid water at low temperatures and the thermodynamic speed of sound c (in the limit of zero frequency) [130]. We find excellent agreement between the experimental data and two classical formulations for the thermal conductivity of liquids, as shown in Fig. 6.3. The first is Bridgman’s 92 230 240 250 260 270 280 290 4.0 4.5 5.0 5.5 Archer and Carter Angell et al. TSEOS c p ( k J k g - 1 K - 1 ) Temperature (K) Figure 6.2: Heat-capacity data of Archer and Carter [6] and Angell et al. [134]. The solid curve shows the prediction of the TSEOS [10]. formula [135] as adapted for polyatomic molecules: λ = 2.8kBv −2/3c, (6.1) where kB is Boltzmann’s constant and v is the molecular volume of the liquid, that is, the molecular mass divided by the mass density ρ. The second formulation is due to Eyring, with a correction from Eucken [136]: λ = 2.8kBv −2/3γ−1/2c, (6.2) where γ is the ratio of the isobaric heat capacity to the isochoric heat capacity, 93 CP/CV . The speed of sound can be related to the adiabatic compressibility κS and the isothermal compressibility κT as follows: c = ( 1 ρκS )1/2 = ( CP CV 1 ρκT )1/2 , (6.3) allowing us to re-write Bridgman’s formula in terms of thermodynamic properties as λ = 2.8kBv −2/3 ( 1 ρκS )1/2 , (6.4) and Eyring’s formula as λ = 2.8kBv −2/3 ( 1 ρκT )1/2 , (6.5) To evaluate these expressions, we have used the version of the TSEOS in [10], with the published values for all parameters. Bridgman’s and Eyring’s formulae differ little in the region where experimental data are available. According to two-structure thermodynamics, both the isothermal and the adiabatic compressibilities show maxima associated with the existence of a virtual liquid-liquid critical point. These maxima are located close to the Widom line, defined as the locus of maximum fluctuations of the order parameter, a con- tinuation of the liquid-liquid transition line into the one-phase region [10,27,28,41]. We note that the magnitudes of these maxima are strongly dependent on the value of the critical pressure, the value of which (possibly ranging from 10 to 40 MPa) is difficult to determine. In conclusion, while the observed behavior of thermal conduc- 94 tivity is anomalous inasmuch as it decreases upon cooling, it tracks the anomalous behavior of the adiabatic and isothermal compressibilities as shown in Fig. 6.3. 220 240 260 280 300 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Taschin et al. Benchikh et al. Bridgman Formula with TSEOS Eyring Formula with TSEOS IAPWS ( W m - 1 K - 1 ) Temperature (K) T W 220 240 260 280 300 1.1 1.3 1.5 TIP5P T W Figure 6.3: Thermal conductivity calculated from thermal diffusivity [36, 130] and the TSEOS [10]. The formulations due to Bridgman (4) and Eyring (5) are both shown, as is the IAPWS correlation for thermal conductivity (which is only guaranteed down to 273 K) [131]. The inset shows the simulation data of Kumar and Stanley for the TIP5P model [133], and the dashed curve on the inset is a quadratic fit to the simulation data. TW refers to the Widom temperature in both the main graph and the inset. Recently, Kumar and Stanley [133] reported evidence of a thermal conductiv- ity minimum in a simulation of the TIP5P model of water. The simulation results indicated that the thermal conductivity of this model at first decreases upon cooling as has been observed experimentally in real water; it then reaches a minimum at approximately 255 K, and increases as the temperature is further decreased [133]. The latter increase might seem at first to conflict with experimental data for real 95 water. However, Kumar and Stanley’s simulation locates the Widom temperature for the TIP5P model at roughly 245 K at atmospheric pressure, while thermo- dynamic equations of state place the Widom temperature for real water close to 228 K [10, 27, 28, 41]. Rescaling the temperatures in the simulation results so that the Widom temperature occurs at 228 K places the predicted minimum at 237 K, several degrees below the lowest-temperature measurements of the thermal conduc- tivity, so there is no real contradiction between simulation data and the predicted behavior of thermal conductivity in real water. Moreover, one can expect the mini- mum of thermal conductivity observed in this model to be smoothed by finite-size effects, as is typical for simulations. As can be seen in Fig. 6.3, both the Bridgman formula (4) and the Eyring formula (5) indicate that thermal conductivty should go through a minimum. We shall return to this topic in more detail in the discussion. Next, we address the possible effects of critical fluctuations on the thermal conductivity of supercooled water. It is well documented that the thermal conduc- tivity of water diverges at its vapor-liquid critical point, and the associated anomaly affects the thermal conductivity noticeably throughout the critical region [137]. We investigate the possibility of such a divergence of the thermal conductivity near the liquid-liquid critical point of H2O, and any effects that this might have on the measurable behavior of the thermal conductivity in supercooled water. 96 6.3 Predictions of Mode-Coupling Theory In the vicinity of a critical point, couplings among the various hydrodynamic modes of a system become increasingly significant as fluctuations in the system become long ranged. This leads to anomalous behavior of the thermal conductivity in the critical region, including a divergence of the thermal conductivity at the vapor-liquid critical point. This divergence has been observed in striking agreement with the mode-coupling theory in many molecular fluids near their respective liquid- vapor critical points [137, 138]. (This mode-coupling theory describes dynamics in the vicinity of a critical point and should not be confused with the mode-coupling theory of the glass transition). Due to a mode-coupling contribution to the thermal diffusivity, which arises in molecular fluids from a coupling between the heat mode and the viscous mode, thermal conductivity can be expected to diverge near any critical point at which the isobaric heat capacity diverges more strongly than the correlation length ξ [138]. Such a strong divergence of the isobaric heat capacity is a feature of the TSEOS [10], as well as other related scaling models [27, 28, 41]; this prompted us to investigate the possibility that a critical enhancement to the thermal conductivity might be experimentally observable. In order to carry out our mode-coupling calculations we once again made use of the TSEOS [10] for thermodynamic properties. This formulation is renormalized by critical fluctuations, and asymptotically close to the critical point it is identical to the scaling models referred to above [27, 28, 41]. For the background value of the thermal conductivity we used the IAPWS formulation [131]. IAPWS provides 97 a formulation for the viscosity of water at atmospheric pressure that is valid for temperatures as low as 253.15 K [139, 140]. Extrapolating below that value is a more subtle task. Data taken by Osipov et al. [141] from 238 K to 273 K show a clear super-Arrhenius dependence on temperature (“fragile” behavior in Angell’s nomenclature [142]). Some researchers [143] have found evidence of Arrhenius tem- perature dependence (“strong” behavior) close to the glass transition, and thus of a fragile-to-strong crossover in water; such a crossover has been observed to occur at 228 K in confined water [144]. Other measurements, however, have found super- Arrhenius behavior at the glass transition [145]. Starr et al. have used Adam-Gibbs theory to estimate the viscosity in the experimentally inaccessible region [146], and this extrapolation includes a fragile-to-strong crossover as well. We have fit a super- Arrhenius equation to the experimental data of Osipov et al. (see Fig. 6.4), and in making the fit we have chosen a hypothetical temperature of structural arrest so that our extrapolation agrees well that of Starr et al. in the portion of the unsta- ble region for which we make predictions. In our calculations we use the IAPWS formulation for viscosity at atmospheric pressure [139, 140] for temperatures above 254 K and our extrapolation for temperatures below 254 K. Mode-coupling theory gives a pair of coupled equations for the critical enhance- ments to viscosity and thermal diffusivity, and these equations should in principle be solved by iteration [147]. However, the viscosity anomaly associated with crit- ical fluctuations is very weak, so we work only to one-loop order in the iteration and simply use the background value of viscosity, which is strongly temperature- dependent (see Fig. 6.4). The anomaly of the thermal diffusivity is additive in 98 230 235 240 245 250 255 260 265 270 275 280 0 4 8 12 16 20 24 28 32 V i s c o s i t y ( m P a s ) Temperature (K) 240 260 280 300 320 0.000 0.002 0.004 ( k g 2 m 2 Figure 6.4: Viscosity data taken by Osipov et al. at atmospheric pressure (dots) [141]. The solid curve shows the IAPWS formulation for viscosity at atmospheric pressure [139,140], while the dashed curve is a fit to the data of a super-Arrhenius or Vogel-Fulcher-Tamman law: η = 0.0885 exp [220/(T − 197)], with T in K. The inset shows the product of viscosity and thermal conductivity. The increase in the viscosity completely dominates the decrease in the thermal conductivity, so it is clear that the decrease in thermal conductivity is not a result of the increase in viscosity. nature, meaning that the thermal diffusivity can be split into a background and a critical contribution [137]: a = ab + ∆a, (6.6) and the thermal conductivity can be treated similarly: λ = ρCPa = λb + ∆λ. (6.7) 99 Even to one-loop order the integral for the thermal diffusivity enhancement cannot be solved exactly except in the asymptotic limit ξ →∞. Our approximation strategy for treating the crossover from critical to mean-field behavior is based on the model put forward by Olchowy and Sengers [148, 149]. It yields an expression of the form ∆λ = ρCP∆a = ρCP RDkBT 6piηξ (Ω− Ω0) , (6.8) where RD is a universal amplitude very close to unity (experiments by Burstyn et al. [150] find RD = 1.02± 0.06). In the limit ξ →∞, we have (Ω− Ω0)→ 1. If we adopt RD = 1, this expression tends to the well-known limit of a Stokes-Einstein law for thermal diffusivity in which the correlation length of the critical fluctuations replaces the hydrodynamic radius of Brownian particles: ∆a = kBT 6piηξ . (6.9) Due to the effects of long-time tails on the hydrodynamic modes, mode-coupling effects do not completely vanish far from criticality [148, 151]. These long-time effects are already present in the background, and so the phenomenological term Ω0 is introduced to subtract these effects from our expression so that it represents only the critical enhancement. Further details of the approximation scheme can be found in the appendix to this Chapter. The path along atmospheric pressure is not the critical isobar, and along this path properties may exhibit finite anomalies but they do not diverge. We find 100 that at atmospheric pressure, the critical enhancement to the thermal conductivity would reach its maximum in the vicinity of the Widom temperature and at a value close to ∆λ = 6 × 10−6 Wm−1K−1 (see Fig. 6.5). This effect is certainly too small to be measurable, either in the experimentally accessible regime or in the region of the phase diagram below the homogeneous nucleation line (predictions for this region appear in Fig. 6.5 as a dashed curve). Figure 6.5 shows experimental data as well as the IAWPS prediction for thermal conductivity [131]. Any effect from the critical enhancement is far too small to be visible on such a graph; it is shown in the inset, magnified by a factor of a million. We note further that error bars on thermal-conductivity measurements in water are typically of the order of 10−2 Wm−1K−1, several orders of magnitude larger than any possible effect induced by critical fluctuations at atmospheric pressure. Even at the critical pressure (for which the TSEOS uses 13.1 MPa [10]), and even if one could somehow carry out measurements in “no man’s land”, any critical enhancement would be undetectable, as it would be confined to to small a range of temperatures. 6.4 Discussion While thermal transport near the vapor-liquid critical point is dominated by a Stokes-Einstein law for thermal diffusivity, near the hypothesized liquid-liquid crit- ical point thermal transport will continue to be governed by the thermodynamic properties. We can identify two immediate reasons for this striking difference be- tween the two critical regions. First of all, mode-coupling theory predicts that the 101 220 240 260 280 300 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Eyring Formula with TSEOS Taschin et al. Benchikh et al. ( W m - 1 ) Temperature (K) 220 240 260 280 300 0 2 4 6 8 10 6 Figure 6.5: The main figure shows thermal conductivity in water, both from experimental data [36,130] and from Eyring’s formula (Eq. (6.5)) evaluated with the TSEOS [10]. The inset shows the critical enhancement in the same units, but magnified by a factor of one million (106). The curves change from solid to dashed at the temperature of homogeneous nucleation. The critical effect is several orders of magnitude smaller than the background and will be completely undetectable. critical enhancement to the thermal conductivity will be inversely proportional to the viscosity. The viscosity of supercooled water increases dramatically as the tem- perature decreases: Osipov et al. [141] found an increase of a factor of 10 between 273 K and 238 K in their experiment. Our extrapolation predicts that the viscosity at the Widom temperature will be yet another order of magnitude larger than at 238 K; this means that the viscosity of water at the Widom temperature for the hypothesized second critical point is between two and three orders of magnitude larger than at the liquid-vapor critical point, which greatly suppresses any critical 102 enhancement to the thermal conductivity. The second physical reason why any measurable anomaly will be confined to such a tiny region of the phase diagram in supercooled water is that while tem- perature is (to a good approximation) the thermal field for the vapor-liquid phase transition, for the hypothesized liquid-liquid phase transition in water it very nearly plays the role of the ordering field [27, 28, 41]. Near the critical point, the ordering field is related to the order parameter according to a power law φ ∼ h1/δ1 , with δ ≈ 4.8 for the Ising-model universality class [37, 38]. Thus, asymptotically close to the critical temperature, small variations in the temperature correspond to large variations in the order parameter, and so for practical purposes a small deviation from the Widom temperature moves the system far from criticality. 6.5 Simulations To further investigate the anomalous behavior of the thermal conductivity as a function of temperature T and pressure P and its possible relation to the speed of sound, we have performed molecular dynamics simulations of the TIP4P/2005 water model. [96] The TIP4P/2005 model is currently the most accurate non-polarizable water model available [152], and yields adequate representations of densities and compressibilities of real supercooled water [153], and of the speed of sound in the cold-stable region [154]. Working with a rigid non-polarizable model enables us to cover long time scales, tens of ns, which are essential to obtain convergent results, particularly at temperatures near the glass transition. The existence of the LLCP 103 180 200 220 240 260 280 300 T (K) 920 940 960 980 1000 1020 1040 1060 ρ (k g m −3 ) 0.1 MPa 70 MPa 120 MPa 180 200 220 240 260 280 300 T (K) 0.0000 0.0004 0.0008 0.0012 0.0016 κ T (1 06 Pa −1 ) 180 200 220 240 260 280 300 T (K) 1000 1500 2000 2500 c (m s− 1 ) Figure 6.6: Density, compressibility, and speed of sound of TIP4P/2005 at 0.1, 70, and 120 MPa as a function of temperature. The symbols indicate the simulated values: squares (with error bars) obtained in this work and circles (without error bars) obtained previously by Abascal and Vega. [77] The curves represent values calculated from the TSEOS. in TIP4P/2005 is the subject of Chapter 4 of this Dissertation. The speed-of-sound minimum that is predicted by two-structure thermody- namics should move to lower temperatures as the pressure increases. If our analysis of the thermal conductivity of supercooled water is correct for the TIP4P/2005 model, then simulations of this model along several isobars should yield the follow- ing results: speed-of-sound minima at temperatures that decrease as the pressure is raised, and thermal-conductivity minima occurring at temperatures that track those of the minima in the speed of sound. With this goal in mind, the thermody- namic properties at three different pressures, namely at 0.1, 70, and 120 MPa, were determined by performing equilibrium molecular dynamics simulations in the NPT ensemble. We employed cubic simulation boxes with full periodic boundary condi- tions. A sample consisting of 878 molecules was simulated at each desired pressure and temperature by using the isotropic Parrinello-Rahman barostat [155, 156] and the Nose´-Hoover thermostat. [108,157] The compressibility for the barostat coupling was set to 5×10−4 MPa−1. The time constants for the thermostat and the barostat were set to 0.2 ps and 1 ps, respectively, while the equations of motion were inte- 104 grated with a time step of 2 fs. The molecular interactions were truncated at 1 nm. Long-range corrections for the pressure P and the energy E were included in our computations and the electrostatic interactions were handled with the particle-mesh Ewald method. To obtain convergent results, our simulations covered times from 0.5 to 0.7 µs. The equations of state were obtained from equilibrium simulations performed in parallel with Gromacs 4.5.5. [105] Figure 6.6 shows the simulated val- ues obtained for the density, isothermal compressibility, and the speed of sound at the three pressures as a function of temperature. Speed of sound was calculated from heat capacities and compressibilities, which were obtained from analysis of the fluctuations. Our equilibrium properties supplement and agree with previous com- putations of Abascal and Vega, [77] also shown in Fig. 6.6 at the pressures considered in this work. We note that TIP4P/2005 in the supercooled regime reaches the diffu- sive regime at short times (10−8 s) compared to our sampling time (7×10−7 s); thus we are confident that we obtained equilibrium properties. We represent the simu- lated thermodynamic properties by the same type of TSEOS that was previously used by Holten et al. to describe the experimental thermodynamic properties of real water, [10] as well as the properties of the mW and ST2 models. [20,21] The curves in Fig. 6.6 represent the values calculated from the TSEOS. The TSEOS generally represents the simulated data to within their accuracy, except for some data points at very low temperatures. Our TSEOS implies a critical temperature Tc = 183 K, in good agreement with a recent estimates of Tc = 185 K [79] and Tc = 182 K [97]. We computed the thermal conductivity of TIP4P/2005 at the same thermo- dynamic states for which the thermodynamic properties were obtained. All simu- 105 lations were performed in a microcanonical NV E ensemble with 500 molecules in cubic boxes with full periodic boundary conditions. The electrostatic interactions were computed by using the particle-particle particle-mesh Ewald (PPPM) method with a 1 nm cutoff for the dispersion interactions. A time step of 1 fs was employed for all the thermal-conductivity simulations. The computations were performed with the parallel code LAMMPS. [158] Equilibrated configurations obtained from the NPT simulations were employed as starting points for the microcanonical sim- ulations. The thermal conductivity was computed with the aid of the Green-Kubo (GK) correlation function: [159] λ = V 3kBT 2 ∫ tm 0 dt 〈Jq(t) · Jq(0)〉 , (6.10) Jq = 1 V [∑ i viei + 1 2 ∑ i 6=j (fij · vi)rij ] . (6.11) In these equations V is the sample volume, Jq the heat flux, ei the energy (kinetic + potential) of atom i, vi the velocity of atom i, and fij the force between atoms i and j. The summations in Eq. (6.11) run over all atoms in the system, and include non-bonded and bonded interactions (see ref. [160]). The computation of the heat flux with the electrostatic interactions has been discussed previously, both for the Ewald-summation approach [161] and for the PPPM approach [160]. We have chosen the GK method, since it is more effective for resolving the ther- mal conductivities of thermodynamic states with similar temperatures. The slow dynamics associated with the supercooled states means that the computation of 106 1 10 100 1000 Time (fs) -0.5 0 0.5 1 < J( t) J (0 )> n 1 2 3 4 5 6 7 8 Figure 6.7: Normalized heat-flux correlation functions, 〈J(t) · J(0)〉n = 〈J(t) · J(0)〉 / 〈J(0) · J(0)〉, for different systems investigated in this work at P=70 MPa. The labels 1 to 8 indicate temperatures: 270.8, 250.5, 240.8, 229.9, 220.7, 211.6 and 191.0 K. the thermal conductivity requires significant sampling for very long times. Since the thermal conductivity exhibits a weak dependence on the temperature near a minimum, the simulations at low temperatures needed trajectories of the order of 80 ns. Only with these long time scales were we able to resolve the minima of the thermal conductivity. Averages obtained from short trajectories, e.g., 2 ns, yielded thermal-conductivity values that were too noisy for us to resolve the presence of a minimum. The heat-flux correlation function in the integrand of Eq (6.10) (Fig. 6.7) exhibits enhanced oscillations at low temperatures requiring short time steps of 1 fs in the evaluation of the integral. With the choice t = 5 ps for the upper 107 180 200 220 240 260 280 T (K) 0.80 0.85 0.90 0.95 1.00 λ (W m −1 K −1 ) 0.1 MPa 70 MPa 120 MPa Figure 6.8: Thermal conductivity of TIP4P/2005 at 0.1, 70, and 120 MPa as a function of temperature. limit, good convergence was found for all the integrations; longer correlation times (up to 10 ps) gave no evidence for decay in the correlation functions. The results of our simulations of the thermal conductivity for the three pressures are shown in Fig. 6.8. We see that the thermal conductivity at each pressure does exhibit a mini- mum as a function of temperature. The temperature Tmin, at which the thermal conductivity exhibits a minimum, decreases with increasing pressure. Within com- putational accuracy, the location of this minimum temperature Tmin is correlated with the temperatures of the maximum of the compressibility and with the mini- mum of the sound velocity, either directly or through the Bridgman equation (6.4), 108 0 20 40 60 80 100 120 140 p (MPa) 190 200 210 220 230 240 T (K ) Temperature of minimum in thermal conductivity Temperature of minimum in speed of sound Temperature of minimum in Bridgman formula Temperature of maximum in isothermal compressibility Figure 6.9: Temperature of minimum thermal conductivity compared with the tempera- ture of the maximum of the isothermal compressibility or the minimum of the sound veloc- ity either directly or through the Bridgman equation (6.4). Extrema the thermodynamic properties and the Bridgman equation were evaluated using our TSEOS; compressibility maxima agree with Ref. [77]. as shown in Fig. 6.9. A simple scale transformation produces a universal curve for the thermal con- ductivity of the water model in the supercooled state. This is shown in Fig. 6.10, where we have plotted the thermal conductivity as a function of T −Tmin, while ac- counting for a small linear dependence of the thermal conductivity on the pressure. This indicates that the depth of the minimum changes only slowly with pressure, in contrast to the anomaly of the sound velocity. It also means that unlike the inverse compressibility, which vanishes at the critical point, the thermal conductivity likely 109 −20 −10 0 10 20 30 40 50 60 T −Tmin (K) 0.75 0.80 0.85 0.90 0.95 1.00 λ − d ·P (W m −1 K −1 ) 0.1 MPa 70 MPa 120 MPa Figure 6.10: Thermal conductivity, corrected for a small linear pressure dependence (d = 0.564 mW m−1 K−1 MPa−1), as a function of T − Tmin. remains finite. Therefore, the thermal conductivity in supercooled water is only partially controlled by thermodynamics. As speed of sound decreases, one should expect other mechanisms of heat transfer, such as particle diffusion, to become more significant. [162] While our results convincingly demonstrate that the anomalous behavior of the thermal conductivity is of a thermodynamic origin, the values obtained for its magnitude from the simulations are significantly larger that the experimental ther- mal conductivities of real water. In Fig. 6.11 we show a comparison between the simulated values of the thermal conductivity (Fig. 6.11a) and the experimental thermal conductivity data of real water [36,130] (Fig. 6.11b) at P = 0.1 MPa. The 110 190 200 210 220 230 240 250 260 270 T (K) 0.5 0.6 0.7 0.8 0.9 λ (W m −1 K −1 ) a 220 230 240 250 260 270 280 290 300 T (K) 0.5 0.6 0.7 0.8 0.9 λ (W m −1 K −1 ) b Taschin et al. (2006) Benchikh et al. (1985) Figure 6.11: Thermal conductivity of TIP4P/2005 at 0.1 MPa (a) and thermal conduc- tivity of real water [36, 130] (b) at 0.1 MPa. The curves represent values calculated from the Bridgman equation (6.4) for TIP4P/2005 (a) and for real water (b). The Bridgman formula in each case was evaluated with the TSEOS; in (a) we used our optimization for TIP4P/2005 as elsewhere in this work, while in (b) we used the parameters for real water as presented in Ref. [10]. simulated thermal conductivities of Kumar and Stanley [133] for TIP5P are even larger (∼ 1.2-1.5 W m−1 K−1) than those found by us for TIP4P/2005. The discrep- ancies between simulated and experimental thermal conductivities are much less at higher temperatures. [163] While the Bridgman equation (6.4) yields a good quan- titative representation of the thermal conductivity of real water, and the Bridgman equation for the model yields values close to the experimental thermal conductivity values in real water, the simulated thermal conductivities of TIP4P are much larger than the values estimated from the Bridgman equation for the model. Hence, it appears that in the supercooled state simulations of thermal conductivity suggest additional heat transport that is not present in real water. A study of the origin of this discrepancy is highly desirable. 111 6.6 Conclusion At high temperatures, density fluctuations associated with the vapor-liquid critical point noticeably affect the thermal conductivity over a fairly broad tem- perature range. Our calculations suggest, however, that this will not be the case for the hypothesized liquid-liquid critical point in supercooled water. The effect of density fluctuations on the liquid-liquid critical point is too small to be measurable. The Stokes-Einstein law that describes thermal transport in the vapor-liquid critical region will not be applicable to the thermal diffusivity in the liquid-liquid critical region. On the other hand, the thermal conductivity and thermal diffusivity of su- percooled water are strongly correlated with the anomalies of the thermodynamic properties associated with the existence of a liquid-liquid transition. The minimum of the thermal conductivity, found in simulations of the TIP5P model by Kumar and Stanley [133], and now confirmed for the TIP4P/2005 model, should also ex- ist in real water. This effect is associated with associated with the maximum of compressibility and with the minimum of the speed of sound. The anomalous behavior of the thermal conductivity in water, then is not a distinct anomaly. It is rather the result of classical dependence anomalous thermo- dynamics. Two-structure thermodynamics gives an elegant account of these phe- nomena. 112 6.7 Appendix: Details of the Mode Coupling Calculations Starting from the mode-coupling integral as given in ref. [149]: ∆a = kBT (2pi)3 ∫ qD 0 dk CP (k) CP (0) k−2sin2(θ) η(k)(1 + a(k)ρ/η(k)) . (6.12) The Olchowy-Sengers approximation [151] reasons that close to the critical point the term aρ/η is negligible. As noted, we use the background value of the viscosity in our approximation, so after carrying out the angular integrals the integral to evaluate is as follows: ∆a = kBT 3pi2η ∫ qD 0 dk CP (k) CP (0) . (6.13) In the TSEOS the isobaric heat capacity can be expressed in terms of the relevant scaling variables and a few parameters of the system as follows: CP = −λ(1 + b∆Pˆ )Tˆ (φ+ 1) + CAP + 1 2 λ2(1 + b∆Pˆ )2Tˆ 2χ1, (6.14) where λ and b are parameters from the TSEOS with the values given in the sup- plement to ref. [10]. CAP is a background term representing the dimensionless heat capacity of pure HDL. The strongest divergence in the isobaric heat capacity is in the strong scaling susceptibility χ1, so for our third approximation we separate the term containing χ1 from the rest of the heat capacity, and treat the remaining terms (the sum of which we shall call A) as having no wave-number dependence. We then 113 use the Ornstein-Zernike approximation for the wave-number dependence of χ1: χ1(k) = χ1(0) 1 + k2ξ2 , (6.15) where ξ is the correlation length characterizing the fluctuations of the order param- eter. Henceforth CP without any explicit wave-number dependence will refer to the hydrodynamic value, CP (k → 0). With that notation, the integral that we must evaluate for the thermal diffusivity takes the form ∆a = kBT 3pi2η ∫ qD 0 dk [ CP − A CP 1 1 + ξ2k2 + A CP ] , (6.16) which yields the modified Stokes-Einstein law ∆a = kBT 6piηξ 2 pi [ CP − A CP arctan(qDξ) + A CP qDξ ] . (6.17) For convenience we define Ω = 2 pi [ CP − A CP arctan(qDξ) + A CP qDξ ] . (6.18) This result is still not entirely satisfactory because in the limit of vanishing correla- tion length, that is, far from criticality, it does not vanish. In fact: lim ξ→0 kBT 6piηξ 2 pi [ CP − A CP arctan(qDξ) + A CP qDξ ] = kBTqD 3pi2η . (6.19) 114 The physical reason for this as summarized by Olchowy and Sengers in refs. [148,151] is that mode-coupling is also responsible for the “long-time-tail effects on transport properties” [148]. These effects are not critical effects and will be observed in the background, so if we want to find the effects due to critical fluctuations we should subtract off this remnant. For this reason we subtract the following term from Ω: Ω0(qDξ) = 2 pi { 1− exp [ − qDξ 1 + (1− A/CP ) (qξ)4 ]} . (6.20) This phenomenological expression has the following limiting behavior: lim x→0 Ω(x) = 2 pi , (6.21) lim x→∞ Ω(x) = 0, (6.22) so that we have lim ξ→0 ∆a = 0. (6.23) The complete expression, using the above definitions, is ∆a = kBT 6piηξ (Ω− Ω0) . (6.24) Because we are interested only in the critical enhancement to the thermal conductivity, for the correlation length we should use only the critical enhancement to the correlation length. We estimate the “background” correlation length ξb by observing the correlation length at a temperature Tref far from any critical point 115 and assuming that the background correlation length in the system is proportional to the space between molecules. So we have ξb(T ) = ξ(Tref) v(Tref)1/3 v(T )1/3. (6.25) For the correlation length in our calculations we use ξc = ξ − ξb, (6.26) where ξ is the correlation length predicted by the TSEOS [10]. For the wave-number cutoff qD, we used a correlation identified by Perkins et al. between the wave-number cutoff and the amplitude ξ0 of the correlation length anomaly [149]: q−1D = 3.683ξ0 − 1.336× 10−10 m. (6.27) 116 Bibliography [1] F. Franks. The Physics and Physical Chemistry of Water. Plenum Press, New York, 1972. [2] P. G. Debenedetti. Supercooled and glassy water. J. Phys.: Condens. Matter, 15:R1669–R1726, 2003. [3] K. Sessen, K. N. Liou, S. Kinne, and M. Griffin. Highly supercooled cirrus cloud water: Confirmation and climatic implications. Science, 227:411–413, 1985. [4] R.W. Bouchard, B.E. Schuetz, L.C. Ferrington, and Kells S.A. Cold hardiness in the adults of two winter stonefly species: Allopcapnia granulata (claassen, 1924) and a. pygmaea (burmeister, 1839) (plecoptera: Capniidae). [5] C. A. Angell, J. Shuppert, and J. C. Tucker. Anomalous properties of su- percooled water. Heat capacity, expansivity, and proton magnetic resonance chemical shift from 0 to −38 ◦C. J. Phys. Chem., 77:3092–3099, 1973. [6] D. G. Archer and R. W. Carter. Thermodynamic properties of the NaCl + H2O system 4. heat capacities of H2O and NaCl(aq) in cold-stable and supercooled states. J. Phys. Chem. B, 104(35):8563–8584, 2000. [7] R. J. Speedy and C. A. Angell. Isothermal compressibility of supercooled water and evidence for a thermodynamic singularity at -45◦C. J. Chem. Phys., 65(3):851–858, 1976. [8] C. Huang, T. M. Weiss, D. Nordlund, K. T. Wikfeldt, L. G. M. Pettersson, and A. Nilsson. Increasing correlation length in bulk supercooled H2O, D2O, and NaCl solution determined from small angle X-ray scattering. J. Chem. Phys., 133:134504, 2010. [9] Peter H. Poole, Francesco Sciortino, Ulrich Essmann, and H. Eugene Stanley. Phase behavior of metastable water. Nature, 360:324–328, November 1992. 117 [10] V. Holten and M. A. Anisimov. Entropy-driven liquid-liquid separation in supercooled water. Sci. Rep., 2:713, 2012. [11] H. Whiting. A New Theory of Cohesion Applied to the Thermodynamics of Liquids and Solids. PhD thesis, Harvard, 1884. [12] W. K. Roentgen. Ueber die constitution des flussigen wassers. Ann. Phys. U. Chim. (Wied), 45:91–97, 1892. [13] A. Taschin, P. Bartolini, R. Eramo, R. Righini, and R. Torre. Evidence of two distinct local structures of water from ambient to supercooled conditions. Nat. Commun., 4:2401, 2013. [14] T. Tokushima, Y. Harada, O. Takahashi, Y. Senba, H. Ohashi, L. G. M. Pet- tersson, A. Nilsson, and Shin. S. High resolution x-ray emission spectroscopy of liquid water: The observation. Chem. Phys. Lett., 460:387–400, 2008. [15] C. Huang, K. T. Wikfeldt, T. Tokushima, D. Nordlund, Y. Haradac, U. Bergmann, M. Niebuhr, T. M. Weiss, Y. Horikawa, M. Leetmaa, M. P. Ljungberg, O. Takahashi, A. Lenz, L. Ojama¨e, A. P. Lyubartsev, S. Shin, L. G. M. Pettersson, and A. Nilsson. The inhomogeneous structure of water at ambient conditions. PNAS, 106:15214–15218, 2009. [16] A. Nilsson and L. G. M. Pettersson. The structural origin of anomalous prop- erties of liquid water. Nat. Comm., 6:8998, 2015. [17] R. S. Singh, J. W. Biddle, P. G. Debenedetti, and M. A. Anisimov. Two- state thermodynamics and the possibility of a liquid-liquid phase transition in supercooled TIP4P/2005 water. J. Chem. Phys., 144:144504, 2016. [18] J. Russo and H. Tanaka. Understanding water’s anomalies with locally favoured structures. Nat. Comm., 5:3556, 2014. [19] G. Bullock and V. Molinero. Low density water is the mother of ice: on the relation between mesostructure, thermodynamics and ice crystallization in solution. Faraday Discussions, 167:371–388, 2013. [20] V. Holten, J. Palmer, P. H. Poole, P. G. Debenedetti, and M. A. Anisimov. Two-state thermodynamics of the ST2 model for supercooled water. J. Chem. Phys., 140:104502, 2014. [21] V. Holten, D. T. Limmer, V. Molinero, and M. A. Anisimov. Nature of the anomalies in the supercooled liquid state of the mW model of water. J. Chem. Phys., 138:174501, 2013. [22] Y. Liu, J. C. Palmer, A. Z. Panagiotopoulos, and P. G. Debenedetti. Liquid- liquid transition in ST2 water. J. Chem. Phys., 137:214505–1–10, 2012. 118 [23] P. H. Poole, I. Bowles, R. K. Saika-Voivod, and F. Sciortino. Free energy surface of ST2 water near the liquid-liquid phase transition. J. Chem Phys, 138:034505–1–7, 2013. [24] O. Mishima, L. D. Calvert, and E. Whalley. An apparently first order tran- sition between two amorphous phases of ice induced by pressure. Nature, 314:76–78, 1985. [25] O. Mishima and H. E. Stanley. Decompression-induced melting of ice IV and the liquid-liquid transition in water. Nature, 392:164–168, 1998. [26] O. Mishima. Polyamorphism in water. Proc. Jpn. Acad. Ser. B, 86:165, 2010. [27] D. A. Fuentevilla and M. A. Anisimov. Scaled equation of state for supercooled water near the liquid-liquid critical point. Phys. Rev. Lett., 97:195702, 2006. [28] C. E. Bertrand and M. A. Anisimov. Peculiar thermodynamics in supercooled water. J. Phys. Chem. B, 115:14099–14111, 2011. [29] V. Holten, J. V. Sengers, and Anisimov M. A. Equation of state for supercooled water at pressures up to 400 MPa. J. Phys. Chem. Data, 43:043101, 2014. [30] R. Speedy. Stability-limit conjecture. an interpretation of the properties of water. J. Phys. Chem., 86:982–991, 1982. [31] P. G. Debenedetti and H. E. Stanley. Supercooled and glassy water. Phys. Today, 56:40–46, 2003. [32] J. S. Rowlinson and F. L. Swinton. Liquids and Liquid Mixtures. Butterworth Scientific, third edition, 1982. [33] S. Magazu, G. Maisano, D. Majolino, F. Mallamace, P. Migliadro, F. Aliotta, and C. Vasi. Relaxation process in deeply supercooled water by mandelstam- brillouin scattering. J. Phys. Chem., 93:942–947, 1989. [34] A. Cunsolo and M. Nardone. Velocity dispersion and viscous relaxation in supercooled water. J. Chem. Phys., 105:3911–3917, 1996. [35] G. Monaco, A. Cunsolo, G. Ruocco, and F. Sette. Viscoelastic behavior of water in the tetrahertz-frequency range: An inelastic X-ray scattering study. Phys. Rev. E, 60:5505–5521, 1999. [36] A. Taschin, P. Bartolini, R. Eramo, and R. Torre. Supercooled water relax- ation dynamics probed with heterodyne transient grating experiments. Phys. Rev. E., 74:031502, 2006. [37] M. E. Fisher. In F. J. W. Hahne, editor, Critical Phenomena, Lecture Notes in Physics, pages 1–139. Springer, 1983. 119 [38] L. D. Landau and E. M. Lifshitz. Statistical Physics. Elsevier, Ltd., third edition, 1980. [39] Hassan Behnejad, Jan V. Sengers, and Mikhail A. Anisimov. Applied Thermo- dynamics of Fluids, chapter Thermodynamic Behavior of Fluids Near Critical Points. Royal Society of Chemistry, London, 2010. [40] M. E. Fisher and G. Orkoulas. The Yang-Yang anomaly in fluid criticality: Experiment and scaling theory. Phys. Rev. Lett., 85. [41] V. Holten, C. E. Bertrand, J. V. Sengers, and M. A. Anisimov. Thermody- namics of supercooled water. J. Chem. Phys., 136:094507, 2012. [42] J. W. Biddle, V. Holten, and M. A. Anisimov. Behavior of supercooled aqueous solutions stemming from hidden liquid-liquid transition in water. J. Chem. Phys., 141:074504, 2014. [43] M. V. Mironenko, G. E. Boitnott, S. A. Grant, and R. S. Sletten. Experimental determinatin of the volumetric properties of NaCl solutions to 253 K. J. Phys. Chem. B, 105:9909–9912, 2001. [44] K. Murata and H. Tanaka. Liquid-liquid transition without macroscopic phase separation in a water-glycerol mixture. Nature Mater., 11:436–443, 2011. [45] Y. Suzuki and O. Mishima. Experimentally proven liquid-liquid critical point of dilute glycerol-water solution at 150 k. J. Chem. Phys., 141:094505, 2014. [46] J. Wang, C. A. Cerdeirin˜a, M. A. Anisimov, and J. V. Sengers. Principle of isomorphism and complete scaling for binary fluids. Phys. Rev. E, 77:031127, 2008. [47] R. B. Griffiths and J. C. Wheeler. Critical points in multicomponent systems. Phys. Rev. A., 2:1047–1064, 1970. [48] W. F. Saam. Thermodynamics of binary systems near the liquid-gas critical point. Phys. Rev. A., 2:1461–1466, 1970. [49] M. A. Anisimov, A. V. Voronel, and E. E. Gorodetskii. Isomorphism of critical phenomena. Sov. Phys. JETP., 33:605–612, 1971. [50] M. A. Anisimov. Critical Phenomena in Liquids and Liquid Crystals. Gordon and Breach Science Publishers, Philadelphia, 1991. [51] M. A. Anisimov, E. E. Gorodetskii, V. D. Kulikov, and J. V. Sengers. Crossover between vapor-liquid and consolute critical phenomena. Phys. Rev. E, 51(2):1199–1215, 1995. 120 [52] M. A. Anisimov, E. E. Gorodetskii, V. D. Kulikov, A. A. Povodyrev, and J. V. Sengers. A general isomorphism approach to thermodynamic and transport properties of binary fluid mixtures near critical points. Physica A, 220:227– 324, 1995. [53] C. A. Angell, P. H. Poole, and M. Hemmati. A new interpretation of liquid- liquid unmxing in classical alkali silicate glasses. In B. Samuneva and Y. Dim- itriev, editors, 12th Conference on Glass and Ceramics. Science Invest Sophia, 1997. [54] J. M. H. Levelt Sengers. Solubility near the solvent’s critical point. J. Super- critical Fluids, 4:215, 1991. [55] A. A. Povodyrev, M. A. Ansimov, J. V. Sengers, and J. M. H. Levelt Sengers. Vapor-liquid equilibria, scaling, and crossover in aqueous solutions of sodium chloride near the critical line. Physica A, 244:298–328, 1997. [56] C. A. Angell. Unpublished, Personal Communication. [57] Anil Kumar. Homogeneous nucleation temperatures in aqueous mixed salt solutions up to elevated pressures. In Proceedings of the International Con- ference on the Properties of Water and Steam XV. Berlin, 2008. [58] H. Kanno and C. A. Angell. Homogeneous nucleation and glass formation in aqueous alkali halide solutions at high pressure. J. Phys. Chem., 81:2639–2643, 1977. [59] O. Mishima. Melting of the precipitated ice IV in LiCl aqueous solution and polyamorphism of water. J. Phys. Chem. B, 115:14064–14067, 2011. [60] O. Mishima. Phase separation in dilute LiCl-H2O solution related to the polyamorphism of liquid water. J. Chem. Phys., 126:244507, 2007. [61] G. Pallares, M. E. M. Azouzi, M. A. Gonzalez, J. L. Aragones, J. L. Abascal, C. Valeriani, and F. Caupin. Proc. Natl. Acad. Science U.S.A., 111:7936–7941, 2014. [62] H. Tanaka. Phase behaviors of supercooled water: Reconciling a critical point of amorphous ices with spinodal instability. J. Chem. Phys., 105:5099–5111, 1996. [63] I. Brovchenko, A. Geiger, and A. Oleinikova. Liquid-liquid phase transitions in supercooled water studied by computer simulations of various water molecules. J. Chem. Phys., 112:044515, 2005. [64] K. Stokely, M. G. Mazza, H. E. Stanley, and G. Franzese. Effect of hydrogen bond cooperativity on the behavior of water. Proc. Nat. Acad. Sci., 107:1301– 1306, 2010. 121 [65] S. L. Meadley and C. A. Angell. Water and its relatives: the stable, super- cooled and particularly the stretched, regimes. Proceedings of the International School of Physics “Enrico Fermi” Course CLXXXVII, pages 18–43, 2014. [66] D. Corradini and P. Gallo. Liquid-liquid coexistence in NaCl aqueous so- lutions: A simulation study of concentration effects. J. Phys. Chem. B, 115:14161–14166, 2011. [67] L. Le and V. Molinero. Nanophase segregation in supercooled aqueous solu- tions and their glasses driven by the polyamorphism of water. J. Phys. Chem. A, 115:5900, 2011. [68] R. Wright. The effect of some electrolytes on the temperature of maximum density of water. J. Chem. Soc., 115:119–126, 1919. [69] M. Despretz. Recherches sur le maximum de densit e´ des dissolutions aqueuses. Annales de Chimie et de Physique, 70:296–310, 1839. [70] M. Despretz. Troisie`me me´moire sur le maximum de densite´ des liquides. Annales de Chimie et de Physique, 73:296–310, 1839. [71] M. F. Cawley, D. McGlynn, and P. A. Mooney. Measurements of the tem- perature of the density maximum of water solutions using a convective flow technique. Int. J. Heat Mass Transfer, 49:1763–1772, 2006. [72] Release on the IAPWS Formulation 2008 for the Thermodynamic Properties of Seawater, 2008. http://www.iapws.org/relguide/Seawater.html. [73] L.-S. Zhao, Z.-X. Cao, and Q. Wang. Glass transition of aqueous solutions involving annealing-induced ice re-crystallization resolves liquid-liquid transi- tion puzzle of water. Sci. Rep., 5:15714, 2015. [74] Zuofeng Zhao and Austen Angell. Apparently first-order liquid-liquid transi- tion with pre-transition density anomaly, in water-rich ideal solutions. Angen- wandte Chemie, 55:1–5, 2016. [75] P. H. Poole, I. Saika-Voivod, and F. Sciortino. J. Phys: Cond. Mat., 17:L431, 2005. [76] Y. Liu, A. Z. Panagiotopoulos, and P. G. Debenedetti. Low-temperature fluid phase behavior of ST2 water. J. Chem. Phys., 131:104508, 2009. [77] J. L. F. Abascal and Carlos Vega. Widom line and the liquid-liquid critical point for the TIP4P/2005 water model. J. Chem. Phys., 133:234502, 2010. [78] T. A. Kesselring, E. Lasearis, G. Franzese, S. V. Buldyrev, H. J. Herrman, and H. E. Stanley. Finite-size scaling investigation of the liquid-liquid critical point in ST2 water and its stability with respect to crystallization. J. Chem. Phys., 138:244506, 2013. 122 [79] T. Yagasaki, M. Matsumoto, and H. Tanaka. Spontaneous liquid-liquid phase separation of water. Phys. Rev. E, 89:020301(R), 2014. [80] J. C. Palmer, F. Martelli, Y. Liu, R. Car, A. Z. Panagiotopoulos, and P. G. Debenedetti. Nature, 510:385, 2014. [81] S. Sastry and C. A. Angell. Nat. Mater., 2:739, 2003. [82] V.V. Vasisht, S. Saw, and S. Sastry. Nat. Phys., 7:549, 2011. [83] F. Smallenburg, L. Filion, and F. Sciortino. Nat. Phys., 10:653, 2014. [84] Francis W. Starr and Francesco Sciortino. “. [85] F. Smallenburg and F. Sciortino. Tuning the liquid-liquid transition by mod- ulating the hydrogen-bond angular flexibility in a model for water. Phys. Rev. Lett., 115:015701, 2015. [86] D. T. Limmer and D. Chandler. The putative liquid-liquid transition is a liquid-solid transition in atomistic models of water. J. Chem. Phys., 135:132503, 2011. [87] D. T. Limmer and D. Chandler. The putative liquid-liquid transition is a liquid-solid transition in atomistic models of water. ii. J. Chem. Phys., 138:2214505, 2013. [88] D. T. Limmer and D. Chandler. Comment on “spontaneous liquid-liquid phase separation of water”. Phys. Rev. E, 91:016301, 2015. [89] N. J. English, P. G. Kusalik, and J. S. Tse. Journal of Chemical Physics, 139:084508, 2013. [90] V. Molinero and E. B. Moore. Water modeled as an intermediate element between carbon and silicon†. J. Phys. Chem. B, 113:4008, 2009. [91] E. B. Moore and V. Molinero. Nature, 479:506, 2011. [92] F. H. Stillinger and T. A. Weber. Phys. Rev. B, 31:5262, 1985. [93] S. D. Overduin and G. N. Patey. Understanding the structure factor and isothermal compressibility of ambient water in terms of local structural envi- ronments. J. Phys. Chem. B, 116:12014, 2012. [94] S. D. Overduin and G. N. Patey. An analysis of fluctuations in supercooled water. J. Chem. Phys., 138:184502, 2013. [95] S. D. Overduin and G. N. Patey. Fluctuations and local ice structure in model supercooled water. J. Chem. Phys., 143:094504, 2015. [96] J. L. F. Abascal and C. Vega. A general purpose model for the condensed phases of water. J. Chem. Phys., 123:234505, 2005. 123 [97] T. Sumi and S. Hideo. Effects of hydrophobic hydration on polymer chains immersed in supercooled water. RSC Adv., 3:12743–12750, 2013. [98] M. W. Mahoney and W. L. Jorgensen. A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J. Chem. Phys., 112(20):8910–8922, 2000. [99] A. K. Soper and M. A. Ricci. Structures of high-density and low-density water. Phys. Rev. Lett., 84:2881, 2000. [100] L. G. M. Pettersson and A. Nilsson. The structure of water; from ambient to deeply supercooled. J. Non-Cryst. Solids, 407:399, 2015. [101] K. T. Wikfeldt, A. Nilsson, and L. G. M. Pettersson. Spatially inhomogeneous bimodal inherent structure of simulated liquid water. Phys. Chem. Chem. Phys., 13:19918–19924, 2011. [102] E. B. Moore and V. Molinero. Growing correlation length in supercooled water. J. Chem. Phys., 130:244505, 2009. [103] M. J. Cuthbertson and P. H. Poole. Mixturelike behavior near a liquid- liquid phase transition in simulations of supercooled water. Phys. Rev. Lett., 106:115706, 2011. [104] F. H. Stillinger and A. Rahman. Improved simulation of liquid water by molecular dynamics. J. Chem. Phys., 60:1545, 1974. [105] D. van der Spoel, E. Lindahl, B. Hess, A. R. van Buuren, E. Apol, P. J. Meulen- hoff, D. P. Tieleman, A. L. T. M. Sijbers, K. A. Feenstra, R. van Drunen, and H. J. C. Berendsen. Gromacs User Manual version 4.5.6, www.gromacs.org, 2010. [106] B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije. Lincs: A linear constraint solver for molecular simulations. J. Comput. Chem., 18:1463–1472, 1997. [107] S. Nose´. A molecular dynamics method for simulations in the canonical en- semble. Mol. Phys., 52:255, 1984. [108] W. G. Hoover. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A, 31(3):1695–1697, 1985. [109] M. Parrinello and A. Rahman. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys., 52:7182, 1981. [110] J. P. Hansen and I. R. McDonald. Theory of Simple Liquids. Academic Press, 2006. 124 [111] D. E. Hare and C. M. Sorensen. The density of supercooled water. II. bulk samples cooled to the homogeneous nucleation limit. J. Chem. Phys., 87:4840, 1987. [112] O. Mishima. Volume of supercooled water under pressure and the liquid-liquid critical point. J. Chem. Phys., 133:144503, 2010. [113] T. Sotani, J. Arabas, H. Kubota, and M. Kijima. High Temp. High Pressures, 332:433, 2000. [114] D. Corradini, M. Rovere, and P. Gallo. A route to explain water anomalies from results on an aqueous solution of salt. J. Chem. Phys, 132:134508, 2010. [115] K. Binder. Proc. Nat. Acad. Sci., 111:9374–9375, 2014. [116] K. T. Wikfeldt, C. Huang, A. Nilsson, and L. G. M. Petterson. Enhanced small-angle scattering connected to widom line in simulations of supercooled water. J. Chem. Phys., 134:214506, 2011. [117] F. Caupin. Escaping the no man’s land: Recent experiments on metastable liquid water. J. Non-Cryst. Solids, 407:441–448, 2015. [118] P. G. Debenedetti and M. C. D’Antonio. On the nature of the tensile instability in metastable liquids and its relationship to density anomalies. J. Chem. Phys., 84, 1985. [119] P. G. Debenedetti and M. C. D’Antonio. On the entropy changes and fluctu- ations occurring near a tensile instability. J. Chem. Phys., 85, 1986. [120] M. C. D’Antonio and P. G. Debenedetti. Loss of tensile strength in liquids without property discontinuities: a thermodynamic analysis. J. Chem. Phys., 86, 1987. [121] S. J. Henderson and R. J. Speedy. A Berthelod-Bourdon tube method for studying water under tension. J. Phys. E: Sci. Instrum., 13:778, 1980. [122] S. J. Henderson and R. Speedy. Temperature of maximum density in water at negative pressure. J. Phys. Chem., 91:3062–3068, 1987. [123] K. Davitt, E. Rolley, F. Caupin, A. Arvengas, and S. Balibar. Equation of state of water under negative pressure. J. Chem. Phys., 133:174507, 2010. [124] M. El Mekki Azouzi, C. Ramboz, J-F. Lenain, and F. Caupin. A coherent picture of water at extreme negative pressure. Nat. Phys., 9:38–41, 2012. [125] A. T. J. Hayward. Negative pressure in liquids: Can it be harnessed to serve man? American Scientist, 59:434–443, 1971. 125 [126] M. A. Gonza´lez, C. Valeriani, F. Caupin, and J. L. F. Abascal. A comprehen- sive scenario of the thermodynamic anomalies of water using the TIP4P/2005 model, 2016. arXiv:1605.05383 [cond-mat.stat-mech]. [127] J. W. Biddle, V. Holten, J. V. Sengers, and M. A. Anisimov. Thermal con- ductivity of supercooled water. Phys. Rev. E, 87:042302, 2013. [128] F. Bresme, J. W. Biddle, J. V. Sengers, and M. A. Anisimov. Communica- tion: Minimum in the thermal conductivity of supercooled water: A computer simulation study. J. Chem. Phys., 140:161104, 2013. [129] F. Mallamace, C. Corsaro, and H. E. Stanley. Relation of water structural relaxation to water anomalies. Proc. Natl. Acad. Sci. USA, 110:4893, 2013. [130] O. Benchikh, D. Fournier, and A. C. Boccara. Photothermal measurement of the thermal conductivity of supercooled water. J. Phys. (Paris), 46:727, 1985. [131] M. L. Huber, R. A. Perkins, D. G. Friend, J. V. Sengers, M. J. Assael, I. N. Metaxa, K Miyagawa, R. Hellman, and E Vogel. New international formulation for the thermal conductivity of H2O. J. Phys. Chem. Ref. Data, 41(3):033102, 2012. This is the official correlation for the thermal conductivity of water from the International Association for the Properties of Water and Steam (IAPWS). [132] I. S. Adrianova, O. Ya. Samoilov, and I. Z. Fisher. Thermal conductivity and structure of water. J. Struct. Chem., 8(5):813–816, 1967. [133] P. Kumar and H. E. Stanley. Thermal conductivity minimum: A new water anomaly. J. Phys. Chem. B, 115:14269–14273, 2011. [134] C. A. Angell, M. Oguni, and W. J. Sichina. Heat capacity of water at extremes of supercooling and superheating. J. Phys. Chem., 86(6):998, 1982. [135] P. W. Bridgman. The thermal conductivities of liquids. Proc. Natl. Acad. Sci. USA, 9:341–345, 1923. [136] Joseph O. Hirschfelder, Charles F. Curtiss, and R. Byron Bird. Molecular Theory of Gases and Liquids. John Wiley and Sons, 1954. [137] J. V. Sengers. Transport properties of fluids near critical points. Int. J. Thermophys., 6(3):203–227, 1985. [138] J Luettmer-Strathmann, J. V. Sengers, and G. A. Olchowy. Non-asymptotic critical behavior of the transport properties of fluids. J. Chem. Phys., 103(17):7482–7501, 1995. [139] M. L. Huber, R. A. Perkins, A. Laesecke, D. G. Friend, J. V. Sengers, M. J. Assael, I. N. Metaxa, E. Vogel, R. Maresˇ, and K. Miyagawa. New international formulation for the viscosity of H2O. J. Phys. Chem. Ref. Data, 38(2):101–125, 2009. 126 [140] J. Pa´tek, J. Hruby´, J. Klomfar, M. Soucˇkova´, and A. H. Harvey. Reference correlations for thermophysical properties of liquid water at 0.1 MPa. J. Phys. Chem. Ref. Data, 38(1):21–29, 2009. [141] Yu. A. Osipov, B. V. Zheleznyi, and N. F. Bondarenko. The shear viscosity of water supercooled to -35 ◦C. Russ. J. Phys. Chem., 51:1264, 1977. [142] C. A. Angell. Formation of glasses from liquids and biopolymers. Science, 267(5206):1924–1935, 1995. [143] K. Ito, C. T. Moynihan, and C. A. Angell. Thermodynamic determination of fragility in liquids and a fragile-to-strong liquid transition in water. Nature, 398:492, 1999. [144] A. Faraone, L. Liu, C.-Y. Mou, C.-W. Yen, and S.-H. Chen. Fragile-to-strong transition in deeply supercooled confined water. J. Chem. Phys., 121:10843– 10846, 2004. [145] R. S. Smith and B. D. Kay. The existence of supercooled liquid water at 150 K. Nature, 398:788–791, 1999. [146] F. W. Starr, C. A. Angell, and H. E. Stanley. Prediction of entropy and dynamic properties of water below the homogeneous nucleation temperature. Physica A, 323:51–66, 2003. [147] K. Kawasaki. Kinetic equations and time correlation functions of critical fluctuations. Ann. Phys. N.Y., 61:1–56, 1970. [148] G. A. Olchowy and J. V. Sengers. A simplified representation for the thermal conductivity of fluids in the critical region. Int. J. Thermophys., 10(2):417– 426, 1989. [149] R. A. Perkins, J. V. Sengers, I. M. Abdulagatov, and M. L. Huber. Simplified model for the critical thermal-conductivity enhancement in molecular fluids. Int. J. Thermophys., 34:191, 2013. [150] H. C. Burstyn, J. V. Sengers, and P. Esfandiari. Stokes-Einstein diffusion of critical fluctuations in a fluid. Phys. Rev. A, 22:282–284, 1980. [151] G. A. Olchowy and J. V. Sengers. Crossover from singular to regular behavior in the transport properties of fluids in the critical region. Phys. Rev. Lett., 61(1):15–18, 1988. [152] C. Vega and J. L. F. Abascal. Simulating water with non-polarizable rigid models: A general perspective. Phys. Chem. Chem. Phys., 13:19663–19688, 2010. [153] J. L. F. Abascal and C. Vega. Note: Equation of state and compressibility of supercooled water: Simulations and experiment. J. Chem. Phys., 134:186101, 2011. 127 [154] I. Shvab and R. J. Sadus. J. Chem. Phys., 139:194505, 2013. [155] M. Parrinello and A. Rahman. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys., 52:7182, 1981. [156] S. Nose´ and M. L. Klein. Constant pressure molecular dynamics for molecular systems. Mol. Phys., 50(5):1055–1076, 1983. [157] S. Nose´. A unified formulation of the constant temperature molecular- dynamics methods. J. Chem. Phys., 81(1):511–519, 1984. [158] S. Plimpton. Fast parallel algorithms for short-range molecular dynamics. J. Comp. Phys., 117, 1995. [159] R. Zwanzig. Time correlation functions and transport coefficients in statistical mechanics. Annu. Rev. Phys. Chem., 16:67, 1966. [160] T. W. Sirk, S. Moore, and E. F. Brown. Characteristics of thermal conductivity in classical water models. J. Chem. Phys., 138:064505, 2013. [161] F. Bresme, B. Kafskjold, and I. Wold. Nonequlibrium molecular dynamics study of heat conduction in ionic systems. J. Phys. Chem., 100:1879, 1996. [162] D. Rozmanov and P. Kusalik. Transport coefficients of the TIP4P-2005 water model. J. Chem. Phys., 136:044507, 2012. [163] F. Ro¨mer, A. Lervik, and F. Bresme. Nonequilibrium molecular dynamics simulations of the thermal conductivity of water: A systematic investigation of the SPC/E and TIP4P/2005 models. J. Chem. Phys., 137:074503, 2012. 128