Proton Disorder in Ice Ih and Inhomogeneous Broadening in Two

Aug 12, 2013 - Joel M. Bowman , Yimin Wang , Hanchao Liu , and John S. Mancini. The Journal of ... Liang Shi , J. L. Skinner , Thomas L. C. Jansen. Ph...
1 downloads 0 Views 1MB Size
Article pubs.acs.org/JPCB

Proton Disorder in Ice Ih and Inhomogeneous Broadening in TwoDimensional Infrared Spectroscopy L. Shi and J. L. Skinner* Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin, Madison, Wisconsin 53706, United States ABSTRACT: It is well-known that in ice Ih the oxygen atoms form a regular hexagonal lattice while the positions of the hydrogen atoms are disordered, called proton disorder in the literature. Various OH (OD) stretch vibrational spectroscopies (e.g., IR, Raman, two-dimensional IR (2DIR), and hole burning) have been used to probe this proton disorder in the past several decades. However, the presence and the magnitude of the inhomogeneous broadening due to this proton disorder in the vibrational spectroscopy is still controversial. In this work, we calculate 2DIR spectroscopy for HOD in D2O ice Ih at 80 K with a mixed quantum/classical approach, and make comparison to a recent 2DIR experiment on the same system. Fair agreement is achieved between theory and experiment, although the calculated 2DIR line shape shows inhomogeneous broadening that was not observed in the experiment. However, the theory reproduces the linear IR for the same system fairly well, and the inhomogeneous broadening from the calculation is consistent with the extrapolation of the experimental IR line-widths in the literature. The effect of this proton disorder on the 2DIR line shape is explored in detail. We also calculate the vibrational three-pulse photon echo peak shift signal, which shows signatures of both low-frequency dynamics and inhomogeneous broadening.

I. INTRODUCTION Ice Ih, the most common form of solid water on earth, has been the subject of great scientific interest for over a century.1 Although the lattice structure (the oxygen atom positions) of ice Ih was determined by X-ray diffraction in the 1920s,2 it took about another 10 years for scientists to propose a complete structure (i.e., Pauling’s structure3), and it was not until the development of neutron diffraction techniques that this structure was confirmed directly.4 The intriguing feature of Pauling’s structure is that the hydrogen atoms are distributed at “random” in the regular oxygen lattice, subject to the Bernal− Fowler ice rules,5 which is termed proton disorder in the recent literature.6 This proton disorder leads to many peculiar properties of ice Ih, such as residual entropy, large static dielectric permittivity, etc.7 OH-stretch vibrational spectroscopy provides an alternative tool to study intermolecular structure in ice due to the sensitivity of vibrational frequencies to local environments. Historically, the similarities between the Raman spectra (OH stretch region) of water vapor, liquid water and ice led Bernal and Fowler in their 1933 paper to speculate that the H2O molecules’ integrity is retained in ice.5 IR and Raman spectra for H2O and D2O ice have numerous peaks and shoulders, and their interpretation has been difficult and controversial.8 We have supported a view9 that the spectral features are primarily due to intermolecular vibrational coupling and inhomogeneous broadening arising from a distribution of local-mode vibrational frequencies, in this case due to proton disorder.8,10,11 Understanding inhomogeneous broadening is particularly important because it provides an experimental window for understanding proton disorder. © 2013 American Chemical Society

In this paper, we focus primarily on this inhomogeneous broadening. To this end, one can eliminate the effects of intermolecular (and intramolecular) vibrational coupling by considering dilute HOD in D2O,12−14 in which case the local OH stretch is essentially a normal mode of the system due to the large frequency mismatch between OH and OD stretches. Unlike the situation for neat H2O or D2O, in this case, experimental OH IR spectra (or OD spectra for dilute HOD in H2O) consist of single peaks. The width of this peak is temperature dependent, ranging from about 125 cm−1 near the melting point to about 65 cm−1 at 86 K.15 Note, however, that one experiment has a line-width below 50 cm−1 at 70 K.16 Raman spectra have been measured down to lower temperatures, in one case finding a line-width as small as 25 cm−1 at 10 K.17 One can then ask: what fraction of the line-width is due to inhomogeneous broadening, and what fraction is due to homogeneous broadening, in this case due to phonons and the excited-state lifetime? Although the nature and extent of proton disorder is essentially temperature independent,3 the observed strong temperature dependence of the line-width shows that it must have at least some homogeneous contribution. On the other hand, for temperatures much less than the Debye temperature (which for ice is approximately 226 K7), the homogeneous contribution from phonons should Special Issue: Michael D. Fayer Festschrift Received: June 13, 2013 Revised: August 9, 2013 Published: August 12, 2013 15536

dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544

The Journal of Physical Chemistry B

Article

become unimportant.18,19 Therefore, at 80 K, one might expect that the line-width would be dominated by an inhomogeneous contribution from proton disorder and a homogeneous contribution from the excited state lifetime. The latter is about 400 fs and roughly temperature independent,20 which by itself would lead to a line-width of 13 cm−1. A convolution of a Lorentzian with this width, and a Gaussian, suggests that the experimental width at 86 K of 65 cm−1 therefore has an inhomogeneous contribution with a width of 58 cm−1. In a previous paper, we have modeled the spectroscopy of the HOD/D2O ice Ih system.15 Our temperature-dependent IR line-widths are in reasonable agreement with experiment. The density of ice depends on temperature (it increases as T decreases),7 so even though the proton disorder is temperature independent, we found (see Figure 12 in ref 15) that the inhomogeneous broadening is weakly temperature dependent, with widths ranging from 59 cm−1 at 245 K to 63 cm−1 at 125 K to 70 cm−1 at 1 K.15 A simple interpolation suggests a theoretical inhomogeneous width of about 65 cm−1 at 80 K due to the proton disorder. Recently, the emergence of ultrafast vibrational spectroscopy renders the experimental extraction of the inhomogeneous width possible at relatively high temperature.21 Graener and coworkers carried out a picosecond IR hole-burning experiment on HOD in D2O ice Ih at 230 K in 1997.22 They found an inhomogeneous width of ∼50 cm−1, attributed to proton disorder. Another cutting-edge nonlinear technique is twodimensional infrared spectroscopy (2DIR), which adds an additional frequency dimension to the linear IR, analogous to two-dimensional NMR.23 If inhomogeneous broadening is dominant, the 2DIR spectrum will usually be elongated along the diagonal. If homogeneous broadening is dominant, the 2DIR spectrum will be more round. If both broadening mechanisms are important, the antidiagonal line-width of the 2DIR spectrum will provide some information about the homogeneous broadening, while the diagonal line-width will reflect both broadening mechanisms.23 A very recent 2DIR study on isotope-diluted ice Ih at 80 K by Hamm and coworkers showed little (less than 25 cm−1) inhomogeneous broadening in their 2DIR spectra.24 To summarize, then, linear IR spectroscopy suggests that the inhomogeneous line-width in HOD/D2O ice Ih at 86 K is about 58 cm−1, hole-burning experiments at 230 K suggest an inhomogeneous width of about 50 cm−1, 2DIR experiments at 80 K are consistent with an inhomogeneous width below 25 cm−1, and theoretical calculations indicate an inhomogeneous width of about 65 cm−1 at 80 K and 59 cm−1 at 245 K. With the success of our theoretical approach in simulating IR and Raman spectra for both the isotope-diluted15 and neat H2O (D2O) ice Ih,8,10,11 and in reproducing the 2DIR spectrum for isotope-diluted liquid water,25 the above controversy over the magnitude of the inhomogeneity due to the proton disorder in vibrational spectroscopy motivates us to extend our theoretical model to nonlinear vibrational spectra for isotope-diluted ice Ih. In this work, we focus on modeling the 2DIR spectrum for HOD in D2O ice Ih (the linear IR spectrum for HOD in H2O ice Ih is complicated by the libration−bend combination band of H2O). In addition, we calculate the vibrational three-pulse photon-echo peak-shift (3PEPS) experiment on isotope-diluted ice Ih, which should be a direct probe of low-frequency dynamics and inhomogeneous broadening. The paper is organized as follows. In section II, we outline our theoretical methods for the spectral calculations. In section

