Intermittent Dynamics, Stochastic Resonance and Dynamical

2221

2009, 113, 2221–2224 Published on Web 01/30/2009

Intermittent Dynamics, Stochastic Resonance and Dynamical Heterogeneity in Supercooled Liquid Water Biman Jana and Biman Bagchi* Solid State and Structural Chemistry Unit, Indian Institute of Science, Bangalore 560012, India ReceiVed: NoVember 4, 2008; ReVised Manuscript ReceiVed: December 11, 2008

Here we find through computer simulations and theoretical analysis that the low temperature thermodynamic anomalies of liquid water arises from the intermittent fluctuation between its high density and low density forms, consisting largely of 5-coordinated and 4-coordinated water molecules, respectively. The fluctuations exhibit strong dynamic heterogeneity (defined by the four point time correlation function), accompanied by a divergence like growth of the dynamic correlation length, of the type encountered in fragile supercooled liquids. The intermittency has been explained by invoking a two state model often employed to understand stochastic resonance, with the relevant periodic perturbation provided here by the fluctuation of the total volume of the system. Understanding the well-known anomalies of water at ambient pressure and temperature has been a long standing goal in physical and chemical sciences.1,2 These anomalies include a density maximum at 4 °C at atmospheric pressure, divergentlike behavior of the linear response functions (CP, κT, and RP) upon cooling, diffusivity maximum upon increasing pressure, and several others.3,4 Several of these anomalies become increasingly prominent as the temperature is lowered below the melting point. This apparent divergence of linear response functions at low temperature has been interpreted in terms of the presence of a liquid-gas spinodal at low temperature and positive pressure.5 Others have interpreted these anomalies as a consequence of the negative slope of the temperature density maximum (TMD) line in the (P, T) plane.6 A recent approach hypothesized a first order liquid-liquid (LL) phase transition separating high density liquid (HDL) (mostly 5-coordinated water molecules) and low density liquid (LDL) (mostly 4-coordinated water molecules) which ends in a critical point.7-9 The presence of a Widom line,8,10 emanating from the LL critical point, may explain the anomalous behavior. Some other studies have revealed the relationship between structural order and the anomalies of liquid water.11,12 In spite of the presence of several explanations, a microscopic origin of these anomalies is yet to be understood. The prime object of this Letter is to provide the same. To understand the molecular origin of these anomalies, we need to identify a microscopic order parameter that shows significant change across the anomalous (T ) 250 K at ambient pressure) region. To this end, we have calculated both the distribution and the average of the distances between a central water molecule (for model and simulation detail see ref 13) and its fifth nearest neighbor.14,15 The distribution and the average fifth neighbor distance show a sharp change between 250 and 260 K. Above 260 K, the distribution has a peak around 3.25 Å. Below 240 K, the peak of the distribution shifts to 3.85 Å. * Corresponding author. Electronic address: [email protected].

10.1021/jp809722w CCC: $40.75

The distribution becomes surprisingly flat between 260 and 250 K, with a mean value now located around 3.5 Å. We have also calculated the distribution of the distances of first four neighbors which shows only a gradual squeezing toward a lower distance with decreasing temperature. These distributions indicate that upon lowering temperature, though the first four neighbors (essential to maintain a tetrahedral geometry) of a water molecule come progressively closer, it is the fifth neighbor that shows maximum change. The temperature at which the peak position of the fifth neighbor distance distribution changes sharply coincides with the region where response functions show anomalous rise (as discussed below). The corresponding free energy surface (with the fifth neighbor distance as an order parameter) is quite flat, like in a critical phenomenon (Landau theory). Therefore, the number of five coordinated water molecules (N5C) can serve as a suitable order parameter to study the anomalous behavior of supercooled water. This order parameter has close similarity with the other structural order parameters (qtet, number of nearest neighbors and number of hydrogen bonded neighbors) previously used to study the anomalous properties of water.9,11 However, the present order parameter provides a more direct understanding of the response function anomalies in supercooled water, as will be discussed below. We have analyzed both the average and the time dependent occupation number of differently coordinated species in the system. We find that, although water molecules with coordination numbers 4 and 5 are abundantly present at all temperatures, they become the most dominant species, to the extent that they together constitute 80% or more of the total population below 270 K. It has been shown by Sciortino et al.16-18 through inherent structure analysis that the presence of 5-coordinated species in the extended network of liquid water should be regarded as normal (that is not as an excitation) constitution of water. Additionally, the isochoric differential X-ray experiments have demonstrated that the water molecules can exist in a  2009 American Chemical Society

