Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Lumin. Author manuscript; available in PMC 2009 January 28.
Published in final edited form as:
J Lumin. 2004 May; 107(1-4): 271–286.
doi:  10.1016/j.jlumin.2003.12.035
PMCID: PMC2632596

Vibrational echo correlation spectroscopy probes of hydrogen bond dynamics in water and methanol


Multidimensional vibrational echo correlation spectroscopy with full phase resolution is used to measure hydrogen bond dynamics in water and methanol. The OD hydroxyl stretches of methanol-OD oligomers in CCl4 and HOD in H2O are studied using the shortest mid-IR pulses (<50 fs, <4 cycles of light) produced to date. The pulses have sufficient spectral bandwidth to span the very broad (>400 cm−1) spectrum of the 0–1 and 1–2 transitions. Hydrogen bond population dynamics are extricated with exceptional detail in MeOD oligomers because the different hydrogen bonded species are spectrally distinct. The experimental results along with detailed calculations indicate the strongest hydrogen bonds are selectively broken through a non-equilibrium relaxation pathway following vibrational relaxation of the hydroxyl stretch. The correlation spectra are also a sensitive probe of the fluctuations in water and provide a stringent test of water models that are widely used in simulations of aqueous systems. The analysis of the 2D band shapes demonstrates that different hydrogen bonded species are subject to distinct (wavelength dependent) ultrafast (~100 fs) local fluctuations and essentially identical slow (0.4 and ~2 ps) structural rearrangements. Observation of wavelength dependent dynamics demonstrates that standard theoretical approaches assuming Gaussian fluctuations cannot adequately describe water dynamics.

Keywords: Infrared, Ultrafast, Vibrational echo, Correlation spectroscopy, Multidimensional, Methanol, Hydrogen bond network, Hydrogen bond dynamics, Hydrogen bond breaking, water

1. Introduction

Hydrogen bonding liquids have generated a great deal of experimental [115] and theoretical [1525] study because of their importance as solvents in chemical and biological systems. The facile evolution of hydrogen bonded networks [26] plays an important role in determining their properties. Examination of the hydroxyl stretch of hydrogen bonding liquids provides information on hydrogen bonding because the stretch frequency is sensitive to the nature and extent of hydrogen bonding [18,19,2629]. Strong evidence for the correlation between the hydroxyl stretch frequency and the hydrogen bond strength in solids is obtained from correlating crystallographic data and spectroscopic data [28]. From the crystallographic data, the length of the hydrogen bond can be determined. A shorter length corresponds to a stronger bond. The same relationship also applies to liquids [18,19,26,27,29].

Linear absorption spectroscopy can provide information on the influence of the number and types of hydrogen bonds on the hydroxyl stretch frequency. However, it cannot be used to examine the structural evolution of hydrogen bonds. Hydrogen bonds are constantly being broken and formed [26]. Stronger hydrogen bonds become weaker and weaker hydrogen bonds become stronger because of the continual structural evolution of the hydrogen bonding networks. The evolution of hydrogen bond strengths causes hydroxyl stretch frequency spectral diffusion [18]. The inhomogeneous broadening of the hydroxyl stretch [27,12,15,17] prevents the spectral evolution from being observed in the absorption line. As a result, the linear absorption spectrum can only provide time independent information.

It is in the context of hydrogen bond dynamics that ultrafast infrared experiments gain their importance [115]. Most ultrafast infrared experiments have utilized pump–probe spectroscopy [1,2,914], which can examine hydroxyl population dynamics. Pump–probe experiments give vibrational lifetimes, hydrogen bond breaking and reformation rates, and some information on spectral diffusion [1,2,1012]. The dynamics of hydrogen bond evolution can be extracted with much greater detail using vibrational echo experiments [3038] because the inhomogeneous contribution to the line shape can be eliminated and the underlying dynamical line shapes can be observed. Because of the complexity of hydrogen bonding systems, analysis of different contributions to dynamical vibrational spectra is enhanced by the use of ultrafast multidimensional methods. Such methods can separate and narrow the dynamic spectral features [35,8,3438], permitting the observation of various contributions to the dynamics. Multidimensional techniques are beginning to be applied to hydrogen bonding systems [35,8]. Multidimensional stimulated vibrational echo methods are akin to multidimensional NMR experiments [39]. However, vibrational stimulated echoes operate on a sufficiently fast time scale (10 s of fs) to directly examine the time evolution of hydrogen bonding networks.

The linear absorption spectra of the hydrogen bonding liquids studied in this work, 10% methanol-OD oligomers in CCl4 solution and 5% HOD in H2O, are displayed in Fig. 1. Examples of the different hydrogen bonded species formed MeOD and H2O are displayed at the top of the figure. A MeOD that is both a hydrogen bond donor and acceptor (called δ), has a very broad hydroxyl stretch band (~150 cm−1) centered at ~2490 cm−1 (see Fig. 1(a)) [22,27,40,41]. (There is a shoulder on the red side of the δ band caused by a Fermi resonance with the overtone of the methyl rocking mode [42].) For a MeOD that is a hydrogen bond donor but not an acceptor (called γ), the frequency shifts to ~2600 cm−1 and the spectrum is significantly narrower than the δ band (~80 cm−1). A MeOD that is a hydrogen bond acceptor but not a donor (called β) absorbs still further to the blue at ~2690 cm−1, and the band is relatively narrow (~20 cm−1) [22,27,40,41]. Water also displays spectral shifts to the blue when the extent of hydrogen bonding decreases [18,19]. However, because water can have many more hydrogen bonds, there is not a simple one-to-one correspondence between the number of hydrogen bonds and the frequency [18,19]. Consequently, methanol oligomers in CCl4 are an important system for studying the dynamics of hydrogen bonded networks [35].

Fig. 1
(a) The linear absorption spectrum of a 10% solution of methanol-OD in CCl4. β’s are hydroxyls that are acceptors but not donors. γ’s are donors but not acceptors. δ’s are both donors and acceptors. These ...

Below, multidimensional stimulated vibrational echoes with full phase information are applied to the study of hydrogen bond dynamics in MeOD oligomers in CCl4 and HOD in H2O. To perform the experiments, we use the shortest mid-IR pulses produced to date (<50 fs or <4 cycles of light) [35]. Because of the very large bandwidth associated with the ultrashort pulses, it is possible to perform experiments on the entire broad hydroxyl stretching 0–1 and 1–2 bands of water and MeOD even though the combined ground state and first excited state bands are ~400 cm−1 wide. The experiments yield time dependent vibrational echo correlation spectra, which provide information about hydrogen bond dynamics, interactions, and structural evolution that cannot be obtained from other experiments.

