Resonant vibrational energy transfer in ice Ih L. Shi, F. Li, and J. L. Skinner Citation: The Journal of Chemical Physics 140, 244503 (2014); doi: 10.1063/1.4883913 View online: http://dx.doi.org/10.1063/1.4883913 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/24?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Resolving anisotropic distributions of correlated vibrational motion in protein hydration water J. Chem. Phys. 141, 22D509 (2014); 10.1063/1.4896073 Vibrational energy relaxation of large-amplitude vibrations in liquids J. Chem. Phys. 137, 024506 (2012); 10.1063/1.4733392 Polarized Raman spectroscopic study on the solvent state of glassy LiCl aqueous solutions and the state of relaxed high-density amorphous ices J. Chem. Phys. 134, 244511 (2011); 10.1063/1.3603965 Water structure, dynamics, and vibrational spectroscopy in sodium bromide solutions J. Chem. Phys. 131, 144511 (2009); 10.1063/1.3242083 Vibrational spectroscopy and dynamics of small anions in ionic liquid solutions J. Chem. Phys. 123, 084504 (2005); 10.1063/1.2000229
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
THE JOURNAL OF CHEMICAL PHYSICS 140, 244503 (2014)
Resonant vibrational energy transfer in ice Ih L. Shi, F. Li, and J. L. Skinner Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin, Madison, Wisconsin 53706, USA
(Received 4 April 2014; accepted 5 June 2014; published online 24 June 2014) Fascinating anisotropy decay experiments have recently been performed on H2 O ice Ih by Timmer and Bakker [R. L. A. Timmer, and H. J. Bakker, J. Phys. Chem. A 114, 4148 (2010)]. The very fast decay (on the order of 100 fs) is indicative of resonant energy transfer between OH stretches on different molecules. Isotope dilution experiments with deuterium show a dramatic dependence on the hydrogen mole fraction, which confirms the energy transfer picture. Timmer and Bakker have interpreted the experiments with a Förster incoherent hopping model, finding that energy transfer within the first solvation shell dominates the relaxation process. We have developed a microscopic theory of vibrational spectroscopy of water and ice, and herein we use this theory to calculate the anisotropy decay in ice as a function of hydrogen mole fraction. We obtain very good agreement with experiment. Interpretation of our results shows that four nearest-neighbor acceptors dominate the energy transfer, and that while the incoherent hopping picture is qualitatively correct, vibrational energy transport is partially coherent on the relevant timescale. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4883913] I. INTRODUCTION
Resonant energy transfer between an excited donor chromophore and an acceptor chromophore has been studied extensively by both experiment and theory for over 50 years.1–8 Most research involves electronic excitation, although due to the recent development of ultrafast infrared spectroscopy, researchers now have the opportunity to examine resonant vibrational energy transfer as well.9–18 Many studies have focused on the OH stretch in neat liquid water,13, 14, 16–18 in which many near-degenerate OH chromophores interact through intramolecular or intermolecular vibrational couplings, and their vibrational frequencies are modulated by other low-frequency translational and rotational modes. If OH chromophores are initially vibrationally excited by interacting with a polarized infrared pulse, the created anisotropy can decay through either molecular rotation or resonant energy transfer to other OH chromophores. The anisotropy decay can be detected experimentally by polarization-resolved pump-probe13, 14, 18 or two-dimensional infrared (2DIR) spectroscopy.16, 17 In neat liquid H2 O, the anisotropy decay occurs within 200 or 300 fs and is dominated by resonant energy transfer. Many theoretical methods have been proposed for understanding resonant energy transfer, which is often characterized as being incoherent or coherent.19–22 The incoherent or hopping picture occurs theoretically when the coupling between chromophores is treated perturbatively, and transport can be described by a master equation in the local site basis.20, 21 If the interaction between chromophores is dipoledipole, for example, this leads to the Förster theory.1, 21, 23, 24 Coherent transport occurs in the weak system-bath coupling limit, and can be described by Redfield25 or Lindblad26, 27 equations. There are also many theoretical works going beyond the Markovian approximation or the standard pertur0021-9606/2014/140(24)/244503/6/$30.00
bative treatment, and some of them treat incoherent/coherent energy transfer in a unified manner.28–35 In interpreting their pump-probe experimental results, Bakker and co-workers adapted the Förster theory to describe the ultrafast resonant energy transfer in liquid water, and found that it fit their experimental data very well.13, 14 Skinner and co-workers have developed a mixed quantum/classical model for OH-stretch spectroscopy in water, and have used it to calculate FTIR, Raman, 2DIR, SFG, and 2DSFG spectra.36–38 The exciton Hamiltonian for the OHstretch vibrations is written in the local-mode basis, and so the diagonal elements are the local-mode (anharmonic) vibrational frequencies, while the off-diagonal elements are the couplings (both intramolecular and intermolecular) between OH-stretch chromophores. The exciton Hamiltonian fluctuates as the classical degrees of freedom evolve in time. Depending on the statistical properties of this fluctuating Hamiltonian, such a model will produce coherent or incoherent transport (or something in between). In fact, application of this model by Yang, Li, and Skinner39 showed that in liquid water the incoherent hopping picture is appropriate, although the interpretation of experiments in terms of the usual Förster theory is not so straightforward. Timmer and Bakker recently extended their anisotropy decay experiments to ice Ih,15 finding a somewhat faster decay on the 100 fs timescale. They also studied what happens when the hydrogen atoms in H2 O ice are replaced by deuterium. That is, when D2 O is added to H2 O, one obtains a mixture of H2 O/HOD/D2 O, and to a good approximation one can model this system by assuming that the probability that each hydrogen atom is actually hydrogen (instead of deuterium) is simply equal to the hydrogen mole fraction. The consequences for energy transfer are dramatic, since the OD stretches are so far off resonance they do not function as energy acceptors. Therefore, as the hydrogen mole fraction
140, 244503-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
244503-2
Shi, Li, and Skinner
J. Chem. Phys. 140, 244503 (2014)
decreases the anisotropy decay slows down dramatically. These beautiful experiments directly confirm the essential picture of energy-transfer-induced anisotropy decay. Timmer and Bakker interpreted their experiments in terms of the Förster theory of incoherent hopping. Within their model the rate of intramolecular hopping was by far the largest (since the distance between the two dipoles is the closest), followed by intermolecular hopping to four near neighbors. They concluded that energy transfer outside the first solvation shell was negligible. Two theoretical papers on this problem, by Poulsen, Nyman, and co-workers actually preceded the experiments.40, 41 They used clusters of 15 water molecules to represent ice Ih, and ran wave-packet dynamics to monitor the survival probability of the excitiation on the initially excited mode.40 This led them to conclude that the resonant OH vibrational energy transfer in ice Ih occurs within 100 fs, in agreement with the recent pump-probe experiment.15 However, they did not directly calculate the anisotropy observable. The goal of this paper is to apply our mixed quantum/classical model to calculate the anisotropy decay in ice Ih, for both neat H2 O and for the isotopic mixtures. This provides a stringent test of the model since there are no free parameters. We have previously used the same model to calculate IR, Raman, and 2DIR spectra of ice Ih.42–46 We find, in fact, that the model performs quite well for the anisotropy decay, for all hydrogen mole fractions. We then analyze our results to examine the relative contribution from intramolecular and intermolecular energy transfer, from first and second solvation shells, and the appropriateness of the incoherent hopping model. The rest of the paper is organized as follows. In Sec. II, we outline our theoretical methodology; in Sec. III, we compare our calculated anisotropy to experiment, and analyze the mechanism, timescales, and transport pathways of OH vibrational energy transfer in ice Ih; and in Sec. IV, we conclude. II. THEORETICAL METHODS
We are interested in the OH stretch in ice Ih of either neat H2 O or H2 O/HOD/D2 O mixtures. In our mixed quantum/classical approach, we partition the nuclear degrees of freedom in ice Ih into two parts, system and bath—the system is treated quantum mechanically and consists of all the OH stretches, and the bath is treated classically and includes the low-frequency phonon and librational modes. We write the Frenkel one-exciton Hamiltonian for the OH stretches in the site basis as ¯ωi (t)|ii| + ¯ωij (t)|ij |, (1) H (t) = i
i,j =i
where |i is the state of the system in which only chromophore i (i.e., the ith OH vibrator) is singly excited, ωi (t) is the fluctuating (due to the classical bath) 0-1 vibrational transition frequency of chromophore i, and ωij (t) is the fluctuating vibrational coupling (intramolecular or intermolecular) between chromophores i and j. The transition frequencies and vibrational couplings depend on the bath degrees of freedom. In principle, they could be calculated on-the-fly through ab initio or DFT electronicstructure methods, but in practice that is not feasible. Instead,
previous work in our group found that the transition frequency calculated from DFT methods on clusters correlates with the electric field on the hydrogen atom along the OH bond, and we parametrized the correlation as a frequency “map.”42, 47 A “map” for the intramolecular coupling (denoted as ωija ) was also determined in a similar way.43–45 For the intermolecular coupling between chromophores i and j, denoted as ωije , we assumed the form of transition dipole interaction, given by ¯ωije =
μi xi μj xj [uˆ i · uˆ j − 3(uˆ i · nˆ ij )(uˆ j · nˆ ij )] rij3
,
(2)
where μi is the dipole derivative of chromophore i, xi is the 0-1 position matrix element for chromophore i, and uˆ i is the unit vector of the transition dipole of chromophore i, which is taken to be along the OH bond. Both the dipole derivative and the position matrix elements have “maps” based on DFT calculations.42, 45 rij is the distance between the point transition dipoles of chromophores i and j, and the position of the point transition dipole has been parametrized to be 0.67 Å from the oxygen atom along the OH bond.43 nˆ ij is the unit vector connecting the two point dipoles from i to j. As for the bath, we treat it classically via the explicit three-body (E3B) rigid water model, recently developed in our group.48, 49 The E3B model has a better freezing point compared to the popular SPC/E and TIP4P models,49 and has been shown to be fairly versatile in simulating the vibrational spectroscopy of liquid water, the air/water interface, water hexamer and ice.37, 42–46, 50–53 Therefore, the bath degrees of freedom evolve classically according to the E3B potential, and the intramolecular bends are frozen. For the OD stretches in the H2 O/HOD/D2 O mixture, even though there are also vibrational couplings between OH and OD stretches, the effect of these couplings is negligible due to the large frequency mismatch between these two stretches. Thus, the OD stretches are not included in the Hamiltonian. The system time-evolution operator U(t) obeys the timedependent Schrödinger equation i (3) U˙ (t) = − H (t)U (t). ¯ Suppose that the initial excitation is localized on site i (i.e., the ith OH chromophore), then the average probability of finding the excitation remaining on site i at time t is given by |Uii (t)|2 , where the brackets denote a time average over the trajectories. Therefore, the average survival probability S(t) is defined as S(t) =
N 1 |Uii (t)|2 , N i=1
(4)
where N is the number of OH chromophores in the system. By monitoring the time evolution of the survival probability, one can characterize the resonant energy transfer in the system. As mentioned in the Introduction, the polarizationresolved vibrational pump-probe spectroscopy can be used to monitor the rotation of the transition dipoles of the vibrational excitations through the anisotropy decay, given by r(t) =
I (t) − I⊥ (t) , I (t) + 2I⊥ (t)
(5)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
Shi, Li, and Skinner
where I (t) is the pump-probe signal when the probe polarization is parallel to the pump polarization, I⊥ (t) is the pump-probe signal when the probe polarization is perpendicular to the pump polarization, and t is the time delay between the pump and probe pulses. The calculation of vibrational anisotropy decay usually involves the doubly-excited vibrational states, as well as the singly-excited states and the ground state. However, Yang, Li, and Skinner39 showed that within some reasonable approximations, the anisotropy decay for frequency-integrated detection can be calculated by 2 j mj (t)2 |Uj i (t)|2 mi (0)2 P2 (uˆ i (0) · uˆ j (t))i , r(t) = 5 j mj (t)2 |Uj i (t)|2 mi (0)2 i (6) where mi ≡ μi xi is the magnitude of the transition dipole of chromophore i, P2 is the second Legendre polynomial, and the brackets with the subscript i indicates a time average over the trajectories and over all initially excited chromophores i. Note that within the assumption that vibrational energy relaxation and energy transfer are independent dynamical processes, the relaxation terms cancel out for the anisotropy, and therefore do not appear in the above equation. In order to evolve the classical bath, we perform a molecular dynamics simulation in the NVT ensemble at 245 K using a modified GROMACS version 3.354 for the E3B model. As before,45 we employ a 432-molecule proton-disordered configuration generated by Hayward and Reimer as our initial configuration of ice Ih,55 and the simulation box is scaled to match the experimental lattice constants.56 Note that for neat H2 O, the simulation is run with H2 O, while for the H2 O/HOD/D2 O mixtures, the simulation is run with D2 O. Other simulation details are the same as those in Ref. 45. For the calculation of the anisotropy for the H2 O/HOD/D2 O mixtures, we assign the D atoms at random to be H atoms with a probability given by the mole fraction of H2 O in D2 O, as is approximately the case in real water.57, 58 Then we perform an average over 80 different random realizations of the H2 O/HOD/D2 O mixture to obtain a result for each concentration.
J. Chem. Phys. 140, 244503 (2014) 0.4 100% 50% 25% 15% 10% 5%
0.3
Anisotropy r(t)
244503-3
0.2
0.1
0 0
0.5
1
1.5
2
t (ps) FIG. 1. Calculated (at 245 K, solid lines) and experimental (at 270 K, symbols)15 pump-probe anisotropy for H2 O/HOD/D2 O ice Ih at various hydrogen mole fractions (as indicated by the colors and the legend).
by Timmer and Bakker,15 the level of the plateau reflects the probability of finding acceptors within the first solvation shell. Assuming that our model captures the physics of the OH stretch resonant energy transfer in ice Ih, we will now attempt to understand this process in more detail, and will focus our attention on the case of neat H2 O hereafter. To this end, we first classify the relevant OH donor-acceptor pairs in the system. In Fig. 2, a six-molecule fragment of ice Ih is shown with a vertical c axis of the hexagonal crystal lattice. The large blue spheres are oxygen atoms, and other color-coded small spheres are hydrogen atoms. For pictorial clarity, we assume
III. RESULTS AND DISCUSSION
Fig. 1 displays the experimental (at 270 K)15 and calculated (at 245 K) (from Eq. (6)) pump-probe anisotropy decay for H2 O/HOD/D2 O ice Ih at various hydrogen mole fractions including 100%. Note that the experimental work by Timmer and Bakker showed that the anisotropy decay is not sensitive to the probe frequency,15 and so our calculation is done for frequency-integrated detection. We choose 245 K for the calculation because that is a few degrees below the melting point of our simulation model.49 As shown in Fig. 1, our calculations are in very good agreement with experiment,15 for all hydrogen mole fractions. Experiment and theory show that the anisotropy decay is very sensitive to hydrogen mole fraction: the higher the mole fraction the faster the decay and the lower the level of the observed “plateau” at around 1 ps. Both of these reflect the fact that for high hydrogen mole fraction there are more nearby OH stretch acceptors, leading to more energy transfer and faster anisotropy decay. As was suggested
FIG. 2. Vibrational resonant energy transfer donor-acceptor pairs in a sixmolecule fragment of ice Ih. The c axis is vertical, and the dotted lines are hydrogen bonds. The large blue spheres are oxygen atoms, and the small spheres of white and other colors are hydrogen atoms. The locations of the vibrational excitation are taken to be on the hydrogen atoms for pictorial clarity. Assume that the initial vibrational excitation is localized on the OH chromophore involving the red hydrogen atom, and the vibrational excitation can transfer to the four green hydrogen atoms through the intermolecular couplings of S class. The excitation can also transfer to the orange hydrogen atom through the intermolecular coupling of WII class, to the yellow hydrogen atoms through the intermolecular coupling of weak classes (only W60, W120, and W180 shown here), and to the purple hydrogen atom through the intramolecular coupling. The classification of the intermolecular pairs is explained in the main text.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
Shi, Li, and Skinner
J. Chem. Phys. 140, 244503 (2014)
the initial excitation locates on the red hydrogen atom, and the water molecule to which the red hydrogen atom belongs is regarded as the central water molecule. The red and purple hydrogen atoms in the central water form an intramolecular donor-acceptor pair, and the magnitude of its coupling is very small as shown in Fig. 14 of Ref. 45. The strongest intermolecular couplings are the four red-green pairs, denoted as S for strong,45 and we expect these four green hydrogens to be the main acceptors of the vibrational energy. The other four yellow hydrogen atoms within the first solvation shell of the central water molecule might also accept the excitation through the weak intermolecular couplings, which can be further classified as W0, W60, W120, and W180 pairs, where W is for weak and the numerals are the dihedral angles formed by the two OH bonds.45 Note that in this particular fragment in Fig. 2, only W60, W120, and W180 pairs are shown. The final important donor-acceptor class is the red-orange pair, which is denoted as WII (W is for weak, and II means the second solvation shell of the central water molecule). Actually, our previous study showed that the vibrational coupling of the WII pair is the second strongest intermolecular coupling in the system.45 Besides the above six classes of donoracceptor pairs, other pairs (e.g., red-white hydrogen pairs in Fig. 2) are negligible since they are very apart and their coupling strengths are very small.45 We thus group them together as the seventh class. We now examine the survival probability S(t) of the OH vibrational excitation, initially localized on one OH chromophore, which is a direct way to investigate the resonant energy transfer problem. As shown in Fig. 3, the survival probability (solid black line, from Eq. (4)) decays on a similar, but slightly longer time scale than the anisotropy decay (as in Fig. 1, 100% hydrogen mole fraction). As mentioned above, the S pairs (red-green hydrogen atoms in Fig. 2) are the main donor-acceptor pairs in the system, which we confirm in Fig. 3 by calculating the survival probability including only the intermolecular couplings of S pairs in the Hamiltonian (dashed black line). Therefore, at short times the excitation will primarily transfer from the initially excited OH chromophore to its four neighbor chromophores through the intermolecular couplings of S pairs. To further investigate the mechanism of the energy transfer in ice Ih, we would like to examine if an incoherent hopping model can reproduce this survival probability. A perturbative Markovian treatment of the Liouville equation for the system density operator leads to a master equation20, 21 P˙i (t) = −
j =i
kij Pi (t) +
kj i Pj (t),
(7)
j =i
where Pi (t) is the probability of finding the excitation on chromophore i at time t, and kij is the hopping rate constant for the energy transfer from chromophore i to j. In the context of our excitonic Hamiltonian (Eq. (1)) and our mixed quantum/classical approach, the hopping rate constant kij can be expressed as21, 23, 39
∞
kij = 2Re 0
t dt ωij (t)ωij (0)ei 0 dτ (ωi (τ )−ωj (τ )) ,
(8)
1
Survival probability S(t)
244503-4
"exact" (Eq. (4)) "exact" (Eq. (4)) with S pairs only master equation (Eq. (7))
0.8
0.6
0.4
0.2
0
0
0.1
0.2
0.3
0.4
0.5
t (ps) FIG. 3. Calculated survival probability of the vibrational excitation. The solid black line is the “exact” result from Eq. (4), and the dashed black line is also from Eq. (4) but with only the intermolecular couplings of S class included in the excitonic Hamiltonian. The red line is the result from the master equation (Eq. (7)) with the hopping rates calculated using Eq. (8).
where ωij is either the intramolecular or intermolecular coupling between chromophores i and j, and the brackets indicate a time average over the trajectory. Note that since the bath is treated classically, we have kij = kji . In principle, we can evaluate the hopping rate for each pair in the system (about 0.4 × 106 pairs in a 432-molecule simulation box), but it is not practical (and not necessary) to do so. Recognizing that rates for two different pairs in the same class will be similar (since their coupling constants will be similar), instead we simply calculate the average rate for each donor-acceptor class, as outlined above. In Fig. 3, the solid red line is the calculated survival probability by solving the master equation numerically, substituting in the average rate constants for pairs in each donor-acceptor class, and using Eq. (4). A significant underestimation of the survival probability is observed compared to the “exact” calculation (solid black line), implying that incoherent hopping is not the right picture for the resonant OH vibrational energy transfer in ice Ih. Therefore, the resonant OH vibrational energy transfer has some coherent features at least up to 100 fs. Even though the incoherent hopping model is not correct here, the hopping rate constants are still helpful for recognizing the important donor-acceptor pairs in the system. Therefore, the average hopping rates for all the seven classes of pairs are shown in Table I. The average hopping rate for S pairs is one order of magnitude larger than those for other pairs, and 6.73 ps−1 agrees with the value (1/140 fs−1 ) estimated by Timmer and Bakker.15 This explains why including only the S pairs is sufficient to get the survival probability semi-quantitatively correct (Fig. 3). The hopping rates of other pairs are all less than 1.00 ps−1 , and the second largest hopping rate is from WII pairs as expected. The hopping rates for W0 and W60 pairs are of the same order of magnitude as that for WII pairs, as is the hopping rate for the intramolecular pair, and other pairs are not efficient energy-transfer partners. Note that these results agree qualitatively with the six-nearestneighbor idea, suggested by Timmer and Bakker, although
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
244503-5
Shi, Li, and Skinner
J. Chem. Phys. 140, 244503 (2014)
TABLE I. Average hopping rates for different donor-acceptor pairs calculated from Eq. (8). The color pairs in the parentheses correspond to the pairs of the color-coded hydrogen atoms in Fig. 2. Average hopping rates (ps−1 )
Donor-acceptor pair class S (red-green) WII (red-orange) Intra (red-purple) W0 (not shown) W60 (red-yellow) W120 (red-yellow) W180 (red-yellow) Other (red-white)
6.73 0.77 0.36 0.64 0.31 0.05 0.02 0.002
their treatment of the intramolecular coupling as arising from dipole-dipole interactions greatly overestimates the hopping rate.15 Also note that a comparison of the coupling matrix elements for the different classes of pairs, as was done in Ref. 45, leads to a similar conclusion about the relative importance of these different classes for energy transfer. Finally, we can discuss the timescale and length scale of the energy transfer in neat H2 O ice Ih. To examine the spread as time evolves of an initially localized excitation, we define a root-mean-square displacement (RMSD) for the excitation, given by59 RMSD(t) =
| rj (t) − r i (0)|2 |Uj i (t)|2 , (9) j
i
0.8
0.8
0.6
0.6 RMSD (nm)
rHH (nm)
where r j (t) is the position vector of the hydrogen atom of chromophore j at time t. The calculated RMSD is plotted in Fig. 4. Up to 50 fs, the RMSD is roughly linear in time, which is a signature of coherent motion.60 Thus, we again conclude
0.4 0.2 0
0.4 0.2
0
1 2 gHH(rHH)
0
3
0
0.1
0.2 0.3 t (ps)
0.4
0.5
0.1
0.2 0.3 t (ps)
0.4
0.5
0.4
r(t)
0.3 0.2 0.1 0 0
FIG. 4. RMSD of the OH vibrational excitation in neat H2 O ice Ih. The black solid line in the top right panel is the RMSD of the excitation calculated from Eq. (9), the top left panel is the H–H radial distribution function, and the bottom right panel is the “exact” anisotropy decay from Eq. (6). The vertical dashed lines indicate the time scale of the excitation spreading, and the horizontal dashed lines indicate the corresponding length scale. The red solid line in the top right panel is the RMSD of the excitation with only the intermolecular couplings of S pairs in the Hamiltonian.
that the vibrational energy transfer is at least partially coherent on the timescale of the anisotropy decay. To the left of the plot is shown the H–H radial distribution function (the intramolecular hydrogen peak is not shown), and below the plot is shown the anisotropy decay. We can then ask, for example, how long does an excitation take to move a distance of about 0.23 nm (the first peak in gHH (r))? From the RMSD we see that it takes about 80 fs. At 80 fs, the anisotropy has decayed to about 0.1. Likewise, at about 180 fs, the excitation has been transported to the second solvation shell of the central molecule (at a distance of about 0.4 nm), and the anisotropy decay has almost vanished. To reiterate the dominance of the S pairs as the most effective energy transfer partners, we calculate the RMSD including coupling for those pairs alone, and the result (red line in the figure) indicates that most of the excitation spreading is due to these pairs. IV. CONCLUSIONS
In this work, we simulate the vibrational pump-probe anisotropy decays for neat H2 O ice Ih and mixtures with deuterium substitution. We have used a mixed quantum/classical approach, where the classical dynamics of the low-frequency modes arose from the E3B water simulation model, and quantum dynamics of the OH-stretch oscillators was described by a one-exciton Hamiltonian in the local-mode basis. The dependence of the Hamiltonian matrix elements on the classical variables has been discussed before, and was used without modification from our earlier work. The comparison between theory and experiment for the time-dependent anisotropy decay, and especially the dependence on the hydrogen mole fraction, was very good. This comparison provides a sensitive test of the vibrational coupling matrix elements in the Hamiltonian, as these lead to the energy transfer, which in turn is responsible for the anisotropy decay. The dramatic slow-down in the anisotropy decay with decreasing hydrogen mole fraction provides very strong evidence, from experiment and theory, about the dominance of vibrational energy transfer in ice Ih. Since vibrational couplings fall off quickly with distance, so do energy transfer rates. This led Timmer and Bakker to argue that both the initial decay of the anisotropy, and the level of the plateau at about 1 ps, are due to the probability of finding hydrogen acceptors in the first solvation shell. We agree with this interpretation. Another interesting issue involves whether the vibrational energy transport is coherent or incoherent on the relevant timescale. We present two pieces of evidence that it is coherent. First, assuming a hopping model we calculate the rate constants and then solve the master equation. We find that the survival probability when calculated with this approach decays much too quickly compared to the “exact” decay. Second, we find that the root-mean-square displacement of the vibrational excitation is linear in time up to 50 fs, which is a signature of coherent transport. These two facts are not unrelated: in a hopping model the root-mean-square displacement grows too fast for short times (like the square root of time). We can contrast this with the situation in liquid water, where an incoherent hopping model is appropriate. We
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
244503-6
Shi, Li, and Skinner
believe there are two reasons for this difference. First, the anisotropy decay, and hence the energy transfer, is slower in liquid water by a factor of 2 or 3. This is because the nicely aligned hydrogen bonds in ice lead to stronger intermolecular vibrational coupling. Second, while transport is certainly coherent in liquid water at very short times, because of the greater disorder in the liquid this timescale is likely significantly shorter than in ice. Therefore, the ratio of the relevant (anisotropy decay) time to the coherence time is likely a factor of 5 or so greater in liquid water, which would explain why the hopping model is appropriate in one situation but not in the other. Finally, there is the question of which donor-acceptor pairs are responsible for the resonant energy transfer. Our analysis indicates that for every OH stretch the energy transfer is dominated by four nearest-neighbor intermolecular acceptors, while the intramolecular energy transfer is much less important. Similar experiments can and should be performed on other fascinating crystalline and amorphous phases of ice; perhaps we will undertake theoretical studies of these as well in the near future. ACKNOWLEDGMENTS
This work was supported by NSF Grant No. CHE1058752. The authors are grateful to Professor Mino Yang for helpful conversations. 1 T.
Förster, Z. Naturforsch. A 4, 321 (1949). R. Selvin, Nat. Struct. Biol. 7, 730 (2000). 3 R. B. Sekar and A. Periasamy, J. Cell. Biol. 160, 629 (2003). 4 T. Ha and P. Tinnefeld, Annu. Rev. Phys. Chem. 63, 595 (2012). 5 G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Manˇ cal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature (London) 446, 782 (2007). 6 G. Panitchayangkoon, D. Hayes, K. A. Fransted, J. R. Caram, E. Harel, J. Wen, R. E. Blankenship, and G. S. Engel, Proc. Natl. Acad. Sci. U.S.A. 107, 12766 (2010). 7 K.-K. Lee, K.-I. Oh, H. Lee, C. Joo, H. Han, and M. Cho, ChemPhysChem 8, 2218 (2007). 8 E. Collini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer, and G. D. Scholes, Nature (London) 463, 644 (2010). 9 V. Botan, E. H. G. Backus, R. Pfister, A. Moretto, M. Crisma, C. Toniolo, P. H. Nguyen, G. Stock, and P. Hamm, Proc. Natl. Acad. Sci. U.S.A. 104, 12749 (2007). 10 E. H. G. Backus, P. H. Nguyen, V. Botan, R. Pfister, A. Moretto, M. Crisma, C. Toniolo, G. Stock, and P. Hamm, J. Phys. Chem. B 112, 9091 (2008). 11 E. H. G. Backus, P. H. Nguyen, V. Botan, A. Moretto, M. Crisma, C. Toniolo, O. Zerbe, G. Stock, and P. Hamm, J. Phys. Chem. B 112, 15487 (2008). 12 M. Schade, A. Moretto, P. M. Donaldson, C. Toniolo, and P. Hamm, Nano Lett. 10, 3057 (2010). 13 S. Woutersen and H. J. Bakker, Nature (London) 402, 507 (1999). 14 L. Piatkowski, K. B. Eisenthal, and H. J. Bakker, Phys. Chem. Chem. Phys. 11, 9033 (2009). 15 R. L. A. Timmer and H. J. Bakker, J. Phys. Chem. A 114, 4148 (2010). 16 M. L. Cowan, B. D. Bruner, N. Huse, J. R. Dwyer, B. Chugh, E. T. J. Nibbering, T. Elsaesser, and R. D. J. Miller, Nature (London) 434, 199 (2005). 2 P.
J. Chem. Phys. 140, 244503 (2014) 17 D.
Kraemer, M. L. Cowan, A. Paarmann, N. Huse, E. T. J. Nibbering, T. Elsaesser, and R. J. D. Miller, Proc. Natl. Acad. Sci. U.S.A. 105, 437 (2008). 18 W. Amir, G. Gallot, and F. Hache, J. Chem. Phys. 121, 7908 (2004). 19 K. Blum, Density Matrix Theory and Applications (Plenum Press, New York, 1980). 20 R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, New York, 2001). 21 M. Yang and G. R. Fleming, Chem. Phys. 282, 163 (2002). 22 A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234110 (2009). 23 S. Jang, J. Chem. Phys. 127, 174710 (2007). 24 D. E. Rosenfeld and M. D. Fayer, J. Chem. Phys. 137, 064109 (2012). 25 A. G. Redfield, IBM J. Res. Develop. 1, 19 (1957). 26 G. Lindblad, Commun. Math. Phys. 48, 119 (1976). 27 V. Gorini, A. Kossokowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976). 28 D. P. S. McCutcheon and A. Nazir, Phys. Rev. B 83, 165101 (2011). 29 D. P. S. McCutcheon and A. Nazir, J. Chem. Phys. 135, 114501 (2011). 30 A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234111 (2009). 31 S. Jang, Y.-C. Cheng, D. R. Reichman, and J. D. Eaves, J. Chem. Phys. 129, 101104 (2008). 32 P. Huo and D. F. Coker, J. Chem. Phys. 133, 184108 (2010). 33 E. N. Zimanyi and R. J. Silbey, Phil. Trans. R. Soc. A 370, 3620 (2012). 34 P. Rebentrost, R. Chakraborty, and A. Aspuru-Guzik, J. Chem. Phys. 131, 184102 (2009). 35 J. Roden, A. Eisfeld, W. Wolff, and W. T. Strunz, Phys. Rev. Lett. 103, 058301 (2009). 36 J. L. Skinner, B. M. Auer, and Y.-S. Lin, Adv. Chem. Phys. 142, 59 (2009). 37 P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Chem. Phys. 135, 044701 (2011). 38 Y. Ni, S. M. Gruenbaum, and J. L. Skinner, Proc. Natl. Acad. Sci. U.S.A. 110, 1992 (2013). 39 M. Yang, F. Li, and J. L. Skinner, J. Chem. Phys. 135, 164505 (2011). 40 J. A. Poulsen, G. Nyman, and S. Nordholm, J. Phys. Chem. A 107, 8420 (2003). 41 C. Backtorp, J. A. Poulsen, and G. Nyman, J. Phys. Chem. A 109, 3105 (2005). 42 F. Li and J. L. Skinner, J. Chem. Phys. 132, 204505 (2010). 43 F. Li and J. L. Skinner, J. Chem. Phys. 133, 244504 (2010). 44 F. Li and J. L. Skinner, J. Chem. Phys. 134, 099901 (2011). 45 L. Shi, S. M. Gruenbaum, and J. L. Skinner, J. Phys. Chem. B 116, 13821 (2012). 46 L. Shi and J. L. Skinner, J. Phys. Chem. B 117, 15536 (2013). 47 S. M. Gruenbaum, C. J. Tainter, L. Shi, Y. Ni, and J. L. Skinner, J. Chem. Theor. Comput. 9, 3109 (2013). 48 R. Kumar and J. L. Skinner, J. Phys. Chem. B 112, 8311 (2008). 49 C. J. Tainter, P. A. Pieniazek, Y.-S. Lin, and J. L. Skinner, J. Chem. Phys. 134, 184501 (2011). 50 P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Am. Chem. Soc. 133, 10360 (2011). 51 I. V. Stiopkin, C. Weeraman, P. A. Pieniazek, F. Y. Shalhout, J. L. Skinner, and A. V. Benderskii, Nature (London) 474, 192 (2011). 52 C. J. Tainter and J. L. Skinner, J. Chem. Phys. 137, 104304 (2012). 53 C. J. Tainter, Y. Ni, L. Shi, and J. L. Skinner, J. Phys. Chem. Lett. 4, 12 (2013). 54 D. van der Spoel, E. Lindahl, B. Hess, A. R. van Buuren, E. Apol, P. J. Meulenhoff, D. P. Tieleman, A. L. T. M. Sijbers, K. A. Feenstra, R. van Drunen et al., Gromacs User Manual Version 3.3, 2005, see www.gromacs.org. 55 J. Hayward and J. Reimers, J. Chem. Phys. 106, 1518 (1997). 56 V. F. Petrenko and R. W. Whitworth, Physics of Ice (Oxford University Press, Oxford, 1999). 57 J. W. Pyper and L. D. Christensen, J. Chem. Phys. 62, 2596 (1975). 58 J. C. Duplan, L. Mahi, and J. L. Brunet, Chem. Phys. Lett. 413, 400 (2005). 59 T. L. C. Jansen, B. M. Auer, M. Yang, and J. L. Skinner, J. Chem. Phys. 132, 224503 (2010). 60 M. Grover and R. Silbey, J. Chem. Phys. 54, 4843 (1971).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
THE JOURNAL OF CHEMICAL PHYSICS 140, 244503 (2014)
Resonant vibrational energy transfer in ice Ih L. Shi, F. Li, and J. L. Skinner Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin, Madison, Wisconsin 53706, USA
(Received 4 April 2014; accepted 5 June 2014; published online 24 June 2014) Fascinating anisotropy decay experiments have recently been performed on H2 O ice Ih by Timmer and Bakker [R. L. A. Timmer, and H. J. Bakker, J. Phys. Chem. A 114, 4148 (2010)]. The very fast decay (on the order of 100 fs) is indicative of resonant energy transfer between OH stretches on different molecules. Isotope dilution experiments with deuterium show a dramatic dependence on the hydrogen mole fraction, which confirms the energy transfer picture. Timmer and Bakker have interpreted the experiments with a Förster incoherent hopping model, finding that energy transfer within the first solvation shell dominates the relaxation process. We have developed a microscopic theory of vibrational spectroscopy of water and ice, and herein we use this theory to calculate the anisotropy decay in ice as a function of hydrogen mole fraction. We obtain very good agreement with experiment. Interpretation of our results shows that four nearest-neighbor acceptors dominate the energy transfer, and that while the incoherent hopping picture is qualitatively correct, vibrational energy transport is partially coherent on the relevant timescale. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4883913] I. INTRODUCTION
Resonant energy transfer between an excited donor chromophore and an acceptor chromophore has been studied extensively by both experiment and theory for over 50 years.1–8 Most research involves electronic excitation, although due to the recent development of ultrafast infrared spectroscopy, researchers now have the opportunity to examine resonant vibrational energy transfer as well.9–18 Many studies have focused on the OH stretch in neat liquid water,13, 14, 16–18 in which many near-degenerate OH chromophores interact through intramolecular or intermolecular vibrational couplings, and their vibrational frequencies are modulated by other low-frequency translational and rotational modes. If OH chromophores are initially vibrationally excited by interacting with a polarized infrared pulse, the created anisotropy can decay through either molecular rotation or resonant energy transfer to other OH chromophores. The anisotropy decay can be detected experimentally by polarization-resolved pump-probe13, 14, 18 or two-dimensional infrared (2DIR) spectroscopy.16, 17 In neat liquid H2 O, the anisotropy decay occurs within 200 or 300 fs and is dominated by resonant energy transfer. Many theoretical methods have been proposed for understanding resonant energy transfer, which is often characterized as being incoherent or coherent.19–22 The incoherent or hopping picture occurs theoretically when the coupling between chromophores is treated perturbatively, and transport can be described by a master equation in the local site basis.20, 21 If the interaction between chromophores is dipoledipole, for example, this leads to the Förster theory.1, 21, 23, 24 Coherent transport occurs in the weak system-bath coupling limit, and can be described by Redfield25 or Lindblad26, 27 equations. There are also many theoretical works going beyond the Markovian approximation or the standard pertur0021-9606/2014/140(24)/244503/6/$30.00
bative treatment, and some of them treat incoherent/coherent energy transfer in a unified manner.28–35 In interpreting their pump-probe experimental results, Bakker and co-workers adapted the Förster theory to describe the ultrafast resonant energy transfer in liquid water, and found that it fit their experimental data very well.13, 14 Skinner and co-workers have developed a mixed quantum/classical model for OH-stretch spectroscopy in water, and have used it to calculate FTIR, Raman, 2DIR, SFG, and 2DSFG spectra.36–38 The exciton Hamiltonian for the OHstretch vibrations is written in the local-mode basis, and so the diagonal elements are the local-mode (anharmonic) vibrational frequencies, while the off-diagonal elements are the couplings (both intramolecular and intermolecular) between OH-stretch chromophores. The exciton Hamiltonian fluctuates as the classical degrees of freedom evolve in time. Depending on the statistical properties of this fluctuating Hamiltonian, such a model will produce coherent or incoherent transport (or something in between). In fact, application of this model by Yang, Li, and Skinner39 showed that in liquid water the incoherent hopping picture is appropriate, although the interpretation of experiments in terms of the usual Förster theory is not so straightforward. Timmer and Bakker recently extended their anisotropy decay experiments to ice Ih,15 finding a somewhat faster decay on the 100 fs timescale. They also studied what happens when the hydrogen atoms in H2 O ice are replaced by deuterium. That is, when D2 O is added to H2 O, one obtains a mixture of H2 O/HOD/D2 O, and to a good approximation one can model this system by assuming that the probability that each hydrogen atom is actually hydrogen (instead of deuterium) is simply equal to the hydrogen mole fraction. The consequences for energy transfer are dramatic, since the OD stretches are so far off resonance they do not function as energy acceptors. Therefore, as the hydrogen mole fraction
140, 244503-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
244503-2
Shi, Li, and Skinner
J. Chem. Phys. 140, 244503 (2014)
decreases the anisotropy decay slows down dramatically. These beautiful experiments directly confirm the essential picture of energy-transfer-induced anisotropy decay. Timmer and Bakker interpreted their experiments in terms of the Förster theory of incoherent hopping. Within their model the rate of intramolecular hopping was by far the largest (since the distance between the two dipoles is the closest), followed by intermolecular hopping to four near neighbors. They concluded that energy transfer outside the first solvation shell was negligible. Two theoretical papers on this problem, by Poulsen, Nyman, and co-workers actually preceded the experiments.40, 41 They used clusters of 15 water molecules to represent ice Ih, and ran wave-packet dynamics to monitor the survival probability of the excitiation on the initially excited mode.40 This led them to conclude that the resonant OH vibrational energy transfer in ice Ih occurs within 100 fs, in agreement with the recent pump-probe experiment.15 However, they did not directly calculate the anisotropy observable. The goal of this paper is to apply our mixed quantum/classical model to calculate the anisotropy decay in ice Ih, for both neat H2 O and for the isotopic mixtures. This provides a stringent test of the model since there are no free parameters. We have previously used the same model to calculate IR, Raman, and 2DIR spectra of ice Ih.42–46 We find, in fact, that the model performs quite well for the anisotropy decay, for all hydrogen mole fractions. We then analyze our results to examine the relative contribution from intramolecular and intermolecular energy transfer, from first and second solvation shells, and the appropriateness of the incoherent hopping model. The rest of the paper is organized as follows. In Sec. II, we outline our theoretical methodology; in Sec. III, we compare our calculated anisotropy to experiment, and analyze the mechanism, timescales, and transport pathways of OH vibrational energy transfer in ice Ih; and in Sec. IV, we conclude. II. THEORETICAL METHODS
We are interested in the OH stretch in ice Ih of either neat H2 O or H2 O/HOD/D2 O mixtures. In our mixed quantum/classical approach, we partition the nuclear degrees of freedom in ice Ih into two parts, system and bath—the system is treated quantum mechanically and consists of all the OH stretches, and the bath is treated classically and includes the low-frequency phonon and librational modes. We write the Frenkel one-exciton Hamiltonian for the OH stretches in the site basis as ¯ωi (t)|ii| + ¯ωij (t)|ij |, (1) H (t) = i
i,j =i
where |i is the state of the system in which only chromophore i (i.e., the ith OH vibrator) is singly excited, ωi (t) is the fluctuating (due to the classical bath) 0-1 vibrational transition frequency of chromophore i, and ωij (t) is the fluctuating vibrational coupling (intramolecular or intermolecular) between chromophores i and j. The transition frequencies and vibrational couplings depend on the bath degrees of freedom. In principle, they could be calculated on-the-fly through ab initio or DFT electronicstructure methods, but in practice that is not feasible. Instead,
previous work in our group found that the transition frequency calculated from DFT methods on clusters correlates with the electric field on the hydrogen atom along the OH bond, and we parametrized the correlation as a frequency “map.”42, 47 A “map” for the intramolecular coupling (denoted as ωija ) was also determined in a similar way.43–45 For the intermolecular coupling between chromophores i and j, denoted as ωije , we assumed the form of transition dipole interaction, given by ¯ωije =
μi xi μj xj [uˆ i · uˆ j − 3(uˆ i · nˆ ij )(uˆ j · nˆ ij )] rij3
,
(2)
where μi is the dipole derivative of chromophore i, xi is the 0-1 position matrix element for chromophore i, and uˆ i is the unit vector of the transition dipole of chromophore i, which is taken to be along the OH bond. Both the dipole derivative and the position matrix elements have “maps” based on DFT calculations.42, 45 rij is the distance between the point transition dipoles of chromophores i and j, and the position of the point transition dipole has been parametrized to be 0.67 Å from the oxygen atom along the OH bond.43 nˆ ij is the unit vector connecting the two point dipoles from i to j. As for the bath, we treat it classically via the explicit three-body (E3B) rigid water model, recently developed in our group.48, 49 The E3B model has a better freezing point compared to the popular SPC/E and TIP4P models,49 and has been shown to be fairly versatile in simulating the vibrational spectroscopy of liquid water, the air/water interface, water hexamer and ice.37, 42–46, 50–53 Therefore, the bath degrees of freedom evolve classically according to the E3B potential, and the intramolecular bends are frozen. For the OD stretches in the H2 O/HOD/D2 O mixture, even though there are also vibrational couplings between OH and OD stretches, the effect of these couplings is negligible due to the large frequency mismatch between these two stretches. Thus, the OD stretches are not included in the Hamiltonian. The system time-evolution operator U(t) obeys the timedependent Schrödinger equation i (3) U˙ (t) = − H (t)U (t). ¯ Suppose that the initial excitation is localized on site i (i.e., the ith OH chromophore), then the average probability of finding the excitation remaining on site i at time t is given by |Uii (t)|2 , where the brackets denote a time average over the trajectories. Therefore, the average survival probability S(t) is defined as S(t) =
N 1 |Uii (t)|2 , N i=1
(4)
where N is the number of OH chromophores in the system. By monitoring the time evolution of the survival probability, one can characterize the resonant energy transfer in the system. As mentioned in the Introduction, the polarizationresolved vibrational pump-probe spectroscopy can be used to monitor the rotation of the transition dipoles of the vibrational excitations through the anisotropy decay, given by r(t) =
I (t) − I⊥ (t) , I (t) + 2I⊥ (t)
(5)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
Shi, Li, and Skinner
where I (t) is the pump-probe signal when the probe polarization is parallel to the pump polarization, I⊥ (t) is the pump-probe signal when the probe polarization is perpendicular to the pump polarization, and t is the time delay between the pump and probe pulses. The calculation of vibrational anisotropy decay usually involves the doubly-excited vibrational states, as well as the singly-excited states and the ground state. However, Yang, Li, and Skinner39 showed that within some reasonable approximations, the anisotropy decay for frequency-integrated detection can be calculated by 2 j mj (t)2 |Uj i (t)|2 mi (0)2 P2 (uˆ i (0) · uˆ j (t))i , r(t) = 5 j mj (t)2 |Uj i (t)|2 mi (0)2 i (6) where mi ≡ μi xi is the magnitude of the transition dipole of chromophore i, P2 is the second Legendre polynomial, and the brackets with the subscript i indicates a time average over the trajectories and over all initially excited chromophores i. Note that within the assumption that vibrational energy relaxation and energy transfer are independent dynamical processes, the relaxation terms cancel out for the anisotropy, and therefore do not appear in the above equation. In order to evolve the classical bath, we perform a molecular dynamics simulation in the NVT ensemble at 245 K using a modified GROMACS version 3.354 for the E3B model. As before,45 we employ a 432-molecule proton-disordered configuration generated by Hayward and Reimer as our initial configuration of ice Ih,55 and the simulation box is scaled to match the experimental lattice constants.56 Note that for neat H2 O, the simulation is run with H2 O, while for the H2 O/HOD/D2 O mixtures, the simulation is run with D2 O. Other simulation details are the same as those in Ref. 45. For the calculation of the anisotropy for the H2 O/HOD/D2 O mixtures, we assign the D atoms at random to be H atoms with a probability given by the mole fraction of H2 O in D2 O, as is approximately the case in real water.57, 58 Then we perform an average over 80 different random realizations of the H2 O/HOD/D2 O mixture to obtain a result for each concentration.
J. Chem. Phys. 140, 244503 (2014) 0.4 100% 50% 25% 15% 10% 5%
0.3
Anisotropy r(t)
244503-3
0.2
0.1
0 0
0.5
1
1.5
2
t (ps) FIG. 1. Calculated (at 245 K, solid lines) and experimental (at 270 K, symbols)15 pump-probe anisotropy for H2 O/HOD/D2 O ice Ih at various hydrogen mole fractions (as indicated by the colors and the legend).
by Timmer and Bakker,15 the level of the plateau reflects the probability of finding acceptors within the first solvation shell. Assuming that our model captures the physics of the OH stretch resonant energy transfer in ice Ih, we will now attempt to understand this process in more detail, and will focus our attention on the case of neat H2 O hereafter. To this end, we first classify the relevant OH donor-acceptor pairs in the system. In Fig. 2, a six-molecule fragment of ice Ih is shown with a vertical c axis of the hexagonal crystal lattice. The large blue spheres are oxygen atoms, and other color-coded small spheres are hydrogen atoms. For pictorial clarity, we assume
III. RESULTS AND DISCUSSION
Fig. 1 displays the experimental (at 270 K)15 and calculated (at 245 K) (from Eq. (6)) pump-probe anisotropy decay for H2 O/HOD/D2 O ice Ih at various hydrogen mole fractions including 100%. Note that the experimental work by Timmer and Bakker showed that the anisotropy decay is not sensitive to the probe frequency,15 and so our calculation is done for frequency-integrated detection. We choose 245 K for the calculation because that is a few degrees below the melting point of our simulation model.49 As shown in Fig. 1, our calculations are in very good agreement with experiment,15 for all hydrogen mole fractions. Experiment and theory show that the anisotropy decay is very sensitive to hydrogen mole fraction: the higher the mole fraction the faster the decay and the lower the level of the observed “plateau” at around 1 ps. Both of these reflect the fact that for high hydrogen mole fraction there are more nearby OH stretch acceptors, leading to more energy transfer and faster anisotropy decay. As was suggested
FIG. 2. Vibrational resonant energy transfer donor-acceptor pairs in a sixmolecule fragment of ice Ih. The c axis is vertical, and the dotted lines are hydrogen bonds. The large blue spheres are oxygen atoms, and the small spheres of white and other colors are hydrogen atoms. The locations of the vibrational excitation are taken to be on the hydrogen atoms for pictorial clarity. Assume that the initial vibrational excitation is localized on the OH chromophore involving the red hydrogen atom, and the vibrational excitation can transfer to the four green hydrogen atoms through the intermolecular couplings of S class. The excitation can also transfer to the orange hydrogen atom through the intermolecular coupling of WII class, to the yellow hydrogen atoms through the intermolecular coupling of weak classes (only W60, W120, and W180 shown here), and to the purple hydrogen atom through the intramolecular coupling. The classification of the intermolecular pairs is explained in the main text.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
Shi, Li, and Skinner
J. Chem. Phys. 140, 244503 (2014)
the initial excitation locates on the red hydrogen atom, and the water molecule to which the red hydrogen atom belongs is regarded as the central water molecule. The red and purple hydrogen atoms in the central water form an intramolecular donor-acceptor pair, and the magnitude of its coupling is very small as shown in Fig. 14 of Ref. 45. The strongest intermolecular couplings are the four red-green pairs, denoted as S for strong,45 and we expect these four green hydrogens to be the main acceptors of the vibrational energy. The other four yellow hydrogen atoms within the first solvation shell of the central water molecule might also accept the excitation through the weak intermolecular couplings, which can be further classified as W0, W60, W120, and W180 pairs, where W is for weak and the numerals are the dihedral angles formed by the two OH bonds.45 Note that in this particular fragment in Fig. 2, only W60, W120, and W180 pairs are shown. The final important donor-acceptor class is the red-orange pair, which is denoted as WII (W is for weak, and II means the second solvation shell of the central water molecule). Actually, our previous study showed that the vibrational coupling of the WII pair is the second strongest intermolecular coupling in the system.45 Besides the above six classes of donoracceptor pairs, other pairs (e.g., red-white hydrogen pairs in Fig. 2) are negligible since they are very apart and their coupling strengths are very small.45 We thus group them together as the seventh class. We now examine the survival probability S(t) of the OH vibrational excitation, initially localized on one OH chromophore, which is a direct way to investigate the resonant energy transfer problem. As shown in Fig. 3, the survival probability (solid black line, from Eq. (4)) decays on a similar, but slightly longer time scale than the anisotropy decay (as in Fig. 1, 100% hydrogen mole fraction). As mentioned above, the S pairs (red-green hydrogen atoms in Fig. 2) are the main donor-acceptor pairs in the system, which we confirm in Fig. 3 by calculating the survival probability including only the intermolecular couplings of S pairs in the Hamiltonian (dashed black line). Therefore, at short times the excitation will primarily transfer from the initially excited OH chromophore to its four neighbor chromophores through the intermolecular couplings of S pairs. To further investigate the mechanism of the energy transfer in ice Ih, we would like to examine if an incoherent hopping model can reproduce this survival probability. A perturbative Markovian treatment of the Liouville equation for the system density operator leads to a master equation20, 21 P˙i (t) = −
j =i
kij Pi (t) +
kj i Pj (t),
(7)
j =i
where Pi (t) is the probability of finding the excitation on chromophore i at time t, and kij is the hopping rate constant for the energy transfer from chromophore i to j. In the context of our excitonic Hamiltonian (Eq. (1)) and our mixed quantum/classical approach, the hopping rate constant kij can be expressed as21, 23, 39
∞
kij = 2Re 0
t dt ωij (t)ωij (0)ei 0 dτ (ωi (τ )−ωj (τ )) ,
(8)
1
Survival probability S(t)
244503-4
"exact" (Eq. (4)) "exact" (Eq. (4)) with S pairs only master equation (Eq. (7))
0.8
0.6
0.4
0.2
0
0
0.1
0.2
0.3
0.4
0.5
t (ps) FIG. 3. Calculated survival probability of the vibrational excitation. The solid black line is the “exact” result from Eq. (4), and the dashed black line is also from Eq. (4) but with only the intermolecular couplings of S class included in the excitonic Hamiltonian. The red line is the result from the master equation (Eq. (7)) with the hopping rates calculated using Eq. (8).
where ωij is either the intramolecular or intermolecular coupling between chromophores i and j, and the brackets indicate a time average over the trajectory. Note that since the bath is treated classically, we have kij = kji . In principle, we can evaluate the hopping rate for each pair in the system (about 0.4 × 106 pairs in a 432-molecule simulation box), but it is not practical (and not necessary) to do so. Recognizing that rates for two different pairs in the same class will be similar (since their coupling constants will be similar), instead we simply calculate the average rate for each donor-acceptor class, as outlined above. In Fig. 3, the solid red line is the calculated survival probability by solving the master equation numerically, substituting in the average rate constants for pairs in each donor-acceptor class, and using Eq. (4). A significant underestimation of the survival probability is observed compared to the “exact” calculation (solid black line), implying that incoherent hopping is not the right picture for the resonant OH vibrational energy transfer in ice Ih. Therefore, the resonant OH vibrational energy transfer has some coherent features at least up to 100 fs. Even though the incoherent hopping model is not correct here, the hopping rate constants are still helpful for recognizing the important donor-acceptor pairs in the system. Therefore, the average hopping rates for all the seven classes of pairs are shown in Table I. The average hopping rate for S pairs is one order of magnitude larger than those for other pairs, and 6.73 ps−1 agrees with the value (1/140 fs−1 ) estimated by Timmer and Bakker.15 This explains why including only the S pairs is sufficient to get the survival probability semi-quantitatively correct (Fig. 3). The hopping rates of other pairs are all less than 1.00 ps−1 , and the second largest hopping rate is from WII pairs as expected. The hopping rates for W0 and W60 pairs are of the same order of magnitude as that for WII pairs, as is the hopping rate for the intramolecular pair, and other pairs are not efficient energy-transfer partners. Note that these results agree qualitatively with the six-nearestneighbor idea, suggested by Timmer and Bakker, although
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
244503-5
Shi, Li, and Skinner
J. Chem. Phys. 140, 244503 (2014)
TABLE I. Average hopping rates for different donor-acceptor pairs calculated from Eq. (8). The color pairs in the parentheses correspond to the pairs of the color-coded hydrogen atoms in Fig. 2. Average hopping rates (ps−1 )
Donor-acceptor pair class S (red-green) WII (red-orange) Intra (red-purple) W0 (not shown) W60 (red-yellow) W120 (red-yellow) W180 (red-yellow) Other (red-white)
6.73 0.77 0.36 0.64 0.31 0.05 0.02 0.002
their treatment of the intramolecular coupling as arising from dipole-dipole interactions greatly overestimates the hopping rate.15 Also note that a comparison of the coupling matrix elements for the different classes of pairs, as was done in Ref. 45, leads to a similar conclusion about the relative importance of these different classes for energy transfer. Finally, we can discuss the timescale and length scale of the energy transfer in neat H2 O ice Ih. To examine the spread as time evolves of an initially localized excitation, we define a root-mean-square displacement (RMSD) for the excitation, given by59 RMSD(t) =
| rj (t) − r i (0)|2 |Uj i (t)|2 , (9) j
i
0.8
0.8
0.6
0.6 RMSD (nm)
rHH (nm)
where r j (t) is the position vector of the hydrogen atom of chromophore j at time t. The calculated RMSD is plotted in Fig. 4. Up to 50 fs, the RMSD is roughly linear in time, which is a signature of coherent motion.60 Thus, we again conclude
0.4 0.2 0
0.4 0.2
0
1 2 gHH(rHH)
0
3
0
0.1
0.2 0.3 t (ps)
0.4
0.5
0.1
0.2 0.3 t (ps)
0.4
0.5
0.4
r(t)
0.3 0.2 0.1 0 0
FIG. 4. RMSD of the OH vibrational excitation in neat H2 O ice Ih. The black solid line in the top right panel is the RMSD of the excitation calculated from Eq. (9), the top left panel is the H–H radial distribution function, and the bottom right panel is the “exact” anisotropy decay from Eq. (6). The vertical dashed lines indicate the time scale of the excitation spreading, and the horizontal dashed lines indicate the corresponding length scale. The red solid line in the top right panel is the RMSD of the excitation with only the intermolecular couplings of S pairs in the Hamiltonian.
that the vibrational energy transfer is at least partially coherent on the timescale of the anisotropy decay. To the left of the plot is shown the H–H radial distribution function (the intramolecular hydrogen peak is not shown), and below the plot is shown the anisotropy decay. We can then ask, for example, how long does an excitation take to move a distance of about 0.23 nm (the first peak in gHH (r))? From the RMSD we see that it takes about 80 fs. At 80 fs, the anisotropy has decayed to about 0.1. Likewise, at about 180 fs, the excitation has been transported to the second solvation shell of the central molecule (at a distance of about 0.4 nm), and the anisotropy decay has almost vanished. To reiterate the dominance of the S pairs as the most effective energy transfer partners, we calculate the RMSD including coupling for those pairs alone, and the result (red line in the figure) indicates that most of the excitation spreading is due to these pairs. IV. CONCLUSIONS
In this work, we simulate the vibrational pump-probe anisotropy decays for neat H2 O ice Ih and mixtures with deuterium substitution. We have used a mixed quantum/classical approach, where the classical dynamics of the low-frequency modes arose from the E3B water simulation model, and quantum dynamics of the OH-stretch oscillators was described by a one-exciton Hamiltonian in the local-mode basis. The dependence of the Hamiltonian matrix elements on the classical variables has been discussed before, and was used without modification from our earlier work. The comparison between theory and experiment for the time-dependent anisotropy decay, and especially the dependence on the hydrogen mole fraction, was very good. This comparison provides a sensitive test of the vibrational coupling matrix elements in the Hamiltonian, as these lead to the energy transfer, which in turn is responsible for the anisotropy decay. The dramatic slow-down in the anisotropy decay with decreasing hydrogen mole fraction provides very strong evidence, from experiment and theory, about the dominance of vibrational energy transfer in ice Ih. Since vibrational couplings fall off quickly with distance, so do energy transfer rates. This led Timmer and Bakker to argue that both the initial decay of the anisotropy, and the level of the plateau at about 1 ps, are due to the probability of finding hydrogen acceptors in the first solvation shell. We agree with this interpretation. Another interesting issue involves whether the vibrational energy transport is coherent or incoherent on the relevant timescale. We present two pieces of evidence that it is coherent. First, assuming a hopping model we calculate the rate constants and then solve the master equation. We find that the survival probability when calculated with this approach decays much too quickly compared to the “exact” decay. Second, we find that the root-mean-square displacement of the vibrational excitation is linear in time up to 50 fs, which is a signature of coherent transport. These two facts are not unrelated: in a hopping model the root-mean-square displacement grows too fast for short times (like the square root of time). We can contrast this with the situation in liquid water, where an incoherent hopping model is appropriate. We
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
244503-6
Shi, Li, and Skinner
believe there are two reasons for this difference. First, the anisotropy decay, and hence the energy transfer, is slower in liquid water by a factor of 2 or 3. This is because the nicely aligned hydrogen bonds in ice lead to stronger intermolecular vibrational coupling. Second, while transport is certainly coherent in liquid water at very short times, because of the greater disorder in the liquid this timescale is likely significantly shorter than in ice. Therefore, the ratio of the relevant (anisotropy decay) time to the coherence time is likely a factor of 5 or so greater in liquid water, which would explain why the hopping model is appropriate in one situation but not in the other. Finally, there is the question of which donor-acceptor pairs are responsible for the resonant energy transfer. Our analysis indicates that for every OH stretch the energy transfer is dominated by four nearest-neighbor intermolecular acceptors, while the intramolecular energy transfer is much less important. Similar experiments can and should be performed on other fascinating crystalline and amorphous phases of ice; perhaps we will undertake theoretical studies of these as well in the near future. ACKNOWLEDGMENTS
This work was supported by NSF Grant No. CHE1058752. The authors are grateful to Professor Mino Yang for helpful conversations. 1 T.
Förster, Z. Naturforsch. A 4, 321 (1949). R. Selvin, Nat. Struct. Biol. 7, 730 (2000). 3 R. B. Sekar and A. Periasamy, J. Cell. Biol. 160, 629 (2003). 4 T. Ha and P. Tinnefeld, Annu. Rev. Phys. Chem. 63, 595 (2012). 5 G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Manˇ cal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature (London) 446, 782 (2007). 6 G. Panitchayangkoon, D. Hayes, K. A. Fransted, J. R. Caram, E. Harel, J. Wen, R. E. Blankenship, and G. S. Engel, Proc. Natl. Acad. Sci. U.S.A. 107, 12766 (2010). 7 K.-K. Lee, K.-I. Oh, H. Lee, C. Joo, H. Han, and M. Cho, ChemPhysChem 8, 2218 (2007). 8 E. Collini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer, and G. D. Scholes, Nature (London) 463, 644 (2010). 9 V. Botan, E. H. G. Backus, R. Pfister, A. Moretto, M. Crisma, C. Toniolo, P. H. Nguyen, G. Stock, and P. Hamm, Proc. Natl. Acad. Sci. U.S.A. 104, 12749 (2007). 10 E. H. G. Backus, P. H. Nguyen, V. Botan, R. Pfister, A. Moretto, M. Crisma, C. Toniolo, G. Stock, and P. Hamm, J. Phys. Chem. B 112, 9091 (2008). 11 E. H. G. Backus, P. H. Nguyen, V. Botan, A. Moretto, M. Crisma, C. Toniolo, O. Zerbe, G. Stock, and P. Hamm, J. Phys. Chem. B 112, 15487 (2008). 12 M. Schade, A. Moretto, P. M. Donaldson, C. Toniolo, and P. Hamm, Nano Lett. 10, 3057 (2010). 13 S. Woutersen and H. J. Bakker, Nature (London) 402, 507 (1999). 14 L. Piatkowski, K. B. Eisenthal, and H. J. Bakker, Phys. Chem. Chem. Phys. 11, 9033 (2009). 15 R. L. A. Timmer and H. J. Bakker, J. Phys. Chem. A 114, 4148 (2010). 16 M. L. Cowan, B. D. Bruner, N. Huse, J. R. Dwyer, B. Chugh, E. T. J. Nibbering, T. Elsaesser, and R. D. J. Miller, Nature (London) 434, 199 (2005). 2 P.
J. Chem. Phys. 140, 244503 (2014) 17 D.
Kraemer, M. L. Cowan, A. Paarmann, N. Huse, E. T. J. Nibbering, T. Elsaesser, and R. J. D. Miller, Proc. Natl. Acad. Sci. U.S.A. 105, 437 (2008). 18 W. Amir, G. Gallot, and F. Hache, J. Chem. Phys. 121, 7908 (2004). 19 K. Blum, Density Matrix Theory and Applications (Plenum Press, New York, 1980). 20 R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, New York, 2001). 21 M. Yang and G. R. Fleming, Chem. Phys. 282, 163 (2002). 22 A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234110 (2009). 23 S. Jang, J. Chem. Phys. 127, 174710 (2007). 24 D. E. Rosenfeld and M. D. Fayer, J. Chem. Phys. 137, 064109 (2012). 25 A. G. Redfield, IBM J. Res. Develop. 1, 19 (1957). 26 G. Lindblad, Commun. Math. Phys. 48, 119 (1976). 27 V. Gorini, A. Kossokowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976). 28 D. P. S. McCutcheon and A. Nazir, Phys. Rev. B 83, 165101 (2011). 29 D. P. S. McCutcheon and A. Nazir, J. Chem. Phys. 135, 114501 (2011). 30 A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234111 (2009). 31 S. Jang, Y.-C. Cheng, D. R. Reichman, and J. D. Eaves, J. Chem. Phys. 129, 101104 (2008). 32 P. Huo and D. F. Coker, J. Chem. Phys. 133, 184108 (2010). 33 E. N. Zimanyi and R. J. Silbey, Phil. Trans. R. Soc. A 370, 3620 (2012). 34 P. Rebentrost, R. Chakraborty, and A. Aspuru-Guzik, J. Chem. Phys. 131, 184102 (2009). 35 J. Roden, A. Eisfeld, W. Wolff, and W. T. Strunz, Phys. Rev. Lett. 103, 058301 (2009). 36 J. L. Skinner, B. M. Auer, and Y.-S. Lin, Adv. Chem. Phys. 142, 59 (2009). 37 P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Chem. Phys. 135, 044701 (2011). 38 Y. Ni, S. M. Gruenbaum, and J. L. Skinner, Proc. Natl. Acad. Sci. U.S.A. 110, 1992 (2013). 39 M. Yang, F. Li, and J. L. Skinner, J. Chem. Phys. 135, 164505 (2011). 40 J. A. Poulsen, G. Nyman, and S. Nordholm, J. Phys. Chem. A 107, 8420 (2003). 41 C. Backtorp, J. A. Poulsen, and G. Nyman, J. Phys. Chem. A 109, 3105 (2005). 42 F. Li and J. L. Skinner, J. Chem. Phys. 132, 204505 (2010). 43 F. Li and J. L. Skinner, J. Chem. Phys. 133, 244504 (2010). 44 F. Li and J. L. Skinner, J. Chem. Phys. 134, 099901 (2011). 45 L. Shi, S. M. Gruenbaum, and J. L. Skinner, J. Phys. Chem. B 116, 13821 (2012). 46 L. Shi and J. L. Skinner, J. Phys. Chem. B 117, 15536 (2013). 47 S. M. Gruenbaum, C. J. Tainter, L. Shi, Y. Ni, and J. L. Skinner, J. Chem. Theor. Comput. 9, 3109 (2013). 48 R. Kumar and J. L. Skinner, J. Phys. Chem. B 112, 8311 (2008). 49 C. J. Tainter, P. A. Pieniazek, Y.-S. Lin, and J. L. Skinner, J. Chem. Phys. 134, 184501 (2011). 50 P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Am. Chem. Soc. 133, 10360 (2011). 51 I. V. Stiopkin, C. Weeraman, P. A. Pieniazek, F. Y. Shalhout, J. L. Skinner, and A. V. Benderskii, Nature (London) 474, 192 (2011). 52 C. J. Tainter and J. L. Skinner, J. Chem. Phys. 137, 104304 (2012). 53 C. J. Tainter, Y. Ni, L. Shi, and J. L. Skinner, J. Phys. Chem. Lett. 4, 12 (2013). 54 D. van der Spoel, E. Lindahl, B. Hess, A. R. van Buuren, E. Apol, P. J. Meulenhoff, D. P. Tieleman, A. L. T. M. Sijbers, K. A. Feenstra, R. van Drunen et al., Gromacs User Manual Version 3.3, 2005, see www.gromacs.org. 55 J. Hayward and J. Reimers, J. Chem. Phys. 106, 1518 (1997). 56 V. F. Petrenko and R. W. Whitworth, Physics of Ice (Oxford University Press, Oxford, 1999). 57 J. W. Pyper and L. D. Christensen, J. Chem. Phys. 62, 2596 (1975). 58 J. C. Duplan, L. Mahi, and J. L. Brunet, Chem. Phys. Lett. 413, 400 (2005). 59 T. L. C. Jansen, B. M. Auer, M. Yang, and J. L. Skinner, J. Chem. Phys. 132, 224503 (2010). 60 M. Grover and R. Silbey, J. Chem. Phys. 54, 4843 (1971).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.63.180.147 On: Sun, 23 Nov 2014 03:59:01