2222 J. Phys. Chem. B, Vol. 113, No. 8, 2009

Letters

Figure 1. Temporal variation of the numbers of 4- and 5-coordinated species in the supercooled water. The variation has been shown for (a) T ) 280 K, (b) T ) 250 K, and (c) T ) 230 K. Note the strong anticorrelated fluctuations between N4C and N5C and intermittency in the fluctuations of both species at T ) 250 K.

distinct configuration with five neighbors in the first coordination shell (within 3.5 Å).19 Parts a-c of Figure 1 display the temporal variation of the number of the 4- and 5-coordinated species (N4C(t) and N5C(t)) in the system at three different temperatures T ) 280, 250, and 230 K, respectively. Both at the high (T ) 280 K) and at the low (T ) 230 K) temperature (away from the anomalous region, T ) 250 K), the fluctuations of N4C(t) and N5C(t) are fairly random and not spatially correlated, as shown in Figure 1a,c, respectively. But, near the anomalous region, however, the fluctuations of N4C(t) and N5C(t) are large, slow oscillatory and, surprisingly, strongly anticorrelated as shown in Figure 1b. The divergent-like growth of the linear response functions (found to peak at 250 K also) may originate from the strongly anticorrelated number fluctuations of 4- and 5-coordinated species, as further discussed below. The plots of the temporal variation of the number of 4- and 5-coordinated species at T ) 250 K (Figure 1b) clearly exhibits an intermittency. The signature of intermittent behavior can be confirmed by calculating the power spectrum of the fluctuation time correlation function.20-25 The power spectrum of the fluctuation of the 5-coordinated species (similarly for the 4-coordinated species also) is presented in Figure 2 for both T ) 300 and 250 K. The linear dependence of the power spectrum for a wide range of frequency starting from low frequency region with slope ∼-1 (in log-log plot) at T ) 250 K confirmed the strong intermittency present in the fluctuation. At T ) 300 K the power spectrum shows a linear dependence only at the higher frequency region, indicating the presence of small time scale intermittency in the fluctuations. This intermittency in the fluctuation of the number of species in the system is responsible for the anomalous behavior of the response functions. Although intermittent dynamics has been observed previously in the energy associated with hydrogen bond network rearrangement at room temperature,20-25 this is the first report of intermittency in N5C(t) at low temperatures. Moreover, 1/f dependence of the power spectrum also indicate the multiple time scale relaxation in the supercooled water which has been confirmed by several experimental studies.26,27

Figure 2. Intermittency and 1/f noise in the number fluctuation of the species in supercooled water. Power spectrum of the time correlation function of the fluctuations in the number of 5-coordinated species in log-log plot at T ) 300 and 250 K. Note the linear dependence over a wide range of frequency with slope ∼-1 (distinctive of 1/f noise) at T ) 250 K, which confirms the large time scale strong intermittency in the fluctuations (see Figure 1b). At T ) 300 K again a straight line with slope ∼-1 appears in the high frequency range, which indicates the presence of intermittency in the short time scale. The dashed line represent a line with slope ) -1 for comparison.

Next, we explain the intermittency by borrowing a model well-known in the area of stochastic resonance.28,29 We consider the motion of a Brownian particle in a bistable potential where two minima are HDL and LDL states,

1 1 V(x) ) - x2 + x4 2 4

(1)

A particle is acted upon by periodic forcing, A0 cos(ωt + φ), and random force η(t). The periodic forcing term appears due to fluctuation of the total volume of the system as two state have different specific volumes (VLDL > VHDL). Now, we can j + δV(t) and then write total volume of the system as V(t) ) V A0 ∝ (〈(δV)2〉)1/2. The random force, η(t), has zero-mean and obeys the fluctuation-dissipation theorem.

〈η(t) η(0)〉 ) 2kBTζδ(t)

(2)

Letters

J. Phys. Chem. B, Vol. 113, No. 8, 2009 2223

Here ζ is the dissipation coefficient. In the absence of the j ) the system will periodic forcing (i.e., at fixed total volume V fluctuate around the local minima (HDL and LDL) with the noise induced hopping between the local minima, which can be described by Kramer’s theory

rK )

1

√2πζ

exp(-βEact) )

1

√2πζ

( β4 )

exp -

(3)