2. Experimental procedures

The ultrashort IR pulses are generated using a Ti:Sapphire regeneratively amplified laser/OPA system. The output of the regenerative amplifier is 26 fs transform limited 23mJ pulses at 1 kHz repetition rate. These are used to pump a substantially modified Spectra Physics short pulse IR OPA. The output of the OPA is compressed to produce <50 fs virtually transform limited IR pulses as measured by collinear autocorrelation. Fig. 2 displays a collinear autocorrelation trace of the shortest IR pulses produced to date. The solid line displays the experimentally measured collinear autocorrelation trace and the dashed line displays a Gaussian fit with a 40 fs full width at half maximum (FWHM). The oscillations in the time trace result from the interference of the pulses with the period of the ~2700 cm−1 optical cycle. The spectrum of the pulse displayed in the inset (solid line) demonstrates that the IR pulse, which is centered at ~2700 cm−1, is transform limited with a spectral FWHM of 368 cm−1, determined by the Gaussian fit to the spectrum (dashed line). For the experiments, the pulses were centered at ~2500 cm−1 and the compression was readjusted to give transform limited pulses in the sample as measured by a sample that gave a purely non-resonant signal. The long-term stability is such that data were collected continuously for as long as 5 days.

Fig. 2
Collinear autocorrelation of shortest mid-IR pulse produced to date. The solid line is the experimental data and the dashed line represents a 40 fs FWHM Gaussian fit to the data. The oscillations result from the ~2700 cm−1 optical cycle of the ...

For the heterodyne detected multidimensional stimulated vibrational echoes, the IR beam is split into 5 beams. Three of the beams are the excitation beams for the stimulated vibrational echo. A fourth beam is the local oscillator (LO) used to heterodyne detect the vibrational echo signal. The fifth beam is used as a probe in pump–probe experiments. All of the beams that pass through the sample are optically identical and are compensated for group velocity dispersion (GVD) simultaneously. The vibrational echo signal combined with the LO or the probe beam is passed through a monochromator. Depending on the experiment, the heterodyne detected echo or the probe beam is detected by a 32 element MCT array. At each monochromator setting, the array detects 32 individual wavelengths. The data were processed using a computer that acted as a digital lock-in amplifier that was synchronized to an optical chopper.

Difference intensity measurements were made for both the vibrational echo and pump–probe experiments. For the first laser pulse, the array measured the local oscillator or probe beam and the corresponding signal field. For the next laser shot, the optical chopper blocked one of the beams to prevent the signal field from being generated and only the local oscillator or probe beam was measured. The difference between these two measurements was then normalized by the spectrum of the local oscillator or probe detected in the second shot to remove the effect of pixel-to-pixel variations in the array. Both the heterodyne detected vibrational echo and the pump–probe signal are quadratic in the intensity of the laser. A single element detector measured the overall intensity of the laser, and was used to do a second normalization of the signals to provide the proper quadratic normalization.

The samples, 10% MeOD in CCl4 or 5% HOD in H2O, were held in a sample cell of CaF2 flats with a spacing of 50 and 6 μm, respectively. The peak absorbance of the samples was <0.2. Such low absorbance is necessary to prevent serious distortions of the pulses as they propagate through the sample.

For the spectrally resolved pump–probe experiments, the change in transmission spectrum was collected as a function of the delay between the pump and probe pulses. The monochromator was stepped across the spectrum so that the arrays of wavelengths measured by the different array blocks matched up to form a continuous spectrum across the entire hydroxyl stretch band. The pump pulse was delayed in variable steps. The spectrum of the pump and probe pulses were identical. The spectrum of the probe was removed from the data as the data were collected. The data were not corrected for the spectrum of the pump; it is sufficiently broad to span the desired spectral region.

The phase-resolved, heterodyne detected, stimulated vibrational echo was spectrally resolved in a monochromator and measured as a function of the frequency variable, ωm, and two time variables, τ and Tw, which are defined as the time between the first and second radiation field–matter interactions and the second and third interactions, respectively. The measured signal is the absolute value squared of the sum of the vibrational echo electric field, S, and the local oscillator electric field, L: |L + S|2 = L2 + 2LS + S2. The L2 term is time-independent and the S2 is negligibly small; hence, neither contributes to the time dependence of the signal. The spectrum of the 2LS term is the ωm frequency axis. As the τ variable is scanned in 2 fs steps, the phase of the echo electric field is scanned relative to the fixed local oscillator electric field, resulting in an interferogram measured as a function of the τ variable for each frequency, ωm, that is detected through the monochromator. The interferogram contains the amplitude, sign, frequency, and phase of the echo electric field as it varies with τ. By numerical Fourier transformation, this interferogram is converted into the frequency variable ωτ. For each ωm frequency measured, a spectrum as a function of ωτ is calculated. In NMR, the ωτ and the ωm axes are generally referred to as the ω1 and the ω3 axes, respectively [39].

The interferogram contains both the absorptive and dispersive components of the vibrational echo signal. To obtain purely absorptive features, two sets of quantum pathways are measured independently by appropriate time ordering of the pulses in the experiment [43,44]. With pulses 1 and 2 at the time origin, pathway 1 or 2 is obtained by scanning pulse 1 or 2 to negative time, respectively. In principle, by adding the Fourier transforms of the interferograms from the two pathways, the dispersive component cancels leaving only the absorptive component [3,43,44]. The 2D vibrational echo correlation spectra are constructed by plotting the amplitude of the absorptive part of the stimulated vibrational echo as a function of both ωm and ωτ.

Lack of perfect knowledge of the timing of the pulses and consideration of chirp on the vibrational echo pulse requires a “phasing” procedure to be used [3]. The projection slice theorem [3,39,43,44] is employed to generate the absorptive 2D correlation spectrum. The projection of the absorptive 2D correlation spectrum onto the ωm axis is equivalent to the IR pump–probe spectrum recorded at the same Tw, as long as all the contributions to the stimulated vibrational echo are absorptive. Consequently, comparison of the projected 2D stimulated vibrational echo spectrum with the pump–probe spectrum permits the correct isolation of the absorptive vibrational echo correlation spectrum from the 2D spectrum obtained from the addition of the two quantum pathways.