III, we compare our calculated linear IR and 2DIR spectra for HOD/D2O ice Ih at 80 K to experiment. We analyze our theoretical spectra, and also calculate the 3PEPS observable for ice Ih. In section IV, we conclude.

II. THEORETICAL SPECTROSCOPY To calculate linear and nonlinear spectra for HOD/D2O ice Ih, we use our mixed quantum/classical approach, where the OH stretch is treated quantum mechanically, OD stretches and all bends are frozen, and molecular rotations and translations are treated classically.15 The expression for the IR absorption spectrum for light with polarization j ̂ is given by15 Ij(ω) ∼ Re

∫0



dt eiωt ⟨μ10j (0)μ10j (t )U10(t )⟩e−t /2T1

(1)

where U10(t) = exp[−i ∫ t0 ω10(τ) dτ], ω10(τ) is the fluctuating 1−0 transition frequency of the OH stretch, μj10 is the jth Cartesian component of the corresponding transition dipole, and the brackets denote a time average for each chromophore and an average over all chromophores. Vibrational relaxation is treated phenomenologically, and the lifetime is T1 = 0.4 ps for HOD in D2O.20 For a polycrystalline sample, one must average over all possible orientations of the polarization j ̂ with respect to the crystal axes. From this, one finds that I(ω) = [Ix(ω) + Iy(ω) + Iz(ω)]/3 for the simulation box frame Cartesian coordinates. The third-order response functions have four polarization indices, and are given by21,26 R1ijkl(t3 , t 2 , t1) = ⟨μ10i (0)μ10j (t1)μ10k (t1 + t 2) * (t1)U10 * (t1 + t 2) μ10l (t1 + t 2 + t3)U10 U10(t1 + t 2 + t3)⟩ R 2ijkl(t3 , t 2 , t1) = R1ijkl(t3 , t 2 , t1) R3ijkl(t3 , t 2 , t1) = −⟨μ10i (0)μ10j (t1)μ21k (t1 + t 2) * (t1)U21 * (t1 + t 2) μ21l (t1 + t 2 + t3)U10 U21(t1 + t 2 + t3)⟩ R 4ijkl(t3 , t 2 , t1) = ⟨μ10i (0)μ10j (t1)μ10k (t1 + t 2) * (t1 + t 2) μ10l (t1 + t 2 + t3)U10(t1)U10 U10(t1 + t 2 + t3)⟩ R 5ijkl(t3 , t 2 , t1) = R 4ijkl(t3 , t 2 , t1) R 6ijkl(t3 , t 2 , t1) = −⟨μ10i (0)μ10j (t1)μ21k (t1 + t 2) * (t1 + t 2) μ21l (t1 + t 2 + t3)U10(t1)U21 U21(t1 + t 2 + t3)⟩

(2)

where t1, t2, and t3 are the three time delays in the four-wavemixing experiments and U21(t) = exp[−i∫ t0 ω21(τ) dτ]. ω21(τ) is the fluctuating 2−1 transition frequency of the OH stretch, and μi21 is the ith Cartesian component of the corresponding 2−1 transition dipole. Herein we will only consider the parallel polarization geometry. For a polycrystalline sample, the average over all orientations of the polarizations with respect to the crystal axes gives27,28 Rn =

∑ R nijkl(δijδkl + δikδjl + δilδjk)/15 ijkl

15537

(3) dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544

The Journal of Physical Chemistry B

Article

for the 2−1 transition (e.g., ω21) are taken from ref 40. To calculate spectra, we average over the putative substitution of each D with an H (and use the maps for the OH spectroscopy).

where each sum is over the three simulation box frame Cartesian coordinates and δij is the Kronecker delta. The 2DIR line shape in the impulsive limit is then given by23,26,29

III. RESULTS AND DISCUSSION In Figure 1, the experimental IR line shape for HOD/D2O ice Ih at 80 K is plotted (black solid line).24 This line shape is more

6

S(ω3 , t 2 , ω1) ∼ Re ∑ Si(ω3 , t 2 , ω1)

(4)

i=1

where S1(ω3 , t 2 , ω1) = S2(ω3 , t 2 , ω1) =

∫0



dt 3

∫0



dt1

exp(iω3t3 − iω1t1) R1(t3 , t 2 , t1)L1(t3 , t 2 , t1) S3(ω3 , t 2 , ω1) =

∫0



dt 3

∫0



dt1 exp(iω3t3 − iω1t1)

R3(t3 , t 2 , t1)L3(t3 , t 2 , t1) S4(ω3 , t 2 , ω1) = S5(ω3 , t 2 , ω1) =

∫0



dt 3

∫0



dt1

exp(iω3t3 + iω1t1) R 4(t3 , t 2 , t1)L4(t3 , t 2 , t1) S6(ω3 , t 2 , ω1) =

