Thermal Conductivity of Supercooled Water: An Equilibrium Molecular

Letter pubs.acs.org/JPCL

Thermal Conductivity of Supercooled Water: An Equilibrium Molecular Dynamics Exploration Niall J. English*,† and John S. Tse*,‡ †

The SEC Strategic Research Cluster and the Centre for Synthesis and Chemical Biology, School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland ‡ Department of Physics and Engineering Physics, University of Saskatchewan, Saskatoon, Saskatchewan S7N 5E2, Canada S Supporting Information *

ABSTRACT: The thermal conductivity of both supercooled and ambienttemperature water at atmospheric pressure has been computed over the 140− 270 K temperature range for three popular water models via equilibrium molecular dynamics in the Green−Kubo setting. No strong temperature dependence of thermal conductivity was observed. The underlying phonon modes contributing to thermal conduction processes have been examined in the present work, and it has been established that (translational) acoustic modes dominate in supercooled water.

SECTION: Liquids; Chemical and Dynamical Processes in Solution

B

decreases continuously with decreasing temperature down to the lowest temperature of measurement, 250 32 and 256 K.33 Thermal diffusivity data at even lower temperatures produce no anomaly.32,33 Thermal conductivity is the product of thermal diffusivity, density, and isobaric heat capacity. Kumar and Stanley have computed TIP5P water thermal conductivity via nonequilibrium (NE-) molecular dynamics (MD), claiming a minimum in the supercooled region, at around 260 K; they conjecture that this mirrors changes in local high- and lowdensity structuring, offering possible support for the two-liquid hypothesis.34 Bresme et al. have also claimed a TIP4P/2005 minimum at around 215−220 K.35 Interestingly, Biddle et al. have developed a two-state thermodynamic model that also predicts a conductivity minimum, with it rising below ∼235 K.36 However, these TIP5P34 and TIP4P/200535 predictions contradict experimental data (continuously decreasing thermal conductivity to 250 K),32,33 while there is no experimental data below the ∼250 K range to compare against two-liquid thermodynamic model predictions.36 Kumar−Stanley TIP5P results34 contrast starkly to those of Mao and Zhang’s careful analysis,37 the latter giving ambient-temperature TIP5P results in much better agreement with experiment, while the present work also discounts minima for TIP4P, TIP4P/2005, and, to some extent, TIP5P.

ulk water has various anomalous characteristics that cannot be explained easily from a simple liquid’s perspective. Recent challenges from X-ray absorption1,2 and emission spectra3 to near-tetrahedral hydrogen bonding explaining water structure show that the liquid is highly distorted with two distinct local structures; a majority of molecules have strong asymmetrical hydrogen bonds, with a minority being tetrahedrally bonded,1 that is, bulk water at ambient conditions is a mixture of two different entities. X-ray scattering density-difference contrasts were attributed to temperature-dependent fluctuations of the two proposed local structure types,4−6 elaborated further to arise from a postulated liquid−liquid transition between high- (HDL) and low-density (LDL) liquid in supercooled water.7−17 The thermodynamics of (supercooled) water presents intriguing anomalies,18 tempting speculation of a liquid−liquid transition7−17 to explain them.18 The “two-liquid hypothesis” has been recast not in terms of coexisting but rather in terms of different local hydrogen-bonding environments19,20 to imply distinct local structures. An apparent fractional Stokes−Einstein relationship in liquid water has been interpreted as a two-liquid picture of water.21 However, these results are highly contested.22−31 The thermal conduction of supercooled water is intriguing; indeed, one might speculate that transport properties of putative two-liquid water display unexpected features in the supercooled region near the hypothesized liquid−liquid point. However, although there is negative dispersion in the speed of sound close to the vapor−liquid critical point that may arguably arise from a thermal conductivity divergence,18 there is no experimental evidence for anomalies in thermal conductivity; it © 2014 American Chemical Society

Received: August 8, 2014 Accepted: October 16, 2014 Published: October 20, 2014 3819

dx.doi.org/10.1021/jz5016179 | J. Phys. Chem. Lett. 2014, 5, 3819−3824

The Journal of Physical Chemistry Letters

Letter

tion)41 is similar to that of Chen et al.,65 who have scrutinized robust, best-practice GK approaches. For the 472-molecule system, Figure 1a shows conductivity essentially invariant between 180 and 270 K (under NPT) for