It is possible to come relatively close to the correct correlation spectrum prior to the “phasing” procedure because the very short pulses permit their time origins to be known within a few fs. The frequency dependent phasing factor used to correct the 2D spectra has the form


Each term in Eq. (1) has a well-defined physical origin. SC, is the correlation spectrum. S1 and S2 are the spectra recorded for pathways 1 and 2, respectively. ΔτLO,E accounts for the lack of perfect knowledge of the time separation of the LO pulse and the vibrational echo pulse; the term involving Q accounts for linear chirp introduced into the echo pulse by the rear window of the sample cell; Δτ1,2 accounts for the lack of perfect knowledge of the time origins of excitation pulses 1 and 2; and the term involving C accounts for the linear chirp caused by propagation of the echo pulse through sample.

The uncertainty in the phasing parameters is minimized by several procedures. The uncertainty in the time variable Δτ1,2 is minimized in a sample that gives a purely non-resonant signal, yet is as similar to the real sample as possible. For the MeOD experiments, a solution of 10% MeOH in CCl4 was used as the non-resonant sample. For the HOD experiments, pure H2O was used. The delay of pulse 1 is scanned relative to the other two fixed pulses and the amplitude of the non-resonant signal is measured. The time origin is set so that the maximum signal occurs at zero delay. Then the delay of pulse 2 is scanned in the same way, and its time origin is set so that the maximum signal also occurs at zero delay. These two steps are iteratively repeated until the maximum signal occurs simultaneously at the zero delay of both pulses. The procedure results in an accuracy of ±2fs for the delay between pulses 1 and 2.

The chirp on the excitation pulses (represented by C) is minimized by measuring a frequency resolved optical grating in the non-resonant sample. Interaction of the first two pulses with the non-resonant sample induces an instantaneous polarization in the sample that only persists while the pulses are in the sample. The third pulse scatters off the resulting polarization grating. By frequency resolving the scattered pulse and scanning the delay of one of the pulses, we directly measure the relative delay of the different frequencies within the pulses. Variable amounts of CaF2 and Ge substrates are inserted into the beam to optimally compress the pulse so that the maximum of the polarization grating signal appears at the zero delay time at all frequencies. After optimally compressing the pulse, the signal maximum varies by <1 fs across the >400 cm−1 bandwidth of the pulses. The window chirp (represented by Q) is always ~50 fs2.

Finally, the uncertainty in the time variable ΔτLO,E is minimized by measuring the spectral interferogram between the vibrational echo signal emitted from the resonant sample and the local oscillator. The local oscillator pulse is intentionally stepped a few ps away from temporal coincidence with the third pulse (which is approximately when the vibrational echo is emitted from the sample). With the time delay between pulses 1 and 2 set to zero, the emitted vibrational echo signal is heterodyne detected with the local oscillator pulse. Inverse Fourier transformation of the resulting spectral interferogram converts the periodic modulation in frequency into the corresponding temporal delay. The local oscillator pulse is then stepped back toward the vibrational echo by the calculated time delay. The uncertainty in this measurement is typically about ±10 fs.

An additional procedure was developed and employed to produce the very high quality correlation spectra shown below. The projection slice theorem reduces a two dimensional entity to a one-dimensional entity. In many instances this is sufficient [43,44]. However, in order to unambiguously assign a correct correlation spectrum it is preferable to use information from both dimensions. Therefore, in the phasing procedure an additional constraint is applied for the ωτ dimension. We can use information from the absolute value correlation spectrum, which is the sum of the absolute value spectra of the pathway 1 (rephasing) and pathway 2 (non-rephasing) discussed above. The absolute value spectrum is independent of the phase factor and peaks at the same frequency along ωτ as does the purely absorptive spectrum. Consequently, the difference in peak positions of a trial phased absorptive spectrum and the absolute value spectrum for each ωm gives an additional criterion on the quality of the correlation spectrum. The correct correlation spectrum is the one that provides the best fit to the pump–probe spectrum and minimizes the difference in peak positions of absorptive and absolute value spectrum. This procedure holds for symmetric line shapes and hence can be applied to the data discussed here, where the dynamical lines can be very well fit to a Gaussian line shape and the small deviations are symmetric. Following phasing, the errors in the time origins are <100 × 10−18 s and the chirp across the entire spectrum is <100 × 10−18 s.

3. Results and discussion

In the discussion that follows, we first focus on the hydrogen bond population dynamics of MeOD. The report focusing on the structural evolution of the hydrogen bond network of MeOD will be presented elsewhere due to space limitations [45]. We then report the structural evolution of the hydrogen bond network in water observed through the spectral diffusion dynamics measured in the hydroxyl stretch band.

3.1. Correlation spectroscopy of MeOD

Fig. 3 displays some of the 2D correlation spectra for MeOD, four contour plots at Tw= 125 fs, 1.2, 1.8, and 5.0 ps. The maximum positive signal has been normalized to unity in each plot. The contours represent equal 10% graduations. Although the amplitudes of the data have been obscured by the normalization, the relative amplitudes between data sets are well defined and will be used in the analysis described below. The ωτ axis is the axis of the first radiation field interaction. The ωm axis is the axis of the emission of the echo pulse. If the frequency of the first interaction and the emission are the same, a peak will appear on the diagonal. If the frequency of the emission is different from the frequency of the first interaction, a peak will appear off diagonal. The positive going band on the diagonal in the 125 fs spectrum corresponds to the 0–1 transition of the δ band with some contribution from the 0–1 γ band on the blue end. The center frequencies of the δ and γ 0–1 transitions are indicated by horizontal and vertical lines. The negative going band below the main diagonal band arises from the 1–2 transition of the δ band. The 1–2 band is off diagonal because the emission is shifted to lower frequency by the vibrational anharmonicity (~150 cm−1).

Fig. 3
Correlation spectra for Tw = 125 fs, 1.2, 1.8, and 5 ps. The contours represent 10% steps. The vertical and horizontal lines mark the center frequencies of the δ and γ 0–1 transitions. The positive going band are the δ ...

As the Tw delay increases, the correlation spectra displayed in Fig. 3 exhibit changes that provide insight into the mechanism of hydrogen bond evolution. By Tw = 1.2 ps, the correlation spectrum has changed dramatically due to vibrational relaxation and hydrogen bond breaking. A good fraction of the initially produced excited state population has decayed to the ground state because the vibrational lifetime of the δ OD stretch is ~0.5 ps [1]. The magnitude of the δ 1–2 off diagonal peak reflects the δ excited state decay. By Tw = 1.8 ps, the δ excited state has completely decayed. Because the lifetime of the γ band (~1 ps) is longer than that of the δ band, the diagonal γ 0–1 band is uncovered as the δ band decays and red shifts. By Tw= 1.8 ps, the diagonal γ peak is almost gone.