∫0



dt 3

∫0



Figure 1. Calculated (from eq 1) and experimental24 infrared (IR) line shapes for HOD/D2O ice Ih at 80 K. Note that the baseline of the original experimental spectrum24 has been subtracted, and the absorbance is converted to the line-shape function. The red dashed line is the line shape for HOD/D2O ice Ih with the approximate linear response function (from eq A6). All line shapes are normalized to have the same peak height.

dt1 exp(iω3t3 + iω1t1)

R 6(t3 , t 2 , t1)L6(t3 , t 2 , t1)

(5)

and L1(t3 , t 2 , t1) = L 2(t3 , t 2 , t1) = L4(t3 , t 2 , t1) = L5(t3 , t 2 , t1) = e−(t3+ 2t2 + t1)/2T1 L3(t3 , t 2 , t1) = L6(t3 , t 2 , t1) = e−(3t3+ 2t2 + t1)/2T1

or less consistent with previous older data,15 although a bit broader (full width at half-maximum (fwhm) is 79 cm−1, compared to 65 cm−1 at 86 K15). The theoretical IR line shape at 80 K is also plotted (red solid line). It is in reasonable agreement with the experiment, although the calculated line shape is slightly red-shifted and narrower (fwhm of 75 cm−1) compared to the experimental line shape (but 10 cm−1 broader than the older experiment). Both line shapes show a single peak with little if any structure. Next we consider the 2DIR spectra (parallel polarization geometry) for HOD in D2O ice Ih at 80 K. The experimental 2DIR spectra were reported in ref 24, and only the 2DIR spectrum with the waiting time t2 = 200 fs is replotted in the right panel of Figure 2 using the raw experimental data.41 The apparent anharmonicity (difference between the 1−0 and 2−1 peak maxima) in this spectrum is about 285 cm−1. Note that, at t2 = 100 fs, the apparent anharmonicity is smaller,24 presumably due to the coherent artifact. The dashed line in the figure is the diagonal. The 1−0 transition peak (red lobe) appears to be primarily homogeneously broadened (it is approximately circular), which led Hamm and co-workers to the conclusion that the proton disorder in ice Ih hardly causes perceivable inhomogeneous broadening.24 The 2−1 transition peak (blue lobe) is much broader than the 1−0 transition peak, which is also observed in the pump−probe42 and hole-burning22 experiments. If we compare the experimental HOD/D2O 2DIR spectra for ice Ih to those for liquid water,43 there are two major differences: (1) the 2DIR spectra for ice barely change with increasing waiting time (at least for the 1−0 transition peak), while the 2DIR spectra for liquid water become circular gradually;44,45 (2) the 2DIR spectrum for ice is much narrower than that for liquid water.

(6)

Note that, as before, the lifetime is taken into account phenomenologically.29,30 In order to implement the mixed quantum/classical approach, we employ the explicit three-body (E3B) rigid water model recently developed31 and reparametrized32 in our group, and a 432-molecule proton-disordered ice Ih configuration generated by Hayward and Reimer.33 A molecular dynamics simulation of D2O ice Ih is then performed at 80 K in the NVT ensemble with a modified GROMACS version 3.334 for the E3B model. The simulation box is scaled to match the lattice constants at 80 K, interpolated from the experimental values at 70 and 85 K.7 The simulation time step is 1 fs. The SETTLE algorithm35 is used for constraining the molecular geometry, and the Berendsen thermostat36 is employed with a 0.5 ps coupling constant to maintain the temperature in the production run. Periodic boundary conditions are applied with the particle-mesh Ewald37,38 method for the electrostatic interactions. The cutoff for the Lennard-Jones potential between oxygen atoms is set to 0.95 nm. The essential step in the mixed quantum/classical approach is to determine the spectroscopic variables (e.g., ω10 and μi10) from the simulations. In principle, these variables depend on all the nuclear degrees of freedom in the system, but previous work in our group found the electric field E at the water hydrogen along the OH bond direction correlates very well with these variables.39 The correlations (also called maps) have been parametrized from DFT calculations on water clusters, extracted from liquid water simulations.8,10,15 The maps used in this work are the same as those in ref 11 except that the maps 15538

dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544

The Journal of Physical Chemistry B

Article

transition frequencies. The dashed lines are the best Gaussian fits for the histograms. The standard deviations are about 26 and 32 cm−1 for the 1−0 and 2−1 transitions, respectively. The corresponding fwhm of the 1−0 distribution is about 62 cm−1. This is a few cm−1 less than the interpolated value of 65 cm−1 from our previous work;15 this small discrepancy presumably arises because the earlier work used our first version of the E3B model.31 Since the width of the theoretical absorption spectrum at 80 K is 75 cm−1, this means that the line shape is dominated by inhomogeneous broadening. The inhomogeneous broadening of the 2−1 peak is larger than that for the 1−0 peak, which is reflected in the wider peak in the theoretical 2DIR spectrum. The reason for the broader 2−1 distribution is that the OH-stretch potential curve in ice is very anharmonic, and the anharmonicity increases with the hydrogen bond strength.40 Since the 1−0 and 2−1 frequencies are highly correlated, it also explains the steeper slope of the 2− 1 transition peak in our theoretical 2DIR. It is worth mentioning that the large breadth of the 2−1 transition peak has been discussed extensively in some experimental studies, and many explanations have been conjectured, such as asymmetric double well potential for the proton between two oxygen atoms,22 vibrational energy relaxation of the second excited state through some intermediate states,42 and the nonadiabatic coupling of the OH stretch with the quantal lowfrequency modes.24 In order to understand the sensitivity of the 2DIR spectrum to the amount of inhomogeneous broadening, it is convenient to derive some approximate expressions for the response functions. Basically, we make the Condon approximation, neglect the rotational dynamics of the chromophores, make the Gaussian approximation for the inhomogeneous broadening, make the cumulant approximation for the frequency fluctuations of each chromophore, and assume that the frequency time-correlation function for each chromophore is the same. The details are described in Appendix A. It turns out that these approximations all work well (compared to the “exact” theoretical results) for ice (although most of them do not work at all well for liquid water). To show this, in Figure 1, we plot the approximate result using eq A6 for the absorption spectrum (red dashed line), and in Figure 4, we plot the approximate result using eq A8 for the 2DIR spectrum for HOD/D2O ice Ih. In both cases, comparison with the exact theoretical results is excellent. Our approximate expressions involve the variances of the time-averaged frequency distribution for the 1−0 transition, σ12, the 2−1 transition, σ32, and the covariance, σ22, as defined in Appendix A. For the sake of argument and illustration, we can artifically decrease σ1, σ2, and σ3 to see the effect on the 1D and 2DIR spectra. We consider three cases: where these values are decreased by 1/3, decreased by 2/3, and set to 0. The results for the 1D spectrum are shown in Figure 5. All three of these new spectra are quite a bit too narrow and not consistent with experiment. The results for the 2DIR spectra are shown in Figure 6. Here we see that the result with the inhomogeneous broadening decreased by 2/3 in panel b is perhaps in best agreement with experiment. Indeed, this would give an inhomogeneous width of about 62/3 = 21 cm−1 for the 1−0 transition, consistent with the claim by Hamm and co-workers that the inhomogeneous broadening is less than 25 cm−1. This leaves the rather uncomfortable situation where our model results are consistent with the experimental 1D but not the 2D spectrum, whereas if we decrease the inhomogeneous