Here β ) (1)/(kBT), and Eact ) 1/4. Next we define the probabilities that the system is either of the local minima as nH/L ) prob(state ) HDL/LDL), and this follows the simple relation nH ) 1 - nL between them. Now the governing rate equation for the system is

dnH dnL )) WL(δV(t))nL - WH(δV(t))nH (4) dt dt ) WL(δV(t)) - (WL(δV(t)) + WH(δV(t)))nH

Figure 3. Observed similarities between the mean squared number fluctuation of the five coordinated species (N5C) and the thermodynamic response function in supercooled water. Temperature dependence of mean square fluctuation 〈(δN5C)2〉, κT from simulation and κT from experiment. The data have been normalized at T ) 250 K. Note the striking similarity in the temperature dependence of both quantities. Note also that the quantities plotted in the figure overlap (within the standard deviation of experiment and simulation) with each other over the entire temperature range. For experimental κT values, the temperature has also been shifted by 18 K for the comparison (κT diverges at T ) 228 K experimental system and at T ) 246 K for model system).

where WH/L(δV(t)) is the volume dependent transition rate out of the HDL or LDL local minima in the presence of periodic forcing. The rate constant (WL/H(δV(t))) can be obtained under the periodic volume fluctuation as WL/H(δV(t)) ) rK exp[(βA0 cos(ωt)]). Equation 4 is now solved to obtain

nH(t) )

1 [n (t ) 0) + g(t) H

∫0t WL(δV(t′)) g(t′) dt′]

(5)

(6)

Figure 4. Time dependence of the nonlinear response function, χ4(t) in supercooled water. χ4(t) has been plotted at various temperatures, starting from 300 K down to 230 K. Note the huge growth of the peak of response function around the T ) 250 K.

At lower values of A0 and higher values of ω, eq 5 shows that nH(t) oscillates around 0.5, showing the local equilibrium situation in both states. This is the scenario at higher temperature away from T ) 250 K. However, as we approach 250 K, A0 becomes larger and ω becomes smaller In this case eq 5 shows that nH(t) oscillates intermittently between 0 and 1 (between two states). Simulations show long-lived periodic fluctuation in the volume of the system. Because the free energy surface and hence the rate of transitions between LDL and HDL forms depend on volume, the present model can explain the intermittency in the number fluctuation of 4- and 5-coordinated species near T ) 250 K (Figure 1b). We have calculated the temporal fluctuation of N5C and N4C using the volume fluctuation as obtained from the simulation at T ) 250 K. In this case, we have assumed a linear dependence of the rate 0 + ((∂WL/H)/(∂V))δV(t). on volume fluctuation as WL/H(δV(t)) ) WL/H 0 WL/H is obtained as the inverse of the corresponding lifetime and ((∂WL/H)/(∂V)) is obtained numerically from the slope of plot of 0 simulated WL/H with variation of the total volume of the system. We find that the calculated temporal fluctuation of N5C and N4C essentially capture the behavior of the fluctuation calculated from simulation (Figure 1b). The remaining differences observed between the fluctuation profiles derived from the model and simulation might be because of the nonlinear dependence of the rate on volume fluctuation. There are several other reports where the two state model of water has been introduced in the literature to explain the thermodynamics anomalies of supercooled water.30,31 The similarity of our model with the previous models is that the two

states of of the system are similar (high density and low density). However, the concept of stochastic resonance due to the modulation of the potential by volume fluctuation has not been introduced before. To further confirm the correlation between the fluctuations in N5C(t) and the anomalies in linear response functions, in Figure 3 we compare the temperature dependence of 〈(δN5C)2〉 with the isothermal compressibility, κT, which is related to the mean-squared volume fluctuation, (κT ) (〈(∆V)2〉)/(VkBT)).3,4 The calculated κT exhibits a maximum precisely at the same temperature where 〈(δN5C)2〉 exhibits a sharp maximum, at T ) 250 K. Moreover, with proper scaling, the values of the two overlap (within the standard deviations) with each other at each temperature over the entire temperature range. In the same figure we have plotted the experimental values of the κT, which shows that the present hypothesis that the species number fluctuation is responsible for the anomalous rise is in agreement with all the data. Similar correlation has also been found between 〈(δN5C)2〉 and 〈(∆E)2〉 (not shown here). Several dynamical features in the temperature dependence of the response functions of supercooled water bear striking similarity to those observed in supercooled liquids near their glass transition temperatures,32 which suggests the possibility of the emergence of a dynamic correlation length in low temperature water. Therefore, we have calculated a nonlinear response function, χ4(t), that measures the dynamical heterogeneity,33-35 which is related to the four-point density correlation function G4 by the following relation