As will be demonstrated below, when the δ OD excited state relaxes, hydrogen bonds break [1,9,46], producing photoproduct γ ODs. The resulting photoproduct γ peak, located above the δ band, is visible by Tw = 1.2 ps. The diagonal δ band progressively contracts to the red as Tw increases. The off diagonal photoproduct γ peak located above the δ 0–1 band grows in with increasing Tw delay. By Tw= 5.0 ps, the only remaining features are the preserved δ and the photoproduct γ bands. Comparison of the center of the preserved δ band at Tw = 5.0 ps and the center of the initially excited δ band (indicated by the horizontal line at 2490 cm−1) demonstrates that the preserved δ band is contracted to the red and shifted off diagonal along the ωm axis. For times longer than ~5 ps, the peaks do not change shape significantly, but slowly decay in magnitude on a time scale of 10’s of ps because of hydrogen bond recombination [1,9,46]. An important aspect of the 2D correlation spectrum is that the photoproduct γ is spectrally separated in the plane from the equilibrium γ band (see Fig. 3, Tw = 1.2 and 1.8 ps), in contrast to 1D spectroscopies.

The Tw dependence of the correlation spectra raises a number of questions. First, why is there a correlation spectrum at all for times long compared to the vibrational lifetime? Normally, the decay of the excited state into the ground state causes the vibrational echo signal to decay. However, there is both an excited state and ground state contribution to the signal [47]. Following vibrational relaxation, a few tens of percent of the initially excited δ’s break a hydrogen bond in ~200 fs, removing them as δ absorbers [1]. Therefore, the decay of the excited state does not completely eliminate the contribution to the signal from the ground state. The ground state signal remains for the hydrogen bond recombination time, which is 10 s of ps [1,9,46]. The detailed analysis of the off diagonal photoproduct γ peak will be discussed in a future publication [45]. Briefly, spectral diffusion dynamics displayed in the photoproduct γ peak demonstrate that the broken oligomers retain a detailed “memory” of the previously intact hydrogen bond network.

The fundamentally important question that will be addressed here is whether vibrational relaxation simply “heats” the hydrogen bond network, producing changes that would be the same as raising the equilibrium temperature, or is there another mechanism (already used in the qualitative discussion above) that gives rise to the time dependence of the bands shown in Fig. 3? To address this question, we have performed comprehensive calculations to model the correlation spectra. The four models will be described in detail in a future publication [5]. Here we give a very brief overview of the modeling procedure and discuss only the two models which gave the best results. The correlation spectra are fit using the sum of two-dimensional Gaussian band shapes for the δ and γ, 0–1 and 1–2 bands. For the Tw= 125 fs correlation spectrum, the widths and positions of the bands are set equal to the equilibrium δ and γ widths and positions as determined from linear spectroscopy of the sample and by determination of the anharmonicity of the hydroxyl stretch (150 cm−1, see Fig. 3, top panel). The Tw= 125 fs correlation spectrum is very well described by this model. The calculations use results of the fits to the 125 fs data and factors that are unique to each model to calculate the 5 ps data. The 1–2 bands do not contribute because 5 ps is much longer than the lifetime.

The experiments also give the ratios of the amplitudes of the 0–1 bands at 125 fs to those at 5 ps. From comparison of the fit to the Tw = 125 fs correlation spectrum with the fit to the 5 ps correlation spectrum using each model, the ratios of the amplitudes are determined for each model.

Model 1: following vibrational relaxation of the initially excited OD stretch, the deposited energy causes hydrogen bond weakening. Weakening of hydrogen bonds shifts their absorption to the blue (higher energy) [18,19,2629]. The δ and γ bleach responsible for the long lived vibrational echo signal is not filled in by the vibrational relaxation because the population is shifted to the blue. The center positions and widths of the preserved δ and γ bleach bands are described by their equilibrium spectra. The weakened δ and γ photoproducts have spectra that are shifted to the blue along the ωm axis. In the calculations, the magnitude of the blue shift is allowed to vary to give the best fit to the 5 ps data (R2 = 1.1, residuals squared). The correlation spectrum calculated from Model 1 is displayed in panel (a) of Fig. 4. The model does a reasonable job of reproducing aspects of the data (compare to lowest panel of Fig. 3). However, it does not have the correct shape. The preserved band (positive going) is elongated to the upper right because of the preserved γ portion of the spectrum. Even more significant, to obtain this best fit, the amplitudes of the bands are unphysical. To get the features in the correct positions along the ωm axis requires a very small shift (5 cm−1) of the weakened band produced by vibrational relaxation. Because of the small shift, the preserved δ bleach band strongly overlaps and nearly cancels the weakened photoproduct band. To obtain the correct ratio of the amplitude of the Tw = 5.0 ps data relative to the 125 fs data, Model 1 requires the magnitude of the preserved bleach δ band to be ~9 times larger than the magnitude of the initial δ bleach, which would require the concentration of initially excited MeOD s increasing by a factor of ~9 between Tw = 125 fs and 5.0 ps. Obviously this result is not possible. Model 1 cannot describe the data. Model 1 is essentially a thermal equilibration model. Detailed measurements of the pump–probe spectrum show that it takes [gt-or-equal, slanted]30 ps for the system to reach thermal equilibrium (at which time the pump–probe spectrum is the same as the temperature difference spectrum), again demonstrating that the mechanism that gives rise to the 5 ps data is not thermal equilibration.

Fig. 4
Simulated correlation spectra calculated from two models of hydrogen bond breaking. Model 2 provides a more quantitative and meaningful description of the data.

Model 2: Vibrational relaxation of initially excited OD stretches selectively breaks hydrogen bonds on the red side of the δ band and produces red shifted γ photoproducts. In the calculations, the spectral region of hydrogen bond breaking is allowed to vary to give the best fit to the 5 ps data (R2 = 0.13). The simulated correlation spectrum calculated from Model 2 is displayed in panel (b) of Fig. 4. This model provides an excellent description of the Tw = 5.0 ps experimental data. The shape of the spectrum is correct, and the R2 value is a factor of ~8 smaller for Model 2 compared to Model 1. In addition to having the correct shape and a better R2, Model 2 predicts a physically reasonable magnitude for the amplitudes of the bands at 5 ps relative to 125 fs. Model 2 gives the preserved δ bleach to ~60% of the initial δ bleach, that is, 60% of the initially excited δ MeOD s break a hydrogen bond, which is within a few tens of percent of the value determined by independent measurements of the quantum yield of hydrogen bond breaking in MeOD oligomers [1].