Figure 2. Calculated (from eqs 2−6) and experimental24 2DIR spectra for HOD/D2O ice Ih at 80 K (t2 = 200 fs). The red lobe is the 1−0 transition peak associated with ground state bleach and stimulated emission processes, and the blue lobe is the 2−1 transition peak associated with the excited state absorption process. The spectra are normalized to have the same 1−0 transition peak height.

In Figure 2, we also show our theoretical 2DIR spectrum (t2 = 200 fs, parallel polarization geometry) for HOD/D2O ice Ih at 80 K (left panel). In contrast to the experiment, both the 1− 0 and 2−1 transition peaks in our theoretical 2DIR spectrum are elongated along the diagonal (actually, along a line steeper than the diagonal for the 2−1 transition), indicative of substantial inhomogeneous broadening. The frequency gap between the 1−0 and 2−1 peaks (apparent anharmonicity of 217 cm−1) is smaller in our calculated 2DIR spectra compared to that in the experiment at this waiting time. Except for these issues, the agreement between the theoretical and experimental 2DIR for HOD/D2O ice Ih is reasonable. Each individual putative OH chromophore has a timeaveraged frequency determined by its proton-disordered environment, and fluctuations about this average frequency due to thermal motion (phonons). The distribution of timeaveraged frequencies represents the inhomogeneous broadening. For the 864 chromophores in our crystal at 80 K, this distribution is shown in Figure 3 for both the 1−0 and 2−1

Figure 3. Histograms of time-averaged transition frequencies ω̅ for HOD/D2O ice Ih at 80 K. The dashed lines are the best Gaussian fits. Red lines are for 1−0 transition frequency, and blue lines are for 2−1 transition frequency. 15539

dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544

The Journal of Physical Chemistry B

Article

broadening by 2/3 our results are consistent with the 2DIR but not the 1D experiments. If we consider, for the moment, reasons why our theoretical inhomogeneous broadening could be too large, the most likely candidate is the accuracy of our frequency map for ice. After all, the map was parametrized for liquid water. However, we have recently compared our map to DFT calculations from ice configurations, and while there is a small more or less uniform frequency shift from the ice results to the map, the slope of the map (versus electric field) and slope of the ice points are about the same, and therefore the amount of inhomogeneous broadening would be unchanged.40 To reiterate, the slope of the map reflects the dispersion in the frequencies for different electric fields, so two maps with the same slope would have the same width in the frequency distribution. Decreasing the inhomogeneous broadening by 2/3 is equivalant to decreasing the slope of the map by 2/3, which is incompatible with the DFT calculations.40 Ultrafast spectroscopy of liquid HOD/D2O has shown that the 3PEPS experiment actually leads to higher resolution dynamical information than does 2DIR (3PEPS shows the oscillation due to the hydrogen bond stretch, where 2DIR does not).46 In fact, at “long” times, the peak shift can be shown to be proportional to the frequency−frequency time-correlation function,47−50 which is the object of primary theoretical interest. For ice, the theoretical phonon contribution to the frequency−frequency time-correlation function (C1(t), see eq A5) is shown in Figure 7. It shows pronounced structure indeed, reflecting the various low-frequency modes of the system. In the vibrational 3PEPS measurement, a series of timeintegrated three-pulse photon-echo measurements with different waiting times are taken. In the impulsive limit, the integrated photon echo intensity is given by51

Figure 4. Calculated 2DIR spectra for HOD/D2O ice Ih at 80 K (t2 = 200 fs): (a) “exact” spectrum with the exact third-order response functions in eq 2, identical to the left panel in Figure 2; (b) approximate spectrum with the approximate third-order response functions in eq A8.

I(t 2 , t1) =

∫0



3

dt3 | ∑ R m(t3 , t 2 , t1)Lm(t3 , t 2 , t1)|2 m=1

(7)

and the peak shift t1*(t2) is defined by50 Figure 5. Absorption line shapes with full inhomogeneous broadening (blue), reduced by 1/3 (green), reduced by 2/3 (red), and set to 0 (black).

∂I(t 2 , t1) ∂t1

=0 t1= t1*(t 2)

(8)

Figure 6. Calculated 2DIR spectra for HOD/D2O ice Ih at 80 K (t2 = 200 fs) with the approximate third-order response functions in eq A8, and inhomogeneous broadening due to the proton disorder (σa, a = 1, 2, 3): (a) set to 0; (b) reduced by 2/3; (c) reduced by 1/3; (d) kept in full, identical to Figure 4b. 15540

dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544

The Journal of Physical Chemistry B

Article