Therefore, we apply equilibrium MD across a wide temperature range (∼180−270 K) for explicit investigation of modes underlying conduction in ambient-pressure water and also for direct comparison to experiment. Wider questions as to water’s putative two-liquid nature are beyond this study’s gambit. TIP4P38 and TIP4P/200539 models were used, with TIP5P40 for selected temperatures; TIP4P/2005 and TIP5P are considered some of the better rigid, nonpolarizable popular water models.39,40 A brief synopsis is given here of the methodology, with further details in the Supporting Information.41 Ewald long-range electrostatics42 was used with tight error control.41 Velocity-Verlet MD42 with RATTLE43 was run with a 1 fs step in the NPT ensemble, although selected runs were performed under NVE at 180 and 270 K for 2048molecule systems to verify consistent results.44 Thermal conductivity can depend on system size;45 therefore, cubic systems containing 472, 1024, and 2048 molecules were simulated, with starting configurations generated and relaxed to 1 bar for 3−5 ns between 180 and 270 K, in 10 K increments, until volume-stabilized and local structure measures did not change. At each temperature, six independent 15−20 ns “production” NPT simulations were performed, giving at least a 100 ns aggregate. TIP5P was run for 2048-molecule systems at 180, 220, 270, and 298 K (with the latter 298 K runs to compare to Mao−Zhang results at 298 K37). No appreciable differences in JACFs (or thermal conductivities) were found from 10 versus 20 ns. Local density remained uniform (and fluctuations by grid-based analysis),46 and the molecular interaction energies’ probability distributions were unimodal. Although supercooled water structural relaxation usually takes place over picoseconds,47−49 increasing to perhaps tens (or hundreds) around (and below) 220 K, relaxation of thermal transport occurs much faster; this arises from energy density fluctuations, largely between nearest-neighboring molecules, itself related to acoustic frequencies.41 Although the homogeneous nucleation temperature of ambient-pressure water is ∼231 K experimentally50,51 (42 K below melting), suggesting the MD equivalent at ∼40−45 K below the respective models’ melting points, self-diffusivities determined from the meansquared displacement42 in the simulations’ latter halves did not yet indicate appreciable coarsening redolent of ice formation. Báez−Clancy ice−liquid distinction criteria were employed,52 where an angular order parameter quantifies the tetrahedrality of bonding for nearest-neighbor molecules; no ice-like structures were identified. The heat flux vector J was evaluated (eq S1, Supporting Information)41,53−55 every 2 fs, and the Green−Kubo (GK) integral of its autocorrelation function (JACF) was used for conductivity (eq S2, Supporting Information);41,56 defining JACF for 50 ps gives sampling of ∼300:1−400:1 for robust estimation.41,57 We have applied thermal conduction analysis to clathrates,45,57−62 water,57 and ices.57,63,64 GK’s essential advantage is that various thermal conduction modes may be gleaned from JACFs’ relevant components, usually characterized by two relaxation times, a rapid initial drop (τi,sh), followed by long decay (τi,lg), superimposed by water’s 30−40 fs librations (optic modes).41,58,59 These may be extracted by fitting JACFs to exponentially decaying functions, with cosine-modulated terms for optic component decay (eq S3, Supporting Information) and Fourier analysis to isolate modes.41,58,59 This approach (eqs S1−S3, Supporting Informa-

Figure 1. Thermal conductivity dependence on temperature for both models, with standard deviation uncertainties indicated: (a) 472, (b) 1024, and (c) 2048 molecules.

both TIP4P and TIP4P/2005 (∼0.77−0.83 and ∼0.88−0.93 W m−1 K−1, respectively). For larger systems, Figure 1b and c shows typical conductivities perhaps 1−2% larger, indicating that the smallest system is insufficiently large to capture all thermal conduction phonon modes; at first glance, the scant difference between results for 1024 and 2048 molecules suggests that between a few hundred and less than a thousand molecules appears reasonable. NVE-based results at 180 and 270 K were very similar to their 2048-molecule NPT counterparts. For larger systems, declining standard deviations are welcome. Upon closer examination, a more subtle “trend” of gently declining conductivity with reducing temperature appears, below the model’s melting point (TIP4P: ∼230 K; 3820

dx.doi.org/10.1021/jz5016179 | J. Phys. Chem. Lett. 2014, 5, 3819−3824

The Journal of Physical Chemistry Letters

Letter

K−1 going from ∼270 down to ∼220 K (with a minimum predicted at ∼250 K), disagree qualitatively, not only being 2.5−3 times larger than experiment but also in trend (i.e., essentially beginning to increase at ∼250 K with dropping temperature). Although Sirk et al.68 found that MP and GK agree (usually within ∼5%) for various models41 (alas, omitting TIP4P, TIP4P/2005, or TIP5P), they used over two thousand molecules, in contrast to only 512;34 moreover, only 40−45 water molecules are in each bin.34 This appears too small for reliable computation to achieve consistent agreement with GK41,68 (and conductivities much closer to experiment). Reflecting on MP,69 systems periodicity along the flux axis has this flux occurring in two directions; effective flux in one direction is halved, by definition.69 However, LAMMPS’s70 MP implementation “tallies” kinetic energy, which does not account for periodicity (or otherwise). After equilibration, if the temperature gradient is nonlinear, energy swapping is likely too frequent, and the system is outside of the linear-response regime; conductivity cannot be inferred accurately. Although overlooking periodicity yields artificial doubling of flux (and a similar conductivity overestimate, notwithstanding postequilibration observations of the gradient’s linearity or otherwise), it can only be conjectured that this is one possible reason for TIP5P’s MP-based 2.5−3-fold overestimate;34 there is insufficient information in ref 34 to conclude this with confidence, although the gradient appears relatively linear. Although Biddle et al.’s two-state model predicts a conductivity minimum,36 validating this is not possible, given that the predictive “departure” of ∼235 K occurs well below the experimental data’s ∼250 K lower limit. Error estimation by linear propagation, in eq 1 of ref 34 (flux equation), versus standard deviations derived from independent simulations here (a statistically rigorous approach),67 produces extremely erratic error bars;34 naturally, this renders robust identification of a minimum (using, e.g., ANOVA)67 very difficult. Unfortunately, there is no such attempt in refs 34 or 35; it is not at all clear if putative minima meet 90−95% confidence for null hypothesis rejection. Certainly, however, there is no conductivity minimum found in the present work. Bearing in mind smaller systems in refs 34 and 35 (512 and 878 molecules, respectively) and the finding here of an artifact of a minimum in the 1024molecule system (absent in the 2048-molecule case), this suggests that systems in refs 34 and 35 are too small for robust estimation, notwithstanding earlier technical comments on MP. One may conjecture that more aggressive subcooling closer to the homogeneous nucleation temperature in refs 34 and 35 (which, at 40−45 K below model melting points, would be ∼10−15 K lower than the “minima” in refs 34 and 35) could lead to partial coarsening toward more ice-like structures at temperatures below the minima. Given that ice’s conductivity is several times that of the liquid,57−59 this would lead naturally to an increase in the simulated value, suggesting a minimum. Although we have confirmed no coarsening here, we are unable to conclude this of refs 34 and 35; interestingly, Bresme et al. note that ∼80 ns was needed to observed any minimum, which was not so for shorter times.35 Examining JACF decay,41 respective short- and longer-term acoustic phonon relaxation times, τi,sh and τi,lg, were ∼44 ± 3 and ∼95 ± 6 fs (TIP4P) and ∼46 ± 3 and ∼106 ± 7 fs (TIP4P/2005), while being supercooled (i.e., below 230 and 250 K, respectively) with little temperature dependence. There was no real system size dependence in these findings. However, above melting temperatures, there was subtle acceleration in

TIP4P/2005: ∼250 K; TIP5P: 274 K),66 although there is little difference between the larger systems in this tentative decline. For the smallest system, conductivities fluctuate more noticeably, with a beguiling “hint” of “minimum” at ∼190 K for the 1024-moecule TIP4P system. However, this is an artifact; the feature is absent in the largest system. Larger system sizes (relative to a few hundred molecules) are required to capture this uncertain trend of declining conductivity with temperature below the melting point; it would appear that well over 1000 molecules are needed to avoid spurious minima. ANOVA67 confirmed null hypothesis (no difference between conductivities across the temperature range) nonrejection. Experimental data, inferred from direct diffusivity measurement (with conductivity as the product of diffusivity, density, and isobaric heat capacity), decline from ∼0.54 to ∼0.48 W m−1 K−1 with cooling from ∼270 to ∼250 K.32,33 For Taschin et al.’s33 diffusivities, uncertainties were estimated to be ±2% for 256−269 K; similar uncertainties are in the data of Benchikh et al. for 250−273 K,32 implying an uncertainty of ±0.01 W m−1 K−1 for 250−273 K, albeit with caution required due to compounding of separate uncertainties in each of the other three quantities. Given this lower uncertainty than that from largest-systems MD (not lower than 0.03 W m−1 K−1), experimental conductivity’s decline down to 23 °C below freezing appears statistically robust, but this drop is more tentative for the 20−25 °C below TIP4P or TIP4P/2005 melting points. TIP4P/2005 conductivities are similar in magnitude to those of Bresme et al., which used 878 molecules35 but with no minimum here; as noted, perhaps well over 1000 molecules is required to avoid spurious minima. Respective 2048-molecule TIP5P results at 180, 220, 270, and 298 K were 0.617 ± 0.029, 0.630 ± 0.033, 0.651 ± 0.032 and 0.668 ± 0.031 W m−1 K−1, which display a faint declining trend with lower temperature redolent of TIP4P-based cases (“echoing” experiment); however, ANOVA does not reject the null hypothesis. The 180−270 K data reflect the supercooled state, below freezing (274 K).66 The 298 K TIP5P result agrees with the 298 K Mao−Zhang37 value of 0.68 ± 0.007 W m−1 K−1 but disagrees with Kumar−Stanley34 data (see below). Despite TIP4P/2005’s more accurate melting point than TIP4P and better performance for (ambientcondition) water for other properties,39,68 TIP4P agrees “better” with experiment (i.e., about 50−70% overestimate of ∼0.82 W m−1 K−1 compared to ∼0.54 to ∼0.48 W m−1 K−1 for 250−273 K) than TIP4P/2005 (∼0.92−0.93 W m−1 K−1, a ∼70−90% overestimate). TIP5P, with an accurate melting point, does so “only” by 20−35%. Crafting a reliable model for supercooled conditions remains elusive.9−31 Naturally, gauging the relative accuracy of equilibrium (GK) and NE-MD is very relevant to comparing the present work to Kumar−Stanley’s findings.34 Such a more general, technical discussion is presented in the Supporting Information,41 drawing in part on Sirk et al.’s68 careful NEMD, Müller− Plathe (MP)69 versus GK comparisons; in summary, unlike “definitive GK” (where system size is an issue but somewhat less serious), several thousand water molecules per system appear necessary for NEMD, several hundred in each subsystem bin, with a thermal gradient of ∼1 K/nm, simulated for tens of nanoseconds, to guarantee good agreement with GK (within ∼5%).68 In contrast to Mao−Zhang 37 and our own TIP5P conductivities, the latter in closest agreement to experiment,32,33 Kumar−Stanley34 MP results, ∼1.2−1.4 W m−1 3821

dx.doi.org/10.1021/jz5016179 | J. Phys. Chem. Lett. 2014, 5, 3819−3824

The Journal of Physical Chemistry Letters

Letter

decay, with TIP4P τi,sh and τi,lg values declining consistently from ∼42.2 ± 2.1 to 36.1 ± 1.7 fs and ∼85.3 ± 5.1 to 66.2 ± 4.3 fs, going from 240 to 270 K; Student’s pairwise t tests suggest the 240 and 270 K results may be statistically different (and certainly lower than supercooled state values). Similar TIP4P/2005 findings applied for 260 and 270 K (i.e., above 250 K “melting”); respective τi,sh and τi,lg values were 38.3 ± 2.3 (260 K) and 37.2 ± 1.6 fs (270 K) and 74.4 ± 5.6 (260 K) and 72.5 ± 4.7 fs (270 K). The 260 and 270 K results are not statistically different, but they are lower than supercooled counterparts. This belies an underlying subtlety not readily evident from conductivity essentially invariant with temperature (cf. Figure 1); in the supercooled state, all others things being equal, it appears that slower relaxation of acoustic modes leads them to become more dominant in conduction. Although speculative at present, we confirm this below. It is instructive to examine the JACF’s full frequency distribution and any contributing phonons, including the rich tapestry of optical modes. The “beating” patterns of “interference” in the time domain, after full acoustic response decay has occurred, are illuminating (cf. latter term, eq S3, Supporting Information).41 A typical example is depicted for the TIP4P case in Figure 2 (for smaller systems, with similar

Figure 3. Power spectra for JACF of TIP4P water at 180, 210, and 270 K. The (librational) optical modes corresponding to the long-time beating patterns in Figure 2 are evident in the ∼600−900 cm−1 range. Also evident are the supercooled (translational) acoustic modes at ∼50 and 250 cm−1 and the rotational acoustic “shoulder” at 470−550 cm−1 not present at 270 K. TIP4P/2005 exhibits similar features. This is for 472 molecules as the larger systems displayed similar features.

2005 exhibited similar features. The ∼50 cm−1 translational acoustic mode arises from hydrogen bond O−O−O bending, while the 250 cm−1 mode originates from O−O stretching.71,72 The shoulder at 470−550 cm−1 arises from rotational motion73 but “softens” as the temperature is increased toward melting (i.e., 180 versus 210 K). These signature phonon “motifs” manifest themselves clearly in thermal conduction, as might be expected from earlier study of energy relaxation in supercooled water.73,74 The natural question arises as to quantitative contribution of each frequency mode to conductivity. Despite a recent phonon theory of liquids, including water,75 phonons are defined for crystals, strictly. As remarked previously for clathrates,58−60 the presence of a JACF spectral feature does not necessarily reveal the contribution that this mode makes (if any) to overall conductivity; this arises from different levels of damping in energy relaxation processes. To gauge the modes’ importance, the cumulative frequency conductivity,58,76 k(ν), obtained by integrating the reverse-transformed low-band-pass-filtered JACF spectrum up to ν,41 is shown in Figure 4. This reveals quantitatively the cumulative conductivity k(ν) arising from contributing modes up to the particular frequency ν. For TIP4P, this is a typical example (with little dependence on

Figure 2. “Beating” patterns of interference evident between “residual” optical modes (with periods of ∼40−48 fs) in TIP4P water’s normalized JACF (cf. the integrand in eq 2) at 180 and 210 K (both supercooled) and at 270 K, compare the latter term in eq S3 (Supporting Information). The acoustic modes have decayed by this time range (48−50 ps), having done so within around the first 5 ps (by examination of ln JACF and data -fitting). The 210 and 270 K data have been displaced upward by intervals of 0.0025 for ease of viewing. Refer also to the corresponding frequencies of the optical modes in the corresponding power spectra (cf. Figure 3) in the ∼700−850 cm−1 range. TIP4P/2005 exhibits similar features. This is for 472 molecules as the larger systems displayed similar features.

features in larger ones), both above and below TIP4P melting. The acoustic modes have decayed by this time (48−50 ps), having done so within the first ∼5 ps. The periods of ∼40−48 fs arise from librational (rotation oscillation); this feature has been observed previously in clathrates.58,59 The mixture of such signature librational vibrations with different periods is evident in the beating of Figure 2 and would give rise to a range of optical mode peaks in JACF power spectra (cosine Fourier transforms) over ∼600−900 cm−1. Figure 3 depicts typical examples of the TIP4P JACF spectra above and below the supercooled temperature range, where optical mode peaks are evident. It is also clear that there are supercooled (translational) acoustic modes at ∼50 and 250 cm−1 and a (rotational acoustic) “shoulder” at 470−550 cm−1 absent at 270 K. TIP4P/

Figure 4. Thermal conductivity evaluated from the GK integral of reverse-transformed low-band-pass-filtered JACF at each given frequency, k(ν), for TIP4P water at 180, 210, and 270 K. TIP4P/ 2005 exhibits similar features. This is for 472 molecules as the larger systems displayed similar features. 3822

dx.doi.org/10.1021/jz5016179 | J. Phys. Chem. Lett. 2014, 5, 3819−3824

The Journal of Physical Chemistry Letters

Letter

Ewald on GK. This material is available free of charge via the Internet at http://pubs.acs.org.

system size), again below and above melting. Frequencies in the GK approach are related to phonon−phonon interactions, which weights lower frequencies more heavily; the difference between GK and individual phonon frequencies is equivalent to that between the mode’s wavelength and its mean free path.76 This initial k(ν) increase arises from dominant contributions of translational acoustic modes, especially in supercooled water; indeed, it can be seen that there are minor shoulders in the vicinity of the respective ∼50 and 250 cm−1 modes, especially for the lower-frequency O−O−O mode, accentuated at deeper levels of supercooling (e.g., 180 K). The greater the extent of supercooling, the greater the degree of relative contribution from acoustic modes, as seen readily for 180 versus 210 K. As mentioned previously, the apparent temperature-invariance of overall conductivity (cf. Figure 1) belies this subtlety of acoustic contributions’ greater importance, hinted by JACF relaxation times (see above). The slight “hump” in the rotational acoustic contribution at ∼470−550 cm−1 emphasizes the more limited role of rotational motion in heat conduction, with softening for less aggressive supercooling. There is more limited contribution from optic modes, despite GK biasing lower frequency ones more heavily, especially for the supercooled case; this is unsurprising given that librational motion is more suppressed at lower temperatures.71−73 This lower optical contribution is not inconsistent with clathrates58,59 and amorphous ices,63,64 where acoustic modes dominate. It was determined, in the context of a rigorous GK formalism with appropriate incorporation of long-range electrostatics, that conductivity of supercooled TIP4P, TIP4P/2005, and TIP5P water exhibits no statistically significant dependence on temperature above 180 K. In fact, conductivity cannot be considered statistically different from the thermodynamically stable liquid state in the 20−40 K range above the respective melting points of the models (cf. Figure 1). The TIP4P and TIP4P/2005 results are in semiquantitative agreement with the experimental data of ∼0.54 to ∼0.48 W m−1 K−1 in cooling from ∼270 to ∼250 K, respectively,32,33 with TIP5P in even better accord; admittedly, the downward trend in conductivity with declining temperature was not observed with full ANOVA statistical rigor; larger systems hinted at this (cf. Figure 1b and c versus a). The lower limit of ∼0.03 W m−1 K−1 for MD uncertainty, compared to ∼0.01 W m−1 K−1 for experimental data, hinders any firm conclusion. What is indeed evident, however, is the absence of a minimum in conductivity and no evidence of rising conductivity with falling temperature below 230 K predicted by Biddle et al.36 Despite TIP5P’s more accurate melting point, these results contradict very directly the Kumar−Stanley34 NEMD-MP conductivity results of ∼1.2−1.4 W m−1 K−1, while the relevance of Biddle et al.’s36 two-state model conductivity predictions relative to experimental data cannot be deduced directly. Further, the underlying phonon modes contributing to heat conduction have been examined, and it has been established that (translational) acoustic modes dominate, with signature motifs of hydrogen bond O−O−O bending and O−O stretch leaving a trace, especially for deeper levels of supercooling.

■

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (N.J.E.). *E-mail: [email protected] (J.S.T.). Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors thank the Ireland Canada University Foundation and Royal Irish Academy. REFERENCES

(1) Myneni, S. C. B.; Luo, Y.; Naslund, L. A.; Ojamae, L.; Ogasawara, H.; Pelmenshikov, A.; Vaterlain, P.; Heske, C.; Pettersson, L. G. M.; Nilsson, A. J. Phys.: Condens. Matter 2002, 14, L213. (2) Wernet, P.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; Näslund, L. Å.; Hirsch, T. K.; Ojamäe, L.; Glatzel, P.; Pettersson, L. G. M.; Nilsson, A. Science 2004, 304, 995. (3) Tokushima, T.; Harada, Y.; Takahashi, O.; Senba, Y.; Ohashi, H.; Pettersson, L. G. M.; Nilsson, A.; Shin, S. Chem. Phys. Lett. 2008, 460, 387. (4) Huang, C.; Wikfeldt, K. T.; Tokushima, T.; Nordlund, D.; Harada, Y.; Bergmann, U.; Niebuhr, M.; Weiss, T. M.; Horikawa, Y.; Leetmaa, M.; Ljungberg, M. P.; Takahashi, O.; Lenz, A.; Ojamäe, L.; Lyubartsev, A. P.; Shin, S.; Pettersson, L. G. M.; Nilsson, A. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 15214. (5) Huang, C.; Weiss, T. M.; Nordlund, D.; Wikfeldt, K. T.; Pettersson, L. G. M.; Nilsson, A. J. Chem. Phys. 2010, 133, 134504. (6) Wikfeldt, K. T.; Huang, C.; Nilsson, A.; Pettersson, L. G. M. J. Chem. Phys. 2011, 134, 214506. (7) Mishima, O. Nature 1996, 384, 546. (8) Mishima, O.; Stanley, H. E. Nature 1998, 396, 329. (9) Debenedetti, P. G.; Stanley, H. E. Phys. Today 2003, 56, 40. (10) Poole, P. H.; Sciortino, F.; Essmann, U.; Stanley, H. E. Nature 1992, 360, 324. (11) Harrington, S.; Zhang, R.; Poole, P. H.; Scriortino, F.; Stanley, H. E. Phys. Rev. Lett. 1997, 78, 2409. (12) Kesselring, T. A.; Franzese, G.; Buldyrev, S. V.; Herrmann, H. J.; Stanley, H. E. Sci. Rep. 2012, 2, 474. (13) Liu, Y.; Palmer, J. C.; Panagiotopoulos, A. Z.; Debenedetti, P. G. J. Chem. Phys. 2012, 137, 214505. (14) Holten, V.; Bertrand, C. E.; Anisimov, M. A.; Sengers, J. V. J. Chem. Phys. 2012, 136 (094), 507. (15) Gallo, P.; Sciortino, F. Phys. Rev. Lett. 2012, 109, 177801. (16) Poole, P. H.; Bowles, R. K.; Saika-Voivod, I.; Sciortino, F. J. Chem. Phys. 2013, 138, 034505. (17) Palmer, J. C.; Martelli, F.; Liu, Y.; Car, R.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Nature 2014, 510, 385. (18) Holten, V.; Bertrand, C. E.; Anisimov, M. A.; Sengers, J. V. J. Chem. Phys. 2012, 136, 094507. (19) Nilsson, A.; Pettersson, L. G. M. Chem. Phys. 2011, 389, 1. (20) Nilsson, A.; Huang, C.; Pettersson, L. G. M. J. Mol. Liq. 2012, 176, 2. (21) Xu, L.; Mallamace, F.; Yan, Z.; Starr, F. W.; Buldyrev, S. V.; Stanley, H. E. Nat. Phys. 2009, 5, 565. (22) Soper, A. K.; Teixeira, J.; Head-Gordon, T. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, E44. (23) Huang, C.; Wikfeldt, K. T.; Tokushima, T.; Nordlund, D.; Harada, Y.; Bergmann, U.; Niebuhr, M.; Weiss, T. M.; Horikawa, Y.; Leetmaa, M.; Ljungberg, M. P.; Takahashi, O.; Lenz, A.; Ojamäe, L.; Lyubartsev, A. P.; Shin, S.; Pettersson, L. G. M.; Nilsson, A. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, E45. (24) Clark, G. N. I.; Hora, G. L.; Teixeira, J.; Soper, A. K.; HeadGordon, T. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 14003.

ASSOCIATED CONTENT

S Supporting Information *

The MD methodology, further notes on relaxation, calculations of the thermal conductivity, GK versus NEMD for water-based systems, and observations of (FFT-based) reciprocal-space 3823

dx.doi.org/10.1021/jz5016179 | J. Phys. Chem. Lett. 2014, 5, 3819−3824

The Journal of Physical Chemistry Letters

Letter

(25) Soper, A. K. Pure Appl. Chem. 2010, 82, 1855. (26) Limmer, D. T.; Chandler, D. J. Chem. Phys. 2011, 135, 134503. (27) Limmer, D. T.; Chandler, D. J. Chem. Phys. 2012, 137, 044509. (28) Limmer, D. T.; Chandler, D. J. Chem. Phys. 2013, 138, 214504. (29) Limmer, D. T.; Chandler, D. Faraday Discuss. 2013, 167, 485− 498. (30) English, N. J.; Kusalik, P. G.; Tse, J. S. J. Chem. Phys. 2013, 139, 084508. (31) Chandler, D. Illusions of phase coexistence: Comments on “Metastable liquid−liquid transition ...” by J. C. Palmer et al., Nature 510, 385 (2014), arXiv:1407.6854. (32) Benchikh, O.; Fournier, D.; Boccara, A. C.; Teixeira, J. J. Phys. (France) 1985, 46, 727−731. (33) Taschin, A.; Bartolini, P.; Ramo, R. E.; Torre, R. Phys. Rev. E 2006, 74, 031502. (34) Kumar, P.; Stanley, H. E. J. Phys. Chem. B 2011, 115, 14269. (35) Bresme, F.; Biddle, J. W.; Sengers, J. V.; Anisimov, M. A. J. Chem. Phys. 2014, 140, 161104. (36) Biddle, J. W.; Holten, V.; Sengers, J. V.; Anisimov, M. A. Phys. Rev. E 2013, 87, 042302. (37) Mao, Y.; Zhang, Y. Chem. Phys. Lett. 2012, 542, 37−41. (38) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (39) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (40) Mahoney, M. W. J. Chem. Phys. 2000, 112, 8910. (41) See the Supporting Information for further details on the MD methodology, structural and thermal relaxation, definition of the heat flux vector and estimation of thermal conductivity, overview of GK versus NEMD approaches, and incorporation of (FFT-based) reciprocal space Ewald contributions into the heat flux vector. (42) Allen, M.P.; Tildesley, D.J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987. (43) Andersen, H. C. J. Comput. Phys. 1983, 52, 24−34. (44) Melchionna, S.; Ciccotti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533−544. (45) English, N. J. Mol. Phys. 2008, 106, 1887−1898. (46) English, N. J.; Tse, J. S. Phys. Rev. Lett. 2011, 106, 037801. (47) Torre, R.; Bartolini, P.; Righini, R. Nature 2004, 428, 296−299. (48) Liu, L. Study of Slow Dynamics in Supercooled Water by Molecular Dynamics and Quasi-Elastic Neutron Scattering. Ph.D. Thesis; MIT: Cambridge, MA, 2005. (49) Saito, S.; Ohmine, I.; Bagchi, B. J. Chem. Phys. 2013, 138, 094503. (50) Cwilong, B. M. Proc. R. Soc. London, Ser. A 1947, 190, 137. (51) Mossop, S. C. Proc. Phys. Soc., London, Sect. B 1955, 68, 193. (52) Báez, L. A.; Clancy, P. Ann. N.Y. Acad. Sci. 1994, 715, 177. (53) Vogelsang, R.; Hoheisal, C. Phys. Rev. A 1987, 35, 3487−3491. (54) Daivis, P. J.; Evans, D. J. Mol. Phys. 1994, 81, 1289−1295. (55) Petravic, J. J. Chem. Phys. 2005, 123, 174503. (56) Evans, D.J.; Morriss, G.P. Statistical mechanics of non-equilibrium liquids; Academic Press: San Diego, CA, 1990. (57) Rosenbaum, E. J.; English, N. J.; Johnson, J. K.; Shaw, D. W.; Warzinski, R. P. J. Phys. Chem. B 2007, 111, 13194−13205. (58) English, N. J.; Tse, J. S. Phys. Rev. Lett. 2009, 103, 015901. (59) English, N. J.; Tse, J. S.; Carey, D. Phys. Rev. B 2009, 80, 134306. (60) English, N. J.; Tse, J. S. Comput. Mater. Sci. 2010, 49, S176− S180. (61) English, N. J.; Tse, J. S. Energies 2010, 3, 1934−1942. (62) English, N. J.; Gorman, P. D.; MacElroy, J. M. D. J. Chem. Phys. 2012, 136, 044501. (63) English, N. J.; Tse, J. S.; Gallagher, R. Phys. Rev. B 2010, 82, 092201. (64) English, N. J.; Tse, J. S. Phys. Rev. B 2011, 83, 184114. (65) Chen, J.; Zhang, G.; Li, B. Phys. Lett. A 2010, 374, 2392−2396. (66) Fernández, R. G.; Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2006, 124, 144506. (67) Hogg, R.V.; Craig, A.; McKean, J.W. Introduction to Mathematical Statistics, 6th ed.; Macmillan: New York, 2004.

(68) Sirk, T. W.; Moore, S.; Brown, E. F. J. Chem. Phys. 2013, 138, 064505. (69) Müller-Plathe, F. J. Chem. Phys. 1997, 106, 6082. (70) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (71) Walrafen, G. E.; Chu, Y. C.; Piermarini, G. J. J. Phys. Chem. 1996, 100, 10363. (72) Garberoglio, G.; Vallauri, R.; Sutmann, G. J. Chem. Phys. 2002, 117, 3278. (73) Yagasaki, T.; Saito, S. J. Chem. Phys. 2011, 135, 244511. (74) Ohara, T. J. Chem. Phys. 1999, 111, 6492. (75) Bolmatov, D.; Brazhkin, V. V.; Trachenko, K. Sci. Rep. 2012, 2, 421. (76) McGaughey, A. J. H.; Kaviany, M. Int. J. Heat Mass Transfer 2004, 47, 1783.

3824

dx.doi.org/10.1021/jz5016179 | J. Phys. Chem. Lett. 2014, 5, 3819−3824