The results presented above demonstrate that the δ band is contracted to the red and the photoproduct γ band is generated because vibrational relaxation selectively breaks the stronger hydrogen bonds (the red side of the line). Breaking of the strongest hydrogen bonds is at first surprising. However, this non-equilibrium mechanism can occur if vibrational relaxation on the low energy side of the line populates different modes than vibrational relaxation on the high energy side of the line, and the modes populated through the low energy relaxation pathway are more effective at breaking hydrogen bonds. Significant support for the proposed origins of the non-equilibrium mechanism comes from the IR pump/Raman probe experiments of Iwaki and Dlott, who observed that the pathway for vibrational relaxation was different on the red and blue sides of the hydroxyl stretching band in OH methanol [48]. Excitation on the red side led to vibrational relaxation that produced substantially more relative population of hydroxyl bends compared to excitation on the blue side that put much more population in the CH bends and other modes [48]. The hydroxyl bends would be expected to be strongly coupled to the hydrogen bond in contrast to CH stretches. This necessary feature of the proposed mechanism is indeed observed experimentally.

Here we have presented the application of vibrational echo correlation spectroscopy to the study of hydrogen bond network dynamics in methanol oligomers. Because of the brevity of this presentation, the analysis focused on the non-equilibrium hydrogen bond breaking dynamics following vibrational relaxation that result in selective breaking of the stronger hydrogen bonds. The data also provide unique information on spectral diffusion that will be presented subsequently to describe hydrogen bond evolution prior to and following hydrogen bond breaking [45].

3.2. Correlation spectroscopy of water

3.2.1. Correlation spectra

As shown above, hydrogen bonds can be broken following vibrational relaxation. Water and as well as other hydrogen bonding liquids are known to break hydrogen bonds following excitation of the hydroxyl stretch [1,3,9,46,49]. We have performed detailed spectrally resolved IR pump–probe experiments on water and find that, in contrast to methanol, broken hydrogen bond photoproducts do not make a measurable contribution to the water spectrum for pump–probe delays <2 ps. Consequently, we only consider water correlation spectra measured at Tw delays <2 ps. Studying the OD stretch in water not only provides dynamical information on water itself, but with a vibrational lifetime of ~1.8 ps, rather than the ~0.7 ps lifetime of the OH stretch in D2O [50], the time window is more than doubled in which we can extract the dynamics of the hydrogen bonded network unperturbed by photoproducts.

Vibrational echo correlation spectra were taken for Tw delays between 100 fs and 1.6 ps. Fig. 5 displays correlation spectra for four of the Tw delays, 100, 400, 800 fs, and 1.6 ps. Two peaks are observed along the ωm axis corresponding to the 0–1 and the 1–2 transitions. The positive going peak, corresponding to the vibrational echo emitted at the 0–1 transition frequency, appears on the diagonal. The negative going peak, corresponding to the vibrational echo emitted at the 1–2 transition frequency, appears red shifted along the ωm axis by the anharmonicity. From a fit we find the anharmonicity of the OD stretch vibration to be 140±20 cm−1.

Fig. 5
The experimental vibrational echo correlation spectra of the OD stretch of HOD in H2O as a function of time Tw. The positive going peak arises from the 0–1 transition. The negative going peak arises from the 1–2 transition. As time progresses ...

Looking in particular at the 0–1 band, it is obvious that the correlation spectra undergo a dramatic change in shape from Tw= 100 fs to 1.6 ps. The width along the ωτ axis (horizontal axis) is the dynamical line width of the 0–1 transition. The width is caused by the very fast fluctuations of the hydrogen-bonded network. As Tw is increased, additional fluctuations come into play that cause the dynamical line to broaden further. Changes in the dynamical line width report the spectral diffusion dynamics [18,51]. Therefore the change in the shape of the correlation spectrum provides information on the hydrogen bond dynamics on different time scales.

3.2.2. MD simulations and comparison to data

The key quantity that connects the vibrational echo correlation spectra and simulations of water dynamics is the time correlation function (TCF) of the hydroxyl stretch frequency. The TCF is defined as C(t) = left angle bracketδω(t)δω(0)right angle bracket, δω is the change in the hydroxyl stretch frequency from its average value. The TCF describes the loss of correlation of the frequency of the ensemble of OD stretch oscillators as time progresses. As the structure of the hydrogen bonding network evolves, the frequency of a given OD oscillator changes. At sufficiently long time, a given oscillator samples all frequencies in the absorption spectrum, and the TCF decays to zero. The time dependence of the vibrational echo correlation spectra provides a good test for the TCF, while the linear absorption spectrum does not.

The methods we use to compare the dynamics predicted from MD simulation with the experimentally observed dynamics are (i) a particular water model simulation is implemented, and (ii) the TCF for a given model of water is generated from a method that involves electronic structure calculations and classical MD simulations [52]. (iii), we calculate the vibrational echo correlation spectra from the TCF using time dependent diagrammatic perturbation theory to obtain the full third order non-linear material response [37,47]. In addition to the TCFs obtained from the water models, we construct a phenomenological tri-exponential TCF and vary its parameters until the calculations resulting from (iii) reproduce the experimental correlation spectra.

The calculations from the two water models and the phenomenological TCF are compared to the data by comparing the dynamical line widths of the corresponding correlation spectra at the center of the 0–1 band (ωm = 2500 cm1). The dynamical line width is the projection of a line through the spectrum at a particular ωm onto the ωτ axis. This method makes it possible to construct a simple plot to compare the data and calculations and, more importantly, look for an ωm dependence of the hydrogen bond network dynamics, as will be discussed below.

TCFs for the OH stretch of HOD in D2O have been calculated previously for the TIP4P model of water [18]. Here calculations for TIP4P and SPC/E water models were performed using a new method [52]. The method is based on ab initio electronic structure calculations of clusters of molecules. First configurations of molecules are generated from a classical molecular dynamics simulation of the HOD/H2O system. Then a representative set of HOD(H2O)n clusters are extracted from the simulation. For each cluster ab initio calculations (using density functional theory) are performed for different values of the OD stretch coordinate, and in doing so, an anharmonic potential curve that leads to the 0–1 OD transition frequency for the HOD molecule in that cluster is generated. From looking at a hundred or so clusters, it was determined that there is a linear correlation between the OD frequency and the component of the electric field from the H2O molecules on the D atom in the direction of the OD bond vector. Assuming this same linear correlation holds for the full liquid, the normalized TCF becomes the normalized electric field time correlation function. The latter can be generated from a completely classical molecular dynamics simulation.