of affairs. As a possible experiment to probe low-frequency dynamics in ice, and also further study inhomogeneous broadening, we encourage experimenters to perform a 3PEPS experiment. To this end, we have presented a theoretical calculation of such. There are several possible reasons for the discrepancy between theoretical and experimental 2DIR spectra. On the theoretical side, we would like to mention the following concerns about our model: (1) The classical treatment of the lattice vibrations might be inadequate at 80 K, as the nuclear quantum effect of hydrogen (deuterium) is pronounced at low temperatures.53−56 However, a model calculation by Sakurai and Tanimura using a reduced-hierarchy equations-of-motion approach shows that 2DIR spectra calculated classically and quantum mechanically are similar except for the frequency shifts of the peaks when the bath modulation is fast,57 consistent with some preliminary results of ours using a harmonic quantum correction factor.58 (2) Both the E3B model and maps are parametrized for liquid water, so their transferability to ice needs further validation. In fact, our recent DFT calculations40 showed that the frequency maps we are using herein work reasonably well for ice, although we estimate that the actual anharmonicity for ice is 254 cm−1 (compared to the value of 217 cm−1 reported herein). (3) We freeze the intramolecular bends so that our model neglects the Fermi resonance between bending overtones and OH stretches, and bend−libration combination bands, which may affect the OH stretch line shape to some extent.59−61 Bakker and co-workers pointed out that the coupling between the 2−1 transition of the OH stretch and the overtone of the H−O−D bending could lead to a very short vibrational lifetime of the second excited state.42 (4) Our spectral calculations assume the impulsive limit, which is a good approximation when the pulse duration is much shorter than the characteristic period of the bath motion.62 The pulse duration in the 2DIR experiments for ice is about 65 fs,24 while the period corresponding to the major Fourier component of the calculated frequency time-correlation functions (about 255 cm−1) is about 130 fs. (5) The initial configuration of ice Ih in our simulation (with an almost zero net dipole) may not be a good representation of the real ice Ih structure. It is possible that we have to average over many initial configurations with different net dipoles to improve our results.63 However, our preliminary results with different initial configurations (some of which have large net dipoles) show very similar 2DIR spectra. Our configuration also does not have defects, which will exist in the experimental samples.7 (6) The experiments were done with 2.5% H2O in D2O,24 so there are still some residual vibrational couplings between OH chromophores. To further validate and improve our model, more systematic work is in progress in our group. We are applying our mixed quantum/classical method to the calculation of the pump− probe anisotropy for ice Ih by Timmer and Bakker.64 The quality of our E3B model in reproducing the low-frequency modes in ice Ih is being examined by calculating the inelastic incoherent neutron scattering and low-frequency IR spectra.65,66 A systematic evaluation of our DFT-based maps for the spectroscopy-related quantities was carried out for liquid water and ice Ih.40 With these, we hope we can provide a more robust model to simulate vibrational spectra for ice Ih in the future. Meanwhile, we hope this work will stimulate more experimental and theoretical work on this ubiquitous but mysterious solid water phase.

Figure 7. (a) Time-correlation function for the dynamic fluctuations; (b) calculated vibrational three-pulse photon echo peak shift (3PEPS) (from eqs 7 and 8) for HOD/D2O ice Ih at 80 K with the approximate third-order response functions in eq A8.

As the approximate third-order response functions (eq A8) perform very well in the calculation of the 2DIR spectrum (Figure 4), we again use them to calculate the 3PEPS signal. In order to get a smooth 3PEPS curve, we increase the resolution of the time-correlation functions to 0.1 fs using the cubic spline interpolation method.52 The calculated 3PEPS signal is shown in Figure 7b, and we see that it does show a similar structure as the time-correlation function itself in panel a. Another interesting feature of the theoretical 3PEPS signal is that it does not decay to 0 as t2 gets large. In fact, as mentioned above, we have previously shown50 that at long times the peak shift is related to the frequency−frequency time-correlation function, but for ice, the full time-correlation function does not decay to zero (as in eq A5). In Appendix B, we show that under certain conditions t1*(∞) = Tpσ12/(σ12 + C1(0))

(9)

where σ1 and C1(0) are defined in Appendix A and Tp is defined in Appendix B. The value of t1*(∞) from Figure 7 is about 7.5 fs, compared to 11.1 fs from above. In any case, since t1*(∞) depends on the amount of inhomogeneous broadening, σ1, it may well be an interesting quantity to measure.

IV. CONCLUDING REMARKS OH (OD) stretch vibrational spectroscopy has been used to study the microscopic structure of ice Ih for almost a century. However, the presence of a spectral signature of the proton disorder is still controversial, even with the power of nonlinear spectroscopy. In this work, we attempt to deepen our understanding of this controversy by calculating 2DIR spectra for HOD/D2O ice Ih with our mixed quantum/classical approach. Our calculated 2DIR line shape shows evident inhomogeneous broadening, which is not observed in the experimental line shape, while our calculated linear IR spectrum is in good agreement with the experiment at the same temperature. We explore the effect of inhomogeneous broadening due to the proton disorder on the 2DIR line shapes. If the amount of inhomogeneous broadening is reduced, the model becomes compatible with the 2DIR spectrum but incompatible with the 1DIR spectrum. Frankly, we are bothered by and do not understand this unsatisfactory state 15541

dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544

The Journal of Physical Chemistry B

Article

RLj (t ) = ⟨μ10j (0)μ10j (t )U10(t )⟩

V. NOTE ADDED AFTER INITIAL SUBMISSION After the manuscript was submitted, we learned that Hamm and co-workers have measured the anisotropic 2DIR spectrum SZZZZ − SZZYY (i.e., the difference between the 2DIR signals measured in the parallel and perpendicular polarization geometries) for HOD in D2O ice Ih at 80 K (unpublished), and it is plotted (with their permission and encouragement) in Figure 8 at t2 = 200 fs. This 2DIR line shape is cleaner than that

(A1)

where the subscript L indicates the linear response function and U10(t) = exp[−i∫ t0 ω10(τ) dτ]. The third-order RFs are already given in eq 2 in the main text. The first two approximations we make are the Condon approximation (i.e., the magnitudes of the transition dipoles are taken to be constants), and the harmonic approximation for their magnitudes (μ21 = (2)1/2μ10). For ice Ih, all water molecules are fully hydrogen bonded to four nearest neighbors. Therefore, the Condon approximation is expected to hold. Regarding the harmonic approximation for the transition dipoles, our DFT calculations show that the anharmonicity for the transition dipole of the OH chromophores is very small throughout the frequency range of interest.40 We can also neglect the orientational dynamics (changes in the directions of the transition dipoles), as the reorientation of the water molecules is hindered due to the rigid hydrogen bond network in ice Ih. With these three approximations, the orientationally averaged linear and third-order RFs become RL(t ) ∼ ⟨U10(t )⟩

(A2)

and * (t1)U10 * (t1 + t 2) R1(t3 , t 2 , t1) = R 2(t3 , t 2 , t1) ∼ ⟨U10 U10(t1 + t 2 + t3)⟩ * (t1)U21 * (t1 + t 2)U21(t1 + t 2 + t3)⟩ R3(t3 , t 2 , t1) ∼ −2⟨U10 * (t1 + t 2) R 4(t3 , t 2 , t1) = R 5(t3 , t 2 , t1) ∼ ⟨U10(t1)U10 U10(t1 + t 2 + t3)⟩