with

g(t) ) exp[

∫0t (WL(δV(t′)) + WH(δV(t′))) dt′]

2224 J. Phys. Chem. B, Vol. 113, No. 8, 2009

χ4(t) )

βV N2

∫ dr1 dr2 dr3 dr4 w(|r1-r2|)

Letters References and Notes

×

w(|r4-r4 |) G4(r1,r2,r3,r4,t) (7)

where w(|r|) is an overlap function that is unity inside a region a and zero otherwise. The four-point density correlation function is defined as

G4(r1,r2,r3,r4,t) ) 〈F(r1,0) F(r2,t) F(r3,0) F(r4,t)〉 〈F(r1,0) F(r2,t)〉〈F(r3,0) F(r4,t)〉 (8) Now, it can be shown that χ4(t) can be written in a more useful form as

χ4(t) )

βV 2 (〈Q (t)〉 - 〈Q(t)〉2) 2 N

(9)

Here Q(t) is the time dependent order parameter that measures the localization of particle between time interval 0 to t inside a region a around its position at time t ) 0. Thus, χ4(t) is dominated by the range of spatial correlation between the localized particles in the fluid.33-36 Figure 4 displays the calculated nonlinear response function, χ4(t), of water molecules (computed with the oxygen atom displacement) for temperatures ranging from 300 K down to 230 K. For every temperature, χ4(t) is nearly zero at short and very long times with a maximum value at some intermediate time t*. Both t* and the amplitude of the peak, χ4(t*), increase sharply while approaching the temperature T ) 250 K (the region where response functions show anomaly for this model). However, after crossing T ) 250 K, the value of χ4(t*) sharply decreases with a further decrease in temperature. This suggests that the spatial correlation between the slow moving water molecules increases as we approach the anomalous region from both the sides. A growing dynamic correlation length is attributed to dynamical heterogeneity often observed in fragile supercooled liquids. Therefore, the huge dynamic heterogeneity observed here (Figure 4) has its origin in the spatially and temporally correlated, intermittent fluctuation in N5C(t) and N4C(t). The results presented here show the characteristics of an underlying self-organization near T ) 250 K. The observed growth of the dynamical correlation length shows that the observed intermittent dynamics is a consequence of a spatially and dynamically correlated process that may arise from a stochastic resonance between the low density and the high density liquid states and is driven by correlated volume fluctuation. In an inverse perspective, we can consider the intermittent fluctuation between N5C (the HDL) and N4C (the LDL) as the cause of the nearly periodic volume fluctuation. The nonlinear interaction between the two, aided by a nearly flat free energy surface of exchange between N5C and N4C, is perhaps the primary reason for the low temperature thermodynamic anomalies of liquid water. Acknowledgment. We thank Dr. S. M. Bhattacharyya, Prof. S. Sastry, Prof. G. Ananthakrishna and Prof. R. Ramaswamy for useful discussions. This work was supported in parts by a grant from DST, India. B.J. thanks CSIR, India, for providing SRF.