The TCFs obtained from MD simulation are exceedingly well described by a sum of three exponential functions, although a small oscillatory feature at early times, which does not change the observables, is not captured. We will henceforth refer to the tri-exponential analytical representation as the TCF. It has the form


The overall amplitude of Eq. (2) was adjusted so that the calculation of the linear absorption line shape was able to reproduce the full width at half maximum of the experimentally measured absorption spectrum [53]. Parameters for the two water model TCFs are summarized in Table 1 along with a phenomenological TCF we obtained by fitting the data as discussed below. All TCFs reproduce the linear line shape with equal accuracy, showing that the linear line shape is not a good test of water dynamics [53].

Table 1
TCF Parameters

Using the above procedures, we calculated the theoretical vibrational echo correlation spectra with the appropriately scaled TCF’s taking into account both the 0–1 and 1–2 transitions of the hydroxyl stretch. Fig. 6 shows the results for four Tws, 100, 400, 800 fs, and 1.6 ps for the SPC/E water model. Comparison of the experimental and MD correlation spectra in Figs. 5 and and66 demonstrate a qualitative, but not a quantitative agreement. The 0–1 peak appears on the diagonal and the 1–2 peak appears shifted off-diagonal by the anharmonicity. The 0–1 peak in the Tw = 100 fs correlation spectrum displays the characteristic elongation along the diagonal that demonstrates the presence of inhomogeneity in the hydroxyl stretch peak. As Tw increases to 1.6 ps, the shape of the 0–1 peak changes, reflecting the substantial increase in the dynamical line width.

Fig. 6
Simulated vibrational echo correlation spectra calculated from the TCF that is derived from the SPC/E model of water.

A quantitative comparison of the data and the results of the water model TCF calculations and the phenomenological TCF calculations are shown in Fig. 7. The main figure plots the dynamical line widths for ωm= 2500 cm−1 extracted from the data and from the calculations. The experimental data appear as diamonds with error bars. The horizontal line at 147 cm−1 is the long time asymptotic line width. This width is narrower than the absorption spectrum (FWHM is 162 cm−1) because the measured dynamical line width is the product of the absorption spectrum and the laser intensity spectrum raised to the power of 1.5. The pulses used in the experiment, which were ~45 fs, had a very broad spectrum, but their spectral width did reduce the measured dynamical line width. Calculations demonstrate that 6 fs pulses are required for the asymptotic value of the measured dynamical line width to be indistinguishable from the absorption spectrum. The top curve (dots) is from the TIP4P model and the next lowest curve (dashes) is from the SPC/E model. The solid line through the data comes from calculations using the fitted phenomenological TCF.

Fig. 7
Dynamical line widths for ωm = 2500 cm1 from the experimental correlation spectra (diamonds with error bars), the TIP4P (dots) and SPC/E (dashes) water models, and the TCF obtained from fitting the data (solid line through the data). ...

In this fit the parameters for a tri-exponential TCF were constrained to be consistent with other complementary experiments and calculations. The magnitude of the TCF is determined by the linear absorption spectrum. The time constant of the fast component was set at 32 fs to be consistent with the fast components of the TCF’s obtained by MD simulation. This time scale is also consistent with recent work by Elsaesser and co-workers who found that a 30 fs component to the TCF provided the best fit to their 2 pulse vibrational echo measurement of water [7]. The time constant for the intermediate component was set equal to the intermediate decay component (400 fs) in a vibrational echo peak shift measurement we performed in the course of these experiments. This can be done since the vibrational echo peak shift has been demonstrated to reproduce the TCF [5456] given a sufficient waiting time, Tw, so that the value of the normalized TCF[double less-than sign]1 [55]. Therefore, in the fitting procedure, three parameters were allowed to vary independently, the pre-exponential factors of two of the components and the time constant of the slowest component.

We find that the dynamical line widths derived from the water models do not reproduce the data within experimental error. The line width from the TIP4P model reaches its asymptotic value already by Tw= 1.6 ps and the SPC/E line width is close to that asymptotic value at that time. However, the data and the line widths obtained from the fit to the data are still removed from the asymptotic value, approaching it relatively slowly.

The time dependence of the model line shapes is correlated to the time dependence of the model TCFs. The inset in Fig. 7 compares model and phenomenological TCFs directly. It is evident that both model TCFs have too much amplitude at short time and the long time components are too fast (see also Table 1).

Of particular interest here is the long time component of the TCF because the time constant and amplitude of this component are in substantial error when calculated from the water models. The SPC/E and TIP4P model give a time constant of the slow component of the TCF of 0.98 and 0.90 ps, respectively, and both have too little amplitude in the slow component. The best least squares fit through the experimental data gives a 1.8 ps time constant. We estimate upper and lower bounds of ±0.5 ps based on a factor of 2 increase in the sum of the squares of the residuals.

Comparison of the phenomenological TCF with those derived from MD simulation demonstrates the TIP4P and SPC/E models of water describe well the fast fluctuations in water but not the slower fluctuations. The TIP4P model essentially misses the slow component. The SPC/E model does better, but over-estimates the magnitude of the intermediate and under-estimates the magnitude of the slow components. In addition, the SPC/E model underestimates the time scale of the slow component by a factor of ~2. Two very different methods were used to calculate the TCF for TIP4P, the one described here [52] and the one used previously [18], but the results are virtually identical. This demonstrates that the failure of the MD simulations to reproduce the data is in the water models themselves, rather than the method used to calculate the OD oscillator frequency fluctuations.

Analysis of the TIP4P model for the OH stretch of the HOD/D2O system indicates that the fast component of the TCF arises from hindered translational motion of the hydrogen bond (stretching of the hydrogen bond length coordinate) [18]. The slower component of the TCF calculated for the HOD/D2O system has been assigned to the influence of hydrogen bond formation and breaking (hydrogen bond equilibration) [18,19,25]. Our analysis demonstrates this time scale (~0.5 ps) is too fast to describe the slow dynamics that we observe in water. Some MD simulations report hydrogen bond breaking times of ~2 ps [24,57]. This time scale corresponds with the slow component of the phenomenological TCF we obtained from the data. Nonetheless, it still is plausible that the slow time scale evolution found in our experiments is due to hydrogen bond equilibration, that is, the breaking and formation of hydrogen bonds.