Figure 8. Experimental anisotropic 2DIR spectrum (SZZZZ − SZZYY) for HOD in D2O ice Ih at 80 K (courtesy of P. Hamm and co-workers).

* (t1 + t 2)U21(t1 + t 2 + t3)⟩ R 6(t3 , t 2 , t1) ∼ −2⟨U10(t1)U21

in the right panel of Figure 2, and now the inhomogeneous broadening (elongation along the diagonal) is evident. The 2DIR spectrum taken in the perpendicular polarization geometry (i.e., SZZYY) is not expected to differ much in shape from that for the parallel 2DIR spectrum (i.e., SZZZZ), as the vibrational anisotropy due to molecular rotation and energy transfer for isotope-diluted ice Ih is very small for short waiting times, but it will suppress the effects of the local heating on the 2DIR spectrum. Thus, it is advantageous to get structural information (e.g., inhomogeneous broadening due to proton disorder) from the anisotropic 2DIR spectrum for ice Ih. Although the new 2DIR spectrum shows evident inhomogeneous broadening, it is still not as pronounced as that in our calculated 2DIR spectrum (left panel in Figure 2). By comparing the diagonal and antidiagonal slices of the new spectrum to those in Figure 6, we find that panel c gives the best agreement with the anisotropic 2DIR spectrum. Therefore, a more reasonable estimate of the inhomogeneous broadening due to the proton disorder from the experimental 2DIR spectrum is about 42 cm−1 at 80 K, more consistent with but still less than the results from the hole-burning experiments at 230 K (50 cm−1), an estimate from the linear IR line-width at 86 K (58 cm−1), and our theoretical calculations at 80 K (62 cm−1).

(A3)

For each chromophore, we can write ω10(τ ) = ⟨ω10⟩ + δω̅10 + δω10(τ )

(A4)

The first term on the r.h.s. is the global average over all chromophores, the second term is the deviation of the timeaveraged frequency for a given chromophore from the global average, and the third term represents the dynamic fluctuations for each chromophore. Thus, the second term results from proton-disorder-induced static inhomogeneity, and the third term comes from phonon fluctuations. The time-averaged frequency for a given chromophore discussed before is the sum of the first two terms. A simple, and it turns out, accurate approximation is to take static fluctuations to be uncorrelated with the dynamic ones, and assume that the dynamic fluctuations are independent of chromophore. In this case, the frequency−frequency time-correlation function is C(τ ) = ⟨(ω10(τ ) − ⟨ω10⟩)(ω10(0) − ⟨ω10⟩)⟩ = ⟨δω̅10 2⟩ + ⟨δω10(τ )δω10(0)⟩ ≡ σ12 + C1(τ )

(A5)

The average in the first term on the r.h.s. is an ensemble average over the chromophores, while the average in the second term is a time average over the fluctuations. At this point, we can make the cumulant expansion, which is valid for the two terms in eq A4 for completely different reasons. The inhomogeneous broadening evidently (see Figure 3) is Gaussian because of the central limit theorem, while the



APPENDIX A In this appendix, we provide the details of the approximations invoked in the model calculations of the linear and third-order response functions (RF) for isotope-diluted ice Ih within our mixed quantum/classical approach. The linear RF is given by15 15542

dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544

The Journal of Physical Chemistry B



time-dependent fluctatuations are Gaussian because the lowfrequency vibrations are approximately harmonic. Therefore, we obtain21 (A6)

where g m (t ) =

∫0

t

dτ(t − τ )Cm(τ )

(A7)

(see below for m = 2, 3). Similarly, 2

2

R1(t3 , t 2 , t1) = R 2(t3 , t 2 , t1) ∼ ei⟨ω10⟩t1− i⟨ω10⟩t3e−(1/2)σ1 (t1− t3) e−g1(t1) + g1(t2) − g1(t3) − g1(t1+ t2) − g1(t2+ t3) + g1(t1+ t2+ t3) 2 2

2 2

2

R3(t3 , t 2 , t1) ∼ − 2ei⟨ω10⟩t1− i⟨ω21⟩t3e−(1/2)(σ1 t1 + σ3 t3 − 2σ2 t1t3) e−g1(t1) + g2(t2) − g3(t3) − g2(t1+ t2) − g2(t2+ t3) + g2(t1+ t2+ t3) 2

2

R 4(t3 , t 2 , t1) = R 5(t3 , t 2 , t1) ∼ e−i⟨ω10⟩t1− i⟨ω10⟩t3e−(1/2)σ1 (t1+ t3) e−g1(t1) − g1(t2) − g1(t3) + g1(t1+ t2) + g1(t2+ t3) − g1(t1+ t2+ t3) 2 2

2 2

2

R 6(t3 , t 2 , t1) ∼ − 2e−i⟨ω10⟩t1− i⟨ω21⟩t3e−(1/2)(σ1 t1 + σ3 t3 + 2σ2 t1t3) e−g1(t1) − g2(t2) − g3(t3) + g2(t1+ t2) + g2(t2+ t3) − g2(t1+ t2+ t3) (A8)

where σ2 = ⟨δω̅ 10δω̅ 21⟩, 2

σ32

= ⟨δω̅ 21 ⟩, 2

C2(t ) = ⟨δω21(t )δω10(0)⟩ C3(t ) = ⟨δω21(t )δω21(0)⟩

(A9)

and g2(t) and g3(t) are defined above.



APPENDIX B In a previous paper,50 we showed that under certain approximations (including C1(t) = C2(t) = C3(t) and σ12 = σ22 = σ32) at long times the peak shift is related to the full frequency−frequency time-correlation function C(t) by t1*(t 2) = TpC(t 2)/C(0)

(B1)

where Tp =

1

ζD(ζ /2)

πC(0) (1 − e−ζ

2

/4

)

2

(B2) 2

ζ = Δ/(C(0))1/2, D(z) = e−z ∫ 0z dx ex is Dawson’s integral, and Δ = ⟨ω10⟩ − ⟨ω21⟩ is the anharmonicity. In the case of ice, C(t) is given by eq A5, so C(0) = σ12 + C1(0). Since C1(∞) = 0, that means that t1*(∞) = Tpσ12/(σ12 + C1(0))

■ ■

REFERENCES