(1) Eisenberg, D.; Kauzmann, W. Structures and Properties of Water; Oxford: London, 1969. (2) Angell, C. A. In Water: A ComprehensiVe Treatise; Frank, F., Ed.; Plenum: New York, 1982; Vol. 7, pp 1-88. (3) Angell, C. A.; Shuppert, J.; Tucker, J. C. J. Phys. Chem. 1973, 77, 3092–3099. (4) Speedy, R. J.; Angell, C. A. J. Chem. Phys. 1976, 65, 851–858. (5) Speedy, R. J. J. Phys. Chem. 1982, 86, 982–991. (6) Sastry, S.; Debenedetti, P. G.; Sciortino, F.; Stanley, H. E. Phys. ReV. E 1996, 53, 6144–6154. (7) Poole, P. H.; Sciortino, F.; Essmann, U.; Stanley, H. E. Nature 1992, 360, 324–328. (8) Stanley, H. E.; Angell, C. A.; Essmann, U.; Hemmati, M.; Poole, P. H.; Sciortino, F. Physica A 1994, 205, 122–139. Stanley, H. E.; Cruz, L.; Harrington, S. T.; Poole, P. H.; Sastry, S.; Sciortino, F.; Starr, F. W.; Zhang, R. Physica A 1997, 236, 19–37. (9) Mishima, O.; Stanley, H. E. Nature 1998, 392, 164–168. (10) Franzese, G.; Stanley, H. E. J. Phys.: Condens. Matter 2007, 19, 205126. (11) Errington, J. R.; Debenedetti, P. G. Nature 2001, 409, 318–321. (12) Mallamace, F.; Corsaro, C.; Broccio, M.; Branca, C.; Gonz, H. E. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 12725–12729. (13) A series of eleven simulations at different temperatures starting from 300 to 220 K have been carried out in isobaric-isothermal ensemble (NPT). The water model used is TIP5P,37,38 which has been examined to reproduce the low temperature properties of water excellently.39,40 A system of 216 water molecules is first equilibrated for several nanoseconds (around 1 ns for 300 K and around 15 ns for 220 K) and then the trajectory is saved for several nanoseconds for the analysis. Both equilibration and subsequent run have been carried out in NPT ensemble. We have used 1.0 ps relaxation time for thermostat and 2.0 ps relaxation time for barostat. The density profile of water37 has been reproduced successfully as well. (14) Yan, Z.; Buldyrev, S. V.; Kumar, P.; Giovambattista, N.; Debenedetti, P. G.; Stanley, H. E. Phys. ReV. E 2007, 76, 051201. (15) Krekelberg, W. P.; Mittal, J.; Ganesan, V.; Truskett, T. M. Phys. ReV. E 2008, 77, 041201. (16) Sciortino, F.; Geiger, A.; Stanley, H. E. Nature 1991, 345, 218– 221. (17) Sciortino, F.; Geiger, A.; Stanley, H. E. J. Chem. Phys. 1992, 96, 3857–3865. (18) Sciortino, F.; Geiger, A.; Stanley, H. E. Phys. ReV. Lett. 1990, 65, 3452–3455. (19) Bosio, L.; Chen, S.-H.; Teixeira, J. Phys. ReV. A 1983, 27, 1468– 1475. (20) Ohmine, I.; Saito, S. Acc. Chem. Res. 1999, 32, 741–749. (21) Sasai, M.; Ohmine, I.; Ramaswamy, R. J. Chem. Phys. 1992, 96, 3045–3053. (22) Mudi, A.; Chakravarty, C.; Milotti, E. J. Chem. Phys. 2006, 125, 074508. (23) Mudi, A.; Chakravarty, C.; Ramaswamy, R. J. Chem. Phys. 2005, 122, 104507. (24) Mudi, A.; Chakravarty, C. J. Phys. Chem. B 2004, 108, 19607– 19613. (25) Sharma, R.; Chakravarty, C.; Milotti, E. J. Phys. Chem. B 2008, 112, 9071–9078. (26) Bellissent-Funel, M. C.; Longeville, S.; Zanotti, J. M.; Chen, S.H. Phys. ReV. Lett. 2000, 85, 3644–3647. (27) Torre, R.; Bertolini, P.; Righini, R. Nature 2004, 428, 296–299. (28) Gammaitoni, L.; Ha¨nggi, P.; Jung, P.; Marchesoni, F. ReV. Mod. Phys. 1998, 70, 223–287. (29) Ray, R.; Sengupta, S. Phys. Lett. 2003, 353, 364–371. (30) Tanaka, H. J. Chem. Phys. 2000, 112, 799–809. (31) Robinson, G. W.; Cho, C. H.; Urquidi, J. J. Chem. Phys. 1999, 111, 698–702. (32) Debenedetti, P. G.; Stillinger, F. H. Nature 2001, 410, 259–267. (33) Morishita, T. Phys. ReV. E 2008, 77, 020501(R).. (34) Dasgupta, C.; Indrani, A. V.; Ramaswamy, S.; Phani, M. K. Europhys. Lett. 1991, 15, 307–312. (35) Sharon, C. G.; Vladimir, N. N.; Thomas, B. S. J. Chem. Phys. 2000, 112, 509–512. (36) Cardenas, M.; Franz, S.; Parisi, G. J. Phys. A 1998, 31, L163L169; J. Chem. Phys. 1999, 110, 1726. (37) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910–8922. (38) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2001, 114, 363– 366. (39) Yamada, M.; Mossa, S.; Stanley, H. E.; Sciortino, F. Phys. ReV. Lett. 2002, 88, 195701. (40) Paschek, D. Phys. ReV. Lett. 2005, 94, 217802.

JP809722W