3.2.3. Wavelength dependent dynamics

Inherent in the TCF analysis presented above was the underlying assumption that the dynamics of water molecules (as interrogated by the OD stretch) are the same at all frequencies within the broad hydroxyl stretching band. A water molecule can make zero to four hydrogen bonds. These species have different transition frequencies, although they are not resolvable within the broad hydroxyl stretch absorption band [16,18]. Furthermore, the strength of a hydrogen bond determines the hydroxyl stretch transition frequency to a great extent [28]. On long time scales, hydrogen bonds will break and form, and strong hydrogen bonds will weaken and weak hydrogen bonds will strengthen. Therefore, it is expected that the dynamics will be independent of frequency as water molecules sample all conformations on long time scales. However, the different hydrogen bonded conformations persist on short time scales. As we show below, hydrogen bond dynamics are frequency dependent, which demonstrates the different hydrogen bonded conformations of water have different dynamics.

The dynamical line widths were determined from the full correlation spectrum at a number of frequencies, ωm, and Tw delays. A procedure, which will be described subsequently [53], was used to determine if the 1–2 band was interfering with the determination of the dynamic line widths of the 0–1 band. It was found that for frequencies [gt-or-equal, slanted]2500 cm−1, the 1–2 band causes no significant perturbation. Some of the results are shown in Fig. 8. The 100 fs data and the line fit to it demonstrate very clearly that the dynamical line widths, and therefore, the hydrogen bond dynamics vary as a function of frequency. As the frequency increases, the dynamical line width increases. The increase in width can be seen in the 100 fs correlation spectrum shown in Fig. 5. The spectrum is clearly broader for higher ωm. However, the 400 fs data and line in Fig. 8 show that by 400 fs, within experimental error, there is no longer frequency dependence to the dynamical line width. The inset in Fig. 8 displays the change in the dynamical line width with ωm as a function of Tw. The points are the slopes of lines like the 100 fs line in the body of Fig. 8. The line through the data yields the time scale for the decay of the inhomogeneity in the dynamics as ~100 fs.

Fig. 8
Dependence of the dynamical line width Γdyn on ωm for two Tws, 100 and 400 fs. For 100 fs there is a distinct dependence of width on ωm. By 400 fs there is no dependence of width on ωm. Inset: the decay of the inhomogeneity ...

The fundamentally important point illustrated by Fig. 8 is that the dynamics of water are different when observed at different wavelengths at short times. Specifically, there is slower spectral diffusion at lower frequencies on short time scales. Recent MD simulations that find the hydroxyl stretch band is composed of differently hydrogen-bonded species may contain the explanation [19,25]. In particular, molecules absorbing at the blue side of the line tend to have weaker hydrogen bonds, and so the local environments would fluctuate faster. If a Gaussian process properly describes the spectral diffusion, then the dynamical line width would be independent of the frequency. Therefore, the observed frequency dependence clearly demonstrates a breakdown of the Gaussian approximation. The wavelength dependent spectral diffusion shows that dynamics in water are more complex than previously assumed and that the development of new theoretical approaches that go beyond the Gaussian approximation/cumulant truncation [17,58] will be necessary to fully understand water.

4. Concluding remarks

The results of ultrafast heterodyne detected multidimensional stimulated vibrational echo correlation spectroscopy with full phase information were employed to study the dynamics of hydrogen bond networks in MeOD and water (H2O). By using extremely short pulses (<50 fs), it is possible to excite the entire very broad hydroxyl-stretching band including the 1–2 transition. The multidimensional correlation spectroscopy makes it possible to explicate the hydrogen bond population dynamics in great detail. Spectral diffusion dynamics of the hydroxyl stretch 0–1 transition of water are elucidated with great sensitivity in a comprehensive manner.

The vibrational echo correlation spectra of MeOD display the hydrogen bond population dynamics in two dimensions, which allows the mechanism of hydrogen bond evolution to be elucidated with exceptional detail. At short time, the correlation spectrum (see Fig. 3) is dominated by the hydrogen bonded OD stretch 0–1 transition (positive going peak) and the 1–2 transition (negative going peak). As Tw is increased the broad diagonal δ band contracts to the red side of the line and shifts off of the diagonal. Detailed calculations used to model the experimental correlation spectra demonstrate that the hydroxyl’s with absorption on the red side of the δ band (strong hydrogen bonds) selectively break hydrogen bonds following vibrational energy relaxation. It is proposed that the strong hydrogen bonds (red side of the hydroxyl stretch band) are selectively broken because the low frequency modes that are excited by vibrational energy relaxation of the hydroxyl stretch differ from those that are excited by vibrational energy relaxation of the higher frequency hydroxyl stretches (weaker hydrogen bonds).

As discussed above, the hydrogen bond network in water can be examined for times <2 ps without complication from hydrogen bond breaking. The vibrational echo correlation spectra of water (see Fig. 5) display clear broadening along the ωτ (horizontal) axis as Tw is increased. A cut across the correlation spectrum through the center of the 0–1 peak (2500 cm−1) reports the dynamical line width for each value of Tw. The dynamical line width of the 0–1 peak is a sensitive spectroscopic probe of the frequency fluctuations in water. From the analysis of the dynamical line widths, it was found that the frequency fluctuations in water occur on multiple time scales out to at least 2 ps. This result is in contrast with dynamics of water calculated using molecular dynamics simulations of two water models, TIP4P and SPC/E. While the SPC/E model is somewhat better, both water models produce dynamical lines that broaden too rapidly. The spectral diffusion dynamics in water were shown to be frequency dependent and are more complex than previously assumed. New theoretical approaches that go beyond the Gaussian approximation/cumulant truncation [17,58] will be necessary to fully understand water.


We would like to thank Professor James L. Skinner, Department of Chemistry, University of Wisconsin, Madison, for the molecular dynamics simulations of water. We would like to thank Professor Peter Hamm, Physikalisch Chemisches Institut, Universität Zürich, for help with the diagrammatic perturbation theory calculations. This work was supported by the AFOSR (F49620-01-1-0018), DOE (DE-FG03-84ER13251), and NSF (DMR-0332692). TS thanks the Emmy Noether program of the DFG for partial support.