(1) Malenkov, G. Liquid Water and Ices: Understanding the Structure and Physical Properties. J. Phys.: Condens. Matter 2009, 21, 283101. (2) Bragg, W. H. The Crystal Structure of Ice. Proc. Phys. Soc., London 1921, 34, 98−103. (3) Pauling, L. The Structure and Entropy of Ice and of Other Crystals with Some Randomness of Atomic Arrangement. J. Am. Chem. Soc. 1935, 57, 2680−2684. (4) Wollan, E. O.; Davidson, W. L.; Shull, C. G. Neutron Diffraction Study of the Structure of Ice. Phys. Rev. 1949, 75, 1348−1352. (5) Bernal, J. D.; Fowler, R. H. A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl Ions. J. Chem. Phys. 1933, 1, 515−548. (6) Singer, S. J.; Knight, C. Hydrogen-Bond Topology and Proton Ordering in Ice and Water Clusters. Adv. Chem. Phys. 2012, 147, 1−74. (7) Petrenko, V. F.; Whitworth, R. W. Physics of Ice; Oxford University Press: Oxford, U.K., 1999. (8) Li, F.; Skinner, J. L. Infrared and Raman Line Shapes for Ice Ih. II. H2O and D2O. J. Chem. Phys. 2010, 133, 244504. (9) Rice, S.; Bergren, M.; Belch, A.; Nielson, N. A Theoretical Analysis of the Hydroxyl Stretching Spectra of Ice Ih, Liquid Water, and Amorphous Solid Water. J. Phys. Chem. 1983, 87, 4295−4308. (10) Li, F.; Skinner, J. L. Erratum: Infrared and Raman Line Shapes for Ice Ih. II. H2O and D2O [J. Chem. Phys. 133, 244504 (2010)]. J. Chem. Phys. 2011, 134, 099901. (11) Shi, L.; Gruenbaum, S. M.; Skinner, J. L. Interpretation of IR and Raman Line Shapes for H2O and D2O Ice Ih. J. Phys. Chem. B 2012, 116, 13821−13830. (12) Hornig, D. F.; White, H. F.; Reding, F. P. The Infrared Spectra of Crystalline H2O, D2O and HDO. Spectrochim. Acta, Part A 1958, 12, 338−349. (13) Haas, C.; Hornig, D. F. Inter- and Intramolecular Potentials and the Spectrum of Ice. J. Chem. Phys. 1960, 32, 1763−1769. (14) Bertie, J. E.; Whalley, E. Infrared Sepctra of Ices Ih and Ic in the Range 4000 to 350 cm−1. J. Chem. Phys. 1964, 40, 1637−1645. (15) Li, F.; Skinner, J. L. Infrared and Raman Line Shapes for Ice Ih. I. Dilute HOD in H2O and D2O. J. Chem. Phys. 2010, 132, 204505. (16) Bergren, M. S.; Schuh, D.; Sceats, M. G.; Rice, S. A. The OH Stretching Region Infrared Spectra of Low Density Amorphous Solid Water and Polycrystalline Ice Ih. J. Chem. Phys. 1978, 69, 3477−3482. (17) Sivakumar, T. C.; Rice, S. A.; Sceats, M. G. Raman Spectroscopic Studies of the OH Stretching Region. J. Chem. Phys. 1978, 69, 3468−3476. (18) Hsu, D.; Skinner, J. L. On the Thermal Broadening of ZeroPhonon Impurity Lines in Absorption and Fluorescence Spectra. J. Chem. Phys. 1984, 81, 1604−1613. (19) Hsu, D.; Skinner, J. L. Nonperturbative Theory of TemperatureDependent Optical Dephasing in Crystals. I. Acoustic or Optical Phonons. J. Chem. Phys. 1984, 81, 5471−5479. (20) Woutersen, S.; Emmerichs, U.; Nienhuys, H.-K.; Bakker, H. J. Anomalous Temperature Dependence of Vibrational Lifetimes in Water and Ice. Phys. Rev. Lett. 1998, 81, 1106−1109. (21) Mukamel, S. Principles of Nonlinear Optical Spectroscopy; Oxford: New York, 1995. (22) Seifert, G.; Weidlich, K.; Graener, H. Picosecond IR HoleBurning Spectroscopy on HDO Ice Ih. Phys. Rev. B 1997, 56, R14231−R14234. (23) Hamm, P.; Zanni, M. Concepts and Methods of 2D Infrared Spectroscopy; Cambridge University Press: Cambridge, U.K., 2011. (24) Perakis, F.; Widmer, S.; Hamm, P. Two-Dimensional Infrared Spectroscopy of Isotope-Diluted Ice Ih. J. Chem. Phys. 2011, 134, 204505. (25) Auer, B. M.; Kumar, R.; Schmidt, J. R.; Skinner, J. L. Hydrogen Bonding and Raman, IR, and 2DIR Spectroscopy of Dilute HOD in Liquid D2O. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 14215−14220. (26) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. Pronounced NonCondon Effects in the Ultrafast Infrared Spectroscopy of Water. J. Chem. Phys. 2005, 123, 044513.

2 2

RL(t ) ∼ e−i⟨ω10⟩t e−(1/2)σ1 t e−g1(t )

Article

(B3)

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.

ACKNOWLEDGMENTS This work was supported by NSF grant CHE-1058752. The authors are grateful to Dr. Fivos Perakis and Prof. Peter Hamm for sending their raw experimental data, and for the new spectrum shown in the “note added” and for helpful correspondence. 15543

dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544

The Journal of Physical Chemistry B

Article