1. Gaffney KJ, Davis PH, Piletic IR, Levinger NE, Fayer MD. J Phys Chem A. 2002;106:12012.
2. Piletic IR, Gaffney KJ, Fayer MD. J Chem Phys. 2003;119:423.
3. Asbury JB, Steinel T, Stromberg C, Gaffney KJ, Piletic IR, Goun A, Fayer MD. Chem Phys Lett. 2003;374:362.
4. Asbury JB, Steinel T, Stromberg C, Gaffney KJ, Piletic IR, Goun A, Fayer MD. Phys Rev Lett. 2003;91:237402. [PubMed]
5. Asbury JB, Steinel T, Stromberg C, Gaffney KJ, Piletic IR, Fayer MD. J Chem Phys. 2003;119:12981.
6. Stenger J, Madsen D, Hamm P, Nibbering ETJ, Elsaesser T. J Phys Chem A. 2002;106:2341.
7. Stenger J, Madsen D, Hamm P, Nibbering ETJ, Elsaesser T. Phys Rev Lett. 2001;87:027401.
8. Yeremenko S, Pshenichnikov MS, Wiersma DA. Chem Phys Lett. 2003;369:107.
9. Graener H, Ye TQ, Laubereau A. J Chem Phys. 1989;90:3413.
10. Laenen R, Rausch C, Laubereau A. Phys Rev Lett. 1998;80:2622.
11. Bratos S, Gale GM, Gallot G, Hache F, Lascoux N, Leicknam JC. Phys Rev E. 2000;61:5211. [PubMed]
12. Gale GM, Gallot G, Hache F, Lascoux N, Bratos S, Leicknam JC. Phys Rev Lett. 1999;82:1068.
13. Bakker HJ, Neinhuys HK, Gallot G, Lascoux N, Gale GM, Leicknam JC, Bratos S. J Chem Phys. 2002;116:2592.
14. Kropman MF, Nienhuys HK, Woutersen S, Bakker HJ. J Phys Chem A. 2001;105:4622.
15. Fecko CJ, Eaves JD, Loparo JJ, Tokmakoff A, Geissler PL. Science. 2003;301:1698. [PubMed]
16. Lawrence CP, Skinner JL. J Chem Phys. 2002;117:5827.
17. Lawrence CP, Skinner JL. J Chem Phys. 2002;117:8847.
18. Lawrence CP, Skinner JL. J Chem Phys. 2003;118:264.
19. Lawrence CP, Skinner JL. Chem Phys Lett. 2003;369:472.
20. Haughney M, Ferrario M, McDonald IR. J Phys Chem. 1987;91:4934.
21. Matsumoto M, Gubbins KE. J Chem Phys. 1990;93:1981.
22. Staib A. J Chem Phys. 1998;108:4554.
23. zum Buschenfelde DM, Staib A. Chem Phys. 1998;236:253.
24. Luzar A. J Chem Phys. 2000;113:10663.
25. Rey R, Møller KB, Hynes JT. J Phys Chem A. 2002;106:11993.
26. Pimentel GC, McClellan AL. The Hydrogen Bond. W. H. Freeman and Co; San Francisco: 1960.
27. Liddel U, Becker ED. Spectrochim Acta. 1957;10:70.
28. Novak A. In: Structure and Bonding. Dunitz JD, editor. Springer; Berlin: 1974. p. 177.
29. Falk M, Ford TA. Canadian J Chem. 1966;44:1699.
30. Zimdars D, Tokmakoff A, Chen S, Greenfield SR, Fayer MD, Smith TI, Schwettman HA. Phys Rev Lett. 1993;70:2718. [PubMed]
31. Tokmakoff A, Fayer MD. J Chem Phys. 1995;102:2810.
32. Rector KD, Engholm JR, Rella CW, Hill JR, Dlott DD, Fayer MD. J Phys Chem A. 1999;103:2381.
33. Hamm P, Lim M, Hochstrasser RM. Phys Rev Lett. 1998;81:5326.
34. Zanni MT, Asplund MC, Hochstrasser RM. J Chem Phys. 2001;114:4579.
35. Merchant KA, Thompson DE, Fayer MD. Phys Rev Lett. 2001;86:3899. [PubMed]
36. Golonzka O, Khalil M, Demirdoven N, Tokmakoff A. Phys Rev Lett. 2001;86:2154. [PubMed]
37. Mukamel S. Ann Rev Phys Chem. 2000;51:691. [PubMed]
38. Venkatramani R, Mukamel S. J Chem Phys. 2002;117:11089.
39. Ernst RR, Bodenhausen G, Wokaun A. Nuclear Magnetic Resonance in One and Two Dimensions. Oxford University Press; Oxford: 1987.
40. Symons MCR, Thomas VK. J Chem Soc Faraday Trans. 1981;1(77):1883.
41. Kristiansson O. J Mol Struct. 1999;477:105.
42. Bertie JE, Zhang SL. J Mol Struct. 1997;413–414:333.
43. Khalil M, Demirdöven N, Tokmakoff A. Phys Rev Lett. 2003;90:047401. [PubMed]
44. Khalil M, Demirdöven N, Tokmakoff A. J Phys Chem A. 2003;107:5258.
45. Asbury JB, Steinel T, Fayer MD. J Phys Chem A. Submitted.
46. Laenen R, Rausch C. J Chem Phys. 1997;106:8974.
47. Mukamel S. Principles of Nonlinear Optical Spectroscopy. Oxford University Press; New York: 1995.
48. Iwaki LK, Dlott DD. J Phys Chem A. 2000;104:9101.
49. Lock AJ, Woutersen S, Bakker HJ. J Phys Chem A. 2001;105:1238.
50. Woutersen S, Emmerichs U, Bakker HJ. Science. 1997;278:658.
51. Berg MA, Rector KD, Fayer MD. J Chem Phys. 2000;113:3233.
52. Corcelli SA, Lawrence CP, Skinner JL. accepted.
53. Asbury JB, Steinel T, Corcelli SA, Lawrence CP, Skinner JL, Fayer MD. J Phys Chem A. in press.
54. Joo TH, Jia YW, Yu JY, Jonas DM, Fleming GR. J Chem Phys. 1996;104:6089.
55. Piryatinski A, Skinner JL. J Phys Chem B. 2002;106:8055.
56. de Boeij W, Pshenichnikov MS, Wiersma DA. Chem Phys Lett. 1996;253:53.
57. Marti J, Padro JA, Guardia E. J Chem Phys. 1996;105:639.
58. Piryatinski A, Lawrence CP, Skinner JL. J Chem Phys. 2003;118:9672.