(27) Hochstrasser, R. M. Two-Dimensional IR-Spectroscopy: Polarization Anisotropy Effects. Chem. Phys. 2001, 266, 273−284. (28) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (29) Schmidt, J. R.; Roberts, S. T.; Loparo, J. J.; Tokmakoff, A.; Fayer, M. D.; Skinner, J. L. Are Water Simulation Models Consistent with Steady-State and Ultrafast Vibrational Spectroscopy Experiments? Chem. Phys. 2007, 341, 143−157. (30) Hamm, P.; Hochstrasser, R. M. Structure and Dynamics of Proteins and Peptides: Femtosecond Two-Dimensional Infrared Spectroscopy. In Ultrafast Infrared and Raman Spectroscopy; Marcel Dekker: New York, 2001; pp 273−348. (31) Kumar, R.; Skinner, J. L. Water Simulation Model with Explicit Three-Molecule Interactions. J. Phys. Chem. B 2008, 112, 8311−8318. (32) Tainter, C. J.; Pieniazek, P. A.; Lin, Y.-S.; Skinner, J. L. Robust Three-Body Water Simulation Model. J. Chem. Phys. 2011, 134, 184501. (33) Hayward, J.; Reimers, J. Unit Cells for the Simulation of Hexagonal Ice. J. Chem. Phys. 1997, 106, 1518−1529. (34) van der Spoel, D.; Lindahl, E.; Hess, B.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Feenstra, K. A.; van Drunen, R.; Berendsen, H. J. C. Gromacs User Manual Version 3.3; www.gromacs.org, 2005. (35) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (36) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (37) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N· log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (38) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (39) Skinner, J. L.; Auer, B. M.; Lin, Y.-S. Vibrational Line Shapes and Spectral Diffusion in Liquid Water. Adv. Chem. Phys. 2009, 142, 59−103. (40) Gruenbaum, S. M.; Tainter, C. J.; Shi, L.; Ni, Y.; Skinner, J. L. Robustness of Frequency, Transition Dipole, and Coupling Maps for Water Vibrational Spectroscopy. J. Chem. Theory Comput. 2013, 9, 3109−3117. (41) Perakis, F.; Hamm, P. Private Communication. (42) Dokter, A.; Bakker, H. Transient Absorption of Vibrationally Excited Ice Ih. J. Chem. Phys. 2008, 128, 024502. (43) Eaves, J. D.; Loparo, J. J.; Fecko, C. J.; Roberts, S. T.; Tokmakoff, A.; Geissler, P. L. Hydrogen Bonds in Liquid Water are Broken Only Fleetingly. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13019−13022. (44) Asbury, J. B.; Steinel, T.; Stromberg, C.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L.; Fayer, M. D. Water Dynamics: Vibrational Echo Correlation Spectroscopy and Comparison to Molecular Dynamics Simulations. J. Phys. Chem. A 2004, 108, 1107−1119. (45) Fecko, C. J.; Loparo, J. J.; Roberts, S. T.; Tokmakoff, A. Local Hydrogen Bonding Dynamics and Collective Reorganization in Water: Ultrafast Infrared Spectroscopy of HOD/D2O. J. Chem. Phys. 2005, 122, 054506. (46) Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Ultrafast Hydrogen-Bond Dynamics in the Infrared Spectroscopy of Water. Science 2003, 301, 1698−1702. (47) de Boeij, W. P.; Pshenichnikov, M. S.; Wiersma, D. A. On the Relation between the Echo-Peak Shift and Brownian-Oscillator Correlation Function. Chem. Phys. Lett. 1996, 253, 53−60. (48) Cho, M.; Yu, J.-Y.; Joo, T.; Nagasawa, Y.; Passino, S. A.; Fleming, G. R. The Integrated Photon Echo and Solvation Dynamics. J. Phys. Chem. 1996, 100, 11944−11953.

(49) Joo, T.; Jia, Y.; Yu, J.; Lang, M. J.; Fleming, G. R. Third-Order Nonlinear Time Domain Probes of Solvation Dynamics. J. Chem. Phys. 1996, 104, 6089−6108. (50) Piryatinski, A.; Skinner, J. L. Determining Vibrational SolvationCorrelation Functions from Three-Pulse Infrared Photon Echoes. J. Phys. Chem. B 2002, 106, 8055−8063. (51) Li, S.; Schmidt, J. R.; Piryatinski, A.; Lawrence, C. P.; Skinner, J. L. Vibrational Spectral Diffusion of Azide in Water. J. Phys. Chem. B 2006, 110, 18933−18938. (52) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran; Cambridge University Press: Cambridge, U.K., 1992. (53) Herrero, C. P.; Ramírez, R. Isotope Effects in Ice Ih: A PathIntegral Simulation. J. Chem. Phys. 2011, 134, 094510. (54) Ramírez, R.; Herrero, C. P. Kinetic Energy of Protons in Ice Ih and Water: A Path Integral Study. Phys. Rev. B 2011, 84, 064130. (55) Waheed, Q.; Edholm, O. Quantum Corrections to Classical Molecular Dynamics Simulations of Water and Ice. J. Chem. Theory Comput. 2011, 7, 2903−2909. (56) Pamuk, B.; Soler, J. M.; Ramírez, R.; Herrero, C. P.; Stephens, P. W.; Allen, P. B.; Fernández-Serra, M.-V. Anomalous Nuclear Quantum Effects in Ice. Phys. Rev. Lett. 2012, 108, 193003. (57) Sakurai, A.; Tanimura, Y. Does ℏ Play a Role in Multidimensional Spectroscopy? Reduced Hierarchy Equations of Motion Approach to Molecular Vibrations. J. Phys. Chem. A 2011, 115, 4009−4022. (58) Lawrence, C. P.; Skinner, J. L. Quantum Corrections in Vibrational and Electronic Condensed Phase Spectroscopy: Line Shapes and Echoes. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6720−6725. (59) Ritzhaupt, G.; Devlin, J. P. Infrared Spectrum of D2O Vibrationally Decoupled in Glassy H2O. J. Chem. Phys. 1977, 67, 4779−4780. (60) Hernandez, J.; Uras, N.; Devlin, J. P. Molecular Bending Mode Frequencies of the Surface and Interior of Crystalline Ice. J. Chem. Phys. 1998, 108, 4525−4529. (61) Liu, H.; Wang, Y.; Bowman, J. M. Quantum Calculations of Intramolecular IR Spectra of Ice Models Using Ab Initio Potential and Dipole Moment Surfaces. J. Phys. Chem. Lett. 2012, 3, 3671−3676. (62) Mukamel, S. Femtosecond Optical Spectroscopy: A Direct Look at Elementary Chemical Events. Annu. Rev. Phys. Chem. 1990, 41, 647−681. (63) Lindberg, G. E.; Wang, F. Efficient Sampling of Ice Structures by Electrostatic Switching. J. Phys. Chem. B 2008, 112, 6436−6441. (64) Timmer, R. L. A.; Bakker, H. J. Vibrational Förster Transfer in Ice Ih. J. Phys. Chem. A 2010, 114, 4148−4155. (65) Li, J.-C.; Ross, D. K. Evidence for Two Kinds of Hydrogen Bond in Ice. Nature 1993, 365, 327−329. (66) Warren, S. G.; Brandt, R. E. Optical Constants of Ice from the Ultraviolet to the Microwave: a Revised Compilation. J. Geophys. Res. 2008, 113, D14220.

15544

dx.doi.org/10.1021/jp405860u | J. Phys. Chem. B 2013, 117, 15536−15544