Molecular Theory of Solute-Pump/Solvent-Probe Spectroscopy and Application to Preferential Solvation Dynamics by Xiang Sun M. A., Brown University, 2010 B. Sc., University of Science and Technology of China, 2008 A Dissertation submitted in partial fulfillment of the requirements for the Degree of Doctor of Philosophy in the Department of Chemistry at Brown University Providence, Rhode Island May 2014 c Copyright 2014 by Xiang Sun This dissertation by Xiang Sun is accepted in its present form by the Department of Chemistry as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Richard M. Stratt, Director Recommended to the Graduate Council Date Jimmie D. Doll, Reader Date Christoph Rose-Petruck, Reader Approved by the Graduate Council Date Peter M. Weber, Dean of the Graduate School iii Curriculum Vitæ Born November 3, 1987, Xinxiang, Henan, China. Education Ph. D., Theoretical Chemistry, Brown University, Providence, RI, 2013. M. A., Chemistry, Brown University, Providence, RI, 2010. B. Sc., Chemical Physics, University of Science and Technology of China, China, 2008. Publications 1. X. Sun, R. M. Stratt, “How a solute-pump/solvent-probe spectroscopy can reveal structural dynamics: Polarizability response spectra as a two-dimensional solvation spectroscopy”, J. Chem. Phys. 139, 044506 (2013). 2. X. Sun, R. M. Stratt, “The molecular underpinnings of a solute-pump/solvent-probe spectroscopy: The theory of polarizability response spectra and an application to preferential solvation”, Phys. Chem. Chem. Phys. 14, 6320 (2012). 3. S. X. Tian, X. Sun, R. Cao, J. Yang, “Thermal Stabilities of the Microhydrated Zwitterionic Glycine: A Kinetics and Dynamics Study”, J. Phys. Chem. A 113, 480 (2009). Presentations 1. Invited talk, School of Chemistry and Materials Science, University of Science and Technology of China, Hefei, China, March 27, 2013. 2. Xiang Sun and Richard M. Stratt, “Theory of solute-pump/solvent-probe spectroscopy and an application to preferential solvation”, Symposium on solvent dynamics at biomolecular interfaces: experiment and theory, the 244th ACS national meeting, Philadelphia, PA, August 19-23, 2012. Honors and Awards Sigma Xi Award for Excellence in Research, Brown University, 2013. Dissertation Fellowship, Brown University, 2012. University Fellowship for First-year Students, Brown University, 2009. Excellent Graduate of University of Science and Technology of China, 2008. MEDY Scholarship, University of Science and Technology of China, 2007. iv Professional Affiliations American Chemical Society Sigma Xi v Dedicated to my parents vi Acknowledgements This dissertation would never have come to fruition without the support of many people. It is with happiness and gratitude that I acknowledge those who were instrumental in all kinds of ways to the success of this study and because of whom my graduate experience has been a great treasure that I will cherish forever. Foremost, I would like to express the utmost appreciation to my advisor, Richard Stratt, for his guidance, encouragement, patience and care during the entire course of my Ph.D. education. I feel fortunate to have such an opportunity of working with him. Rich is truly a role model as well as a great mentor for young scientists like myself. I have been amazed by his broad knowledge, enlightening insights and systematic investigating styles. Rich showed me how to think in a logical manor, how to grasp the point of an idea, how to ask right questions and how to present in a clear and enjoyable way. In addition, he has always been there when I have questions or problems. This work would not have been possible without his help, and I have to mention that he kindly read the thesis draft line by line and offered valuable revision suggestions including the use of articles. It gives me immense pleasure to declare gratitude to my thesis committee, Jimmie Doll and Christroph Rose-Petruck for reading through my thesis and providing thoughtful comments. I really appreciate that Jimmie considerately took part in the telephone meeting of my defense when he is off sick. I am extremely grateful to Jimmie for his constant interest in my research and constructive suggestions for all the years I spend here at Brown. I also would like to thank Christoph for his invaluable questions and scientific merits. I vii would like to extend my gratitude especially to Matthew Zimmt for serving on my original research propositional committee. He has always been very supportive, enthusiastic and willing to help me to learn new things. I gratefully acknowledge Norbert Scherer, David Blank and their co-workers for imple- menting the idea of solute-pump/solvent-probe spectroscopy experimentally, and especially I would like to thank Margaret Hershberger for sharing her experience and knowledge of the spectroscopy and much more beyond. I am also deeply grateful to Graham Fleming, Nancy Levinger, Richard Saykally, Kevin Kubarych and David Birch for their brilliant comments and helpful remarks from experimentalists’ perspective. My sincere gratitude and warm appreciation goes to distinguished theorist Branka Ladanyi who collaborated with our group long termly and recently hosted my visit to the beautiful Colorado State University. I was impressed by the progress that we made on our molecular liquid simulation in such a limited time, and I believe that is due to the fact that Branka always put herself in my place and hence gave me spot-on advice. I must say a big thank you to Guoqing Zhang for inviting Rich and me to deliver lectures at my Alma Mater, University of Science and Technology of China, early this year. I appreciate very much his lasting inspiration, encouragement and friendship in the past decade. My acknowledgement will never be complete without the special mention of our group seniors who have introduced me the basic skill set for theoretical chemists. I am especially indebted to Crystal Nguyen for her helps on my first molecular dynamics simulation and furthermore the insights on preferential solvation and Baofeng Zhang for bouncing ideas with me. I also would like to thank our group members, Xiao Liang, Qingqing Ma, Beth Shimmyo, Joe Isaacson, Dan Jacobson, Vale Cofer-Shabica and Gengqi Wang for their intriguing discussions on research and stuff. It has been a great privilege to spend best several years in the Department of Chemistry at Brown University. Its members will always remain dear to me. Besides those I have already mentioned above, I would like to extend my appreciation to Margaret Doll for viii maintaining our computing resources. I also want to thank James Baird for teaching me quantum mechanics and lending me books. I benefited from the Edward Mason collection of books in room GC246 as well. I am truly grateful to our librarian Lee Pedersen who helped me finding books and journal articles that I needed, particularly for purchasing a copy of Shaul Mukamel’s nonlinear spectroscopy book on my request. Life would not have been as colorful without the many wonderful friends I met in this pretty little state of Rhode Island. Special thanks also go to my lifelong friends scattered all over the world. No amount of word is enough to express my regards to them who make me understand the truth of life. I feel blessed to have such friends. With you, I share some lyrics from song “I was here” by Beyonc´e: I was here. I lived, I loved. Above and beyond all, my heartfelt gratitude goes to my parents for their much unwavering supports and understanding in every possible way. Their endless love and everything made all this possible. Xiang Sun December 2013 Providence, RI ix Contents List of Tables xiii List of Figures xviii List of Symbols xxii List of Abbreviations xxiv 1 Introduction 1 1.1 Liquid Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Theoretical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Instantaneous-normal-mode theory . . . . . . . . . . . . . . . . . 3 1.2.2 Time correlation function . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.3 Statistical mechanics and spectroscopy . . . . . . . . . . . . . . . 7 1.2.4 Langevin equation and generalized Langevin equation . . . . . . . 9 1.2.5 Linear response theory . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.6 Computer simulations . . . . . . . . . . . . . . . . . . . . . . . . 19 1.3 Experimental Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.3.1 Time-dependent fluorescence spectroscopy . . . . . . . . . . . . . 23 1.3.2 Transient absorption . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.3.3 Scattering experiments . . . . . . . . . . . . . . . . . . . . . . . . 26 1.3.4 Coherent nonlinear spectroscopy . . . . . . . . . . . . . . . . . . . 29 1.4 Design of Solute-Pump/Solvent-Probe Spectroscopy . . . . . . . . . . . . 42 2 Linear Response Theory 48 2.1 Preparations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.2 Traditional Linear Response Theory . . . . . . . . . . . . . . . . . . . . . 52 2.2.1 Average on excited state . . . . . . . . . . . . . . . . . . . . . . . 53 x 2.2.2 Average on ground state (Schr¨odinger picture) . . . . . . . . . . . 55 2.2.3 Average on ground state (Heisenberg picture) . . . . . . . . . . . . 58 2.2.4 Application to four-wave-mixing light scattering . . . . . . . . . . 60 2.3 Gaussian Statistics and Linear Response Theory . . . . . . . . . . . . . . . 62 2.3.1 Gaussian distribution . . . . . . . . . . . . . . . . . . . . . . . . . 63 2.3.2 Generating function . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.3.3 Linear response . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 2.4 Application to Solute-Pump/Solvent-Probe Spectroscopy . . . . . . . . . . 67 3 Preferential Solvation Dynamics and Optical Kerr Effect Spectroscopy 77 3.1 Preferential Solvation Dynamics . . . . . . . . . . . . . . . . . . . . . . . 77 3.2 Molecular Dynamics Simulation . . . . . . . . . . . . . . . . . . . . . . . 80 3.3 Optical Kerr Effect (OKE) Spectroscopy . . . . . . . . . . . . . . . . . . . 86 3.3.1 OKE experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.3.2 Many-body polarizability (dipole-induced-dipole model) . . . . . . 93 3.3.3 OKE spectra for our model . . . . . . . . . . . . . . . . . . . . . . 100 3.3.4 Difference spectra . . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.4 Projection Operator Analysis of Molecular Contributions . . . . . . . . . . 112 3.5 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 4 Nonequilibrium Two-Dimensional Solute-Pump/Solvent-Probe Spectroscopy 139 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 4.2 Instantaneous-Normal-Mode (INM) Theory . . . . . . . . . . . . . . . . . 141 4.2.1 Mass-weighted INM . . . . . . . . . . . . . . . . . . . . . . . . . 142 4.2.2 INM density of states in liquid argon . . . . . . . . . . . . . . . . 145 4.2.3 INM influence spectrum . . . . . . . . . . . . . . . . . . . . . . . 154 4.2.4 INM approach to OKE spectra . . . . . . . . . . . . . . . . . . . . 156 4.2.5 INM solvation influence spectrum . . . . . . . . . . . . . . . . . . 167 4.3 Evaluation of Solute-Pump/Solvent-Probe Response . . . . . . . . . . . . . 167 4.3.1 Hybrid INM/MD method . . . . . . . . . . . . . . . . . . . . . . . 169 4.3.2 SHO-like approximation . . . . . . . . . . . . . . . . . . . . . . . 172 4.3.3 Non-INM numerical evaluation . . . . . . . . . . . . . . . . . . . 177 4.4 Short-time Expansion of Solute-Pump/Solvent-Probe Response . . . . . . . 179 4.4.1 Exact expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 4.4.2 INM approximation . . . . . . . . . . . . . . . . . . . . . . . . . 182 xi 4.4.3 SHO-like approximation . . . . . . . . . . . . . . . . . . . . . . . 186 4.5 2D Solute-Pump/Solvent-Probe Spectra for Our Model . . . . . . . . . . . 187 4.6 Structural Information in Preferential Solvation . . . . . . . . . . . . . . . 192 4.6.1 First-shell population dynamics . . . . . . . . . . . . . . . . . . . 192 4.6.2 Bond-angle distribution . . . . . . . . . . . . . . . . . . . . . . . 195 4.7 Discussions on Analytical Approximations for 2D Spectroscopic Responses 219 5 Concluding Remarks 224 5.1 Future Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230 Appendix A Poisson Bracket 234 Appendix B Derivation of Linear Contribution to exp(Aˆ + B) ˆ 237 Appendix C Rotational Average 240 Appendix D Simulation Details 246 Appendix E Spatial Derivative of the Many-body Polarizability 253 Appendix F Snapshot of Initial OKE Response 261 Bibliography 266 Index 281 xii List of Tables 1.1 Light–matter interaction Hamiltonian Hint = −Vˆ · F(r,t) with different choices of external field F(r,t) and coupling molecular property Vˆ . . . . . 41 3.1 Molecular polarizability and dipole of common solvents. . . . . . . . . . . 79 3.2 Averaged first shell population and 1/e relaxation time (ps). . . . . . . . . . 86 3.3 Fitting results for the long-time diffusive portion of the OKE response function in 10% S and 50% S mixtures on the equilibrium ground and excited state. The OKE response function is calculated by averaging 107 infinite order DID polarizabilities, sampled every 5 time steps. . . . . . . . 103 3.4 Fitting results for the long-time diffusive portion of the OKE response function in 10% S and 50% S mixtures on the equilibrium ground and excited state. The OKE response function is calculated by averaging 108 first order DID polarizabilities, sampled every 10 time steps. . . . . . . . . 105 3.5 Summary of single snapshot analysis of kB T R(0) and its largest term(s) (unit in σ 5 · (ε /m)1/2 ) of each category with “u” for the solute, “v” for a solvent. s and w stand for strong solvent and weak solvent with sub- scripts j, k, l to distinguish different solvent atoms. The same ground state equilibrium configuration is used in both 10% S and 50% S, and different excited equilibrium configurations are used in 10% S and 50% S. Details are included in Appendix F. . . . . . . . . . . . . . . . . . . . . . . . . . . 130 4.1 Compositions of first-shell/first-shell, first-shell/second-shell and second- shell/second-shell contributions to f cc and hcp lattice bond-angles cos θ =0.5 and 0.87 with strong solvent percentage p = 0.1, 0.5 and 1. . . . . . . . . . 215 xiii List of Figures 1.1 Instantaneous-normal-mode spectrum (density of states) of liquid argon, simulated at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80. . . . . . . 5 1.2 Schematic illustration of time-dependent fluorescence spectroscopy (time- resolved fluorescence Stokes shift). . . . . . . . . . . . . . . . . . . . . . . 16 1.3 Experimental Sν (t) and simulated C∆E (t) solvation response functions for Coumarin 343 in water. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4 Rotational energy relaxation of CN in liquid Ar. . . . . . . . . . . . . . . . 18 1.5 Time scales relevant to physical, chemical and biological processes. . . . . 22 1.6 Molecular structure of coumarin 153. . . . . . . . . . . . . . . . . . . . . 23 1.7 Solvation response functions of coumarin 153 in acetonitrile/benzene mix- ture measured by time-dependent fluorescence spectroscopy. . . . . . . . . 25 1.8 Double-sided Feynman diagrams for time-domain four-wave mixing in a two-level system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.9 (a) Comparison of experimental and calculated 3-pulse photon echo peak shift for tricarbocyanine dye IR144 in acetonitrile at 297K; (b) Comparison of the instantaneous normal mode solvation spectrum of acetonitrile, and the rescaled inertial spectral density calculated with parameters fitting experimental result. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 1.10 Illustration of solute-pump/solvent-probe spectroscopy for preferential sol- vation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 1.11 Experimental sequence of laser pulses in resonant-pump polarizability- response spectroscopy (RP-PORS). . . . . . . . . . . . . . . . . . . . . . . 44 1.12 (a) Two-dimensional anisotropic resonant-pump polarizability-response spec- trum (RP-PORS) response function for coumarin 153 in acetonitrile; (b) Anisotropic transient solvation polarizability measurement for coumarin 153 in acetonitrile, which is the derivative of RP-PORS spectrum with respect to solvation time T . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.1 Preparation of nonequilibrium state in response to an external perturbation. 53 2.2 Sequence of polarizability interactions in four-wave mixing. . . . . . . . . 61 2.3 Sequence of excitations in solute-pump/solvent-probe spectroscopy. . . . . 69 2.4 Schematic illustration of pulse sequence and liquid structural evolution in solute-pump/solvent-probe spectroscopy. . . . . . . . . . . . . . . . . . . . 72 xiv 3.1 Solute–solvent radial distribution function in atomic liquid with ground- state solute, 10% S system with excited-state solute, and 50% S system with excited-state solute. . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.2 Velocity autocorrelation function for Lennard–Jones liquid at density ρσ 3 = 0.8 and temperature kB T /ε = 1.00. . . . . . . . . . . . . . . . . . . . . . . 83 3.3 Solvation responses for the atomic liquid preferential solvation model in a range of solvent compositions. Shown here are both the equilibrium solva- tion correlation function C(t) and the nonequilibrium solvation relaxation function S(t), with the detailed view of subpicosecond relaxations in the inset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.4 Time evolution of the first shell population in 10% S solvent mixture after solute excitation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.5 Schematic illustration of optical Kerr effect (OKE) spectroscopy. . . . . . . 88 3.6 Optical Kerr effect spectra for 2,4,6-trifluoropyridine at 306 K. . . . . . . . 92 3.7 Anisotropic and isotropic OKE responses in 10% S and 50% S systems on equilibrated ground and excited states. . . . . . . . . . . . . . . . . . . . . 101 3.8 Time-domain OKE response functions for two different solvent mixtures, 10% S and 50% S systems in the equilibrated ground state and excited state. 102 3.9 Frequency-domain OKE spectra for 10% S and 50% S systems. For each system, shown are the results with equilibrated ground-state (g) solute, excited-state (e) solute and the difference between the two (e-g). . . . . . . 104 3.10 Comparison of the OKE reduced spectral densities for 10% S and 50% S solvent mixtures with ground-state solute (g), excited-state solute (e) and their difference spectrum (e-g). . . . . . . . . . . . . . . . . . . . . . . . 106 3.11 Comparison of the long waiting time limit anisotropic RP-PORS spectra for 10% S and 50% S solvent mixtures. . . . . . . . . . . . . . . . . . . . 107 3.12 Comparison of OKE spectra in 10% S system, (a) infinite order DID approximation and full solute polarizability, (b) infinite order DID approx- imation with the solute’s polarizability set to zero and (c) first order DID approximation with full solute polarizability. . . . . . . . . . . . . . . . . . 108 3.13 Comparison of OKE spectra in 50% S system, (a) infinite order DID approximation and full solute polarizability, (b) infinite order DID approx- imation with the solute’s polarizability set to zero and (c) first order DID approximation with full solute polarizability. . . . . . . . . . . . . . . . . . 109 3.14 Comparison of different polarizability contributions to the long waiting limit anisotropic RP-PORS spectra for 10% S and 50% S solvent mixtures. 111 3.15 Polarizability-velocity correlation function G(t) projected to first shell and everything else contributions. . . . . . . . . . . . . . . . . . . . . . . . . . 117 3.16 Polarizability-velocity correlation function G(t) projected to first shell, second shell and everything else contributions. . . . . . . . . . . . . . . . . 118 3.17 Projection of the OKE responses for 10% S and 50% S solvent mixtures into first-shell, outer-shells and cross contributions. . . . . . . . . . . . . . 121 3.18 Projection of the OKE spectral densities for 10% S and 50% S solvent mixtures into first-shell, outer-shells and cross contributions. . . . . . . . . 122 xv 3.19 Projection of the e–g difference in spectral densities for 10% S and 50% S solvent mixtures into first-shell, outer-shells and cross contributions. . . . . 123 3.20 Projection of the OKE reduced spectral densities for 10% S and 50% S solvent mixtures into first-shell and outer-shells. . . . . . . . . . . . . . . 124 3.21 Projection of the long waiting time limit anisotropic solute-pump/solvent- probe spectra for 10% S and 50% S solvent mixtures into first-shell and outer-shells contributions. . . . . . . . . . . . . . . . . . . . . . . . . . . 125 3.22 Projection of the OKE spectral densities for 10% S and 50% S solvent mixtures into first-shell, second-shell and outer-shells contributions. . . . . 127 3.23 Projection of the e–g difference in spectral densities for 10% S and 50% S solvent mixtures into first-shell, second-shell and outer-shells contributions. 128 3.24 Illustration of how interaction-induced polarizabilities are affected by num- ber of spectroscopic bright solvent, leading to sign changes in the solute- pump/solvent-probe spectra. . . . . . . . . . . . . . . . . . . . . . . . . . 135 3.25 Rotational averaging of a solute–solvent–solvent equilateral triangle. . . . . 137 4.1 Instantaneous-normal-mode spectrum (density of states) of liquid Ar at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80. . . . . . . . . . . . . . 147 4.2 Instantaneous-normal-mode spectrum (density of states) of liquid Ar at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80 with and without zero frequencies removed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 4.3 Distribution of eigenvalues of the dynamical matrix for liquid Ar at tem- perature kB T /ε = 1.00 and density ρσ 3 = 0.80. . . . . . . . . . . . . . . . 149 4.4 Comparison of the instantaneous-normal-mode spectra of liquid Ar using frequency histogram and eigenvalue histogram. . . . . . . . . . . . . . . . 150 4.5 Instantaneous-normal-mode spectra (density of states) of (a) liquid argon, 10% S on the ground and the excited state, and (b) 10% S vs 50% S on the ground and the excited state. . . . . . . . . . . . . . . . . . . . . . . . . . 153 4.6 Instantaneous-normal-mode polarizability influence spectra of (a) 10% S and (b) 50% S on the ground state and on the excited state at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80. . . . . . . . . . . . . . . . . . . . 157 4.7 Instantaneous-normal-mode OKE influence spectra of (a) 10% S and (b) 50% S on the ground state and on the excited state at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80. . . . . . . . . . . . . . . . . . . . . . . . . . 162 4.8 Comparison of instantaneous-normal-mode and exact molecular-dynamics predictions of OKE spectra for the 10% S solvent mixture. . . . . . . . . . 163 4.9 Comparison of e-g difference of the OKE spectra using instantaneous- normal-mode approximation and exact molecular dynamics results for 10% S and 50% S solvent mixtures. . . . . . . . . . . . . . . . . . . . . . . . . 164 4.10 Comparison of instantaneous-normal-mode and exact molecular-dynamics predictions of OKE spectra for the 10% S solvent mixture. The upper panel is the OKE spectrum of the system with an excited-state solute and the lower panel is the difference spectrum between the excited- and ground- state-solute OKE spectra. . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 xvi 4.11 Comparison of instantaneous-normal-mode and exact molecular-dynamics predictions of e–g difference of OKE spectra for the 10% S solvent mixture with an atomic solute. Different system sizes and numbers of configura- tions for averaging are compared. . . . . . . . . . . . . . . . . . . . . . . . 166 4.12 Instantaneous-normal-mode polarizability influence spectrum vs. solvation influence spectrum of 10% S solvent mixture with an excited-state solute. . 168 4.13 Comparison of response functions using SHO-like approximation at T = 0 ps, T = 50 ps, long T with the equilibrium excited state OKE response function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 4.14 Comparison of full DID and first order DID response functions using SHO- like approximation at T = 0 ps. . . . . . . . . . . . . . . . . . . . . . . . 177 4.15 Two-dimensional solute-pump/solvent-probe spectra for our preferential solvation system 10% S solvent mixtures. . . . . . . . . . . . . . . . . . . 189 4.16 Two-dimensional solute-pump/solvent-probe spectra for our preferential solvation system 50% S solvent mixtures. . . . . . . . . . . . . . . . . . . 190 4.17 Comparison of the fluctuations of solute-pump/solvent-probe responses at a finite waiting time T with that of the equilibrium e–g difference in INM OKE influence spectra for our preferential solvation systems 10% S and 50% S. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 4.18 The principal result of this chapter: two-dimensional solute-pump/solvent- probe spectra for our 10% S preferential solvation system. The 2D spectra show how the solution’s intermolecular vibrational spectrum evolves with increasing delay time T after the solute excitation. . . . . . . . . . . . . . . 196 4.19 Most revealing analysis of this chapter: comparison of the relaxation profile observed by solute-pump/solvent-probe spectra with structure- and potential-energy-sensitive equilibrium solvation correlation functions for our 10% S preferential solvation system. . . . . . . . . . . . . . . . . . . . 197 4.20 Evolution of strong solvent population distributions in the first solvation shell and in the second solvation shell for both 10% S solvent mixture and 50% S solvent mixture following the resonant solute excitation. . . . . . . . 198 4.21 Comparison of the structure- and potential-energy-sensitive equilibrium solvation correlation functions for our 50% S preferential solvation system. 199 4.22 Two-dimensional solute-pump/solvent-probe spectra for our polarizability- switched 10% S preferential solvation system. . . . . . . . . . . . . . . . . 200 4.23 Illustration of the bond angle between strong solvent pair and the solute. . . 201 4.24 T -dependent angular distributions of strong solvent (a) in the first shell, (b) in the second shell and (c) in the first and second shells of 10% S solvent mixture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 4.25 T -dependent angular distributions of strong solvent (a) in the first shell and (b) in the first and second shell of 50% S solvent mixture. . . . . . . . . . . 204 4.26 Comparison of first-shell population relaxation, solvation energy relaxation and the bond-angle relaxations in the 10% S system. . . . . . . . . . . . . 205 4.27 Illustration of the f cc and hcp close packing. . . . . . . . . . . . . . . . . 206 xvii 4.28 Angular distributions for liquid-equivalent f cc lattices for atoms (a) in the first shell, (b) in the second shell and (c) in the first and second shells taken together. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208 4.29 Angular distributions for liquid-equivalent hcp lattices for atoms (a) in the first shell, (b) in the second shell and (c) in the first and second shells taken together. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209 4.30 Angular distributions for liquid-equivalent f cc lattices for atoms (a) in the first shell, (b) in the second shell and (c) in the first and second shells with varying bright solvent percentages p = 0.1, 0.5 and 1. . . . . . . . . . . . . 211 4.31 Angular distributions for liquid-equivalent hcp lattices for atoms (a) in the first shell, (b) in the second shell and (c) in the first and second shells with varying bright solvent percentages p = 0.1, 0.5 and 1. . . . . . . . . . . . . 212 4.32 Comparison of T -dependent angular distributions in 10% S solvent mixture and angular distributions of liquid-equivalent f cc and hcp lattices in the first shell and in the second shell. . . . . . . . . . . . . . . . . . . . . . . . 214 4.33 Comparison of T -dependent angular distributions in 10% S solvent mixture and angular distributions of liquid-equivalent f cc and hcp lattices in the first and second shells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 215 4.34 Enlarged view of comparison of T -dependent angular distributions in 10% S mixture and angular distributions of liquid-equivalent f cc and hcp lat- tices in the first shell, in the second shell and in the first and second shells. . 216 4.35 Comparison of T -dependent angular distributions in 50%S mixture and angular distribution of liquid-equivalent f cc and hcp lattices in the first shell and in the first and second shells. . . . . . . . . . . . . . . . . . . . . 217 4.36 Schematic illustration of the solvent structural change induced by the solute excitation in a preferential solvation system. . . . . . . . . . . . . . . . . . 218 C.1 Euler angle φ , θ , χ relating the space-fixed XYZ and molecule-fixed xyz frames. φ = 0 ∼ 2π , θ = 0 ∼ π , χ = 0 ∼ 2π . . . . . . . . . . . . . . . . . 241 D.1 Order parameter of the ground-state simulated Lennard-Jones liquid (N = 256). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 248 E.1 Comparison of numerical and analytical OKE response functions for the ground-state 50% S mixture. . . . . . . . . . . . . . . . . . . . . . . . . . 259 E.2 Comparison of numerical and analytical OKE response functions for the ground-state 50% S mixture. Both first order and infinite order DID approximations used to calculate the many-body polarizability and the spatial derivative of the many-body polarizability are compared in the same plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 260 xviii List of Symbols Symbol Meaning i, j, k, . . . molecule indices µ, ν, γ, δ , . . . Cartesian coordinate indices x, y, z A, B, C, . . . vectors, matrices, or tensors N number of particles R coordinate, liquid configuration r single molecular/atomic position P momentum or polarization v single molecular/atomic velocity X phase space point = (R, P) A(X(t)) dynamical variable ρ (X,t) phase space distribution hAi ensemble average of A A nonequilibrium average or rotational average Aˆ operator A˙ time derivative of A {A, B} classical Poisson bracket [A, B] quantum commutator δA fluctuation of A σ (A) standard deviation of A β (kB T )−1 H Hamiltonian V potential energy µ molecular dipole moment α isolated molecular polarizability xix α isotropic molecular polarizability γ anisotropy of molecular polarizability π effective molecular polarizability Π many-body polarizability Π element of many-body polarizability m mass of atom or molecule ε energy well depth of Lennard-Jones potential σ length parameter of Lennard-Jones potential ˆ Ω molecular orientational unit vector Pˆ X projection operator Z mass-weighted coordinate D(R) dynamical matrix D(ω ) density of states or Fourier transform of OKE response function R(t, . . .) response function G(t) spectroscopic “velocity” autocorrelation function C(t) equilibrium solvation relaxation or correlation function S(t) nonequilibrium solvation relaxation ∆R(0, T, T + t) solute-pump/solvent-probe response ∆R(ω , T ) solute-pump/solvent-probe spectra S(T ; ω ) solute-pump/solvent-probe spectral relaxation qα the α -th instantaneous normal mode (INM) pα the conjugate momentum of the α -th INM ρA (ω ) influence spectrum cα coupling coefficient with respect to the α -th INM χ (t1 ,t2) fundamental Poisson bracket E(r,t) electric field χ (n) nth-order electric susceptibility k wave vector ω light frequency or vibrational frequency ρ (t) density matrix ρ (cos θ ) bond-angle distribution ρ (k) translational order parameter xx g(r) radial distribution function L Liouville operator F Fourier transform ∆V potential energy gap = Ve −Vg F(t) external perturbing field x n-dimensional Gaussian variable δt MD time step g ground state of the solute e excited state of the solute u solute v solvent S strong solvent W weak solvent xxi List of Abbreviations Abbreviation Meaning Introduced page 2D IR two-dimensional infrared spectroscopy 41 3PEPS 3-pulse photon echo peak shift 38 C153 coumarin 153 23 DID dipole-induced-dipole 93 DMSO dimethyl sulfoxide 79 DRSK DeVane-Ridley-Space-Keyes approximation 222 FWHM full width at half maximum 231 FWM four-wave mixing 30 GLE generalized Langevin equation 11 INM instantaneous normal mode 3 LJ Lennard-Jones 80 LO local oscillator 89 LRT linear response theory 9 MC Monte Carlo 20 MD molecular dynamics 20 OHD optical heterodyne detection 36 OKE optical Kerr effect spectroscopy 30 RaPTORS resonant-pump third-order Raman spectroscopy 228 RIKES Raman-induced Kerr effect spectroscopy 37 RP-PORS resonant-pump polarizablity response spectroscopy 44 RSD reduced spectral density 91 SD spectral density 7 SHO simple harmonic oscillator 172 TCF time-correlation function 6 xxii TDF time-dependent fluorescence 6 TG transient grating 30 TRFSS time-resolved fluorescence Stokes shift 15 TSP transient solvation polarizability spectrum 46 xxiii Abstract of “Molecular Theory of Solute-Pump/Solvent-Probe Spectroscopy and Ap- plication to Preferential Solvation Dynamics” by Xiang Sun, Ph.D., Brown University, May 2014 Common resonant spectroscopic methods used to study dynamics in solutions, such as time-dependent fluorescence, share a feature of probing the solute directly, or more pre- cisely the solute-solvent interaction energy. One has to infer how the solvents move from the solute’s time-dependent information. By contrast, nonresonant light scattering experiments report on the dynamics of a liquid as a whole, but cannot concentrate on dynamics of any local portion of the solution. A recently demonstrated two-dimensional solute-pump/solvent-probe spectroscopy, a combination of the two approaches mentioned above, enables us to follow the nonequilibrium dynamics of solvents after the solute’s electronic excitation. This dissertation is a theoretical attempt at understanding the molec- ular information behind this kind of spectroscopy. After developing the general lin- ear response theory for these spectra using classical statistical mechanics, I apply the resulting formalism to a preferential solvation model system consisting of an atomic solute dissolved in an atomic-liquid mixture. In the experimentally interesting limit of long solute-pump/solvent-probe delays, the spectra become the differences in light- scattering spectra between solutions with equilibrated ground- and excited-state solutes. The drastically distinctive spectra for various solvents in this limit suggest how changing liquid structure affects intermolecular liquid dynamics and how local a portion of the solvent dynamics can be accessed by the spectra. For the more general nonequilibrium case of the spectra with finite solute-pump/solvent-probe delays, a practical hybrid calculation method combining instantaneous-normal-mode theory with molecular dynamics shows a great advantage in dealing with two-dimensional spectroscopies especially with separated time scales. The full two-dimensional spectra can serve as a solvation spectroscopy capable of distinguishing the structural and energetic solvation dynamics. Calculations of our preferential solvation model indicate that the spectra indeed display the same relaxation profile as the local solvent population changes, which is measurably different from the solute-solvent interaction energetic relaxation measured by time-dependent fluorescence. Thus the two-dimensional spectroscopy effectively singles out structural dynamics of local solvents around the solute. Chapter 1 Introduction 1.1 Liquid Dynamics From glassware in the laboratory to warm and wet interior of living cells, a tremendous number of important chemical reactions take place in liquid phases. Like all condensed- matter phases, the structure and dynamics of a liquid are drastically different from the gas phase — in the gas phase, each molecule is nearly isolated and moves freely most of the time and only collides with another molecule much more infrequently than in a condensed phase. Unlike other condensed phases, liquids possess the highest translational and rotational symmetries,1 i.e. homogeneity and isotropy. In other words, a liquid is so random that any direction or any location is equally important. These non-broken symmetries might make one believe that liquids are simple at first glance. As a matter of fact, liquid dynamics could be complicated especially from a microscopic perspective. Liquids are highly disordered and dense, so we cannot investigate liquid dynamics with the established approaches that are used to treat gases and solids. For example, molecular motions in crystals are only small oscillations about their equilibrium positions. Oscillations in solids can be described as combinations of collective motions of the entire crystal. These collective motions are called phonons and are independent from each other: 1 they can be just as simple as individual molecules in the gas phase. Furthermore, the structure of a liquid can vary over a broad range of temporal and spatial scales, which makes the problem of liquid dynamics more intricate. Molecules in liquids constantly interact with one another through translation, rotation and vibration. These elementary motions are necessary for a chemical reaction to occur. The ultimate question that one would ask is exactly how these microscopic processes take place. In particular, how is the dynamics of these degrees of freedom affected by the interactions between molecules, the liquid structures and external perturbations? For ex- ample in a solution, there are two typical fundamental processes arising between a molecule and its surroundings: molecules exchange energy and rearrange their positions.2 The particular form of energy exchange relevant to how a vibrationally hot species dissipates its excess energy to other molecules or to other degrees of freedom within itself is called vibrational relaxation. An essential concept in reaction dynamics named Intramolecular Vibrational-energy Redistribution (IVR)3, 4 describes energy transfers from the initially populated vibrational state to other vibrational states inside a polyatomic molecule. The second of the two fundamental processes is known as solvation, during which a group of ambient solvent molecules reorganize their positions in order to accommodate a newly excited solute. Solvation is primarily related to the energetics of solute-solvent interaction and the intermolecular motions involving both solute and solvent. We will get back to the topic of solvation later in this thesis, but now let us have an overview of theoretical and experimental methods that have been applied to the field of liquid dynamics. 1.2 Theoretical Methods Owing to the complexity of liquid dynamics on different spatial and temporal scales, there is no universal theory that can solve all the problems related to the dynamics of liquid. Yet a number of theories have flourished in the last few decades intended to tackle 2 different aspects of liquid dynamics. An excellent example is the Instantaneous-Normal- Mode (INM) theory, which is developed in the early 1990s by Stratt5–9 and Keyes10, 11 independently. This theory offers an accurate picture of the short-time dynamics of liquid at the molecular level. The notion of INM originates from the phonon concept of solid- state matter.9 Phonons are well defined in solid — independent collective oscillations of the entire crystal — that require the intermolecular forces are harmonic, in other words, atoms can be regarded as being held together by springs. But the spring analogy does not apply to liquid in that the arrangement of atoms or molecules in liquid is not stable. In a stable arrangement like solid, only small oscillations of atoms around some fixed positions are allowed. Atoms in liquid diffuse and the whole liquid flow, so in conventional sense, phonons do not exist in liquid, at least persistently. Nevertheless when looking at a time scale that is so short during which atoms do not move far from its original positions, comparable to the dislocations of atoms in solid, the normal-mode concept is still applicable in dealing with liquid dynamics. 1.2.1 Instantaneous-normal-mode theory The instantaneous-normal-mode approach is based on the normal-mode analysis on an instantaneous liquid configuration R0 (a 3N-dimensional coordinate vector in an N-atom system, for example).9 Suppose one would like to know what is the configuration at a short time t later, Rt . We can expand the total potential energy V (R) as a function of the config- urational difference between two times to the second order (harmonic approximation). 1 V (Rt ) = V (R0 ) + ∇V (R0 ) · (Rt − R0 ) + (Rt − R0 ) · ∇∇V (R0 ) · (Rt − R0 ), (1.1) 2 ∂ where gradient operator ∇ ≡ ∂R. Although the first derivative of the potential energy does not vanish like in the solid case, one can still perform the normal-mode analysis by diagonalizing the 3N × 3N dynamical Hessian matrix ∇∇V (R0 ). 3 U(R0 ) · ∇∇V (R0 ) · UT (R0 ) = m · diag{ω12 , . . . , ω3N 2 }. (1.2) In the equation above, m is the mass of an atom, U(R0 ) is the unitary matrix that diagonal- izes the dynamical Hessian matrix and ωα (α = 1, . . ., 3N) are frequencies corresponding to every INM. Now the instantaneous normal modes are qα (t; R0) = [U(R0 ) · (Rt − R0 )]α , (α = 1, . . ., 3N). (1.3) Among 3N instantaneous normal modes, each mode corresponds to an independent collective harmonic motion. The independency warrants a limpid physical picture of liquid dynamics, a collection of uncoupled harmonic oscillators. By rewriting the above relation, it gets clearer that the dynamics of a liquid at short times is determined by these harmonic modes and their initial conditions which include the initial instantaneous forces and initial velocities. In other words, the time evolution of the whole liquid can be predicted based on the dynamical information at time 0 alone, as long as the time interval is short enough. It is worth noting that most of the harmonic modes barely finish a single cycle within such a short time, so full oscillations especially low-frequency ones are rare in liquid. Rt = R0 + UT (R0 ) · q(t; R0). (1.4) The instantaneous-normal-mode spectrum is defined as the distribution of INM fre- quencies. For example, in an atomic liquid with total N atoms, the density of states * + 1 3N 3N α∑ D(ω ) = δ (ω − ωα ) . (1.5) =1 Figure 1.1 shows the density of states of liquid argon simulated by molecular dynamics, in which the imaginary frequencies are plotted on the negative frequency axis. These imaginary frequencies corresponds to unstable instantaneous normal modes. And it is these 4 modes with imaginary frequencies that differentiate liquids from solids — no imaginary frequency in solids. An imaginary-frequency INM implies a negative local curvature on the potential energy surface,5 or a local maximum along the direction of this INM. By contrast, a real-frequency INM indicates a positive local curvature on the potential energy surface (consider a local minimum). 0.020 Liquid Ar 0.015 k T/ =1.00 ) (cm) B 3 =0.80 0.010 D( 0.005 0.000 -50 0 50 100 150 -1 /2 c (cm ) Figure 1.1 Instantaneous–normal–mode spectrum (density of states) of liquid argon. Molecular dynamics simulation with 256 atoms is performed at temperature kB T /ε = 1.00, density ρσ 3 = 0.80. 20,000 instantaneous liquid configurations are sampled. The imaginary frequencies are plotted on the negative frequency axis. Three translational modes with zero frequencies are removed. The real power of the instantaneous-normal-mode theory is not merely an extrapolation of liquid dynamics from any instant to an immediate future, but a new statistical mechanical point of view to inspect the dynamics in liquids. An instantaneous-normal-mode analysis on a single liquid configuration will not tell too much about a macroscopic system, but what is more practical is the ensemble average through sampling an equilibrium distribu- tion. And the ensemble average connects the most microscopic events with macroscopic measurable quantities. Besides the density of states spectrum, instantaneous-normal-mode 5 theory enables one to calculate ensemble-averaged time correlation functions (TCFs). 1.2.2 Time correlation function Time correlation functions play a central role in nonequilibrium statistical mechanics, as the partition function does in equilibrium statistical mechanics.12 Many properties of systems out of equilibrium are determined by various time correlation functions. For example, diffusion, viscosity and thermal conductivity can be described by time correlation functions. Suppose we look at some time-dependent quantity A(t), and the time correlation function of this quantity is hA(0)A(t)i, (h· · · i denotes the ensemble average) R dX A(X)A(X(t))e−β H(X) hA(0)A(t)i ≡ R , (1.6) dXe−β H(X) where X = (R, P) is the phase space point, R is the position coordinate and P is the con- jugate momentum. The time-dependent quantity A is often termed a dynamical variable13 that is a function of the phase space, A(X). A dynamical variable A with initial condition at phase space point X is written as A = A(0) = A(X) = A(X(0)). And A(t) = A(X(t)) denotes the value of the dynamical variable at time t, evolved from the initial state X. The R integration dX is over the whole of phase space. In addition, β = 1/kB T and H(X) is the Hamiltonian of the system equal to the total energy, which governs the time evolution of the whole system. The equations of motion of a classical system are the Hamilton’s equations:14(the overdot denotes time derivative d/dt throughout the thesis.) ˙ = ∂H , R and P˙ = − ∂H . (1.7) ∂P ∂R The probability density of finding the system at state around phase point X is e−β H(X) ρeq (X) = R . (1.8) dXe−β H(X) 6 Then Equation 1.6 can be rewritten as Z hA(0)A(t)i = dX A(X)A(X(t))ρeq(X). (1.9) The most general assumption that serves as the foundation of statistical mechanics is ergodicity. This ergodic hypothesis states that the ensemble average is equivalent to the time average, e.g. for a time correlation function,12 Z τ 1 hA(0)A(t)i = lim ds A(s)A(t + s). (1.10) τ →∞ τ 0 Time correlation functions with the form hA(0)A(t)i are called autocorrelation function. Another form of time correlation function is called a cross-correlation function with a gen- eral expression hA(0)B(t)i (sometimes when discussing general properties of correlation functions, this form is preferred.). Physically, time correlation functions describe how a given property of a system affects the same property (or other property) at another time under statistical averaging, in other words, how much is left of the memory of a property at time 0 that still influences the same property (or other property) at time t. Because the correlation function is an ensemble average, this memory tends to be averaged out by interactions. A classic example is the velocity correlation function which determines the diffusion coefficient, Z ∞ 1 D= dthv(0) · v(t)i. (1.11) 3 0 1.2.3 Statistical mechanics and spectroscopy The Fourier transform converts the time correlation function C(t) = hA(0)A(t)i from time domain to frequency domain. The frequency-domain counterpart of the response function is generally called spectra density (SD) or power spectrum, Z ∞ C(ω ) = dt e−iω t C(t). (1.12) −∞ 7 The above relation is essential because it gives us a unified view of dynamics and its spec- trum with a time–frequency correspondence. For example, in infrared spectroscopy (IR), the absorption lineshape is the Fourier transform of dipole-dipole correlation function,15 Z ∞ 1 I(ω ) = dt e−iω t hM(0) · M(t)i, (1.13) 2π −∞ where M is the total dipole moment of the system. It would be wise to keep the time-domain/frequency-domain Fourier transformation in mind whenever looking at a time correlation function or a spectrum. Useful conversion between time and frequency is that 1 ps (=10−12 s) time scale corresponds to a frequency of 33 cm−1 (or 1 terahertz, THz) and 0.33 ps corresponds to 100 cm−1 . Most of many- body dynamics in liquid falls into the picosecond and sub-picosecond regime, which in frequency domain means a few wave numbers to hundreds of wave numbers. Now we go back to the topic of the instantaneous-normal-node theory. The theory enables one to calculate the so-called influence spectrum9, 16, 17 that is defined as follows,   ρA (ω ) = ∑ c2α δ (ω − ωα ) , (1.14) α where the sensitivity of variable A with respect to α -th INM is ∂A cα = . (1.15) ∂ qα The influence spectrum is the Fourier transform of a “velocity” correlation function GAA (t) = 2 ˙ A(t)i hA(0) ˙ = − dtd 2 hA(0)A(t)i, Z ˙ A(t)i hA(0) ˙ = kB T ρA (ω ) cos ω tdω . (1.16) The time correlation function hA(0)A(t)i can be derived by double time-integration of the 8 “velocity” correlation function GAA (t). So far, the solvation energy autocorrelation,16, 18 friction–force autocorrelation,17, 19 polarizability autocorrelation20 have been studied with the instantaneous-normal-mode approach. We will discuss the instantaneous-normal-mode theory and its application to nonlinear spectroscopy in detail in Chapter 4. 1.2.4 Langevin equation and generalized Langevin equation The fluctuation–dissipation theorem21, 22 has many a profound implication in statistical mechanics, one of which is the linear response theory (LRT). The connotations of the theorem can be understood with a Brownian motion model. The equation of motion for the random motion of a heavy particle like a dust particle immersed in a liquid bath, the Brownian motion, is called the Langevin equation. Generally speaking, Brownian motion refers to random motion of a particle or a mode in a macroscopic dynamical system caused by a large number of small particles or degrees of freedom. It is a basic model connecting macroscopic measurable random motion with thermal motion of microscopic degrees of freedom, which has wider applications besides heavy dust in liquid, for example, molecular simulations and nonlinear spectroscopy. The Langevin equation incorporates both frictional force and random noise. In one dimension, dv m = −ζ v + F(t), (1.17) dt where v is the velocity of the particle and ζ is the friction coefficient. The negative sign in the first term on the right side of the equation is determined by the fact that the friction drags the particle in the opposite direction of its velocity. In this equation, the frictional force is only dependent on the present velocity, which has no memory on the history of velocity. This amnesic process is often called Markovian. The last term F(t) is the random 9 noise (fluctuating force) satisfying the condition hF(t)i = 0, hF(t1 )F(t2 )i = 2Bδ (t1 − t2 ), (1.18) where constant B is the strength of the noise and δ (t) is the Dirac delta function. The condition indicates that the fluctuating force has no correlation between any distinct time instants, or is “white” — the Fourier transform of the correlation function of the noise is independent of frequency. Although the random term makes the Langevin equation a stochastic differential equation, the first-order linear inhomogeneous differential equation has a general formal solution Z t −ζ t/m 1 ′ v(t) = v(0)e + dt ′F(t ′)e−ζ (t−t )/m . (1.19) m 0 The velocity has an exponential decay of the initial velocity (first term) but the random noise produces extra velocity (second term). If we consider the mean squared value of velocity, B hv(t)2i = hv(0)2 ie−2ζ t/m + (1 − e−2ζ t/m ). (1.20) ζm In the long time limit, the first pure frictional term dies out, whereas the random noise keeps the system alive. Eventually, the mean squared value of velocity approaches the equipartition value kB T /m. Thus, B = ζ kB T. (1.21) The above result relates the fluctuating random noise B with the friction ζ (energy dissi- pation), and is a version of the fluctuation–dissipation theorem.21 This theorem is one of most fundamental cornerstones in nonequilibrium statistical mechanics, and its intrinsic connection to the linear response theory will be discussed in the next section. The mean squared displacement of Brownian particle at long times h∆x(t)2i → 2kBT t/ζ 10 gives the diffusion coefficient kB T D= , (1.22) ζ which is called the Stokes–Einstein formula.23 The Langevin equation can be only applied to a Markovian process, like the Brownian motion of all kinds, but far-reaching generalizations of the Langevin equation to non- Markovian cases give us more freedom describing complicated yet realistic systems. By adding a time-dependent friction kernel ζ (t), so-called memory function, we have the generalized Langevin equation (GLE),24, 25 d2 x(t) Z t ∂ W (x) m 2 =− − dt ′ ζ (t ′) · x(t ˙ − t ′ ) + F(t), (1.23) dt ∂x 0 where x(t) is the special coordinate or any degrees of freedom you are interested in, W (x) is potential of mean force on x coordinate (averaged over all the bath degrees of freedom), and F(t) is again the random noise. The convolution of the memory function ζ (t) and the velocity x(t) ˙ associates the evolution of x˙ with its earlier history. A GLE can be derived from a harmonic bath Hamiltonian by simply solving Hamilton’s equations of motion.26, 27 The corresponding non-Markovian version of the fluctuation–dissipation theorem reads ζ (t) = β hF(0)F(t)i. (1.24) The above equation explains the intrinsic equivalence between the friction in liquid and the fluctuations of random forces. The generalized Langevin equation is useful for studying relaxation processes. One can deduce the rate of vibrational energy relaxation from the microscopic friction by the Landau–Teller relation,28 1 ζR (ω ) = . (1.25) T1 µ The ζR (ω ) is the real part of the Fourier transform of the friction kernel, ζR (ω ) = 11 R∞ 0 dt cos ω t ζ (t), and T1 is the energy relaxation time for vibrational mode at frequency ω and reduced mass µ . The fluctuating force in the correlation function in particular is the solvent force on the vibrational mode of the solute. In a absorption spectroscopy, for instance, the homogeneous lineshape broadening results from the lifetime of a vibrational population relaxation T1 and the pure dephasing time T2∗ that is the time scale for losing phase information within a given vibrational state due to fluctuations of the environment. The pure dephasing rate can be calculated using the equation below Z ∞ 1 = dthδ ω (0)δ ω (t)i, (1.26) T2∗ 0 where δ ω (t) is the fluctuation of vibrational frequency.29 Then the following relation gives us the total dephasing time T2 that is directly related to homogeneous line width, 1 1 1 = + ∗. (1.27) T2 2T1 T2 The key term in generalized Langevin equation is the memory function.30 For some systems, there is no exact form for the memory function, so we need an approximated form. For example, the instantaneous-normal-mode theory provides a route to calculate the memory function as well as a microscopic perspective on what specific molecular mecha- nisms contribute to the vibrational friction. The solvent force correlation function that is equivalent to the friction kernel has a corresponding INM influence spectrum for vibrational relaxation.19 The influence spectrum is the solvent INM frequency distribution weighted by the sensitivity of the solvent force to each solvent mode. Theoretical calculations enables one to isolate typical liquid configurations and project out of the influence spectrum the contributions of a selected subset of solvent molecules or degrees of freedom. Simulation of simple liquids showed that the most effective modes during vibrational relaxation involve solute/nearest-single-solvent pairs, especially at high frequencies. This idea is called the instantaneous–pair (IP) theory.31 12 The choice of dynamical variable in the generalized Langevin equation is not limited to some peculiar particle positions and velocities; the generalized Langevin equation can be applied to any dynamical variable.32, 33 And there is more than one derivation of the generalized Langevin equation, such as that using the Zwanzig–Mori projection operator formalism.34–38 1.2.5 Linear response theory The most essential and fundamental connection between equilibrium and nonequilibrium statistical mechanics is the linear response theory. The fundamental idea of the linear re- sponse theory is first enunciated in 1931 by Onsager’s regression hypothesis: the relaxation of a macroscopic nonequilibrium disturbance will obey the same laws as the regression of spontaneous microscopic fluctuations in an equilibrium system.39, 40 This hypothesis is now viewed as an important consequence of the fluctuation–dissipation theorem proved by Callen and Welton in 1951.41 The fluctuation–dissipation theorem has many forms, two of which are shown in the previous two subsections, and so does the linear response theory. The significance of the linear response theory is that one will not have to actually perturb the system to get the response, instead all the information needed to predict the response is embedded in the fluctuation of the interested quantity in an equilibrated system. To put the linear response theory in a statistical mechanical language,42 if we observe a dynamical variable A in a system that is subject to an external perturbation F(t), with the perturbation Hamiltonian H ′ = −A · F(t) that is linear in the perturbation, then in the linear regime the nonequilibrium average value of A(t) is can be written as Z t A(t) = dt ′ R(t,t ′) · F(t ′), (1.28) −∞ where R(t,t ′) is the response function or the generalized susceptibility which describes how much the influence on the system at time t is left, due to an external perturbation at time t ′ . 13 The integration runs from time −∞ to time t, incorporating all the influence of the external perturbation till the observation time t. Generally, the response function has the following properties, (i) R(t,t ′) = R(t − t ′ ), (stationarity of the unperturbed system) (1.29) (ii) R(t − t ′ ) = 0, when t − t ′ < 0. (causality) (1.30) In linear response, the response function is independent of the external field as the response function is a property of the system. It will be shown in the next chapter that the response function has the following form, where θ (t) is the Heaviside step function. d R(t) = −β θ (t) hδ A(0)δ A(t)i, (1.31) dt where β = 1/kB T , and δ A(t) ≡ A(t) − hAi is the fluctuation of dynamical variable A at time t. For convenience, we choose the external perturbation F(t) = F · θ (−t). One can get that the nonequilibrium average by combining Equation 1.28 and Equation 1.31. ∆A(t) = A(t) − hAi = β Fhδ A(0)δ A(t)i, (1.32) another form of the fluctuation–dissipation theorem: the nonequilibrium quantity ∆A(t), the amount of “dissipation” in A that occurs as equilibrium is approached, has a natural connection to itself’s equilibrium “fluctuation” hδ A(0)δ A(t)i. We define the nonequilib- rium relaxation profile by normalizing ∆A(t), A(t) − hAi S(t) = . (1.33) A(0) − hAi This S(t) represents how fast the nonequilibrium quantity approaches its equilibrium value 14 starting from S(0) = 1 and ending up with S(∞) = 0. We can also define the normalized correlation function of the fluctuation in A as hδ A(0)δ A(t)i C(t) = . (1.34) hδ A(0)2 i This normalized correlation function describes how soon the system loses the memory of the fluctuation in A at an equilibrium condition. C(0) = 1 corresponds to no memory loss and C(∞) = 0 corresponds to a complete memory loss. The Onsager regression hypothesis is then expressed in the following equation, S(t) = C(t). (1.35) The above relation is also the criterion for testing if a system or a dynamical variable obeys the linear response theory. Time-dependent fluorescence spectroscopy (TDF), (which has the alternative name the time-resolved fluorescence Stokes shift (TRFSS)) has become a standard approach to investigate solvation dynamics (see Figure 1.2).2 In this measurement, the observable is some time-dependent characteristic fluorescence emission frequency, for example, the peak of fluorescence emission spectrum at time t. The dynamical variable is conventionally chosen as the solute-solvent interaction energy gap between ground and excited states, ∆E that is equivalent to the fluorescence frequency since ∆E = hν . The corresponding nonequilibrium solvation response S(t) and the energy-gap time correlation function C(t) are defined as ν (t) − ν (∞) hδ ∆E(t)δ ∆E(0)i Sν (t) = , and C∆E (t) = . (1.36) ν (0) − ν (∞) hδ ∆E(0)2 i Figure 1.3 shows remarkable agreement between experimental Stokes shift response Sν (t) and simulated energy-gap time correlation function C∆E (t) for coumarin 343 in water 15 Free energy excited state ground state Solvent coordinate Figure 1.2 Schematic illustration of time-dependent fluorescence spectroscopy (time- resolved fluorescence Stokes shift). Shown here is how an electronic transition in solute can be used to study solvation dynamics. In a polar solvent, a dipole is created in the solute with an ultrafast laser pulse, which brings the solute from the ground state to electronic excited state. The electronic transition is so fast compared with the nuclear motions of solvent that the initially prepared configuration of the excited state is still the same as in equilibrium ground state. As the solvation process happens the solvent rearranges to accommodate the change in the solute and to lower the solvation free energy. The relaxation can be monitored by detecting the emission from the solute as a function of time after excitation. A time- dependent red shift of the emission spectrum as a result of solvation is shown schematically. 16 solution.43 This agreement validates the linear response theory in the aqueous coumarin solution. The physical content behind is that, if the perturbation is not too large the relaxation of the nonequilibrium system perturbed by electronic transition is intrinsically the same with the relaxation of the spontaneous fluctuations in an unperturbed system at equilibrium. The solvation responses reveal a common bimodal pattern for solvation dynamics, initial ultrafast inertial (mostly librational motions) response followed by a slow diffusive component. Figure 1.3 Experimental and simulated solvation response functions for Coumarin 343 in water. The experimental nonequilibrium response (Sν (t): “expt.”) is the fluorescence Stokes shift function measured in 10−4 M Coumarin 343 sodium salt aqueous solution. The curve marked “∆q” is a classical molecular dynamics simulation result C∆E (t) using a charge distribution difference between ground and excited states calculated by semi- empirical quantum chemical methods, with SPC/E water potential. Also shown is a simulation for a neutral atomic solute with the Lennard-Jones parameters of the water oxygen atom (S0 ). [From R. Jimenez, G. R. Fleming, P. V. Kumar, M. Maroncelli, Nature 369, 471 (1994).] The traditional linear response theory uses the same approximation as in the fluctuation– 17 dissipation theorem which assumes the external disturbance is small. But in some cases the disturbance is not small at all. As long as the fluctuation obeys Gaussian statistics, the linear response is still valid.2, 44, 45 The linear response theory with Gaussian statistics will be derived in Chapter 2. The near universality of Gaussian statistics for microscopic quantities is guaranteed by the central limit theorem, because the “microscopic” quantity is usually a statistical average over the order of 1023 events. For most cases, the linear response theory can successfully predict the solvation relaxation.2, 46, 47 Figure 1.4 Rotational energy relaxation of CN in liquid Ar. The equilibrium linear- response prediction CE (t) (black curve) is compared with the nonequilibrium response functions SE (t) (colored curves) calculated for six different choices of the initial CN rotational energy EROT (0). The lower panel shows the first picosecond in greater detail, that the initial relaxation (about 100 fs) is given exactly by linear response and that the onset of deviations from linear response is marked by distinguishable wiggles. [From G. Tao, R. M. Stratt, J. Chem. Phys. 125, 114501 (2006).] However, the breakdown of linear response, in solvation processes for example, often accompanies a large initial disturbance of the the solute that pushes the system far from 18 equilibrium, and the subsequent relaxation involves a dramatic rearrangement of liquid structure.48–50 To make the breakdown happen, a sufficient separation between the solute time scale and that of the solvent geometry evolution is required.51, 52 Moskun et al. reported a linear-response failure in a rotational excited diatomic molecule solution.53 Experimentally, a highly excited CN rotor is generated by photodissociating ICN with an ultrafast deep ultraviolet laser pulse, and found that nearly free rotation lasts for tens of rota- tional periods (several picoseconds) in solution. Simulation shows that the nonequilibrium rotational relaxation is slower than the equilibrium time correlation function (Figure 1.4). The nonlinear response originates from the fact that the rotationally hot CN solute kicks a nearby solvent out of the innermost solvation shell creating a bubble in solution which makes the solvent structural relaxation time scale much slower than that of the solute.51 Among a few cases of linear-response breakdown, the nonequilibrium relaxation could be faster than the equilibrium linear-response prediction, with examples like solvated electron in methanol by Turi et al.54 and solute’s dipole flip in methanol by Ladanyi and co- workers.55, 56 The electron solvation calculations shows that a rapid decrease in the size of the solute alters the solvent structure, specifically different hydrogen bonding patterns between equilibrium initial and final states. The solvent collapses inwards to accommodate the solute’s transformation, which is easier than a swollen solute, thus faster. 1.2.6 Computer simulations Computer simulation is nothing but a well-controlled experiment on computer, which serves as a bridge connecting experiment and theory. One can test models by comparing simulation with real experimental results; one can test theories by comparing simulation with theoretical predictions.57 Computer simulation provides detailed microscopic infor- mation of molecular systems, as well as macroscopic measurable quantities of experimental interest. The flexibility to artificially choose models and selectively turn on or turn off certain contributions enables one to scrutinize molecular mechanisms and unearth hidden 19 details behind macroscopic measurements. For example, the diffusion coefficient can be determined by integrating velocity autocorrelation function as in Equation 1.11. The two main families of simulation technique are molecular dynamics (MD) and Monte Carlo (MC). Molecular dynamics simulation provides a direct avenue to real-time motions of in- dividual molecules in many-body systems.58 The trajectory, a sequence of positions and momenta of all the particles, is numerically evaluated by integrating classical equations of motion in a step-by-step fashion. The particles interact with each other via inter-/intra- molecular potentials, which is also called the force field. Once the initial condition (posi- tions and momenta of all particles in the system) is given, the trajectory is deterministically set. A widely used integration method for propagating the classical equations of motion is the velocity Verlet algorithm.59 Denote the time-dependent position, velocity and acceleration of a particle as r(t), v(t) and a(t) respectively. The time step is δ t, for instance several femtoseconds. The propagation of the position and velocity from time t to time t + δ t is expressed as follows. 1 a(t) = − ∇V (r(t)), (1.37) m 1 r(t + δ t) = r(t) + v(t)δ t + a(t)δ t 2, (1.38) 2 1 v(t + δ t) = v(t) + [a(t) + a(t + δ t)]δ t. (1.39) 2 Conventional molecular dynamics simulation is essentially a phase-space sampling technique from the microcanonical ensemble — constant NVE (number of particles, vol- ume, energy). The conservation of energy is automatically guaranteed by the nature of classical equations of motion. Although MD in other ensembles, like canonical (constant NVT) and isothermal-isobaric (constant NPT) can be achieved using artificial thermostat 20 and barostat constraints on system properties with correct equilibrium ensemble averages, some elements of the dynamics might no longer be correct. The most natural canonical-ensemble sampling method is Monte Carlo60 which is a random walk through the configuration space with probability consistent with thermal equilibrium ensemble distribution, in which the temperature is embedded. For example, in Metropolis MC algorithm,61 the probability accepting a random chosen trial move from R to R′ is the minimum value between 1 and the canonical Boltzmann factor ratio of the initial and final configuration. h ′ i P(R′ |R) = min 1, e−β [V (R )−V (R)] , (1.40) which assumes moving downwards on the energy surface is always acceptable, but climb- ing up the energy surface has to be determined by comparing a uniform distributed random number with the exponential Boltzmann probability. Monte Carlo simulation does not provide true time-dependent dynamical information, the sequence of configurations has no chronological connection, but it has a great number of applications to molecular systems such as statistical mechanics of rare events62 and molecular electronic structures.63 1.3 Experimental Methods The development of laser technology over the past decades has made time-domain obser- vation of molecular structure and interactions in condensed system a reality. A widely- used Ti:Sapphire laser oscillator generates pulse with only a few femtoseconds (1 fs = 10−15 s) duration, which is shorter than the time scale of most chemical reactions and molecular motions. Figure 1.5 shows time scales for different kinds of molecular motions and fundamental processes in physics, chemistry and biology.64 To date, the shortest laser temporal resolution reaches attosecond regime (1 as = 10−18 s) that has been used to investigate electron motions.65–67 21 Figure 1.5 Time scales relevant to physical, chemical and biological processes. The fundamental molecular vibrational motion is on femtosecond time scale. [From A. Zewail, J. Phys. Chem. A 104, 5660 (2000).] 22 Graham Fleming pioneered the field of ultrafast laser spectroscopy68 and Ahmed Ze- wail won the Nobel prize in chemistry in 1999 for building up the field of femtochemistry.64 1.3.1 Time-dependent fluorescence spectroscopy As mentioned in the previous section (Figure 1.2 and Figure 1.3), the time-dependent fluorescence spectroscopy has long been a main tool in understanding the dynamics of solvation,.2, 69–80 A newly excited dye solute with an abrupt change in charge distribution causes the solvent to reposition to lower the solute–solvent interaction energy, leading to a decreasing fluorescence frequency emitted by the solute. By tracking the time evolution of the fluorescence frequency shift (the Stokes shift), one gets the solvation relaxation profile Sν (t). The solvation relaxation reports on how fast the solution find its new optimized organization and the fluorescence frequency is equivalent to the solute–solvent interaction energy which is mostly affected by the interaction between the solute and the solvent in the first shell.81 This time-dependent-fluorescence idea has been applied to analyze solvation dynamics in simple liquids like chromphore solutions82–84 and in complex systems such as proteins,74 DNAs75–78 and reverse micelles.79, 80 However, this solute–solvent interaction energetic probe cannot directly tell the structure of the solvent — the dye solute serves as the sole reporter. One can only infer how the solvent structure changes based on the energetic information of the solute–solvent interaction. Figure 1.6 Molecular structure of coumarin 153. Levinger and co-workers studied the solvation response for coumarin 153 (C153) solute (Figure 1.6) in acetonitrile–benzene mixtures using time-dependent fluorescence 23 technique.85 Acetonitrile has a dipole around 4 Debye, and benzene is nonpolar but possesses quadrupole. The solute C153 has a large dipole creation (about 8 Debye) upon photoexcitation, thus the excited C153 favors the dipolar acetonitrile. What they found interesting is that the solvent mixtures relax more slowly than either pure solvent, with the 5% acetonitrile mixture the slowest as shown in Figure 1.7. This phenomenon is called preferential solvation,86 the solvation process with multiple solvents when the solvating abilities of each solvent are different. The slow-down of solvent mixture has its roots in the time-consuming exchanging process involving both solvents that is absent in a pure-solvent case.85–98 Starting from chapter 3, we will discuss the molecular mechanism of preferential solvation thoroughly and use it as a model system for applying solute-pump/solvent-probe spectroscopy. 1.3.2 Transient absorption The transient absorption measurement is implemented by a pump–probe configuration of ultrafast laser pulses.99 The pump–probe configuration refers to a general description of the set-up of many ultrafast spectroscopies. The time-zero-defining pump laser pulse initiates chemical reactions and the time-delayed probe laser pulse is used to take a snapshot of the molecular behavior (by measuring absorption for example). To study molecular motions, the relative timing of ultrafast laser pulses has to be controlled accurately. The time delay between synchronized pump and probe pulses is varied by diverting the probe pulse through an adjustable optical path length. The numerous choices for wavelengths of pump and probe pulses and specific probing technique lead to different kinds of pump–probe spectroscopies that are sensitive to distinct properties of matter. Transient IR absorption experiment excels at observing the evolution of molecular vibrations under the influence of ambient chemical environment in real time.100–106 This technique has been successfully applied to study water dynamics.107–109 For instance, Fayer and co-workers studied hydrogen bond breaking dynamics through tracking transient 24 S(t) S(t) Figure 1.7 Solvation response functions of coumarin 153 in acetonitrile–benzene binary mixtures measured by time-dependent fluorescence spectroscopy. The mole fraction of benzene is shown in the legend. Upper panel (a) depicts the longer-time-scale solvation relaxation; the lower panel (b) shows the initial solvation response in the first 2.5 ps. [From B. M. Luther, J. R. Kimmel, N. E. Levinger, J. Chem. Phys. 116, 3370 (2002).] 25 absorption of the hydroxyl OD stretching mode of HOD in water at room temperature.109 In their pump–probe setup, a pump pulse excites a selected vibrational mode(s) followed by a time-delayed probe pulse that detects the instantaneous IR absorption. Ultrafast X-ray absorption spectroscopy (XAS)110–118 is implemented in a laser-pump/X- ray-probe geometry. It measures the difference in X-ray absorption between the laser- excited sample and the unexcited sample. Bressler and Chergui111, 117 reviewed the time- resolved X-ray absorption spectroscopy and its ability to reveal ultrafast molecular geo- metrical structures as well as electronic structures in liquids with a typical time resolution of 50–100 ps. Recently, the temporal resolution has been achieved in the sub-picosecond regime119 which makes this technique promising in studying ultrafast chemical dynamics in condensed phases. 1.3.3 Scattering experiments Scattering experiments offer accurate detection of microscopic structures by impinging elementary particles/waves on matter and watching the scattered particles/waves at all directions. The scattered particle could be a photon, electron or neutron. Photons and electrons are scattered mostly by electrons in the material, whereas neutrons are scattered primarily by nuclei. A popular photon scattering technique is X-ray diffraction (XRD)120 which can probe structures at the angstrom scale in bulk samples, primarily solid materials, and is particularly sensitive to lattice dynamics, melting and phase transitions. In light of the million-fold larger scattering cross sections for electrons than X-rays, time-resolved electron scattering has been utilized to study transient structures in dilute gas phase and thin materials.115 Neutron scattering is able to provide complete microscopic information on time-dependent nuclear positions of all atoms, but the neutron source requires a expensive nuclear reactor.121, 122 X-ray diffraction provides three-dimensional atomic-scale static structure of lattices with impressive examples like unraveling highly complex structures of proteins and DNA. 26 Thus usually, X-ray diffraction from solid matter is also called X-ray crystallography. The periodic structures of atoms enhances the diffraction intensity at particular angles, the Bragg angles. Although the well-defined crystalline planes are absent in disordered samples, X-ray diffraction of liquids (diffuse scattering) gives structural information with the form of radial distribution functions.123, 124 In X-ray scattering, the key quantity is the momentum transfer q = k − k0 , where k0 and k are the incident and the diffracted X-ray wave vectors. The magnitude of momentum transfer obeys q = 4π sin θ /λ , where λ is the wavelength of the X-rays and 2θ is the scattering angle. The amplitude of the diffracted electric field is proportional to Z E(q) = ρe (r)e−iq·rdr = ∑ fn (q)e−iq·rn , (1.41) n in which ρe (r) is the electron density and fn (q) is the atomic form factor that is defined as R Fourier transform of the electron density of atom n, fn (q) = ρn (r)e−iq·rdr. So the scattering intensity 2 I(q) = ∑ En (q) = ∑ fn (q) fm (q)e−iq(rn−rm ) (1.42) n n,m and in isotropic medium, the orientational average of the scattering intensity leads to sin qrnm I(q) = ∑ fn2 (q) + ∑ fn (q) fm (q) , (1.43) n m |e> k2 k3 t2 |g> k3 -k2 t1 |g> 2. So Gaussian distribution is fully described by the covariance matrix that can be diagonalized by an orthogonal matrix U (UT = U−1 ) by M = UT Λ U, where the diagonal matrix Λ = UMUT = diag(λ1 , · · · , λn ). Define the new basis q = Ux, we then have xT Mx = (xT UT )(UMUT )(Ux) = qT Λ q = ∑ λα q2α . (2.61) α 63 We test the normalization condition for the distribution, Z ∞ r Z ∞ r Z ∞ det M −xT Mx det M ∂ (x1 , . . . , xn ) −qT Λ q dxρ (x) = dxe = dqe −∞ πn −∞ πn −∞ ∂ (q1 , . . . , qn ) r Z ∞ Z ∞ det M n 2 = ··· | det UT |dq1 · · · dqn e− ∑α λα qα . (2.62) πn −∞ −∞ Using orthogonal matrix property det UT = ±1, we have r r r det M n det M n Z ∞ Z ∞ −λα q2α π π n α∏ π n α∏ dxρ (x) = dqα e = = 1. (2.63) −∞ =1 −∞ =1 λα The average of the square of qα can be calculated as follows, r Z ∞ r  Z ∞  λα 2 λα ∂ −λα q2α hqα2 i = q2α e−λα qα dqα = · − e dqα π −∞ π ∂ λα −∞ r  r  λα ∂ π 1 = − = . π ∂ λα λα 2λ α Arranging the above equation, the relation between the eigenvalues and the average of corresponding basis square is 1 = 2hqα2 i. (2.64) λα 2.3.2 Generating function T The generating function (or the characteristic function) is defined as G ≡ heF x i, with co- efficient vector F = (a, b, c, . . .)T . The generating function for the n-dimensional Gaussian variable can be derived in the following two ways. 64 (1) Derivation of basic Gaussian identity by diagonalizing the covariance matrix. T G = heF x i (2.65) T UT )(Ux) T = he(F i = hef q i (define f = UF) r det M n Z ∞ 2 dqα e−λα qα + fα qα π n α∏ = =1 −∞ r   2 fα det M n −λα q2α − λfα qα + Z ∞ 2 fα α 2 π n α∏ 4λα = dqα e e 4λα =1 −∞ n 2 fα 1 2 2 = ∏ e 4λα = e 2 ∑α fα hqα i. α =1 Since ∑ fα hqα2 i fα = ∑ fα δαβ hqβ qγ iδγα fγ = f(UT U)hqqT i(UT U)f = FT hxxT iF, α αβ γ the generating function can be expressed as 1 T hxxT iF G = e2F . (2.66) (2) The other way to derive the generating function does not involve diagonalizing the covariance matrix. Z ∞ Z ∞ FT x −xT Mx+FT x T Mx G = he i= dxe dxe−x . (2.67) −∞ −∞ 65 To complete the square term, we need to obtain a vector c such that xT Mx − FT x = (x + c)T M(x + c) − cT Mc, −FT x = cT Mx + xT Mc = cT Mx + (Mc)T x, −FT = 2cT M. ⇒ F = −2Mc. 1 ⇒ c = − M−1 F. 2 1 1 cT Mc = FT M−1 MM−1 F = FT M−1 F. 4 4 So, " # T 1 T −1 1 G = heF x i = e 4F M F = exp ∑ Fµ (M−1 )µν Fν . 4 µν Then, the derivative of the generating function   ∂G G −1 −1 ∂ Fµ = 4 ∑(M )µν Fν + ∑ Fν (M )µν . ν ν ∂ 2 G 1 T So, hxx iµν = = M−1 . (2.68) ∂ Fµ ∂ Fν F=0 2 µν Thus, the generating function T 1 T −1 1 T hxxT iF G = heF x i = e 4 F M F = e2F . (2.69) 2.3.3 Linear response From the form of exact result Equation 2.14, we choose x = (x, y), and x = δ A(t)e, y = δ ∆V and the coefficient F = (a, β ). Assuming x and y are Gaussian variables, using 66 Equation 2.69, hxeβ y i ∂ lnheax+β y i ∂ ln G = = (2.70) heβ y i ∂a ∂ a a=0 a=0    ∂ 1  hx2 i hxyi   a  = (a β )    ∂a 2 hxyi hy2 i β a=0 1 ∂  2 2  = a hx i + 2aβ hxyi + β 2 hy2 i a=0 2 ∂a = β hxyi, (2.71) Substituting x and y by δ A(t)e and δ ∆V respectively, we get δ A(t) = β hδ A(t)eδ ∆V (0)i. (2.72) This is the identical to the traditional linear response result. But assuming Gaussian statistics does not require the perturbation δ ∆V to be small, and is thus more general. 2.4 Application to Solute-Pump/Solvent-Probe Spectroscopy So far, we have talked about four approaches to the linear response theory, two of which are particularly useful in the treatment of the solute-pump/solvent-probe spectra. The first one is the linear response with averaging on the excited state (Equation 2.17). This approach is useful since the solute excitation from the solute-pump is a resonant perturbation, although the linearization with respect to the energy gap is optional. The second one is the linear response with averaging on the ground state in the Heisenberg picture. This approach is useful because the solvent excitation from the four-wave-mixing pump is a nonresonant perturbation and it is natural to discuss the evolution of an observable rather than that of a density. 67 In two-dimensional time-domain solute-pump/solvent-probe spectroscopy, the nonequi- librium state is prepared by one resonant excitation and one nonresonant excitation. At time zero, the equilibrium ground-state system (g) undergoes an electronic transition to its excited state (e) by a resonant solute-pump, and the Hamiltonian changes from Hg to He = Hg + ∆V . The excited system evolves on the excited-state potential surface (e) for time T (this e surface is the “ground” state with respect to nonresonant electric field in the previous Heisenberg picture treatment) until a pair of solvent-pump pulses perturb the system with the interaction Hamiltonian 1 H ′ (τ ) = −B(X)F(τ ) = − 2∑ Eµ (τ )Πµν (X)Eν (τ ), (2.73) µν then the system starts to propagate on excited state e′ with Hamiltonian He′ = He +H ′ (τ ) for another short time t (this e′ surface is the “excited” state with respect to nonresonant electric field in the previous Heisenberg picture treatment). Finally at time T + t, the many-body polarizability Πγδ (T + t) is observed. This excitation process is depicted in Figure 2.3. For the sake of generality and simplicity in notation, we still use A(T +t) = Πγδ (T +t) as the observable, and B(X) = Πµν (X) as the coupled variable with the perturbing field F(τ ) = 12 Eµ (τ )Eν (τ ). The experiment measures the nonequilibrium average Z   δ A(T + t) = dX0 ρg (X0 ) A(X(T + t))|e,e′ − hAie (2.74) Z 1 = dX0 ρe (X0 )eβ δ ∆V (X0 ) δ A(X(T + t))|e,e′ . (2.75) heβ δ ∆V ie The notation δ A(X(T + t))|e,e′ is essential, which stands for the instantaneous value of fluctuation in dynamical variable A with preparation starting from initial condition X0 , evolving on e surface for time T and evolving one e′ surface for time t. The two subscripts e, e′ indicate which surface the system evolve on (or which Hamiltonian governs the propagation) during time periods T and t. This notation applies to the phase space point as 68 X0 X1 X(T+t) e’ e g 0 T T+t Figure 2.3 Sequence of excitations in solute-pump/solvent-probe spectroscopy. The molecular system is prepared in an equilibrated ground state (g), which gets electronically excited by the resonant solute-pump pulse at time 0, and promoted to the excited state (e). Propagating on the excited-state surface e for time T , the system is then perturbed by the nonresonant solvent-pump pulse that brings the system to excited state e′ . After evolving for a short time t on e′ surface, the many-body polarizability is observed at time T + t. The phase space points X0 , X1 and X(T + t) correspond to resonant solute-excitation time 0, nonresonant solvent-excitation time T , and observation time T + t. 69 well. As shown in Figure 2.3, starting from phase space point X0 , the system evolves for time T on e surface ending up with phase space point X1 = X(X0 ; T )|e . Then starting from phase space point X1 the system evolves on e′ surface for time t ending up with phase space point X(T + t) = X(X1 ;t)|e′ = X(X0 ; T,t)|e,e′ . Inside the parenthesis the first argument is the initial phase space point, and separated by a semicolon are propagation durations during which the system evolves on the surface marked by the subscripts in the order written. We define Liouville operators Le = −{He , }, (2.76)   ∂ Le′ = Le + + L1 , (2.77) ∂τ L1 = F(τ ){B(X), }. (2.78) The instantaneous value of dynamical variable A with the above mentioned preparation can be expressed as h i tLe′ T Le ′ A(X(t))|e,e′ = e e A(X0 ) = etLe A(X1 ) = A(X1 ;t)|e′ (2.79) ∂ = et (Le + ∂ τ )+L1 A(X1 )  Z t  t (Le + ∂∂τ ) ′ (t−t ′ )(Le + ∂∂τ ) t ′ (Le + ∂∂τ ) ≈ e + dt e L1 e A(X1 ) 0 Z t = A(t)e + dt ′F(t − t ′ ){B(X1;t − t ′ ), A(X1;t)}. (2.80) 0 In above derivation, the exponential of sum of operators, Equation 2.46 is used. The notation A(X1 ;t) means the value of dynamical variable A starting at phase space point 70 X1 and propagating on e surface for time t. The nonequilibrium average observed is Z  1 β δ ∆V (X0 ) δ A(T + t) = dX0 ρe (X0 )e δ A(T + t)e+ heβ δ ∆V ie Z t  ′ ′ ′ dt F(t − t ){B(X1;t − t ), δ A(X1;t)} (2.81) 0 heβ δ ∆V (0) δ A(T + t)eie Z 1 = β δ ∆V + β δ ∆V dX0 ρe (X0 )eβ δ ∆V (X0 ) he ie he ie Z t · dt ′F(T + t − t ′ ){B(X0; T + t − t ′ ), δ A(X0 ; T + t)} (2.82) 0 ≈β hδ ∆V (0)δ A(T + t)e ie Z t 1 + dt ′ F(T + t − t ′ )heβ δ ∆V (0) {B(T + t − t ′ ), δ A(T + t)}ie. heβ δ ∆V i e 0 (2.83) Now we substitute the physical quantities of solute-pump/solvent-probe spectroscopy into the above result, A(T +t) = Πγδ (X; T +t), B(X) = Πµν (X) and F(τ ) = 21 Eµ (τ )Eν (τ ). The sequence of pulses are shown in Figure 2.4. Z t 1 δ Πγδ (T + t) =β hδ ∆V (0)δ Πγδ (T + t)eie + dt ′Eµ (T + t − t ′ )Eν (T + t − t ′ ) 2 0 1 heβ δ ∆V (0) {Πµν (T + t − t ′ ), δ Πγδ (T + t)}ie. (2.84) heβ δ ∆V i e The first term above can be selected out by choosing particular phase matching condition in the signal detection. Thus, we have the solute-pump/solvent-probe response function that 71 time 0 T T+t u u u u Figure 2.4 Sequence of events in solute-pump/solvent-probe spectroscopy. The molecular system is initially in equilibrium with a ground-state solute (u with open black circle), which gets electronically excited by a resonant solute-pump pulse at time 0, and promoted to the excited state (u with filled black circle). This excitation changes the solute–solvent interactions and leads to structural reorganization in the surrounding solvent (colored circles). After waiting for time T , the experiment measures the ultrafast dynamics by hitting the sample with a pair of nonresonant laser pulses and scattering off the sample with a third nonresonant pulse at a shorter time t later. So the time T reflects the progress of the liquid structural change, and the time duration t corresponds to the ultrafast intermolecular vibrations of the liquid (black arrows) characteristic of the liquid structure at each time T . 72 is the most important result in this chapter, 1 R(0, T, T + t) = heβ δ ∆V (0) {Πµν (T ), δ Πγδ (T + t)}ie (2.85) heβ δ ∆V i e 1 = heβ δ ∆V (0) {δ Πµν (T ), δ Πγδ (T + t)}ie (2.86) heβ δ ∆V ie 1 = heβ δ ∆V (0) {Πµν (T ), Πγδ (T + t)}ie, (2.87) heβ δ ∆V ie since {B, A} = {δ B, A} = {δ B, δ A}. The anisotropic response function is practically calculated using the rotational average R(0, T, T + t) = Rµν µν (0, T, T + t) 1 = heβ δ ∆V (0) {Πµν (t1), Πµν (t2)}rotational ie , (2.88) heβ δ ∆V i e but for convenience, the rotational average will not be shown in next chapters. Until now, we have derived only the solute-pump-on response function (Equation 2.87). However, the solute-pump/solvent-probe experiment looks at the pump-on/pump- off difference, 1 ∆Rµνγδ (0, T, T + t) = heβ δ ∆V (0) {Πµν (T ), Πγδ (T + t)}ie − h{Πµν (0), Πγδ (t)}ig, heβ δ ∆V ie (2.89) which is the difference between the solute-pump/solvent-probe response and the four- wave-mixing response equilibrated with a ground-state solute. With optical-Kerr-like xzxz polarized laser configuration, the RP-PORS response function is ∆R(0, T, T + t) = ∆Rxzxz (0, T, T + t) 1 = heβ δ ∆V (0) {Πxz (T ), Πxz(T + t)}ie − h{Πxz (0), Πxz(t)}ig. (2.90) heβ δ ∆V i e 73 To understand this finding, consider the situation when the waiting time T is longer than the typical structural relaxation time in the liquid. Then the RP-PORS response function approaches 1 ∆R(0, T, T + t)|T →∞ = heβ δ ∆V (0) ie h{Πxz (T ), Πxz(T + t)}ie − h{Πxz (0), Πxz(t)}ig heβ δ ∆V i e =h{Πxz (T ), Πxz (T + t)}ie − h{Πxz (0), Πxz(t)}ig (2.91) d =−β [hΠxz(0)Πxz (t)ie − hΠxz (0)Πxz(t)ig] , (2.92) dt which is simply the difference between four-wave-mixing response functions on two distinct electronic states. At long waiting time limit, the pump-on response is the ordinary OKE response in equilibrium with an excited-state solute. We will concentrate on this large T limit in chapter 3. By contrast, at the limit of T = 0, the pump-on response reflects the light scattering that samples the initial configuration on the equilibrium ground state, with the subsequent dynamics propagating on the excited state. § Z R(0, 0,t) = ˙ dXo ρg (X0 ){Π(0), Π(t)e} = h{Π(0), Π(t)e}ig = β hΠ(0)Π(t) e ig . (2.93) The corresponding solute-pump/solvent-probe response function is ˙ ∆R(0, 0,t) = β [hΠ(0)Π(t) ˙ eig − hΠ(0)Π(t)i g], (2.94) which highlights the pump-on response is the correlation between the ground-state polar- ˙ = dΠ/dt and the excited-state polarizability Π(t)e. In other words, the izability velocity Π solute excitation does not change the instantaneous value of polarizability velocity; what § Throughout this thesis, Π stands for any element of the many-body polarizability (usually the off- diagonal elements when discussing the anisotropic four-wave-mixing response), and Π denotes the second- rank many-body polarizability tensor. 74 gets measured by the four-wave-mixing probe is the instantaneous polarizability velocity correlated with the time-delayed polarizability evolving on the final solute electronic state. When the waiting time T is in an immediate range, the solute-pump/solvent-probe response function (Equation 2.87) cannot be simplified into some ordinary two-time or three-time correlation function as in ordinary four-wave-mixing response (Equation 2.52). The Poisson bracket in the solute-pump/solvent-probe response is part of the reason, since for functions with form hA(0){B(t1),C(t2)}i, there is no relation like h{A(0), B(t)}i = ˙ β hA(0)B(t)i that can transform the solute-pump/solvent-probe response function to some ordinary correlation functions. A similar situation happens in the fifth order Raman spectra196–217 where the initial excitation is a nonresonant light scattering event instead of a resonant solute-pump. The Poisson brackets exist in these 5th order experiments due to the fact that they are measuring the sensitivity of a four-wave-mixing light scattering to a initial perturbation, rather than measuring the light scattering itself. A Poisson bracket describes the sensitivity of a system’s trajectories because it reports on how the value of some dynamical variable at time t is affected by variations in another dynamical variable at time 0. In principle, one can resort to nested molecular dynamics to evaluate a Poisson bracket. For example, in our case, the Poisson brackets {Π(T ), Π(T + t)} can be expressed as follows, since the many-body polarizability depends only on the liquid configuration R. 3N      ∂Π ∂Π ∂ Rk (T + t) {Π(T ), Π(T + t)} = {Π[R(T )], Π[R(T + t)]} = ∑ ∂ Ri ∂ Rk ∂ Pi (T ) , i,k=1 T T +t (2.95) where the derivative ∂ R(T + t)/∂ P(T ) can be numerically evaluated by performing a series of perturbed trajectories (T → T + t) and averaging the changes in positions at time T + t with respect to some initial momenta P(T ) that has a finite difference from that in the unperturbed trajectory. This numerical approach has been applied to calculate nonlinear spectroscopic responses in some well-behaved systems,196, 197, 234, 235 but the 75 Poisson bracket is essentially fluctuating and practically difficult to compute, because of the intrinsically chaotic nature of classical many-body systems.236–238 The correlation function of a Poisson bracket with another fluctuation would be even more noisy. An approach that one can take to avoid the issue of chaotic Poisson brackets is modeling the real experiment through carrying out (finite field) nonequilibrium simulations with explicit external perturbations.206, 210, 239 To increase the computational efficiency in evaluating 5th order Raman spectra, for instance, another approach has been proposed which combines equilibrium MD with nonequilibrium finite field trajectories.211, 216, 217 However, these approaches utilizing nonequilibrium trajectories abnegate the real power of linear response theory that enables one to predict a system’s response based on equilibrium averages without actually perturbing the system. Besides computational advances in evaluating the Poisson bracket, analytical attempts using generalized-Langevin-equation language240, 241 and the mode-coupling theory207, 212, 242–245 have been made. Moreover, quantum mechanics motivates valuable approximations into classical nonlinear spectro- scopic responses.201–205 We will get back to the topic of evaluating Poisson brackets in chapter 4 where we will show a practical approach to calculate Poisson bracket taking advantage of the instantaneous-normal-mode approximation. In this chapter, both traditional linear response theory and the linear response theory with Gaussian statistics are derived in a general way using classical statistical mechanics. The connection between the linear response theory and nonlinear spectroscopic responses has been made, and the solute-pump/solvent-probe response function is derived using the linear-response-theory formalism. Brief discussion about the difficulty in evaluating Poisson bracket is given. In the next chapter, we will focus on the long waiting time T case of the solute-pump/solvent-probe spectroscopy. 76 Chapter 3 Preferential Solvation Dynamics and Optical Kerr Effect Spectroscopy 3.1 Preferential Solvation Dynamics The concept of solvation dynamics is introduced in chapter 1, a structural rearrangement process of solvent molecules to accommodate a change in the solute such as a charge redistribution triggered by a laser-induced electronic excitation. The driving force of the structural reorganization originates from the tendency of the solute–solvent composite system to lower its free energy. The questions relevant to the solvation process are: How fast is the solvation relaxation? What molecular features determine the relaxation rate? More interestingly, what is the molecular mechanism during the solvation pro- cess? For example, how would the presence of solvent alter the reaction pathway in chemical dynamics. How can one characterize the solvation dynamics theoretically and experimentally? One of the most representative approach to solvation dynamics is time- dependent fluorescence study on chromophore solutions.2, 83, 141 By watching characteristic emission fluorescence frequency of the excited chromophore solute, one can construct a solute–solvent interaction energy relaxation profile, which serves as the measure for the 77 progression of solvation.246, 247 In dipolar solvents, the solvation relaxation has been found to be contributed mostly by reorientational motions of the solvent.81 The solvation profiles exhibit bimodal behavior — subpicosecond inertial motion of solvent and picosecond diffusive motion.83 The ultrafast solvation time is mostly determined by the dynamics of the first solvation shell.81, 248 The reason why only a local portion of molecules governs the solvation process is that the observable, potential energy of a single solute is only sensitive to the ambient solvent near the solute. Following an electronic excitation of the solute, such as creating or increasing the solute’s dipole, a net number of solvent molecules translate into the first solvation shell leading to a more crowded arrangement (“electrostriction”). However, in solvent mixtures, another translational relaxation pathway emerges as a result of the new choices of replacing one kind of the solvents with another in the first shell (“redistribution”). This redistribution process is the root of the preferential solvation.85–98 In preferential solvation, multiple solvents with different solvating abilities are mixed, for example, dipolar acetonitrile and nonpolar benzene85 (Figure 1.7). An interesting phenomenon in preferential solvation is that the solvent mixture displays a much slower solvation than either pure solvent case, which is now believed to be a consequence of the slow redistribution process.87 Moreover, the solvation rate depends on the solvent composition with observations such as a small percentage of the favored solvent exhibits the slowest solvation.85, 87, 220 Although the redistribution of different species of solvents is always present in liquid mixtures, a significant slowdown of the solvation is not guaranteed. In order to see a conspicuous preferential solvation, the dynamical properties of the solvents need to be different enough such that one solvent frustrates another solvent’s relaxation thus impeding the overall solvation.90, 92, 95 Solvent mixtures are often used in organic chemistry when a pure solvent cannot dissolve all the reactants. Table 3.1 shows some common solvents’ average polarizabilities and dipoles. To date, solvation properties of many binary solvent mixtures have been inves- tigated microscopically including dimethylsulfoxide (DMSO)/water,89, 91, 92, 249–252 ben- 78 zene/acetonitrile (ACN),85, 95 water/alcohol,89, 253 alkane/alcohol,94, 254 cyclohexane/ben- zene,161 tetrahydrofuran (THF)/water,255 THF/ACN, ethanol,256 N, N-dimethylformamide (DMF)/water, alcohol,257 DMSO/benzene,258, 259 and CO2 /ACN, methanol, cyclohexane.260, 261 Solvent ˚ 3) Polarizability (A Dipole (D) Water 1.47 a 1.88 c Dimethylsulfoxide (DMSO) 7.97 b 3.96 d N, N-Dimethylformamide (DMF) 7.81 d 3.82 d Ethanol 5.11 c 1.44 c Methanol 3.31 c 1.70 c Acetonitrile (ACN) 4.51 c 3.91 c Tetrahydrofuran (THF) 8.0 e 1.75 d p-Dioxane 8.6 d 0 Acetone 6.33 d 2.88 d Dichloromethane (DCM) 6.48 d 1.60 d Chloroform 8.5 c 1.04 c Cyclohexane 11.0 d 0 Hexane 11.9 d 0 Benzene 10.6 c 0 Table 3.1 Molecular polarizability and dipole of common solvents. (Conversion of polarizability unit, 1 a.u. = 0.148185 A˚ 3 .) (a) W. F. Murphy, J. Chem. Phys. 67, 5877 (1977); (b) K. J. Miller, J. Am. Chem. Soc. 112, 8533 (1990); (c) C. G. Gray and K. E. Gubbins, Theory of Molecular Fluids, Vol.1: Fundamentals (Oxford, New York, 1984); (d) W. M. Haynes, Ed., CRC Handbook of Chemistry and Physics (CRC, Boca Raton, 2013), 93rd ed.; (e) H. Reis, A. Grzybowski, and M. G. Papadopoulos, J. Phys. Chem. A 109, 10106 (2005). In this thesis, the idea of solute-pump/solvent-probe spectroscopy will be applied to a model problem showing preferential solvation. In the rest of this chapter, I will begin with introducing an atomic liquid model for molecular dynamics simulation and show the solvation responses with respect to different solvent compositions. Then the following section describes the solute-pump/solvent-probe spectra at the infinite waiting time limit for this preferential solvation problem (which gives the difference in optical Kerr effect spectra between equilibrated excited-state system and ground-state system). This result will enable us to distinguish how the ultrafast intermolecular dynamics varies with solute’s electronic states. Molecular-level analysis of the spectra reveals local solvent dynamics. 79 3.2 Molecular Dynamics Simulation The major reason for the slow preferential solvation is the solvent redistribution, which is mostly translational motions.85–88, 93–95, 254 It is not only straightforward but also wise to choose an atomic liquid model to simulate a preferential solvation process over some realistic molecular liquid models. It turns out that an atomic solute with mixture of two different atomic solvents is sufficient to reproduce the whole dynamical phenomenology of preferential solvation that we are interested in. The atomic liquid model we are going to use was first proposed by Sakurai and Yoshimori.97 It consists of a single atomic solute (u) dissolved in a mixture of strongly solvating (S) and weakly solvating (W) solvents. The solute could be in its ground electronic state (u-g) or its excited electronic state (u-e). All atoms have identical mass (m) and diameter (σ ). The interaction between every pair of atoms are described by the Lennard–Jones (LJ) potentials,   σ 12  σ 6 u(rab ) = 4εab − , (a, b = u-g, u-e, S,W). (3.1) rab rab εab = ε , (all interactions except those involving u-e) εu-e,S = 3ε , εu-e,W = 1.5ε . The dynamics of the system with a ground-state solute is no different from that of a pure Lennard–Jones liquid, whereas when the solute is promoted to the excited state, the solute–solvent attraction increases and the S solvent has a deeper energy well depth than the W solvent. The physical timescale is defined by choosing argon LJ parameters m = 39.948 amu, σ = 3.405 A,˚ ε /kB = 119.8 K, and the corresponding time unit τLJ = p mσ 2 /ε = 2.16 ps. Molecular dynamics simulations for the system are performed with thermal conditions density ρσ 3 = 0.8 and temperature kB T /ε = 1.00 ± 0.03. In simulations with total atom number N = 256, there are one solute and 255 solvent 80 4.0 ground 3.5 excited 10% S 3.0 excited 50% S 2.5 g(r) 2.0 1.5 1.0 0.5 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 r / Figure 3.1 Solute–solvent radial distribution function in atomic liquid with ground-state solute, 10% S system with excited-state solute, and 50% S system with excited-state solute. Total number of atoms is 256 including one solute. Radial distribution functions are averaged over 106 liquid configurations. 81 atoms. In a 10% S system, there are 26 S solvent atoms and 229 W solvent atoms. Molecu- lar dynamics is simulated with cubic periodic boundary conditions using the velocity Verlet propagation algorithm.59 The time step δ t is 0.0025 τLJ = 5.4 fs, and ensemble averages are computed by sampling the liquid configurations every 10 time steps. Equilibration is achieved by performing standard procedure whose details are shown in Appendix D. Figure 3.1 shows the solute–solvent distribution function g(r) for the ground-state system and excited-state systems with 10% S and 50% S. The radial distribution function is defined as g(r) = n(r)/(ρ · 4π r2 dr), where n(r) is the number of solvent atoms at a distance between r and r + dr from the solute, and ρ is the average number density. Shown in Figure 3.1 are typical equilibrium liquid-state radial distributions. The location of first minimum of the solute–solvent radial distribution defines the boundary of the first ˚ Similarly the second shell cutoff radius solvation shell, r1st min = 1.556 σ = 5.298 A. ˚ It is worth noting that the minimum locations do not vary with r2nd min = 2.548σ = 8.676A. the solute’s electronic state nor the solvent composition. Another check for the equilibrium in ground state is the velocity autocorrelation function (Figure 3.2), whose definition is as follows. hv(0) · v(t)i Cvv (t) = . (3.2) hv2 i The negative velocity correlation in the region 0.3 < t < 0.7 indicates that liquid atoms have more probability moving in the opposite directions to those they used to move along at a time t earlier. After 1.25 ps, atoms almost lose all the memory of their velocities at time 0 corresponding to Cvv → 0. The solute–solvent interaction energy is the excited-state/ground-state energy gap that can be expressed as follows h i ∆V = ∑ uu-e,v (ruv ) − uu-g,v (ruv ) , ruv = |ru − rv |. (3.3) v∈S,W In our case, the normalized nonequilibrium solvation relaxation profile S(t) and the nor- 82 1.0 k T/ = 1.00 0.8 B m) = 0.800 0.6 (t) ( 0.4 vv C 0.2 0.0 0.0 0.5 1.0 1.5 t (ps) Figure 3.2 Normalized velocity autocorrelation function for Lennard–Jones liquid at density ρσ 3 = 0.8 and temperature kB T /ε = 1.00. Total number of atoms N = 256. Velocity correlation function is averaged over 2 × 106 liquid configurations. malized equilibrium correlation function C(t) defined in chapter 1 can be written as ∆V (t) − h∆V i hδ ∆V (0)δ ∆V (t)i S(t) = , C(t) = . (3.4) ∆V (0) − h∆V i hδ ∆V (0)2i Figure 3.3 shows the simulated equilibrium C(t) and the nonequilibrium S(t) for this model in a range of solvent compositions including pure solvent cases. An evident solvent composition dependence is observed and the slowest relaxation takes place in the case most dilute in S, 10% S system, which is consistent with the experimental observations in Figure 1.7. Thus this atomic liquid mixture model imitates the phenomenology of preferential solvation successfully. The excellent agreement between C(t) and S(t) also validates the linear response theory in this model system. Moreover, the pure W solvent (0% S) and the pure S solvent (100% S) exhibit identical solvation responses for time longer than a picosecond. 83 1.0 1.0 Solvation response 0.8 0.5 0.6 0.0 0.0 0.5 1.0 0% S 0.4 Solid: C(t) 10% S Dashed: S(t) 50% S 80% S 0.2 100% S 0.0 0 50 100 150 t (ps) Figure 3.3 Solvation dynamics for the atomic liquid preferential solvation model in a range of solvent compositions. Shown here are both the equilibrium solvation correlation function C(t) and the nonequilibrium solvation relaxation function S(t), with the detailed view of subpicosecond relaxations in the inset. The C(t) and S(t) profiles are averaged over 5 × 106 and 105 liquid configurations respectively. 84 Figure 3.4 Time evolution of the first shell population in 10% S solvent mixture after solute excitation. The same atomic liquid model is used as described in the text. Total number of solvent in the first shell is the black curve, the number of S solvent in the first shell is the red curve and the number of the W solvent in the first shell is the blue curve. [From C. N. Nguyen, R. M. Stratt, J. Chem. Phys. 133, 124503 (2010).] The two stages of preferential solvation, electrostriction and solvent redistribution, can be identified in Figure 3.4 where the time evolution of the first shell population in 10% S mixture is plotted against the time after switching the solute’s electronic state. The 1/e time for the relaxation of the total number of solvents in the first solvation shell is less than 1 ps, corresponding to the solvent compression (“electrostriction”), while the 1/e time for S solvent population relaxation is 12 ps corresponding to the solvent exchanging process (“redistribution”). A recent geodesic (most efficient) pathway262–264 investigation to the same atomic liquid model performed by Nguyen and Stratt demonstrates that the slow solvent exchanging process exhibits rather a sequential mechanism than a concerted mechanism.87 The number of S and W solvent population and their 1/e time constants are summarized in Table 3.2. Two systems are particularly interesting: the 10% S solvent mixture has the largest S solvent percentage change (1.3 S → 3.0 S) and the 50% S solvent mixture has the largest S solvent number change (+3.2 S). Hence these two systems are 85 going to be investigated using the solute-pump/solvent-probe spectroscopy, and the long waiting time T limit spectra will be discussed in rest of this chapter. 1st-shell eq. population time scale g e e-g S W Tot 10% S 1.3 S + 11 W 3.0 S + 9.7 W +1.7 S - 1.3 W 12 15 0.7 50% S 6.2 S + 6.1 W 9.4 S + 3.5 W +3.2 S - 2.6 W 10 12.8 0.8 80% S 10 S + 2.4 W 12 S + 1.1 W +2.0 S - 1.3 W 5.5 8.7 0.8 Table 3.2 Averaged first shell population and 1/e relaxation time (ps). S and W solvent populations on equilibrium ground state (g), equilibrium excited state (e) along with the difference between the two (e-g) are shown on the left side. The 1/e times for population relaxation of S, W and total solvent in the first shell are shown on the right side. [From C. N. Nguyen, Ph.D. thesis, Brown University (2011), chapter 2.] 3.3 Optical Kerr Effect (OKE) Spectroscopy In chapter 2, the solute-pump/solvent-probe response function at the long waiting time T limit is proven to be the difference between the excited-state optical Kerr effect (OKE) response and the ground-state OKE response (Equation 2.92). As mentioned in chapter 1, OKE spectroscopy20, 142, 145–159 is a four-wave-mixing technique with a pump–probe laser configuration. Essentially, OKE spectroscopy is a kind of third-order time-domain nonresonant coherent Fourier transform Raman spectroscopy, and is an excellent tool for the interrogation of ultrafast dynamics in liquids. In this section, a more detailed description of the OKE experiment is given followed by our theoretical route. 3.3.1 OKE experiment The experimental implementation of optical Kerr measurement is straightforward by using a pump–probe geometry of lasers. Figure 3.5 shows the pump–probe setup of the ex- periment. The simplest way to measure an OKE signal is called homodyne detection.142 The pump pulse polarized at 45◦ with respect to the vertical direction perturbs the sample and induces a transient birefringence. Then the probe pulse polarized vertically (0◦ ) 86 arrives at the sample at a delayed time to read out the effect of the perturbation. After the sample, the probe pulse hits an analyzing polarizer (analyzer) before it gets detected. The analyzer is set to allow horizontally (90◦ ) polarized light to pass, so if there were no nonlinear optical effect (for instance the index of refraction were still isotropic), the probe pulse would not be able to pass the analyzer and thus no signal could be detected. In fact, the transient birefringence created by the pump pulse provokes depolarization of the probe pulse, allowing a signal to leak through the analyzer. By scanning the delay time between the pump pulse and the probe pulse, time-dependent optical Kerr effect can be measured. This polarization spectroscopy configuration is sensitive to the depolarized (3) response Rxzxz (t) that is proportional to the difference between two elements of the third- (3) (3) (3) order response, Rzzzz (t) − Rxxzz (t).265 In an isotropic medium, this Rxzxz (t) is known as the anisotropic response, 1 h (3) (3) i (3) Raniso (t) = Rzzzz (t) − Rxxzz(t) = Rxzxz (t). (3.5) 2 In the homodyne detection, the intensity of the signal is measured by the square-law (3) 2 photodiode detector, IOKE (t) ∝ Rxzxz (t) . By setting the angle between the pump pulse polarizer and the probe pulse polarizer to the magic angle (m = 54.7◦ which is the root of the second order Legendre polynomial P2 (cos θ ) = 0), one can measure the isotropic response152 1 h (3) (3) i (3) Riso (t) = Rzzzz (t) + 2Rxxzz(t) = Rzzmm (t). (3.6) 3 Partitioning the many-body polarizability into the isotropic (scalar) and the anisotropic 87 (traceless) parts Π = Πiso + Πaniso , the isotropic and anisotropic responses152, 266 β d aniso Raniso (t) = − Π (t) Π aniso (0) , (3.7) 15 dt d iso Riso (t) = −β Π (t) Π iso (0) . (3.8) dt pump detector probe λ/4 sample delay delay Figure 3.5 Schematic illustration of optical Kerr effect (OKE) spectroscopy. The upper panel shows the pump–probe polarization spectroscopy implementation of OKE experiment. The pump pulse has a polarization at 45◦ and the delayed probe pulse is polarized vertically. Before the detector, there is an analyzer polarizer that allows horizontally polarized light. In optical heterodyne detection (OHD), a quarter-wave plate (λ /4) is inserted after the polarizer in the probe beam. The bottom panel depicts a transient molecular structural change induced by the pump pulse — transient birefringence that is a result of a net alignment of the molecules along some specific direction defined by the pump pluse. In isotropic media, there are 21 nonzero elements out of total 34 = 81 elements of third- (3) order response tensor Rµνγδ , but only 3 elements are independent Rzzzz = Rxzzx + Rxxzz + Rxzxz .130 As in third-order Raman measurement, Rxzxz = Rxzzx ,266 only two independent response functions are enough to completely describe the third-order nonlinear optical behavior. Conventionally one can choose independent pairs of responses to characterize an isotropic fluid, anisotropic/isotropic (Raniso & Riso ) or polarized/depolarized (Rzzzz & Rxzxz ). 88 The traditional OKE spectroscopy refers to the anisotropic or the depolarized spectra. The isotropic spectra can be measured through the spatial mask technique developed by Tokmakoff and co-workers.152 However, the signal-to-noise ratio in isotropic OKE response measurement is not as good as that in anisotropic studies. Recent development of OKE spectroscopy can be found in Hunt, Jaye and Meech.154 Applying the rotational average, the response functions of third-order four-wave-mixing spectroscopy are given by267 (where pair product PP and trace Tr are defined in Equations 2.55, 2.56) (3) dD 1  1  E Raniso (t) = Rxzxz (t) = −β PP Π(t), Π(0) − Tr Π(t) · Tr Π(0) , (3.9) dt 10 30 (3) d D1  E Riso (t) = Rzzmm (t) = −β Tr Π(t) · Tr Π(0) , (3.10) dt 9 (3) dD 2  1  E Rzzzz (t) = −β PP Π(t), Π(0) + Tr Π(t) · Tr Π(0) . (3.11) dt 15 15 Optical heterodyne detection (OHD) Heterodyne detection overcomes the disadvantages in the homodyne detection including a generally weak OKE signal that is quadratic in the response function. In optical heterodyne detection (OHD),145, 147, 149, 150, 154 a large-amplitude phase-controlled electric field called a local oscillator (LO) is mixed with the signal field. So the total measured intensity has three terms, I(t) ∝ |ELO + EOKE (t)|2 = ILO + 2Re [E∗LO EOKE (t)] + IOKE (t). (3.12) The first term ILO is the largest, but does not depends on delay time t, thus this constant term can be filtered out by a lock-in detector. The second term is the heterodyne signal (3) which is linear in the OKE signal and therefore linear in the response function Rxzxz . The last term is the weak homodyne signal that can be removed by taking the difference of two 89 heterodyne measurements with opposite-signed local oscillators. The implementation of heterodyne detection is accomplished simply142, 154 by inserting a quarter-wave plate after the probe-pulse polarizer (Figure 3.5). The quarter-wave plate is aligned with one of its fast or slow axes along the polarization of the probe pulse. Then if one rotates the probe-pulse polarizer by a small angle like 1◦ , a 90◦ out-of-phase field with respect to the probe is generated and can pass through the analyzing polarizer, producing the local oscillator. What is measured by the detector is the intensity of the total electric field consisting the signal field and the local oscillator. If one rotates the probe-pulse polarizer by the same amount but in the opposite direction, the cross term (the heterodyne signal) changes its sign but the squared OKE intensity does not. So by taking the difference of the heterodyne detected signals with probe-pulse polarizer oriented at ±1◦ , the OHD- OKE signal SOKE (t) can be obtained. Since the laser pulses are not infinitely sharp, appreciable dynamics do occur on the time scale of the pulse duration. McMorrow and Lotshaw145, 268 developed the Fourier- transform deconvolution technique that removes the effects of finite pulse duration from the OKE signal. Given that the local oscillator is out-of-phase with respect to the probe ELO = iε E pr (ε is small), the OHD-OKE signal SOKE (t) can be shown to be the convolution of the OKE response function with the second order instrument correlation function269 Z ∞ (3) (3) SOKE (t) ∝ dt ′ G(2) (t ′)Rxzxz (t − t ′ ) = G(2) (t) ∗ Rxzxz(t), (3.13) −∞ where ∗ demotes a convolution, and the instrument response G(2) (t) is identical to the second-harmonic-generation intensity cross-correlation (convolution) between pump and probe pulses, G(2) (t) = I pump (t) ∗ I probe (t). Using the property of convolutions F [G(2) (t) ∗ (3) (3) Rxzxz (t)] = F [G(2) (t)] × F [Rxzxz (t)], the Fourier transform of the third-order response D(ω ) is simply (3) F [SOKE (t)] D(ω ) = F [Rxzxz (t)] = . (3.14) F [G(2) (t)] 90 In fact, the total response function contains contributions from the electronic and nuclear polarization, where the electronic response is a delta function peaked at time zero, bδ (t), so its Fourier transform is a real constant. Therefore, the imaginary portion of D(ω ) which is often called the spectral density (SD) is determined solely by nuclear response (dynamics). The time-domain nuclear response can be retrieved by performing inverse Fourier transform, Rnucl (t) = θ (t)F −1Im[D(ω )]. Z ∞ (3) SD(ω ) = Im[D(ω )] = dt Rxzxz (t) sin ω t. (3.15) 0 In most cases, we are interested in the ultrafast intermolecular dynamics, so we could re- move the long-time orientational diffusion contribution from the spectral density Im[D(ω )] to obtain the reduced spectral density (RSD), Im[D′ (ω )].142, 145, 154 The long-time diffusive contribution of the time-domain response (for example t > 2–3 ps) can be fitted by a sum of exponentials, R f it (t) = A1 e−t/τ1 + A2 e−t/τ2 , (A1 , A2 > 0). (3.16) After subtracted the diffusive portion, the response function is then Fourier transformed to calculate the RSD.   (3) Rdi f f (t) ≈ 1 − e−t/τrise R f it (t), (3.17) Z ∞ (3) (3) RSD(ω ) = Im[D′ (ω )] = F [Rxzxz (t) − Rdi f f (t)] = dt R′ (t) sin ω t. (3.18) 0 In the diffusive portion of response function, the rise time τrise can be evaluated by R R τrise = (2hν i)−1 with hν i = hω /2π i = dω 2ωπ SD(ω )/ dω SD(ω ), giving a typical value of several hundred fs. Its precise value does not have a substantial effect on the shape of the RSD.153 Experimental OKE spectra of 2,4,6-trifluoropyridine at 306 K are shown in Figure 3.6. The spectral density and the reduced spectral density are the upper and lower panels 91 (a) SD (b) RSD Figure 3.6 Optical Kerr effect spectra for 2,4,6-trifluoropyridine at 306 K. Upper panel (a) shows the spectral density (SD) and the lower panel (b) shows the reduced spectral density (RSD) (that is the spectral density removed the contribution from orientational diffusion). Insets in both panels are time-domain response functions which are the Fourier transform of the spectral density and the reduced spectral density, respectively. The peak near 230 and 270 cm−1 results from intramolecular vibrations. [From Q. Zhong, J. T. Fourkas, J. Phys. Chem. B 112, 15529 (2008).] 92 respectively in which the low-frequency (0 – 150 cm−1 ) spectra arise from intermolecular motions and the two peaks with frequency higher than 200 cm−1 represent Raman-active intramolecular vibrations. Note that the RSD has less low-frequency diffusive contribution, as can be also seen in the time-domain response functions in the insets. The persistent oscillations in these responses result from the intramolecular vibrations. 3.3.2 Many-body polarizability (dipole-induced-dipole model) Many-body polarizability is a central quantity in Raman spectroscopies with examples like OKE and transient grating spectroscopy. For a single molecule, the polarizability enters simply: the induced dipole moment is defined as µ ind = α · E, where α is the isolated molecular polarizability and E is the applied electric field. The isolated molecular polariz- ability α tensor has three principal axes components α1 , α2 , α3 and can be partitioned into the spherical isotropic part 1 1 α = (α1 + α2 + α3 ) = Trα (3.19) 3 3 and the traceless anisotropic part γ = α − α 1. Note that the polarizability anisotropy is defined as270 1 γ 2 = [(α1 − α2 )2 + (α2 − α3 )2 + (α3 − α1 )2 ]. (3.20) 2 For instance, for a cylindrically symmetrical molecule (only one symmetric axis, e.g. linear or symmetric top molecule), α1 = α2 = α⊥ , α3 = αk . Thus, the isotropic polarizability α = 31 (αk + 2α⊥ ) and the anisotropy γ = αk − α⊥ . The isolated molecular polarizability tensor for a symmetrical top is given by  γ ˆ iΩ ˆ i, αi = α − 1+γ Ω (3.21) 3 which is determined by its isotropic(α ) and anisotropic(γ ) polarizability and by orienta- 93 ˆ prescribing the principle axis of the molecule.20 tional unit vector Ω For isotropic atoms, the anisotropic part γ = 0 since there is no orientational degrees of freedom, therefore α = α 1. As for a group of N molecules, the many-body polarizability Π is not simply the sum of all the isolated molecular polarizabilities since the dipoles on other molecules can create an extra induced electric field. In the absence of an external electric field, the induced electric field at molecule i arises from the dipole moments of other molecules in the liquid, N Eind (ri ) = ∑ Ti j · µ j , (3.22) j6=i where Ti j is the dipole–dipole interaction tensor between molecules i and j   1 3ˆrrˆ − 1 Ti j ≡ T(ri j ) = ∇∇ = (3.23) ri j r3   2 2  x − r /3 xy xz  3   = 5 yx y2 − r2 /3 yz . r     zx zy z2 − r2 /3 r = ri j = ri − r j = (x, y, z), r = |r|, rˆ = r/r, (3.24) where 1 is the unit tensor. On applying an external electric field, the total induced dipole moment at molecule i is given by " # N µ i = α i · Etotal (ri ) = α i · Eext + ∑ Ti j · µ j . (3.25) j6=i This is the so-called dipole-induced-dipole (DID) model.271 The effective molecular polarizability of a single molecule i, π i , (i = 1, . . ., N) is defined as follows µ i = π i · Eext . (3.26) 94 Substituting Equation 3.26 into Equation 3.25, then gives the expression for this effective polarizability ∗ " # N π i = α i · 1 + ∑ Ti j · π j . (3.28) j6=i Finally the full many-body polarizability of the entire liquid is the sum of all the effective molecular polarizabilities N Π = ∑ π i. (3.29) i=1 Calculation of many-body polarizability To calculate the effective molecular polarizability in Equation 3.28, there are multiple methods we can use corresponding to various orders in the DID approximation. (1) First order DID approximation. N π i = α i + α i · ∑ Ti j · α j . (3.30) j6=i (2) Second order DID approximation. N N N π i = α i + α i · ∑ Ti j · α j + α i · ∑ Ti j · α j · ∑ T jk · α k . (3.31) j6=i j6=i k6= j (3) Infinite order DID approximation using iteration method.272 " # N (n+1) (n) πi = αi · 1+ ∑ Ti j · π j . (3.32) j6=i The effective molecular polarizabilities are evaluated by solving the equation above it- eratively. The initial guess of π i is just the isolated molecular polarizabilities α i . The ∗ The DID approximation can also be written as the following infinite series, N N N π i = α i + α i · ∑ Ti j · α j + α i · ∑ Ti j · α j · ∑ T jk · α k + · · · (3.27) j6=i j6=i k6= j 95 convergence condition is the sum of the absolute value of difference between every element in present and previous iteration steps less than 1 × 10−6 σ 3 or 3.95×10−6 A˚ 3 . (n+1) (n) ∑ Πµν − Πµν < 1 × 10−6 σ 3 . (3.33) µ ,ν =x,y,z (4) Infinite order DID approximation (Applequist method)273 by inverting a large matrix.271 Equation 3.25 may be rearranged to give the external electric field N α −1 i µ i − ∑ Ti j · µ j = Ei . (3.34) j6=i This equation is a system of N 3 × 3 matrix equations       α −1 1−T12 · · · −T1N  µ1   E1        −T21 α −1 ··· T2N  µ2   E2   2      . .. ..  = ..    .. , (3.35)  . ..     . . . .  .   .       −TN1 −TN2 ··· α −1 N µN EN or briefly ˜ Bµ˜ = E, =⇒ µ˜ = B−1 E˜ ≡ AE, ˜ (3.36) where µ˜ and E˜ are 3N × 1 vectors, A and B are a pair of reciprocal 3N × 3N matrices with 3 × 3 block matrix elements Ai j and Bi j = α −1 i δi j − Ti j (1 − δi j ). µν A3(i−1)+µ , 3( j−1)+ν = Ai j , (3.37) µν B3(i−1)+µ , 3( j−1)+ν = Bi j , (3.38) in which µ and ν denote Cartesian components. 96 The equivalent N matrix equations of Equation 3.36 are N  N  N µi = ∑ Ai j E j = ∑ Ai j E =⇒ π i = ∑ Ai j . (3.39) j=1 j=1 j=1 Assuming uniform field E j = E, then the total dipole moment is  N N  µ = ∑µi = ∑ ∑ Ai j E. (3.40) i i=1 j=1 Thus, the many-body polarizability is N N N Π = ∑ πi = ∑ ∑ Ai j . (3.41) i=1 i=1 j=1 Polarizability parameters for our model In order to calculate the many-body polarizability for our atomic liquid mixture model, we have to assign the isolated polarizabilities for each species, i.e. the solute, strong- solvating S solvent, and weak-solvating W solvent. The commonly used coumarin solutes in solvation studies usually carry the largest polarizabilities in the solution; one such example is coumarin 153 (C153) whose structure is shown in Figure 1.6. The ground-state and excited-state isotropic polarizabilities of C153 are 23.4 A˚ 3 and 35.7 A˚ 3 respectively.222 It is a common case that the S and W solvents have different polarizabilities, sometimes they are significantly distinct from one another. In widely-studied preferential solvation systems, for example the DMSO/water binary mixture,89, 91, 92, 251 αDMSO = 7.97 A˚ 3274 and αwater = 1.47 A˚ 3 .275 In this particular case, the ratio of solvents’ polarizabilities is about 5. In light of the fact that the most obvious preferential solvation occurs when only a small fraction of S solvent is present, it is reasonable to study cases in which the polarizability of the S solvent is much larger than that of the W solvent, because it is easier to distinguish the polarizability change due to a few large-polarizability solvent immersed in a large number of small-polarizability solvent background than the case with a few small-polarizability 97 molecules against an abundant large-polarizability solvent background. Consequently, in order to get a maximized signal, we have chosen the solute and solvent polarizabilities for our atomic liquid model, αu = 0.2 σ 3 = 7.90 A˚ 3 , (3.42) αS = 0.101 σ 3 = 3.99 A˚ 3 , (3.43) αW = 0.0186 σ 3 = 0.73 A˚ 3 . (3.44) The polarizability of the S solvent has a factor of 5.5 difference with respect to the W solvent, the same ratio as DMSO/water. The largest polarizability of the solute is kept constant on resonant solute excitation and is twice as large as the strong-solvent value. The reason for choosing a constant solute polarizability is to avoid complicating the inceptive calculations of solute-pump/solvent-probe spectroscopy. Changing the polarizability of the solute is definitely worth consider in future studies.222 The reader may have noticed that the polarizabilities of solvents chosen is half the size of the experimental values of DMSO/water. The reason for this choice is to circumvent the so-called polarizability catastrophe which refers to the divergence of the many-body polarizability with infinite order DID approximation, when the distance between two interacting sites approaches (4αi α j )1/6 .276 This issue of polarizability divergence is solved by Thole through modifying the dipole–dipole tensor.277, 278 When more complex polarizable systems are concerned, the Thole model is a good candidate.279–281 With the polarizability parameters above, we tested the four evaluation methods in a 10% S solvent mixture plus an excited-state solute, with total number of particles N = 256. The following many-body polarizabilities (unit σ 3 ) are averaged over 200 liquid configurations. In the infinite-order iterative calculation, the numbers of iteractions till convergence are typically 9 – 14 with absolute tolerance 10−5 . The zeroth order many- body polarizability is Π 0th = (∑i αi ) · 1 = 7.0854 · 1. 98 (1) First order DID:    7.073 ± 0.033 −0.005183 ± 0.027377 0.02739 ± 0.02424    Π 1st =   −0.005183 ± 0.027377 7.088 ± 0.029 0.01504 ± 0.01641     0.02739 ± 0.02424 0.01504 ± 0.01641 7.096 ± 0.037 (2) Second order DID:    7.130 ± 0.029 −0.002191 ± 0.029720 0.03451 ± 0.02616    Π 2nd =   −0.002191 ± 0.029720 7.151 ± 0.032 0.008034 ± 0.020589     0.03451 ± 0.02616 0.008034 ± 0.020589 7.162 ± 0.030 (3) Infinite order DID – iteration:    7.130 ± 0.032 −0.001822 ± 0.032146 0.03753 ± 0.02820    Π inf =   −0.001822 ± 0.032146 7.153 ± 0.035 0.01003 ± 0.02237     0.03753 ± 0.02820 0.01003 ± 0.02237 7.164 ± 0.033 (4) Infinite order DID – matrix inversion:    7.130 ± 0.032 −0.001822 ± 0.032146 0.03753 ± 0.02820    Πinf =   −0.001822 ± 0.032146 7.153 ± 0.035 0.01003 ± 0.02237     0.03753 ± 0.02820 0.01003 ± 0.02237 7.164 ± 0.033 From above results, we can conclude that the first order DID is qualitatively correct compared with the infinite order DID, and the second order is quantitatively correct compared with the infinite order DID. Two infinite order DID methods give the same results. However, the iterative approach is less time-consuming compared with the matrix inversion approach. Furthermore, in the iterative approach, controllable maximum cycles of iteration make it easy to calculate arbitrary order DID approximated polarizabilities. We 99 will choose the first order DID approximation for computational efficiency and iterative infinite-order method for accuracy in evaluation of polarizabilities. It is worth noting that the fluctuations of the off-diagonal elements make a difference when we looks at the difference in the OKE spectra with a ground- and excited-state solute, even though the DID approximated many-body polarizability has almost identical values as the zeroth-order many-body polarizability. 3.3.3 OKE spectra for our model The time-domain OKE response functions for 10% S and 50% S systems (total number of atoms N = 256), each equilibrated in the presence of ground- and excited-state solute, are shown in Figure3.7 and Figure 3.8. These typical response functions are computed using Equation 3.9 and 3.10 taking advantage of tensor invariants and the time derivatives are calculated by central-difference finite-difference methods. Figure 3.7 also includes the isotropic responses which are much weaker than the anisotropic responses making subsequent spectra analysis more difficult, thus we will focus on the anisotropic OKE spectra later in this thesis. For the same anisotropic OKE responses, Figure 3.7 is a linear plot and Figure 3.8 is a log plot. From both figures, we can observe that there is a sharp peak in the ultrafast region (1–2 ps) corresponding to the vibrational intermolecular dynamics and after 2–3 ps there is mostly long-time diffusion. Since we are mainly interested in the ultrafast dynamics in liquids, it is useful to remove the long-time diffusive portion of the response function and obtain the intermolecular vibrational and librational part of the response. By then performing Fourier transformation on the response without the diffusive portion, one obtains the reduced spectral density, the conventional OKE spectra. The frequency-domain OKE spectra are shown in Figure 3.9. Both the full spectral den- sity and the reduced spectral density are derived from the imaginary parts of corresponding time-domain responses. (Equation 3.15 and 3.18). The diffusive portion fitting results are shown in Table 3.3. (Table 3.4 is the first order equivalent fitting results.) The diffusion- 100 ) 0.0035 1/2 m) 0.0030 10%S (g) 10% S (e) 0.0025 aniso (Response funciton 0.0020 iso 0.0015 0.0010 0.0005 0.0000 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 t (ps) t (ps) (a) ) 1/2 0.030 m) 50%S (g) 50%S (e) 0.025 5 aniso ( 0.020 iso Response funciton 0.015 0.010 0.005 0.000 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 t (ps) t (ps) (b) Figure 3.7 Anisotropic and isotropic OKE responses in 10% S and 50% S systems on equi- librated ground (g) and excited (e) states. Using exact infinite-order DID approximation for polarizability, the results are averaged over 107 configurations sampled every 5 time steps. 101 0.01 10% S g e ) 1/2 1E-3 m 5 1E-4 ( Response function 0.1 50% S g 0.01 e 1E-3 1E-4 0 2 4 6 8 10 12 14 t (ps) Figure 3.8 Time-domain OKE response functions for two different solvent mixtures, 10% S and 50% S systems in the equilibrated ground state (g) and excited state (e). Exact infinite-order DID approximation are used to calculate many-body polarizability, and the results are averaged over 107 configurations sampled every 5 time steps. Note there is a order of magnitude difference in scales between the panels. 102 removal procedure only alters the spectra in the low frequency region (<30 cm−1 ) and the important high-frequency spectra are not changed. The difference between the ground- and excited-state OKE responses is more distinguishable in the frequency domain than in the time domain. In spite of a large solvent background which is presumably not be affected by the solute excitation, the e–g differences are 5%–10% of the full excited- and ground-state spectra, an amount which is significant and reproducible. These substantial e–g differences suggest the capability of OKE spectra to single out the effect of the solute excitation on solvent dynamics. The most striking feature is that the e–g difference has opposite signs in 10% S and 50% S systems, an occurrence whose molecular origin will be discussed soon in the last section of this chapter. A1 τ1 A2 τ2 hν i τrise adj.R2 10% S (g) 5.7743E-4 1.1426 3.8276E-4 4.2377 35.424 0.4708 0.99939 10% S (e) 5.1617E-4 1.6498 3.6480E-4 5.7469 37.0284 0.4504 0.99969 50% S (g) 0.00207 1.6605 0.00151 4.6133 38.8421 0.4294 0.99696 50% S (e) 0.00266 2.2324 6.8201E-4 7.5971 39.2645 0.4248 0.99502 Table 3.3 Fitting results for the long-time diffusive portion of OKE response function in 10% S and 50% S mixtures on the equilibrium ground and excited state. The OKE response function is calculated by averaging 107 infinite order DID polarizabilities, sampled every 5 time steps. Fit starts from 2 ps for 10% S and from 3 ps for 50% S. As in Equation 3.16 and 3.17, A1 and A2 have unit σ 5 /(ε m)1/2 , τ1 , τ2 and τrise have unit ps, and average frequency hν i based on the spectral density has unit cm−1 . When fitting the time-domain response, the data points are dt = 0.054 ps apart and consequently in the frequency-domain the spectral points are dν = 2.41294 cm−1 apart. The rise time is defined as τrise = (2hν i)−1. The adj.R2 is the indicator of goodness of the fit. We begin our analysis of the solute-pump/solvent-probe spectra of the preferential solvation model by comparing the reduced spectral densities of 10% S and 50% S with their e–g differences on the same scale as shown in Figure 3.10. Here the first order DID approximation is utilized to improve the computational efficiency which still is able to capture the general behavior of the non-diffusive OKE response. The most direct observation about Figure 3.10 is that the magnitude of OKE spectra for 50% S solvent mixture is much larger than the 10% S mixture with both ground- and excited-state solutes. 103 6 g 5 10% S / ) e 6 4 e-g - 4 3 (10 2 RSD Anisotropic Spectral Density 1 SD 0 35 g 50% S 30 e 25 e-g 20 15 10 RSD 5 SD 0 -5 0 50 100 150 200 -1 c (cm ) Figure 3.9 Frequency-domain OKE spectra for 10% S and 50% S systems (Fourier transformed from Figure 3.8). For each system, shown are the results with equilibrated ground-state (g) solute, excited-state (e) solute and the difference between the two (e-g). For each case, there are full spectral density (SD, dashed lines) and reduced spectral density (RSD, solid lines) with long-time diffusive behavior removed. Exact infinite-order DID approximation are used to calculate many-body polarizability, and the results are averaged over 107 configurations sampled every 5 time steps. Note that the scales in two panels are different. 104 A1 τ1 A2 τ2 hν i τrise adj. R2 10% S (g) 3.8239E-4 2.1184 1.7869E-4 5.4728 31.8039 0.5244 0.99993 10% S (e) 3.8454E-4 2.0506 2.6838E-4 6.0915 32.6986 0.5101 0.99995 50% S (g) 0.00292 1.6523 0.00136 5.0666 34.4044 0.4848 0.99992 50% S (e) 0.0028 1.8203 0.00114 6.0463 34.6765 0.4810 0.99987 Table 3.4 Fitting results for the long-time diffusive portion of OKE response function in 10% S and 50% S mixtures on the equilibrium ground and excited state. The OKE response function is calculated by averaging 108 first order DID polarizabilities, sampled every 10 time steps. Fit starts from 3 ps for all cases. As in Equation 3.16 and 3.17, A1 and A2 have unit σ 5 /(ε m)1/2 , τ1 , τ2 and τrise have unit ps, and average frequency hν i based on the spectral density has unit cm−1 . When fitting the time-domain response, the data points are dt = 0.054 ps apart and consequently in the frequency-domain the spectral points are dν = 2.41294 cm−1 apart. The rise time is defined as τrise = (2hν i)−1. The adj.R2 is the indicator of goodness of the fit. The areas under the 50% S RSD curves differ from the areas of 10% S RSDs by a factor of 10. Consequently the e–g difference of 50% S system is larger than that of 10% S system. But the magnitude difference of full OKE spectra has little dynamical importance since there are 5 times more of the spectroscopically visible S solvents in the 50% S case than in the 10% S case, which predominantly composed of an insipid background. 3.3.4 Difference spectra The e–g difference, on the other hand, tells the dynamical change after the solute excitation. Figure 3.11 highlights these essential e–g differences in 10% S and 50% S mixtures. These are the primary results of this chapter: the long waiting time T limits of the solute- pump/solvent-probe spectra. The most astonishing trait of these difference spectra is that the signs of the spectra in the low-frequency region (< 70 cm−1 ) are reversed when switching the solvent mixture, yet the relative magnitudes of these difference spectra are also worth noting. The absolute area of the e–g difference in 50% S is only 3 times as large as that of 10% S although the ratio of full spectra area is about 10. So the fraction for the e–g difference is about 3 times higher in 10% S than in 50% S, implying that the solute excitation requires more dramatic structural rearrangement for each S solvent in the 10% S 105 / ) 30 6 50% S - 4 25 g (10 20 e Anisotropic Spectral Density e-g 15 10 5 10% S 10% S 0 50% S -5 0 50 100 150 200 -1 c (cm ) Figure 3.10 Comparison of the OKE reduced spectral densities for 10% S and 50% S solvent mixtures with ground-state solute (g), excited-state solute (e) and their difference spectrum (e-g). First-order DID approximation are used to calculate many-body polariz- ability, and the results are averaged over 108 liquid configurations. For each case, fitting parameters of the long-time diffusive portion are shown in Table 3.4. 106 case, which is also consistent with the previous solvation dynamics studies that show that the dilute S case is the slowest.85, 87, 93, 98 A second feature of the e–g difference is that 10% S has a flat shaped difference spectra with peak position around 60 cm−1 , and 50% S has a sharper peak with position around 25 cm−1 . The shapes of spectra indicate that in 10% S system, the intermolecular vibrational modes are excited more uniformly than those in 50% S system. Furthermore, above 75 cm−1 both systems have positive e–g differences, and the e–g difference of 10% S is even larger than that of 50% S. As high frequency is dominated by pair motions, does it mean 10% S gets more pair motions excited than 50% S? If so, what are these pairs? e/g difference spectrum / ) 1.0 6 RSD - 4 (10 0.5 SD Anisotropic Spectral Density 0.0 10% S -0.5 50% S -1.0 0 50 100 150 -1 c (cm ) Figure 3.11 Comparison of the long waiting time limit anisotropic RP-PORS spectra for 10% S and 50% S solvent mixtures. The equilibrated excited-state-solute/ground- state-solute RP-PORS spectra correspond to the e–g differences in the OKE reduced spectral density (RSD, solid lines), and the e–g differences in full spectral density (SD, dashed lines) are also shown. First-order DID approximation are used to calculate many- body polarizability, and the results are averaged over 108 liquid configurations. The e–g differences in RSD plotted here are the same curves as in Figure 3.10. To understand this result, one can selectively turn on/off certain molecular effects 107 0.0006 / ) g 6 SD ( 0.0005 RSD e Anisotropic Spectral Density e-g 0.0004 10%S = 0.2, 0.0003 U =0.101, S 0.0002 =0.0186 W 7 # Conf = 1x10 0.0001 Sample steps = 5 0.0000 0 50 100 150 200 -1 c (cm ) (a) 0.0006 / ) 6 SD g ( 0.0005 RSD e Anisotropic Spectral Density e-g 0.0004 10%S = 0 U 0.0003 =0.101 S 0.0002 =0.0186 W 7 # Conf = 1x10 Sample steps = 5 0.0001 0.0000 0 50 100 150 200 -1 c (cm ) (b) 0.0006 / ) SD g 6 0.0005 RSD e ( e-g Anisotropic Spectral Density 0.0004 First order DID 10%S = 0.2 0.0003 U =0.101 S 0.0002 =0.0186 W 7 # Conf = 1x10 0.0001 Sample steps = 5 0.0000 0 50 100 150 200 -1 c (cm ) (c) Figure 3.12 Comparison of OKE spectra in 10% S system, (a) infinite order DID approximation and full solute polarizability, (b) infinite order DID approximation with the solute’s polarizability set to zero and (c) first order DID approximation with full solute polarizability. Number of configurations to average is 1 × 107 and configurations are sampled every 5 steps. 108 0.0035 ) SD g / 6 0.0030 e ( RSD Anisotropic Spectral Density 0.0025 e-g 50%S 0.0020 = 0.2, U 0.0015 =0.101, S =0.0186 0.0010 W 7 # Conf = 1x10 0.0005 Sample steps = 5 0.0000 0 50 100 150 200 -1 c (cm ) (a) 0.0035 / ) SD g 6 0.0030 ( RSD e Anisotropic Spectral Density 0.0025 e-g 50%S 0.0020 = 0 U 0.0015 =0.101 S =0.0186 0.0010 W 7 # Conf = 1x10 0.0005 Sample steps = 5 0.0000 0 50 100 150 200 -1 c (cm ) (b) 0.0035 ) SD g / 6 0.0030 RSD e ( e-g Anisotropic Spectral Density 0.0025 First order DID 50%S 0.0020 = 0.2 U 0.0015 =0.101 S =0.0186 0.0010 W 7 # Conf = 1x10 0.0005 Sample steps = 5 0.0000 0 50 100 150 200 -1 c (cm ) (c) Figure 3.13 Comparison of OKE spectra in 50% S system, (a) infinite order DID approximation and full solute polarizability, (b) infinite order DID approximation with the solute’s polarizability set to zero and (c) first order DID approximation with full solute polarizability. Number of configurations to average is 1 × 107 and configurations are sampled every 5 steps. 109 in order to isolate the roles of different contributions, in our case, to the many-body polarizability. First of all, one can compare the first order DID approximation with the exact infinite order results in order to see if higher order DID terms have a crucial influence on the difference spectra. Moreover, one can set the solute’s polarizability as zero making the solute spectroscopically invisible so as to see if the solute–solvent terms are important to the e–g difference at all. Accordingly, Figure 3.12, 3.13 offer evidence of which contribution to the many-body polarizability is predominant in the OKE spectra. The effects are difficult to see on this scale, but Figure 3.14 is the summary of the above mentioned figures plotting only the e–g differences in all the cases. One can observe that the first order DID and the infinite order DID spectra are qualitatively equivalent except for the 10% S spectrum in the low-frequency region, although the higher order DID terms contribute almost a half of the result. Now we are convinced that the first order DID approximation warrants a qualitative prediction to the solute-pump/solvent-probe spectra. Moreover, when turning off the solute’s polarizability αu =0, almost all the high-frequency portions get demolished. We thereby conclude that the physically most important polarizability terms detected by the spectra are those involving the solute/S-solvent pairs. From Equation 3.30, the dominant polarizability terms are 2αu αS ∑ T(r0 j ), r0 j = r(solute) − r(solvent- j). (3.45) j∈S Then exactly which group of the S solvents is responsible for the measured difference? Is it a local portion of the solvent near the solute? To answer these questions, one can take advantage of the projection operator technique that will be discussed in the upcoming section. 110 0.4 0.3 10% S / ) 0.2 6 - 4 0.1 (10 Anisotropic Spectral Density 0.0 total 1.5 st 1 order DID 1.0 =0 (full DID) u 0.5 50% S 0.0 -0.5 -1.0 -1.5 0 50 100 150 -1 c (cm ) Figure 3.14 Comparison of different polarizability contributions to the long waiting limit anisotropic RP-PORS spectra for 10% S and 50% S solvent mixtures. “Total” uses the infinite order DID approximation for many-body polarizability; “1st order DID” curve neglects all higher order DID terms; and “αu =0” uses the infinite order DID approximation but with a zero-polarizability solute. The results shown are reduced spectral density differences averaged over 107 , 108 , 107 liquid configurations respectively. 111 3.4 Projection Operator Analysis of Molecular Contribu- tions The projection operator discussed here, e.g. Pˆ X , acts on a summation of the atom coordi- nates and results in a selected sum of certain degrees of freedom ( j ∈ X ). There must be a complementary projection operator Pˆ Y that has no common degree of freedom as in Pˆ X . Moreover, projection operator is idempotent — it can be applied multiple times without changing the result of applying just once. N N Pˆ X ∑ = ∑, Pˆ Y ∑=∑ (definition), (3.46) j=1 j∈X j=1 j6∈X Pˆ X + Pˆ Y =1ˆ (complementariness), (3.47) Pˆ X Pˆ Y =0 (orthogonality), (3.48) Pˆ 2X =Pˆ X (idempotent). (3.49) For a correlation function C(t) = hA(0)B(t)i, if variables A(t) and B(t) contain a sum ∑ j , one can use projection operator to isolate certain types of contributions. For example, if projected into two sets X and Y, the choice of the partitioning sets could be based on regional separation like the first solvation shell vs. outer shells or based on different degrees of freedom like translational vs. rotational motions. C(t) = hA(0)B(t)i = h(Pˆ X + Pˆ Y )A(0) (Pˆ X + Pˆ Y )B(t)i (3.50) = hPˆ X A(0)Pˆ X B(t)i + hPˆ Y A(0)Pˆ Y B(t)i +hPˆ X A(0)Pˆ Y B(t)i + hPˆ Y A(0)Pˆ X B(t)i = CXX (t) +CYY (t) +CXY,YX (t). (3.51) 112 Partitioning into two complementary sets is the simplest case. One certainly can do more, such as partitioning the whole set into three subsets, X, Y, Z. Similarly, we have 1ˆ = Pˆ X + Pˆ Y + Pˆ Z , (3.52) C(t) = hA(0)B(t)i = h(Pˆ X + Pˆ Y + Pˆ Z )A(0) (Pˆ X + Pˆ Y + Pˆ Z )B(t)i (3.53) = hPˆ X A(0)Pˆ X B(t)i + hPˆ Y A(0)Pˆ Y B(t)i + hPˆ Z A(0)Pˆ Z B(t)i +hPˆ X A(0)Pˆ Y B(t)i + hPˆ Y A(0)Pˆ X B(t)i +hPˆ X A(0)Pˆ Z B(t)i + hPˆ Z A(0)Pˆ X B(t)i +hPˆ Y A(0)Pˆ Z B(t)i + hPˆ Z A(0)Pˆ Y B(t)i (3.54) = CXX (t) +CYY (t) +CZZ (t) +CXY,YX (t) +CXZ,ZX (t) +CY Z,ZY (t). (3.55) The first three “self” terms are resulted from the same subset of degrees of freedom, and the last three terms are called cross terms. Polarizability-velocity correlation function A spectroscopic “velocity” autocorrelation function, G(t) is defined as the autocorrelation function of time derivative of variable A(t) (chapter 1) and it can be also partitioned into different contributions using projection operators. ˙ A(t)i ˙ d ˙ d2 GAA (t) = hA(0) = hA(0)A(t)i = − 2 hA(0)A(t)i (3.56) dt dt = hA˙ X (0)A˙ X (t)i + hA˙ Y (0)A˙Y (t)i +hA˙ X (0)A˙ Y (t)i + hA˙ Y (0)A˙ X (t)i = GXX (t) + GYY (t) + GXY,Y X (t), (3.57) 113 where N ∂ A(t) ∂ r j ∂ A(t) A˙ X (t) = Pˆ X A(t) ˙ = Pˆ X ∑ · = ∑ · v j (t), (3.58) j=1 ∂ r j ∂ t j∈X ∂ r j N ∂ A(t) ∂ r j ∂ A(t) A˙Y (t) = Pˆ Y A(t) ˙ = Pˆ Y ∑ · = ∑ · v j (t). (3.59) j=1 ∂ r j ∂t j6∈X ∂ r j In OKE spectroscopy, the quantity A(t) is the many-body polarizability Π (t), so the ˙ (0) ⊗ Π corresponding “velocity” autocorrelation function is G(t) = hΠ ˙ (t)i, which is also called the polarizability-velocity correlation function. The time derivatives of many-body polarizability of certain sets of degrees of freedom X and Y are ˙ X (t) = Pˆ X Π Π ˙ (t) = ∑ ∇ j Π (t) · v j(t), (3.60) j∈X ˙ Y (t) = Pˆ Y Π Π ˙ (t) = ∑ ∇ j Π (t) · v j(t). (3.61) j∈Y Appendix E shows the detailed derivation of the spatial derivative of many-body polariz- ability with respect to a particular molecule k’s coordinate, ∇k Π . ∂ ∂ ∇k Π = ∑ ∇k π i or (Π)αβ = ∑ (πi )αβ . (3.62) i ∂ rk µ i ∂ rk µ The spatial derivative of the effective molecular polarizability of molecule i in component form is (Equation E.7) " x,y,z 3 r 2 (r ∂ i j i j µ δαγ + ri jα δµγ + ri jγ δµα ) − 15 ri j µ ri jα ri jγ (πi )αβ = αi · ∑ ∑ ∂ rk µ j6=i γ ri7j # 3 ri jα ri jγ − ri2j δαγ ∂ ·(δki − δk j ) · (π j )γβ + · (π j )γβ . (3.63) ri5j ∂ rk µ To select the local portion of the OKE spectra, one can partition contributions into (1) first shell and everything else or (2) first shell, second shell and everything else. We expect 114 that the innermost first-solvation shell would be the largest contribution to the change in many-body polarizability, thus the e–g difference spectra. For the polarizability-velocity correlation function G(t), we have (1) First choice of partition: X includes all particles in the first solvation shell, and Y is everything else. G(t) = GXX (t) + GYY (t) + GXY,YX (t) (3.64) ˙ X (0) ⊗ Π = hΠ ˙ X (t)i + hΠ ˙ Y (0) ⊗ Π ˙ Y (t)i ˙ X (0) ⊗ Π +hΠ ˙ Y (t)i + hΠ ˙ Y (0) ⊗ Π ˙ X (t)i, (3.65) where the first/outer shell partition is based on the instantaneous (t = 0) liquid configuration — whether a given solvent is inside or outside the cutoff radius of the first shell that is the first minimum of the solute-solvent radial distribution function, 1.556 σ (Figure 3.1). Figure 3.15 shows all components of the polarizability-velocity correlation G(t) for partition of 1st shell and all else. The first order DID approximation is used to calculate the spatial derivatives of many-body polarizability. Comparing Figure 3.15(a) and Figure 3.15(b), we can observe that in 10% S the pure first-shell term (XX) changes relatively more dramatically when the solute gets excited, than the pure first-shell contribution does in 50% S. As for the pure outer shell component (YY), the absolute magnitude keeps almost the same when the solute gets excited for both systems and noticeably larger than the pure first- shell component. Consequently, we wonder if the second-shell term is the contribution to the pure outer shell component. This is the reason we test the second shell partition below. Figure 3.15(c) are the e–g differences in these projections. This result is of only qualitative significance. From this figure one can observe that the outer-shell component in 50% S has a large oscillation (which is not trustworthy) with amplitude about 1.0 σ 4 /m which results from different oscillation frequencies of the excited- and ground-state GYY (t); all components in 10% S have a positive value in the first hundred femtoseconds. The cross 115 terms are negligible compared with “self” terms. It seems that the larger magnitude of 50% S projections does not guarantee the same effect on e–g differences. (2) Second choice of partition: X=1st shell, Y=2nd shell and Z=else. The correspond- ing projected polarizability-velocity correlation function is G(t) = GXX (t) + GYY (t) + GZZ (t) + GXY,YX (t) + GXZ,ZX (t) + GYZ,ZY (t). (3.66) Figure 3.16 shows the projected components under the first-shell/second-shell/else partition for both systems, where 10% S has more obvious changes in all “self” terms before and after the solute’s excitation, whereas 50% S system exhibits no significant change. Figure 3.16(a) and 3.16(c) tells us that the first shell and second shell both have an increase in polarizability on solute excitation, consistent with the fact that more atoms come into these closest shells, and so the outer shells lose some polarizabilities. Figure 3.16(c) and 3.16(d) shows the e–g differences of G(t), from which we can observe that all components in 10% S have comparable contributions and only the outer-shell contribution has a significant fluctuation in 50% S. To summarize these two choices of projection partitioning, 10% S has more obvious changes in G(t) in early times, and 50% S seems doesn’t change much except for the outer- shell term oscillating with a big amplitude all the time (questionably). These observations in time-domain polarizability-velocity correlation function cannot lead us to an answer to the question about the signs of e–g difference OKE spectra. Probably the frequency-domain projections can tell us more information. Projection Analysis of the OKE Spectra The relation between OKE response function and polarizability-velocity correlation func- R ˙ (0) ⊗ Π (t)i = G(t)dt. But numerical integration has more tion is as follows R(t) = hΠ instability than directly calculation of OKE response function using the following projec- 116 0.6 XX g YY e XYYX 0.4 10% S X = first shell 1st order DID G(t) Y = else # conf=2500 0.2 Sample steps=2 0.0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 t (ps) (a) 8 XX g YY XYYX e 6 50% S 4 X = first shell 1st order DID G(t) Y = else # conf=2500 Sample steps=2 2 0 -2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 t (ps) (b) (c) Figure 3.15 Polarizability-velocity correlation function G(t) projected to first shell (X) and everything else (Y) contributions. Shown are the first-shell, outer-shells and cross components of G(t) (unit σ 4 /m). (a)10% S solvent mixture, (b) 50% S solvent mixture, (c) e–g differences of both 10% S and 50% S mixtures. Results are computed by averaging 2500 liquid configurations. 117 0.6 8 XX g XX g YY YY e 6 e ZZ ZZ 0.4 X = 1st shell 10% S 4 G(t) Y = 2nd shell 1st order DID X = 1st shell 50% S G(t) Z = else Y = 2nd shell 1st order DID 0.2 2 Z = else 0 0.0 -2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 t (ps) t (ps) (a) (b) 0.2 2 XX (e-g) X = 1st shell XX (e-g) X = 1st shell YY (e-g) Y = 2nd shell YY (e-g) Y = 2nd shell 0.1 ZZ (e-g) Z = else 1 ZZ (e-g) Z = else G(t) G(t) 0.0 0 10% S 50% S -0.1 -1 1st order DID 1st order DID -0.2 -2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 t (ps) t (ps) (c) (d) Figure 3.16 Polarizability-velocity correlation function G(t) projected to first shell (X), second shell (Y) and else (Z) contributions. Shown are pure first-shell, second-shell and outer-shells components of G(t) (unit σ 4 /m). (a) 10% S solvent mixture, (b) 50% S solvent mixture, (c) e–g differences in 10% S solvent mixture and (d) e–g differences in 50% S solvent mixture. Results are computed by averaging 2500 liquid configurations. 118 tion method. We can apply projection operators on both a variable and tis time derivative. ˙ X (t) = Pˆ X Π Π ˙ (t) = ∑ ∇ j Π (t) · v j(t), (3.67) j∈X N Π X (t) = Pˆ X ∑ π j (t) = ∑ π j (t). (3.68) j=1 j∈X (1) Choosing partition X=1st shell and Y=else like in the previous section, we have R(t) = RXX (t) + RYY (t) + RXY,YX (t), (3.69) ˙ X (0) ⊗ Π X (t)i, RXX (t) = hΠ (3.70) ˙ Y (0) ⊗ ΠY (t)i, RYY (t) = hΠ (3.71) ˙ X (0) ⊗ ΠY (t)i + hΠ RXY,Y X (t) = hΠ ˙ Y (0) ⊗ Π X (t)i. (3.72) Figure 3.17 shows these time-domain OKE response functions with partitioning into the first shell and the outer shells. Performing a Fourier transformation on these responses, one gets the frequency-domain spectral densities which is plotted in Figure 3.18. From these frequency-domain spectra, the outer shells compose the largest contribution to the OKE spectral densities and the cross terms are negligible. The first-shell contributions are weighted differently in the two solvent mixtures: in 10% S, the first-shell contribution is modest, about a third of the largest outer-shells contribution; in 50% S, the first-shell contribution is comparable to the cross term which is less than 10% of the largest outer- shells contribution. Key result here is that the change of first-shell contribution due to the g → e transition is noticeably larger than the other two contributions in 10% S mixture, whereas in 50% S mixture all three contributions contribute roughly the same amount in the frequency domain. The differences between the excited-state and ground-state spectral densities are shown in Figure 3.19. In both mixtures, the contributions of cross terms to the e– 119 g difference are not negligible any more, and are likely to counteract the outer-shells contributions. The frequency-domain results most relevant to ultrafast intermolecular dynamics though are the reduced spectral densities. Figure 3.20 shows the projection of OKE reduced spectral densities into the first shell and outer shells partition for 10% S and 50% S mixtures on equilibrated ground and excited states. Taking the e–g difference, we have the most important result of this section, Figure 3.21 presenting the projection of long waiting time limit solute-pump/solvent-probe spectra into first shell and outer shells. For the 10% S solvent mixture, the full e–g difference must be determined by the first shell contribution, at least for frequencies greater than 30 cm−1 , since the first-shell projection has the same positive sign as the full response (Figure 3.11). The outer-shells projection, on the other hand, is negative and is partially canceled by the cross term, which is easy to see in e–g differences of projected spectral densities in Figure 3.19. By contrast, for 50% S solvent mixture, the outer shells seem to be more important in that the outer shells have the largest contribution to the overall negative response (30 cm−1 – 70 cm−1 ); the positive cross term partially cancels the first shell’s negative contribution. (2) Choosing partition X=1st shell, Y=2nd shell, Z=everything else, we have R(t) = RXX (t) + RYY (t) + RZZ (t) + RXY,YX (t) + RXZ,ZX (t) + RY Z,ZY (t). (3.73) Figure 3.22 shows the OKE spectral densities with partitioning into 1st shell, 2nd shell and else. One can observe that the “self” terms (RXX , RYY , RZZ ) are less noisy than the cross terms (RXY,Y X , RXZ,ZX , RY Z,ZY ) in both 10% S and 50% S solvent mixtures. Figure 3.23 shows the e–g difference spectra of these three “self” terms. These results are noisy due to much less averaging than the first-shell/outer-shells partitioning, but it is still useful in providing some qualitative hints. The first shell contribution reproduces the signs in Figure 3.21. The second shell contribution in 10% S is almost zero and in 50% S it has negative e–g 120 ) 0.0018 1/2 m) 0.0016 first shell (g) first shell (e) 0.0014 outer shells (g) ( 0.0012 outer shells (e) Response funciton cross (g) 0.0010 cross (e) 0.0008 0.0006 0.0004 10% S 0.0002 0.0000 0 2 4 6 8 10 12 14 t (ps) (a) ) 1/2 0.020 first shell (g) m) first shell (e) outer shells (g) 0.015 outer shells (e) ( cross (g) Response funciton cross (e) 0.010 0.005 50% S 0.000 0 2 4 6 8 10 12 14 t (ps) (b) Figure 3.17 Projection of the OKE responses for 10% S and 50% S solvent mixtures into first-shell, outer-shells and cross contributions. The ground-state (g, solid lines) and excited-state (e, dashed lines) curves are computed with first-order DID polarizabilities and averaged over 5 × 105 and 2.5 × 106 liquid configurations for 10% S and 50% S solvent mixtures respectively. 121 / ) 0.00025 first shell (g) 6 first shell (e) (Anisotropic Spectral Density outer shells (g) 0.00020 outer shells (e) cross (g) 0.00015 cross (e) 0.00010 10% S 0.00005 0.00000 0 50 100 150 200 -1 c (cm ) (a) 0.0025 / ) first shell (g) 6 first shell (e) ( 0.0020 Anisotropic Spectral Density outer shells (g) outer shells (e) 0.0015 cross (g) cross (e) 0.0010 50% S 0.0005 0.0000 0 50 100 150 200 -1 c (cm ) (b) Figure 3.18 Projection of the OKE spectral densities for 10% S and 50% S solvent mixtures into first-shell, outer-shells and cross contributions. The ground-state (g, solid lines) and excited-state (e, dashed lines) spectral densities are computed with first-order DID polarizabilities and averaged over 5 × 105 and 2.5 × 106 liquid configurations for 10% S and 50% S solvent mixtures respectively. 122 / ) 0.4 6 0.3 - 4 10% S (10 0.2 Anisotropic Spectral Density 0.1 0.0 -0.1 first shell -0.2 outer shells -0.3 cross -0.4 0 50 100 150 200 -1 c (cm ) (a) / ) 1.5 6 - 4 1.0 50% S (10 Anisotropic Spectral Density 0.5 0.0 -0.5 first shell outer shells -1.0 cross -1.5 0 50 100 150 200 -1 c (cm ) (b) Figure 3.19 Projection of the e–g difference in spectral densities for 10% S and 50% S solvent mixtures into first-shell, outer-shells and cross contributions. The e–g differences are computed with first-order DID polarizabilities and averaged over 5 × 105 and 2.5 × 106 liquid configurations for 10% S and 50% S solvent mixtures respectively. 123 / ) -4 2.0x10 6 first shell (g) ( Anisotropic Spectral Density first shell (e) -4 1.5x10 outer shells (g) outer shells (e) -4 1.0x10 -5 10% S 5.0x10 0.0 0 50 100 150 200 -1 c (cm ) (a) / ) 6 -3 2.0x10 first shell (g) ( Anisotropic Spectral Density first shell (e) outer shells (g) -3 1.5x10 outer shells (e) -3 1.0x10 50% S -4 5.0x10 0.0 0 50 100 150 200 -1 c (cm ) (b) Figure 3.20 Projection of the OKE reduced spectral densities for 10% S and 50% S solvent mixtures into first-shell and outer-shells. The cross terms are not shown. The equilibrium ground (g) and excited (e) state OKE spectra are computed with first-order DID polarizabilities and averaged over 5 × 105 and 2.5 × 106 liquid configurations for 10% S and 50% S solvent mixtures respectively. 124 0.4 0.3 10% S / ) 0.2 0.1 6 0.0 - 4 (10 -0.1 Anisotropic Spectral Density -0.2 first shell -0.3 outer shells -0.4 1.5 50% S 1.0 0.5 0.0 -0.5 first shell -1.0 outer shells -1.5 0 50 100 150 200 -1 c (cm ) Figure 3.21 Projection of the long waiting time limit anisotropic solute-pump/solvent- probe spectra for 10% S and 50% S solvent mixtures into first-shell and outer-shells contributions. Cross terms are not shown. The long waiting time limit results are the e–g differences in the partitioned reduced spectral densities that are computed with first-order DID polarizabilities and averaged over 5 × 105 and 2.5 × 106 liquid configurations for 10% S and 50% S solvent mixtures, respectively. 125 difference below 25 cm−1 and positive e–g difference above 25 cm−1 , which is similar with the first shell behavior in Figure 3.21. Back to Figure 3.23, the outer-shells contribution (ZZ) to the e–g difference seems to have more negative data points than positive ones; these negative points serves as the major contribution for the negative total e–g difference spectra in frequency region 30 cm−1 – 70 cm−1 . Thus, for 10% S mixture, the second shell contribution is almost zero and for 50% S mixture, the second shell resembles the first shell projection for frequency higher than 30 cm−1 and the outer shells determine the overall negative response for these high frequencies. Finally after all the projection analysis we can conclude that in 10% S solvent mixture the first shell determines the total e–g difference, whereas in 50% S solvent mixture, all the solvation shells participate in determining the overall e–g difference. Snapshot of the initial response From previous analysis we know that the first order DID and solute(u)/S-solvent(v) inter- action are important for the e–g difference. Thus, (uvuv) and (uvvv) terms defined below are the leading ingredients in the e–g difference spectrum. ˙ u j α j )t ⊗ (α u Tuk α k )t , j, k ∈ S (α u T (uvuv) terms (3.74) 1 2 ˙ u j α j )t ⊗ (α k Tkl α l )t , j, k, l ∈ S (α u T (uvvv) terms (3.75) 1 2 We also know that in the 10% S the first shell is the most visible region in the difference spectra, but in 50% S system both the first shell and the outer shells have significant contributions to the overall e–g difference. We now step back and ask what are those most important terms contributing to the instantaneous value of response at time 0, R(0), 126 0.0005 0.0003 ) ) / XX_10 (g) / XYYX_10 (g) 6 6 XX_10 (e) 0.0002 ( XYYX_10 (e) ( 0.0004 Anisotropic Spectral Density Anisotropic Spectral Density XX_50 (g) XYYX_50 (g) XX_50 (e) XYYX_50 (e) 0.0001 0.0003 X= first shell Y= second shell 0.0000 Z= else 0.0002 -0.0001 X= first shell 1st order DID Y= second shell 0.0001 #conf=5000 Z= else -0.0002 sample step=10 0.0000 -0.0003 0 50 100 150 200 0 50 100 150 200 -1 -1 c (cm ) c (cm ) (a) (b) 0.0005 0.0004 ) YY_10 (g) ) / XZZX_10 (g) 6 0.0003 / YY_10 (e) 6 XZZX_10 (e) ( 0.0004 ( YY_50 (g) Anisotropic Spectral Density 0.0002 XZZX_50 (g) Anisotropic Spectral Density YY_50 (e) XZZX_50 (e) 0.0003 X= first shell 0.0001 Y= second shell 0.0000 Z= else 0.0002 -0.0001 X= first shell 1st order DID Y= second shell -0.0002 #conf=5000 Z= else 0.0001 sample step=10 -0.0003 0.0000 -0.0004 0 50 100 150 200 0 50 100 150 200 -1 -1 c (cm ) c (cm ) (c) (d) 0.0025 ) ) ZZ_10 (g) 0.0008 / YZZY_10 (g) / 6 6 ZZ_10 (e) YZZY_10 (e) ( ( 0.0020 Anisotropic Spectral Density ZZ_50 (g) Anisotropic Spectral Density 0.0006 YZZY_50 (g) ZZ_50 (e) YZZY_50 (e) 0.0015 X= first shell 0.0004 X= first shell Y= second shell Y= second shell Z= else 0.0002 Z= else 0.0010 1st order DID 0.0000 0.0005 #conf=5000 sample step=10 -0.0002 0.0000 0 50 100 150 200 0 50 100 150 200 -1 -1 c (cm ) c (cm ) (e) (f) Figure 3.22 Projection of the OKE spectral densities for 10% S and 50% S solvent mixtures into first-shell (X), second-shell (Y) and outer-shells (Z) contributions. “Self” terms are shown on the left column and cross terms are shown on the right column. The spectral densities are computed with first-order DID polarizabilities and averaged over 5000 liquid configurations. 127 XX e/g difference spectrum 0.0003 ) / 6 0.0002 ( 10% S Anisotropic Spectral Density 50% S 0.0001 0.0000 -0.0001 X= first shell Y= second shell Z= else -0.0002 -0.0003 0 50 100 150 -1 c (cm ) (a) YY e/g difference spectrum ) 0.0003 / 6 ( 0.0002 10% S Anisotropic Spectral Density 50% S 0.0001 0.0000 -0.0001 X= first shell Y= second shell Z= else -0.0002 -0.0003 0 50 100 150 -1 c (cm ) (b) ZZ e/g difference spectrum 0.0010 ) / 6 ( 10% S Anisotropic Spectral Density 0.0005 50% S 0.0000 X= first shell Y= second shell -0.0005 Z= else -0.0010 0 50 100 150 -1 c (cm ) (c) Figure 3.23 Projection of the e–g difference in spectral densities for 10% S and 50% S solvent mixtures into first-shell (a), second-shell (b) and outer-shells (c) contributions. The e–g differences are computed with first-order DID polarizabilities and averaged over 5000 liquid configurations. 128 for a single phase space point in the first shell to the first order DID. d ˙ (0) ⊗ Π (0) R(0) = − β Π (t) ⊗ Π (0) = −β Π (3.76) dt t=0   β ˙ (0), Π (0)) − Tr(Π˙ (0)) · Tr(Π =− 3PP(Π Π(0)) . (3.77) 30 To the first order DID, the many-body polarizability and its time derivative are Π (t) = ∑ α i + ∑ α iTi j α j , (3.78) i i6= j ˙ (t) = Π ∑ α i·Ti j α j . (3.79) i6= j ˙ (0)⊗Π So all the contributing terms to kB T R(0) = Π Π(0) are (u is the solute, vi is solvent, and the sum of all particles in the first shell)   kB T R(0) = α u ∑ T˙ uv j α v j + ∑ α vi ∑ ˙ vi v j α v j T vj vi v j 6=vi   ⊗ ∑ α vi + α u ∑ Tuvk α vk + ∑ α vk ∑ Tvk vl α vl (3.80) vi vk vk vl 6=vk     ˙ uv j + ∑ ≡ ∑Π ∑ ˙ vi v j ⊗ Π ∑ αi + ∑ Πuvk + ∑ ∑ Πvk vl (3.81) vj vi v j 6=vi i vk vk vl 6=vk = ∑ ∑ Πuvuv ( j; k) j k + ∑ ∑ ∑ Πuvvv ( j; k, l) j k l6=k + ∑ ∑ ∑ Πvvuv (i, j; k) i j6=i k + ∑ ∑ ∑ ∑ Πvvvv (i, j; k, l), (3.82) i j6=i k l6=k 129 h i ˙ i j ⊗ ∑k αk 1 = αi α j ∑k αk 1 ˙ 1 ˙ where Π 10 PP(T, 1) − 30 Tr(T)Tr(1) = 0 (see Equation 3.88). For instance, a typical number of particles in the first shell is 12 , so (uvuv) has about 122 terms, (uvvv)/(vvuv) has 123 terms and (vvvv) has 124 terms. Listing all these terms in the indices order is tedious and pointless. Instead, sorting these terms based on their absolute value may give us some inspiration. So the largest 100 absolute values of each term tagged with their solvent indices for the four categories are sorted out. The largest 20 terms in each category for 10% S and 50% S on ground and excited states are listed in Appendix F. Table 3.5 is the summary of the dominant terms of each system. ground excited 10% S kB T R(0) = −4 × 10−3 kB T R(0) = −8.5 × 10−4 (uvuv) us j us j −9 × 10−4 dominant us j us j & usk usk −2 × 10−4 dominant (uvvv) us j wk s j 4 × 10−5 us j wk sl −5 × 10−5 (vvuv) s j wk us j 5 × 10 −5 s j wk usl −5 × 10−5 (vvvv) s j wk s j wk −7 × 10−6 s j wk s j sl −1 × 10−5 50% S kB T R(0) = 2.6 × 10−3 kB T R(0) = 2.8 × 10−3 (uvuv) us j us j −9 × 10 −4 us j us j 1.3 × 10−3 usk usk 6 × 10−4 (uvvv) us j sk sl −5 × 10−4 us j sk sl 6 × 10−4 (vvuv) s j sk us j 3 × 10 −4 s j sk usl 3 × 10−4 (vvvv) s j sk s j sk −2 × 10−4 s j sk s j sl −3 × 10−4 Table 3.5 Summary of single snapshot analysis of kB T R(0) and its largest term(s) (unit in σ 5 · (ε /m)1/2 ) of each category with “u” for the solute, “v” for a solvent. s and w stand for strong solvent and weak solvent with subscripts j, k, l to distinguish different solvent atoms. The same ground state equilibrium configuration is used in both 10% S and 50% S, and different excited equilibrium configurations are used in 10% S and 50% S. Details are included in Appendix F. We can see that kB T ∆R(0) of g → e transition are positive for both 10% S and 50% S, but the overall signs of kB T R(0) are negative for 10% S and positive for 50% S. The dominant term’s sign matches kB T R(0)’s sign (a dominant term should be about a order of magnitude larger than the second largest contribution). More obviously in 10% S system, the leading contributions to kB T R(0) are solute/single-S-solvent terms (us j us j ). First order 130 DID us j us j term can be expressed analytically:     αu α j T˙ u j ⊗ αu α j Tu j = αu2 α 2j T˙ u j ⊗ Tu j , ( j ∈ S). (3.83) Denoting the unit vector as rˆ = (ex , ey , ez ), we can express the dipole-dipole tensor as follows. ∂ ∂ 1 3rµ rν − r2 δµν Tµν (t) = = ∂ rµ ∂ rν (x2 + y2 + z2 )1/2 r5 1 = (3eµ eν − δµν ), (3.84) r3 where δµν is the Kronecker delta. The pair product and the trace of the dipole-dipole tensor are 1 Tr(T(t)) = ∑ Tµ µ (t) = r3 ∑(3e2µ − 1) = 0. (3.85) µ µ 1 2 PP(T(0), T(0)) = ∑ Tµν (0)Tµν (0) = r6 ∑ (9eµ eν eµ eν − 6eµ eν δµν + δµν ) µ ,ν µ ,ν 6 = . (3.86) r6 In addition, the time derivative of the dipole-dipole tensor is then d 3rµ rν − δµν  3 15 3 T˙µν (t) = 5 3 = 5 (˙rµ rν + rµ r˙ν ) − 6 rµ rν r˙ + 4 δµν r˙. (3.87) dt r r r r r So the trace of the derivative of dipole-dipole tensor is ˙ =PP(T, ˙ 1) = ∑ T˙µ µ = r˙ Tr(T) (6 − 15 + 9) = 0. (3.88) µ r4 131 The pair product of the derivative of dipole-dipole tensor is h i 3r r − r2 δ  ˙ 3 15 3 µ ν µν PP(T(0), T(0)) = ∑ (˙rµ rν + rµ r˙ν ) − 6 rµ rν r˙ + 4 δµν r˙ · (3.89) µ ,ν r5 r r r5  9 45 9 =∑ 10 (˙rµ rµ rν2 + r˙ν rν r2µ ) − 11 r2µ rν2 r˙ + 9 rµ2 r˙ µ ,ν r r r  3 15 2 3 2 − 8 (2˙rµ rµ ) + 9 rµ r˙ − 7 δµν r˙ r r r r˙ =(2 × 9 − 45 + 9 − 6 + 15 − 9) r7 r˙ = − 18 , (3.90) r7 dr 1 1 1 where r˙ = ∑ dt γ rγ2 = q · ∑ 2rγ r˙γ = ∑ rγ r˙γ . 2 ∑ r2 γ r γ (3.91) γ γ The instantaneous response function at time 0 is then   β ˙ 9β 2 2 r˙ R(0) = − PP(Π (0), Π (0)) = α α (3.92) 10 5 u S r7 r=solute/S-solvent distance The above equation can be rewritten in a more illuminating form.   d 3 2 21 R(0) = −β α α . (3.93) dt 10 u S r6 This equation highlights the most important contributions to the magnitude of the initial response: αu2 αS2 reflects the importance of solute/single-S-solvent pairs, and fast decaying 1/r6 emphasizes the contribution from the innermost solvation shell. 3.5 Discussions In this chapter, we have seen the long waiting time limit solute-pump/solvent-probe spectra for a preferential solvation model. In this limit, the solute-pump/solvent-probe spectra 132 is just the excited-state/ground-state OKE difference. The traditional OKE spectroscopy measures the entire liquid polarizability fluctuation, which cannot tell us about any local dynamical information. The OKE e–g difference spectra, on the other hand, emphasize the change of dynamics focusing on the nearby region around the solute. A striking feature of the e–g difference spectra is for 10% S system, the e/g difference is positive and flat; for 50% S system, the e–g difference is negative in low frequency region (ω /2π c < 75cm−1 ). Another interesting feature is above 75cm−1 , both systems have positive e–g difference, and 10% S is even larger than 50% S. To understand these, we turn off the solute’s polarizability (set as 0) and find that the e–g difference is almost zero; then we compare first order DID approximation with the exact infinite order DID approximation. The e–g difference remains even if the magnitude is smaller. So the e–g difference must primarily come from the solute–solvent interaction up to the first order DID. Projection operators are used to selectively isolate the influence of certain physical degrees of freedom on a physical variable. The physical variable in our case is the many- body polarizability. Applying projection operator to correlation functions leads to different contributions made up with “self” terms and cross terms. The partition of degrees of freedom we choose here is (1) first shell and else, (2) first shell, second shell and else. The correlation functions we calculated are polarizability-velocity correlation function, G(t) and the OKE response function, R(t). It turns out that all components using both partitions of G(t) or R(t) of 10% S has a notable change after a transition from g to e. In frequency domain, projected e–g difference spectra show that first shell contribution qualitatively agrees with full e–g difference for 10% S solvent mixture. As for 50% S solvent mixture. the role of outer shells seems to be more substantial. In order to have a closer look at the first shell polarizability interaction, we take a snapshot of first shell contributing terms to an instantaneous response R(0) in the first order DID approximation. In the analyzed 10% S configuration, the dominant term in kB T R(0) determines magnitude of kB T R(0). The dominant term in 10% S system, if any, tends to be 133 a (us j us j ) term, or instantaneous dipole-induced dipole interaction between the solute and the same strong solvent. The 50% S system is more likely to have multiple important terms or a dominant term involving the solute and two distinct solvents. In brief, the solute/1st-shell-S-solvent term to the first order is essential to the OKE e–g difference spectra, especially for the 10% S mixture. The big question that still puzzles us is why the signs of e–g difference are different in 10% S and 50% mixtures. To interpret the sign of e–g difference spectra in Figure 3.11, one have to reconsider exactly how many S solvents locate in the immediate vicinity of the solute in both solvent mixtures. According to the first shell solvent populations in 10% S and 50% S (Table 3.2), the average number of S solvent increases drastically from 1.3 S to 3.0 S in the 10% S mixture, and the average number of S solvent increases moderately from 6.2 S to 9.4 S in the 50% mixture. It is easy to understand that in 10% S, more than 100 percent rise in the population of bright S solvent leads to a gain in the spectra. In the 50% S, the number of S solvent only increases by 50%, not to mention that even on the ground state, the S solvent populates more than a half if the first shell. The relative population in the first shell is a key quantity. A typical total number in the first shell is 12–13 solvents, so for the 10% S mixture, the spectroscopic bright S solvent is a minority with the solute on whatever electronic state. But for 50% S mixture, the S solvent is about 3 times more population than the S solvent in the 10% S mixture, and is always a majority, which indicates that the first shell is crowded with spectroscopically bright atoms. Since we know the solute/S-solvent terms dominate the OKE response function, we have a general expression as follows * + ˙ hΠ(0)Π(t)i = ∑ π˙ j (0)πk (t) j,k * + ≈4αu2 αS2 ˙ 0 j (0)) ∑ T(r0k (t)) ∑ T(r . (3.94) j∈S k∈S In the dilute S solvent case (10% S), the major contribution for the many-body polarizability 134 probably comes from two-body terms involving the solute and a single S solvent ( j = k terms); yet in a dense S solvent case (50% S), the major contribution for the response could come from three-body terms ( j 6= k terms). A schematic illustration in Figure 3.24 helps us to think about the signs of the responses. It depicts typical arrangements of instantaneous dipoles in 10% S and 50% S mixtures. Although there is no permanent dipole in either the solute or the solvent in our model, instantaneous dipoles can exist in any atom. In the dilute S solvent case (10% S), the instantaneous dipole on an S solvent (µ j ) always tend to align with the instantaneous dipole on the solute (µ 0 ) in order to minimize the energy. By contrast, in the dense S solvent case (50% S), with a second S solvent, the minimum-energy arrangement can lead to instantaneous dipoles on two solvents pointing to opposite directions such that µ j · µ k < 0. As a result of this partial cancellation, the many-body polarizability could become lower when the S solvent number increases in the 50% S case. Consequently, the e–g differences have opposite signs. 10% S 50% S Figure 3.24 Illustration of how interaction-induced polarizabilities are affected by number of spectroscopic bright solvent, leading to sign changes in the solute-pump/solvent-probe spectra. When an instantaneous dipole is created in the solute (black circle), instantaneous dipoles induced in the most polarizable solvent (denotes as S, red circles) tend to have an arrangement that leads to a minimum energy. In the case of dilute S solvent (10% S), the dipole–dipole interaction results in aligned instantaneous dipoles on the solute and on the single S solvent, which results in a reinforced many-body polarizability. By contrast, in the case of dense S solvent (50% S), three-body solute-solvent-solvent triplets begin to be more important, inducing a diminished many-body polarizability since the unaligned instantaneous dipoles on solvents are partially canceled. 135 The following simple calculation is illustrative. As the average of the OKE response at time 0 vanishes, hR(0)i = 0 since h˙ri = 0, we consider the integrated OKE correlation function C(0) in the first order DID approximation. Again, the biggest contribution comes from the solute/S-solvent terms, so the integrated OKE correlation function is approximately * + C(0) = hΠxz(0)Πxz (0)i ≈ 4αu2 αS2 ∑ Txz (r0 j )Txz (r0k ) (3.95) j,k∈S * +  xz   xz  = 36αu2 αS2 ∑ (3.96) j,k∈S r5 r0 j r5 r0k * + 36αu2 αS2 = ∑ Xˆ0 j Zˆ0 j Xˆ0k Zˆ0k . (3.97) r6 j,k∈S In the 10% S case, consider only the solute and a single strong solvent j with distance r, we apply rotational averaging (using the notation defined in Appendix C) by integrating over all possible θ , φ for the solute–solvent displacement vector rˆ 0 j = (Xˆ0 j , Yˆ0 j , Zˆ0 j )T = (sin θ cos φ , sin θ sin φ , cos θ )T . (3.98) The rotational average in this case is 2 Xˆ0 j Zˆ0 j = (sin θ cos φ · cos θ )2 Z 2π Z π 1 1 = dφ · dθ sin θ (sin θ cos φ · cos θ )2 = . (3.99) 4π 0 0 15 Therefore,   36αu2 αS2 1 C(0) = > 0. (3.100) r6 15 In the 50% S case, consider there is a solute-(solvent j)-(solvent k) equilateral triangle 136 in the first shell. We fix the angle between rˆ 0 j and rˆ 0k as 60◦ and perform the rotational averaging. (Figure 3.25) Z z r0j r0j r0k 60o 60o X x O r0k O Figure 3.25 Rotational averaging of a solute–solvent–solvent equilateral triangle. Still, rˆ 0 j has two dimensional degrees of freedom, θ , φ . We can choose a fix point √ rˆ 0k = (xˆ0k , yˆ0k , zˆ0k )T = ( 3/2, 0, 1/2)T that has a 60◦ angle with respect to z-axis in the molecule-frame xyz, then transform it into laboratory frame XYZ,      ˆ  X0k   cφ cθ cχ − sφ sχ −cφ cθ sχ − sφ cχ cφ sθ   xˆ0k        Yˆ  =  sφ cθ cχ + cφ sχ −sφ cθ sχ + cφ cχ sφ sθ   yˆ . (3.101)  0k     0k       Zˆ 0k −sθ cχ sθ s χ cθ zˆ0k rˆ 0k has three degrees of freedom: θ , φ associated with solvent j and χ associated with the rotation along the rˆ 0 j direction with a fixed 60◦ angle. The rotational average in this case 137 is "√ 3 Xˆ0 j Zˆ0 j Xˆ0k Zˆ0k =sin θ cos φ · cos θ · (cos φ cos θ cos χ − sin φ sin χ ) 2  " √ # 1 3 1 + cos φ sin θ · − sin θ cos χ + cos θ 2 2 2 1 =− . (3.102) 120 Therefore,   36αu2 αS2 1 C(0) = − < 0. (3.103) r6 120 Unlike the dilute 10% S system, in which the instantaneous dipole gets reinforced by just the solute/S-solvent pair, the instantaneous dipole in the more concentrated 50% S solution gets frustrated by the non-collinear solute/S-solvent/S-solvent arrangement within the first shell, resulting in a negative e-g difference distinct from that of the more dilute 10% S case. † 2 † Note: the polarization choice of YZYZ is equivalent to XZXZ by symmetry, checking Yˆ0 j Zˆ 0 j = 1/15, and Yˆ0 j Zˆ0 jYˆ0k Zˆ0k = −1/120. 138 Chapter 4 Nonequilibrium Two-Dimensional Solute-Pump/Solvent-Probe Spectroscopy 4.1 Introduction In the last chapter, we have seen the long waiting time limit (T → ∞) of the solute- pump/solvent-probe spectra in our model atomic liquid mixtures; the response degenerates into the difference between the optical Kerr effect (OKE) spectra for a solution with an electronically excited solute and that with a ground-state solute. Even at this large T limit, the solute-pump/solvent-probe spectroscopy proves its capability at unveiling microscopic features of preferential solvation dynamics. The sign of the difference changes when the mole fraction of the solvent components is adjusted. Perhaps coincidentally, the sign change occurs in experimental long-T measurements as well, when different pure dipolar solvents are examined, displaying a sensitivity of the spectra to the specifics of the solvent.226 Figure 1.12 shows the first experimentally measured full two-dimensional (2D) solute-pump/solvent-probe spectra scanning both solvation time T and ultrafast dy- 139 namical time t, the resonant-pump polarizability response spectra (RP-PORS).227 This RP- PORS experiment is designed to detect the evolution of the solvent dynamics as structural relaxation proceeds during solvation. However, the 2D spectra still lack a molecular interpretation. Finding this interpretation is the primary aim of this chapter. Most commonly used experimental avenues for probing solvation dynamics have been relying on the solute’s energy gap as the indicator of the solvent’s dynamics. For example, as discussed in chapter 1, traditional time-dependent fluorescence spectroscopy tells us about the relaxation of solvents around a newly excited dye solute by measuring the solute’s red-shifting fluorescence emission frequency, which is equivalent to measuring the energy gap between the excited-state-solute and ground-state-solute potentials.2, 83, 246 Moreover, alternative photon-echo route into solvation dynamics such as 3-pulse photon echo peak shift (3PEPS) tracks the same solute–solvent interaction energy following an electronic excitation of the solute.141, 172, 176, 177, 282 How much does the evolution of the energy gap disclose the full dynamics of solvent structural rearrangement? Since different solvent structures may share common solute– solvent interaction energies, the time scale for solvent structural relaxation may or may not be the same as that for the energetic relaxation. For example, consider the situation that the solute-solvent interaction energy is dominated by solvents in the first shell, yet the solvent molecules in the outer shells needed more time to reorganize even after the first- shell equilibrium solute-solvent organization is achieved. Any energy-gap spectroscopy would lead to a potential-energy relaxation that is faster than the actual full structural relaxation. On the contrary, if considerable rearrangement in the outer shells were needed before the first shell starts its own repositioning, the energetic spectroscopy would display a potential-energy relaxation that is slower than most of the structural relaxation. The solute-pump/solvent-probe spectroscopy, on the other hand, grants a possibility of directly monitoring the solvent structure following a resonant solute excitation, and thus could serve as a 2D solvation spectroscopy. This spectroscopy is essentially a nonequi- 140 librium four-wave-mixing light scattering measurement. A resonant electronic excitation in the solute triggers the solvation, and the four-wave-mixing spectra are measured at a time T later. The nonresonant four-wave-mixing light scattering employs the many-body polarizability as the observable that is sensitive to the intermolecular distances between all the spectroscopically bright species, not just the solute-solvent distances, and thus may be more responsive to solvent structure than energy-gap spectroscopies. As briefly mentioned in chapter 2, the practical difficulty of calculating the full solute- pump/solvent-probe spectra lies in the evaluation of the classical Poisson bracket. In the following sections of this chapter, a hybrid method using a combination of molecular dynamics and instantaneous-normal-mode theory to evaluate the 2D response will be described. The full 2D spectra for our previously studied preferential solvation model will be shown, as well as a molecular analysis on what information can be extracted from the them. Finally we offer some general discussions on the utility of the hybrid method to two-dimensional spectroscopic responses. 4.2 Instantaneous-Normal-Mode (INM) Theory The idea behind instantaneous normal modes is stated in chapter 1. The INM theory works well in predicting short-time dynamics in liquids and also can provide a practical approach to evaluating the classical Poisson bracket,20 which is a computational challenge in our calculation of the solute-pump/solvent-probe spectra. The independent harmonic nature of the INMs makes the Poisson bracket easy to evaluate. In the current section, we will familiarize ourselves with the notations and then discuss the connection of INM theory to the OKE spectra through the INM influence spectra. The INM predictions on OKE spectra for the previously studied preferential solvation model will be demonstrated to have a good agreement with the exact MD result. 141 4.2.1 Mass-weighted INM The instantaneous-normal-mode theory treats the microscopic short-time liquid9 and clus- ter6 dynamics as a set of independent simple harmonic motions. For example, for an atomic liquid with N atoms, the classical Hamiltonian can be written as N 1 ˙ = H(R, R) ∑ 2 m j r˙2j +V (R), (4.1) j=1 where R is the complete set of 3N coordinates R = {r1 , . . . , rN }. The coordinate dimen- sionality is 5N (including 2N rotational coordinates) for rigid linear molecules16 or 6N (including 3N rotational coordinates) for rigid nonlinear molecules.20, 283 In the mass- weighted approach of INM theory that provides a general description on systems with multiple atom types, we define the mass-weighted coordinates as √ Z = {z j µ ; j = 1, . . ., N; µ = x, y, z}, z jµ = m j r jµ . (4.2) If denote diagonal mass matrix M = diag{m1 , m1 , m1 , m2 , m2 , m2 , . . . , mN , mN , mN }, in the case of pure atomic liquid, M = m1, with 1 the 3N × 3N identity matrix and m the atom mass. Then the mass-weighted coordinates can be written as Z = M1/2 · R. (4.3) Assuming that the time t is short enough, the Hamiltonian at time t can be expand about initial configuration R0 ,7 1˙ ˙ 1 H≈ Zt · Zt +V (Z0 ) − F(R0 ) · (Zt − Z0 ) + (Zt − Z0 ) · D(R0 ) · (Zt − Z0 ), (4.4) 2 2 where except for the first kinetic energy term, the right-hand side is the potential energy. 142 The instantaneous force vector is defined as ∂ V 1 ∂ V F j µ (R0 ) = − =− √ , (4.5) ∂ z j µ R0 m j ∂ r j µ R 0 and the dynamical (instantaneous Hessian) matrix is defined as ∂ 2V 1 ∂ 2V D j µ ,kν (R0 ) = = √ . (4.6) ∂ z j µ ∂ zkν R0 m j mk ∂ r j µ ∂ rkν R 0 If the orthogonal matrix U(R0 ) can diagonalize the real symmetric dynamical matrix,    ω12     ω22  UDUT =    , (4.7)  ..   .    2 ω3N where ωα are frequencies of the α -th INM qα (t; R0). The mutually independent INMs are the eigenvectors as follows qα (t; R0) =[U(R0) · (Zt − Z0 )]α , (α = 1, . . . , 3N), (4.8) whose conjugate momenta are ˙ t − Z˙ 0 )]α , pα (t; R0) = [U(R0 ) · (Z (α = 1, . . ., 3N). (4.9) The INMs compose the optimum basis for describing the short-time dynamics, and the 143 transformation matrix to the Cartesian coordinate is ∂ qα 1 ∂ qα Uα , j µ (R0 ) = =√ , (4.10) ∂ z jµ m j ∂ r jµ ∂ zk ν √ ∂ rk ν UTkν ,β (R0 ) = = mk . (4.11) ∂ qβ ∂ qβ The instantaneous forces associated with each mode α are fα =[U(R0 ) · F(R0)]α . (4.12) Now the potential energy can be written in the INM basis as  3N  1 2 2 V (Rt ) ≈V (R0 ) + ∑ − fα qα (t) + ωα qα (t) , (4.13) α =1 2 fα vα (0) where qα (t) = 2 (1 − cos ωα t) + sin ωα t, (4.14) ωα ωα fα vα (t) =q˙α (t) = vα (0) cos ωα t + sin ωα t. (4.15) ωα So we have the dynamics of a liquid that is governed by INMs at short times, 1 r j µ (t) =r j µ (0) + √ ∑ UTjµ ,α qα (t). (4.16) mj α In other words, one would only need to know the probability distribution of the set of fα , ωα , and vα (0) so as to know the complete dynamics of a liquid in a short time scale. The instantaneous normal mode spectrum is the density of states of a liquid (or the phonon spectrum), i.e. the distribution of instantaneous normal mode frequencies averaged over liquid configurations, * + 1 3N 3N α∑ D(ω ) = δ (ω − ωα ) . (4.17) =1 144 The density of states is normalized such that the integral Z dω D(ω ) = 1, (4.18) and for us the unit of D(ω ) is 1/[ω ] = 1/cm−1 = cm. 4.2.2 INM density of states in liquid argon In chapter 1, the INM spectrum of liquid argon is shown (Figure 1.1). Here we show some related details starting from looking closer at the structure of the dynamical matrix D. If our system has N atoms with identical mass m, interactions between which are described by pair potential u(r), then the 3×3 dynamical matrix for a pair of atoms i and j is defined as Dµ ,ν (i, j) = Diµ , jν , and in tensor form,5   1 1  ∑k6=i t(rik ), i = j D(i, j) = ∇i ∇ j V (R(0)) = · , (4.19) m m  −t(ri j ), i = 6 j t(r) =[u′ (r)/r]1 + [u′′(r) − u′ (r)/r]ˆr rˆ . (4.20) In practice, liquid argon is modeled by molecular dynamics simulation in which there are 256 argon atoms under thermal conditions of temperature kB T /ε = 1.00 and density ρσ 3 = 0.80. The potential energy is the sum of pair potentials between each pair of atoms, and the pair potential used here is the Lennard-Jones potential, h i V (R) = ∑ u(ri j ) = ∑ 4εi j (σ /ri j )12 − (σ /ri j )6 . (4.21) i< j i< j 145 The dynamical matrix is then 1 ∂ 2V (t) [D(t)]iµ , jν = (4.22) m ∂ riµ (t)∂ r jν (t) " ! 1 σ 12 σ6 =δi j ∑ 4εik 168 · 16 − 48 · 10 (riµ − rkµ )(riν − rkν ) m k6=i rik rik ! # σ 12 σ6 + −12 · 14 + 6 · 8 δµν rik rik " ! 1 σ 12 σ6 − (1 − δi j ) 4εi j 168 · 16 − 48 · 10 (riµ − r j µ )(riν − r jν ) m ri j ri j ! # σ 12 σ6 + −12 · 14 + 6 · 8 δµν . (4.23) ri j ri j Periodic boundary condition and minimum image convention are utilized in calculating forces and dynamical matrices.284 There are two fundamentally equivalent ways to calculate the INM spectrum; the first one is to bin the INM frequencies into a histogram (the discrete method to evaluate the average of a delta function) and the other way is applying the same approach to the dynamical matrix eigenvalues and then convert the distribution of eigenvalues into that of frequencies. We denote the eigenvalues of dynamical matrix as λα = ωα2 . From the normalization condition of the distribution of eigenvalues, we can derive the relation between the frequency distribution and the eigenvalue distribution: Z 1= d λ ρ (λ ) Z = d ω 2ω · ρ ( λ = ω 2 ) Z = d ω D(ω ), ∴ D(ω ) =2ω · ρ (λ = ω 2 ). (4.24) 146 The INM spectrum using the frequency histogram (200 bins) is shown in Figure 4.1, the comparison with zero frequency removed INM spectrum is shown in Figure 4.2 and the distribution of eigenvalues is shown in Figure 4.3. The comparison of frequency-binning and eigenvalue-binning is shown in Figure 4.4. 0.020 Liquid argon k T/ =1.00 B 0.015 3 =0.80 ) (cm) 0.010 D( 0.005 0.000 -50 0 50 100 150 -1 /2 c (cm ) Figure 4.1 Instantaneous-normal-mode spectrum (density of states) of liquid Ar at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80 averaged 20,000 configurations. The imaginary frequencies are plotted on the negative frequency axis. Note that there is a small peak at zero frequency, corresponding to the three translational modes. Notice that in simulation, the LJ reduced unit is used with the argon reference tempera- ture ε /kB =119.8 K, length σ =3.405 A˚ and mass mAr =40.0 amu.284 Any quantity using the reduced unit is indicated by a star superscript. m∗ =m/mAr = 1, (4.25) ω ∗ =ωτLJ , (4.26) ω ∗2 λ λ∗ = 2 = . (4.27) 1/τLJ ε /(mAr σ 2 ) 147 0.020 zero freq removed total 0.015 ) (cm) 0.010 D( 0.005 0.000 -50 0 50 100 150 -1 /2 c (cm ) Figure 4.2 Instantaneous-normal-mode spectrum (density of states) of liquid Ar with (black, solid) and without zero frequencies removed (red, dashed). The simulation is performed at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80 and 20,000 configurations are averaged. The imaginary frequencies are plotted on the negative frequency axis. When the three zero-frequency translational modes are not counted in histogram, the small peak at zero frequency vanishes leading to a continuous smooth INM spectrum. 148 0.006 0.005 0.004 m 0.003 0.002 0.001 0.000 0 1000 2000 3000 Figure 4.3 Distribution of eigenvalues of the dynamical matrix for liquid Ar at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80 averaged 20,000 configurations. Eigenvalue is in Lennard-Jones reduced units. Finally we usually convert frequency unit to wavenumbers (cm−1 ) as is commonly done in spectroscopies and the conversion rule is as follows ω ω∗ 3.33565 × 10−11(cm−1 /Hz) ν˜ (cm−1 ) = = = ω∗ · −12 = ω ∗ × 2.4578(cm−1 ). 2π c 2π cτLJ 2π · 2.1560 × 10 (s) (4.28) Check simulations The matrix diagonalization is performed by a Fortran library LAPACK.285 To make sure the simulation is correct and the INM indexing is right, one has to check some properties of the INMs. Those checks are listed below. (1) The sum of elements of any row or column of the dynamical matrix is zero. ∑k D( j, k) = ∑k D(k, j) = 0. I checked sum of one row is zero (< 10−13 ). (2) Three zero eigenvalues corresponding to three zero-frequency translational modes. 149 0.020 frequency statistics eigenvalue statistics 0.015 Liquid argon ) (cm) k T/ =1.00 B 3 0.010 =0.80 D( 0.005 0.000 -50 0 50 100 150 -1 /2 c (cm ) Figure 4.4 Comparison of the instantaneous-normal-mode spectra of liquid Ar using frequency histogram and eigenvalue histogram, at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80 averaged 20,000 configurations. Both spectra are conceptually equivalent. Although eigenvalue approach gives a zero density of states at zero frequency, it has less data points at the low frequency region than the high frequency region, not like the frequency histogram approach that bins the frequency evenly. 150 Checked one configuration and the 9 eigenvalues around zero are as follows -0.804426857204 -0.363150031495 -0.0872967319105 -1.993958527e-16 1.14488423151e-13 1.38218188184e-13 0.0440473625128 0.349046510601 0.493732811159 (3) The averaged sum of diagonal elements of the dynamical matrix is equal to the averaged sum of all eigenvalues of the dynamical matrix. The simulation gives two identical sums (normalized, divided by 3N), 279.505020419 (ε /mσ 2 ). (4) Einstein frequency, ωE . The square of Einstein frequency is equal to the average of the diagonal element of dynamical matrix.6 * + * + N Z 1 1 1 ωE2 ∇2j u(r jk ) drg(r)∇2 u(r). 3N ∑ ∑ = D jµ , jµ = = ρ (4.29) jµ 3Nm j=1 3m From the result of check (3), we know that average of the diagonal element is 279.505 ε /mσ 2 . And we can use the formula at the last of the above equation where g(r) is the radial distribution function and also shown in Figure 3.1. In evaluating g(r), 700 bins are used and I assume g(r) is 1 for r is larger than the half of the simulation box. Then, we can calculate the integral as follows Z ∞ 1 ωE2 = ρ 4π dr · r2 g(r)∇2u(r) = 270.740 (ε /mσ 2 ). (4.30) 3m 0 These results are of about 3% difference, and can be considered the same. One could 151 increase the number of bins for g(r) and average more configurations to have a better sampling around 0.9–1.3 σ region where the integrand has a significant value. Another origin of error comes from the assumption for g(r)=1 for large r values. (5) Check the order of eigenvectors by performing matrix multiplication of UDUT (= Λ ). The first 3 × 3 submatrix of UDUT is -125.052838464 -1.01167119861e-13 8.58009376167e-14 -6.5685077413e-14 -116.405159786 3.22481841578e-14 8.81489552982e-14 2.72363749121e-14 -113.388715178 where the diagonal elements are the first three eigenvalues (first three diagonal elements of Λ ). Figure 4.5 summarizes the density of states for liquid argon and atomic solvent mixtures (10% S and 50% S) and under thermodynamic condition of temperature kB T /ε = 1.00 and density ρσ 3 = 0.80. They have the identical density of states distributions. The real part is dominant in area (real:imaginary=2.9:1) but the imaginary part has higher peak height. Part of the reason of there are imaginary modes in an atomic liquid is that the translational motion has a diffusive portion that could contribute to the imaginary modes and in our atomic liquid model there is no orientational motion (especially librations) that usually has a relatively large contribution to real-frequency INMs. The imaginary frequencies correspond to the unstable INMs meaning the negative curvature of the potential energy surface5 (curvature is proportional to the square of the frequency) such that the motion will not posses any oscillatory character. We will focus on stable real modes when discuss the influence spectrum later in this chapter. The non-diffusive part of the imaginary modes is not completely known and how imaginary modes contribute to a spectral observable is still an open question. But we are sure that some local negative curvature do contribute to a well-defined global oscillation as illustrated by the Morse oscillator with a energy higher than the energy at the curvature turning point x0 at which u′′ (x0 ) = 0. 152 0.020 10% S, excited 10% S, ground 0.015 Argon ) (cm) 0.010 D( 0.005 0.000 -50 0 50 100 150 200 -1 /2 c (cm ) (a) 0.020 10% S, ground 10% S, excited 0.015 50% S, ground ) (cm) 50% S, excited 0.010 D( 0.005 0.000 -50 0 50 100 150 200 -1 /2 c (cm ) (b) Figure 4.5 Instantaneous-normal-mode spectra (density of states) of (a) liquid argon, 10% S on the ground and the excited state, and (b) 10% S vs 50% S on the ground and the excited state. The molecular dynamics simulations are performed at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80 and averaged 20,000 configurations (sampling 1 configuration every 10 time steps, and time step is δ t = 0.0025τLJ ). The imaginary frequencies are plotted on the negative frequency axis. All density of states curves overlap on this scale. 153 4.2.3 INM influence spectrum The INM influence spectrum17, 18 connects the INMs with the correlation function of the physical quantity we are interested in. For any dynamical variable A (or called generalized force), the Taylor expansion at some initial configuration in powers of the INMs is 3N     ∂A 1 ∂ 2A A(t) = A(0) + ∑ ∂ qα qα (t) + ∑ 2 α ,β ∂ qα ∂ qβ qα (t)qβ (t) + · · · . (4.31) α =1 t=0 t=0 In the linear INM theory,81 neglecting terms after the first order in q, we have the time derivative of A as follows     ˙ =∑ ∂A dv j µ ∂A A(t) =∑ q˙α (t) = ∑ cα q˙α (t), (4.32) jµ ∂ r jµ t=0 dt α ∂ qα t=0 α where the coupling coefficient is   ∂A ∂ A ∂ r j µ cα = =∑ · . (4.33) ∂ qα t=0 jµ ∂ r j µ ∂ qα t=0 The generalized force velocity autocorrelation function17, 81 is then * +   ˙ A(0)i GAA (t) =hA(t) ˙ = ∑ cα q˙α (0)cβ q˙β (t) = ∑ cα2 q˙α (0)q˙α (t) α ,β α   Z   =kB T ∑ cα2 cos ωα t = kB T dω cos ω t ∑ cα2 δ (ω − ωα ) α α Z =kB T dω cos ω t ρA (ω ). (4.34) As any autocorrelation function is even, the cosine transformation is equivalent to R the Fourier transformation, i.e. GAA (t) = kB T dω eiω t ρA (ω ). According to Wiener- ˙ Khintchine theorem,12 ρA (ω ) is the spectral density (or power spectrum) of A(t). † Now † The R property of delta function dx f (x)δ (x − a) = f (a) has been used in above derivation as well as the change of order of integration with respect to ω and the phase space h·i and the equipartition theorem18 154 we define the influence spectrum ρA (ω ),   ρA (ω ) = ∑ c2α δ (ω − ωα ) , (4.35) α which is a weighted INM density of states, and the weighting factor is square of the derivative of the generalized force with respect to the INM (the coupling coefficient cα ). There are two normalization methods for the influence spectrum. (1) the area of the influence spectrum is normalized to the sum of cα2 with all frequencies; (2) the area of the influence spectrum for real frequencies is normalized to the sum of cα2 with only real frequencies. Z * + (1) dωρA (ω ) = ∑ cα2 , (4.36) all all ωα Z * + (2) dωρA (ω ) = ∑ cα2 . (4.37) real real ωα The above formulas are these two methods correspondingly and we use the all-frequency normalization only when the comparison of real and imaginary of the influence spectra is needed (Figure 4.6) and we use the real-frequency normalization for OKE related spectra (Figure 4.7, 4.9) and our solute-pump/solvent-probe spectra. It is worth noting that both methods give almost the same normalization factors, whose relative difference is only about 2 × 10−5 . The generalized force autocorrelation function CAA (t) can be calculated from the velocity autocorrelation function GAA (t) defined by Equation 4.34,286 using the following hq˙α q˙β i = kB T δαβ . 155 relation. hδ A(t)δ A(0)i CAA (t) = , (4.38) hδ A(0)δ A(0)i d 2CAA (t) GAA (t) 2 =− , (4.39) dt hδ A(0)δ A(0)i Z kB T CAA (t) = 1 − d ωρA (ω )(1 − cos ω t)/ω 2 . (4.40) h(δ A)2 i 4.2.4 INM approach to OKE spectra The OKE response function is the time correlation function of the many-body polarizabil- ity, ˙ xz (0)Πxz(t)i, Rxzxz (t) = h{Πxz (0), Πxz(t)}i = β hΠ Z ∞ d ˙ xz (0)Π ˙ xz(t)i = Rxzxz (t) = β hΠ dωρxzxz (ω ) cos ω t, (4.41) dt 0 by analogy to Equation 4.34. So in OKE spectroscopy, the dynamical variable A is Πxz and the polarizability influence spectrum can be expressed as ∂ Πxz ρxzxz (ω ) = (Πxz,α )2 δ (ω − ωα ) , where Πxz,α ≡ . (4.42) ∂ qα 1 In practice, the rotational invariance in isotropic system is applied, hAµν Bµν irot = 10 PP(A, B)− 1 30 Tr(A)Tr(B), thus we have * 2 +     ∂ Πxz 1 ∂Π ∂Π 1 ∂Π 2 = PP , − Tr . (4.43) ∂ qα 10 ∂ qα ∂ qα 30 ∂ qα rot Figure 4.6 shows the polarizability influence spectra for 10% S and 50% S systems. Atomic polarizabilities are chosen as in the previous chapter, αU =0.2σ 3 , αS =0.101σ 3 , and αW =0.0186σ 3. The first order DID approximated many-body polarizabilities are 156 0.0008 10% S, excited 10% S, ground 0.0006 -1 m cm 0.0004 4 0.0002 0.0000 -50 0 50 100 150 200 -1 /2 c (cm ) (a) 0.007 50% S, excited 0.006 50% S, ground -1 0.005 m cm 0.004 0.003 4 0.002 0.001 0.000 -50 0 50 100 150 200 -1 /2 c (cm ) (b) Figure 4.6 Instantaneous-normal-mode polarizability influence spectra of (a) 10% S and (b) 50% S on the ground state (black curve) and on the excited state (red curve) at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80 and averaged 100,000 configurations (sampling 1 configuration every 10 time steps, and time step is δ t = 0.0025τLJ ). The imaginary frequencies are plotted on the negative frequency axis. Normalization of the area of the influence spectrum ρ (ω ) of all frequencies to the sum of c2α with all frequencies is utilized. 157 used. From this figure, we observe that the shape of influence spectra are similar and the magnitude of excited state is higher than ground state for 10% S and the magnitude of excited state is lower than ground state for 50% S which is consistent with our previous MD results in chapter 3. However, the polarizability influence spectrum is not the OKE spectrum itself. The INM prediction for the OKE spectrum is called the OKE influence spectrum that is shown in the following equation. Note that the OKE influence spectrum is called the “modified” influence spectrum in literature.20 Z ∞ π ρxzxz (ω ) Im[Rxzxz (ω )] = dtRxzxz (t) sin ω t = . (4.44) 0 2 ω To derive the above relation, we start by express the Poisson bracket in the INM basis. ∂ Πxz (0) ∂ Πxz (t) {Πxz (0), Πxz(t)} = ∑ · (4.45) jµ ∂ r j µ (0) ∂ p j µ (0) ∂ Πxz ∂ qα (0) ∂ Πxz ∂ qβ (t) =∑ ∑ · (4.46) j µ α ,β ∂ qα 0 ∂ r j µ (0) ∂ qβ t ∂ p j µ (0) ∂ Πxz ∂ Πxz =∑ {qα (0), qβ (t)} . (4.47) α ,β ∂ qα 0 ∂ qβ t The INM fundamental Poisson bracket198 below can be derived from Equation 4.14 and 4.15, ∂ qα (t1 ) ∂ qβ (t2 ) χ α ,β (t1 ,t2) ≡{qα (t1 ), qβ (t2 )} = ∑ jµ ∂ r j µ (t1 ) ∂ p j µ (t1 ) sin ωα (t2 − t1 ) = δαβ = χα (t2 − t1 )δαβ , (4.48) ωα Z t α sin ωα t where χα (t) = dτ Cvv (τ ) = , (4.49) 0 ωα α hvα (0)vα (τ )i Cvv (τ ) = = cos ωα τ . (4.50) hv2α i 158 We then have the OKE response function in INM basis: * + ∂ Πxz ∂ Πxz R(t) = h{Πxz(0), Πxz(t)}i = {qα (0), qβ (t)} , (4.51) ∂ qα 0 ∂ qβ t In the linear INM theory, we assume the coupling coefficient is constant in such a short time, ∂ Πxz ∂ Πxz ≈ = Πxz,α . (4.52) ∂ qα t ∂ qα 0 So the OKE response function becomes   2 RINM (t) = ∑(Πxz,α ) χα (t) . (4.53) α The OKE spectral density is the imaginary part of the Fourier transform (sine transform) of the response function, Z ∞ Im[RINM (ω )] = dt R(t) sin ω t 0 Z ∞   2 = 0 dt ∑(Πxz,α ) χα (t) sin ω t α  Z ∞  2 sin ωα t = ∑(Πxz,α ) 0 dt ωα sin ω t . α Using following identity Z ∞ π sin ω t sin ω ′t dt = [δ (ω − ω ′ ) − δ (ω + ω ′ )]. (4.54) 0 2 To prove above relation, one needs to use the following property of delta function Z ∞ Z ∞ 1 1 δ (ω ) = dt cos ω t = dte±iω t . (4.55) 2π −∞ 2π −∞ 159 Thus, the OKE spectrum using INM approximation is given by   π 1 2 Im[RINM (ω )] = 2 ∑(Πxz,α ) ωα [δ (ω − ωα ) − δ (ω + ωα )] (4.56) α π ρxzxz (ω ) = . (4.57) 2 ω Check simulations The sum of the squares of spatial derivative of the polarizability are equivalent using either the INM qα basis or Cartesian coordinates r j µ basis as a result of the transformation matrix being unitary:  2  2 ∂ Πxz ∂ Πxz ∑ qα =∑ r jµ . (4.58) α jµ For one configuration of 10%S system, the sum is 0.056834 using both methods, and the sum is 0.049215 if only real modes are included in the INM-basis method which is about 13% less than the full INM result. OKE spectra: comparison between INM and MD methods Figure 4.7 shows the OKE influence spectra for 10% S and 50% S systems and the relative peak heights of the excited vs. ground is in the same order with the polarizability influence that is also consistent with previous results. Comparing with the molecular-dynamics results for the OKE spectra (Figure 4.8), the INM approximation works well for frequency larger than 25 cm−1 in 10% S system with ground- and excited-state solute. But what we are most interested in is the difference of excited-state and ground-state OKE influence spectra, so we show the comparisons with the previous molecular-dynamics OKE spectra in Figure 4.9 and 4.10. ‡ Figure 4.9 and 4.10 are the major findings of this section. They indicate a good ‡ The OKE influence spectrum π2 ρ (ωω ) has a unit of σ 4 /(m ∗ cm−2 ) which can be transformed to the LJ −1 p reduced unit σ 6 /ε by multiplying 2.45782 since 1τLJ = 2.4578cm−1 and τLJ = mσ 2 /ε . 160 agreement between INM and MD methods for predicting OKE e-g difference spectra in the high-frequency region. First of all, the signs of e–g difference for both 10% S and 50% S in Figure 4.9 are the same as the MD results. Furthermore, the magnitudes of the difference spectra using the INM approach is quite similar with that of the exact MD results. Figure 4.10 compares full OKE spectra with e–g difference calculated with INM and exact MD methods as well as corresponding results of spectra densities. Highlighted is the region where the INM prediction works well in 10% S system — for frequencies above 25 cm−1 . The top panel compares the INM prediction to the OKE spectra for the 10% S with an excited-state solute with the exact MD results (showing both reduced spectral densities and spectral densities). The overall OKE spectra shapes are in a quantitative agreement, where the spectral densities agree better with the INM predictions than the reduced spectral densities. The bottom panel shows the solute-pump/solvent-probe spectrum at the large T limit, which reflects the fact that the INM theory can capture the distinction between the excited- and ground-state dynamics, although INM prediction for the difference spectrum is less quantitative than for the full OKE spectra. In short, the INM theory can provide us a qualitative and quasi-quantitative estimate of the OKE spectra for the most interesting high frequencies (larger than 25 cm−1 ) and it will be demonstrated to be expandable to the nonequilibrium 2D solute-pump/solvent-probe spectra with finite waiting time T . Figure 4.11 provides comparisons between INM predictions of e–g difference using different system sizes and different averaging with MD results. One can observe that the system size has some influence on the overall magnitude of the e–g difference spectra: the larger system, the larger the e–g difference, because the first-order DID approximated po- larizability increases with system size. But the system-size dependence is not pronounced. One also can observe that more averaging gives a bit lower e–g difference. However, all mentioned variations will not affect the qualitative accuracy of the INM approximation. 161 0.00008 -2 10% S, excited m cm 0.00007 10% S, ground 0.00006 0.00005 4 0.00004 0.00003 0.00002 0.00001 2 0.00000 0 50 100 150 200 -1 /2 c (cm ) (a) 0.0007 -2 50% S, excited m cm 0.0006 50% S, ground 0.0005 4 0.0004 0.0003 0.0002 0.0001 2 0.0000 0 50 100 150 200 -1 /2 c (cm ) (b) Figure 4.7 Instantaneous-normal-mode OKE influence spectra of (a) 10% S and (b) 50% S on the ground state (black curve) and on the excited state (red curve) at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80 and averaged 100,000 configurations (sampling 1 configuration every 10 time steps, and time step is δ t = 0.0025τLJ ). Normalization of the area of the influence spectrum ρ (ω ) of real frequencies to the sum of cα2 with real frequencies is utilized. 162 40 MD / ) INM 30 6 - 5 20 (10 10 g R 0 -50 0 50 100 150 200 -1 /2 c (cm ) (a) 40 MD / ) INM 30 6 - 5 20 (10 e 10 R 0 -50 0 50 100 150 200 -1 /2 c (cm ) (b) Figure 4.8 Comparison of instantaneous-normal-mode (INM) and exact molecular- dynamics (MD) predictions of optical Kerr effect spectra for the 10% S solvent mixture (255 solvents in total) with an atomic solute. The upper panel is the OKE spectrum of the system with a ground-state solute and the lower panel is the OKE spectrum with an excited-state solute. The imaginary INM frequencies are plotted on the negative frequency axis. Molecular dynamics results are reduced spectral densities (solid black curves) that are retrieved from Figure 3.10, calculated by averaging over 108 liquid configurations, and the INM results are calculated by averaging over 1.6 × 106 liquid configurations. The first order DID approximation is utilized in evaluating the many-body polarizability. 163 0.000020 10% S Reduced Spectral Density 0.000015 INM MD 0.000010 0.000005 0.000000 0 50 100 150 200 -1 /2 c (cm ) (a) 0.00002 0.00000 Reduced Spectral Density 50% S -0.00002 INM MD -0.00004 -0.00006 0 50 100 150 200 -1 /2 c (cm ) (b) Figure 4.9 Comparison of e-g difference of the OKE spectra using instantaneous-normal- mode approximation (black) and exact molecular dynamics results (red) for (a) 10% S and (b) 50% S at temperature kB T /ε = 1.00 and density ρσ 3 = 0.80. Molecular dynamics results are reduced spectral densities that are retrieved from Figure 3.10, calculated by averaging over 108 liquid configurations, and the INM results are calculated by averaging over 105 configurations. The first order DID approximation is utilized in evaluating the many-body polarizability. 164 Figure 4.10 Comparison of instantaneous-normal-mode (INM) and exact molecular- dynamics (MD) predictions of optical Kerr effect spectra for the 10% S solvent mixture (255 solvents in total) with an atomic solute. The upper panel is the OKE spectrum of the system with an excited-state solute and the lower panel is the difference spectrum between the excited- and ground-state-solute OKE spectra. Molecular dynamics results including reduced spectral densities (RSD, solid black curves) that are retrieved from Figure 3.10 and spectral densities (SD, dashed black curves) are calculated by averaging over 108 liquid configurations, and the INM results are calculated by averaging over 1.6 × 106 liquid configurations. The first order DID approximation is utilized in evaluating many-body polarizability. The shaded region for frequency less than 25 cm−1 indicates that the INM methods are not suitable for this domain. 165 3.5 MD, N=256, conf=1e8 3.0 INM, N=256, conf=1e5 INM, N=256, conf=1.6e6 / ) 2.5 INM, N=108, conf=1.6e6 6 INM, N=108, conf=1.6e7 2.0 - 5 (10 1.5 g 1.0 -R e R 0.5 0.0 -50 0 50 100 150 200 -1 /2 c (cm ) Figure 4.11 Comparison of instantaneous-normal-mode (INM) and exact molecular- dynamics (MD) predictions of e–g difference of optical Kerr effect spectra for the 10% S solvent mixture with an atomic solute. System sizes (N=108, 256) and numbers of configurations for averaging are indicated in the legend. The imaginary INM frequencies are plotted on the negative frequency axis. The first order DID approximation is utilized in evaluating the many-body polarizability. 166 4.2.5 INM solvation influence spectrum In studying the solvation dynamics, the physical variable A is the solute-solvent interaction energy gap ∆V = Ve −Vg and the solvation influence spectrum is defined as286 * 2 + ∂ ∆V ρsolv (ω ) = δ (ω − ωα ) . (4.59) ∂ qα In a Lennard-Jones system, the coupling coefficients are ∂ ∆V ∂ ∆V ∂ r j µ =∑ , (4.60) ∂ qα j µ ∂ r j µ ∂ qα  " #  σ 12 σ6 ∂ ∆V    ∑ 4∆ε ( j) −12 r14 + 6 r8 (r0µ − r jµ ), (i = 0) j6=0 0j =  0 j (4.61) ∂ ri µ  σ 12 σ6   −4∆ε (i) −12 14 + 6 8 (r0µ − riµ ), (i = 6 0) r0i r0i e g where ∆ε ( j) = εsolute-solvent( j) − εsolute-solvent( j) . (4.62) Figure 4.12 shows how different the influence could be when different physical variable A is used in the coupling factor. Here the polarizability influence and the solvation influence are compared and clearly the solvation influence spectrum has more high frequency contri- bution and much smaller imaginary frequency contribution than the polarizability influence spectrum. 4.3 Evaluation of Solute-Pump/Solvent-Probe Response Although in the large T limit, the response function of a solute-pump/solvent-probe spec- troscopy can be easily calculated by taking the difference of ordinary four-wave-mixing responses (such as optical Kerr effect spectra) on two equilibrated states, h i ˙ ˙ ∆R(0, T, T + t)|T→∞ = β Πxz (0)Πxz (t) e − Πxz (0)Πxz(t) g , (4.63) 167 0.016 0.014 solvation OKE 0.012 0.010 0.008 0.006 0.004 0.002 0.000 -50 0 50 100 150 200 -1 /2 c (cm ) Figure 4.12 Instantaneous-normal-mode polarizability influence spectrum (OKE, red curve) vs. solvation influence spectrum (solvation, black curve) of 10% S solvent mixture with an excited-state solute at thermal condition of temperature kB T /ε = 1.00 and density ρσ 3 = 0.80. Both the polarizability influence spectrum and the solvation influence spectrum are calculated by averaging 20,000 configurations. The imaginary frequencies are plotted on the negative frequency axis. The total areas of both influence spectra are normalized to 1. 168 for finite T , the solute-pump/solvent-probe response function is given by Equation 2.90 which we restate here as follows. D E eβ δ ∆V (0) {Πxz (T ), Πxz (T + t)} ∆R(0, T, T + t) = e − β Π˙ xz (0)Πxz(t) . (4.64) eβ δ ∆V (0) e g As mentioned earlier, a practical problem in evaluating the response is the Poisson bracket, which is a measurement of the interference of two classical trajectories start- ing their evolution from initial conditions with a small difference.287 The difficulty for calculating the Poisson bracket is that for each time a new configuration generated, the whole trajectory of the Poisson bracket has to be updated, which is governed by a different equation of motion than the Newton’s law in molecular dynamics simulation.197 Another concern is the chaotic nature of classical dynamics gives rise to the divergence issue of the stability matrix for some systems236 (the 2nd order stability matrix is the fundamental Poisson bracket),237 and the ensemble averaging can tame the divergence only for a very short time scale.197, 234 The instantaneous-normal-mode theory, on the other hand, provides an analytical approach to evaluate classical Poisson bracket. The INM treatment to the OKE spectra has been discussed in the previous section. In this section, we will focus on the solute-pump/solvent-probe spectra starting with a hybrid approach incorporating the INM approximation and molecular dynamics followed by another analytical attempt. For completeness, the direct numerical approach to evaluate Poisson bracket will be included at the end of this section. 4.3.1 Hybrid INM/MD method There are two time scales associated with the 2D solute-pump/solvent-probe spectra, the solvation time T and the ultrafast dynamical time t. In most of the situations in solvation, these two time scales tend to be widely separated. The solvation time T needs to be long enough to distinguish an overall liquid geometrical reorganization, whereas the t time only 169 needs to be long enough to reflect the intermolecular vibrations. The ultrafast time scale measured by optical Kerr effect spectra lies in sub-picosecond regime corresponding to vibrational frequencies larger than 30 cm−1 . The INM theory gives a reasonably reliable prediction to the ultrafast dynamics corre- sponding to any such high-frequency intermolecular motions.20, 177 We also demonstrated the applicability of INM theory to our preferential solvation model in the previous section. Here we present a hybrid approach, in which the ultrafast dynamics is treated with INM approximation but the long-time solvation is still simulated with the exact molecular dy- namics. Like the INM treatment to OKE spectra, the {Πxz (T ), Πxz(T +t)} Poisson bracket in the solute-pump/solvent-probe response (Equation 4.64) can be evaluated analytically using dynamical information at time T from the primary molecular dynamical trajectory. We perform the INM analysis on each instantaneous configuration RT , obtaining the INM qα (t; RT ), (α = 1, . . . , 3N). Equation 2.95 can be evaluated in that INM basis, 3N   ∂ qβ (t; RT ) {Πxz (T ), Πxz(T + t)} = ∑ Πxz,α (T )Πxz,β (T + t) (4.65) α ,β =1 ∂ pα (0; RT ) 3N   ∂ qβ (t; RT ) ≈ ∑ Πxz,α (T )Πxz,β (T ) , (4.66) α ,β =1 ∂ pα (0; RT ) where the sensitivity of the polarizability with respect to the INM α is ∂ Πxz 1 3N ∂ Πxz Πxz,α (T ) ≡ =√ ∑ · UT (RT ). (4.67) ∂ qα T m i=1 ∂ ri T i,α Taking advantage of the facts that the INM modes are mutually orthogonal and each mode α has a harmonic frequency ωα (RT ), we have the INM fundamental Poisson bracket ∂ qβ (t; RT ) sin[ωα (RT )t] = δαβ . (4.68) ∂ pα (0; RT ) ωα (RT ) In Equation 4.66, the observable is linearized. This approximation (which is an example of 170 the linear INM theory we referred to earlier) will be discussed in more detail in section 4.7. The Poisson bracket we need is then simply 3N sin[ωα (T )t] {Πxz (T ), Πxz (T + t)} = ∑ Π2xz,α (T ) ωα (T ) , (4.69) α =1 where for convenience, ωα (T ) ≡ ωα (RT ).198 As a result, one avoids calculating any secondary (T → T + t) trajectories altogether. The response function that is a three-time ensemble average ∆R(0,T ,T + t) now can be evaluated based on two-time dynamical information (at time 0 and time T ). Substituting this result in Equation 4.64, we have the time-domain solute-pump/solvent-probe response function * + 3N sin ωα (T )t eβ ∆V (0) ∑ Π2xz,α (T ) α =1 ωα (T ) ∆R(0, T, T + t) = e −β Π˙ xz (0)Πxz(t) . (4.70) eβ ∆V (0) e g Taking imaginary part of the Fourier transform of the response with respect to the ultrafast time scale t, Z ∞ ∆R(ω , T ) = dt sin ω t∆R(0, T, T + t), (4.71) 0 we have the solute-pump/solvent-probe spectra * + 3N eβ ∆V (0) ∑ Π2xz,α (T )δ (ω − ωα (T )) * 3N + π α =1 π ∆R(ω , T ) = 2ω eβ ∆V (0) e e − 2ω ∑ Π2xz,α δ (ω − ωα ) α =1 g (4.72) This equation is the most important result of this chapter. The total solute-pump/solvent- probe spectrum consists of a T -dependent part and a T -independent part. The T -dependent part (the first term in above equation) is the nonequilibrium four-wave-mixing spectrum explicitly dependent on the time delay between the solute excitation and the solvent-probe 171 light scattering. The T -independent part (the last term in above equation) is the OKE spectra of solution with a ground-state solute. We can check the long T limit behavior of the nonequilibrium part: when T → ∞, the exponential eβ ∆V (0) gets uncorrelated with the Poisson bracket, giving rise to an ordinary OKE spectra for the solution with an excited- state solute. This long T limit agrees with our exact expression without INM approximation in Equation 2.92. Before showing the calculated 2D spectra, let us see another analytical attempt to evaluate the Poisson bracket. 4.3.2 SHO-like approximation From the previous discussions, the central question in calculating solute-pump/solvent- probe spectra is how to evaluate the correlation function of the following form, R(0,t1,t2) = hA(0){B(t1),C(t2)}i. (4.73) Although the dynamical variable A corresponding to the fluctuation of potential energy gap δ ∆V in Equation 4.64, is only dependent on liquid configuration, for generality, we still treat this variable as a function of coordinates R = {r j µ } and momenta P = {p j µ }, where the label j is the individual atom, and µ = x, y, z is every degree of freedom of an individual atom. A = A(R, P), A(t) = A(R(t), P(t)). (4.74) Since polarizability is a function of positions, we assume B and C are functions of coordinates R alone, B = B(R), C = C(R). (4.75) 172 The Poisson bracket of B and C is then ∂ B (R(t1)) ∂ C (R(t2)) ∂ C (R(t2 )) ∂ B (R(t1 )) {B(t1 ),C(t2)} = · − · (4.76) ∂ R(t1 ) ∂ P(t1 ) ∂ R(t1 ) ∂ P(t1 )   ∂ B (R(t1)) since =0 ∂ P(t1 )   ∂ B ∂ C ∂ rkν (t2 ) = ∑ · (4.77) j µ ,kν ∂ r j µ t1 ∂ rkν t2 ∂ p j µ (t1 ) =∇B(t1) · χ (t1,t2 ) · ∇C(t2), (4.78) where the fundamental Poisson bracket is defined as ∂ rkν (t2) χ j µ ,kν (t1 ,t2) = {r j µ (t1), rkν (t2 )} = , (4.79) ∂ p j µ (t1) ∂ denoting ∇ ≡ ∂R. So the exact T -dependent response function can be rewritten as follows D E eβ ∆V (0) ∇Π(T ) · χ (T, T + t) · ∇Π(T + t) e R(0, T, T + t) = β ∆V (4.80) e e * + ∂ Π xz ∂ Π xz eβ ∆V (0) ∑ χ j µ ,kν (T, T + t) j µ ,kν ∂ r j µ T ∂ r k ν T +t e = β ∆V . (4.81) e e For a 1-dimensional simple harmonic oscillator (SHO), the time evolution of its position is x(t) = x0 cos ω t + mpω0 sin ω t, thus the fundamental Poisson bracket is a function of (t2 − t1 ) and independent of phase space point x(t1 ) or p(t1). ∂ x(t2 ) sin ω (t2 − t1 ) χ (t1 ,t2) = = . (4.82) ∂ p(t1 ) mω The SHO-like approximation to the fundamental Poisson bracket is the following: 173 replacing the fundamental Poisson bracket by its ensemble average and assuming it is diagonal. χ j µ ,kν (t1,t2 ) ≈δ j µ ,kν h{r j µ (t1), rkν (t2)}i (4.83) =δ j µ ,kν h{r j µ (0), r j µ (τ )}i (τ = t2 − t1 )   ∂ r j µ (τ ) =δ j µ ,kν ∂ p j µ (0) ! p2j µ (0) 1 Z Z −β 2m +V ∂ r j µ (τ ) =δ j µ ,kν dr j µ (0) dp j µ (0)e Q ∂ p j µ (0) Z p jµ (0)=+∞ 1 −β H =δ j µ ,kν dr j µ (0)e r j µ (τ ) Q p jµ (0)=−∞ !   −β p2j µ (0) 2m +V Z Z 1 β − δ j µ ,kν dp j µ (0) dr j µ (0) r j µ (τ ) − · 2p j µ (0) e Q 2m Z Z 1β =δ j µ ,kν dp j µ (0) dr j µ (0) r j µ (τ )p j µ (0)e−β H Qm =δ j µ ,kν β hv j µ (0)r j µ (τ )i (4.84) Z τ 1 =δ j µ ,kν Cvv (t)dt. (4.85) m 0 where the normalized velocity autocorrelation functions in one-dimensional and 3N-dimensional cases are respectively hv j µ (0)v j µ (t)i Cvv (t) = , hv2j µ i = kB T /m (4.86) hv2j µ (0)i hv(0) · v(t)i Cvv (t) = , hv2 i = 3NkB T /m (4.87) hv2 (0)i 1 Rτ sin ω t For 1d-SHO, Cvv (t) = cos ω t and χ (τ ) = m 0 Cvv (t)dt = mω . Thus, we get the SHO-like approximated evaluation of the solute-pump/solvent-probe 174 response function,  Z t2 −t1  1  R(0,t1,t2 ) ≈ Cvv (t)dt hA(0)∇B(t1) · ∇C(t2)i hAi , (4.88) m 0 where A = eβ ∆V , B = C = Πxz as in Equation 4.81. Check properties: (1) when t1 = t2 , the response should vanish, and I checked for the SHO-like approxi- mation R(0,t1,t1 ) = 0. (2) if B = C, switching t1 and t2 should result in a sign change of the response, and I R −τ Rτ Rτ checked R(0,t1,t2) = −R(0,t2,t1 ), since 0 Cvv (t)dt = − 0 Cvv (−t ′ )dt ′ = − 0 Cvv (t ′)dt ′. (3) checked that the response function is exact for 1d-SHO. Figure 4.13 shows the comparison of the SHO-like approximation with the exact equilibrium OKE response. At long waiting time T , the response function should converge to the excited-state equilibrium OKE response function. We observe that T = 50 ps and “long T ” curves overlap with each other, suggesting that 50 ps is long enough for the system to get almost relaxed. Also, we can observe that the large T results including both T =50 ps and “long T ” responses are a little higher in magnitude than short T result (T = 0 ps). This observation is in agreement with the fact that in 10% S system, the excited-state OKE response is higher than the ground-state OKE response. For a short T , the geometric structure of system has not changed significantly thus should be close to the equilibrium ground state geometry. Figure 4.13 (b) shows the SHO-like approximation works well in the long T limit for very short time (t < 0.2 ps) comparing T =50 ps and “long T ” results with exact result. But in the longer dynamical time t, SHO-like approach could not give a good estimate of the response (Figure 4.13(a)). Figure 4.14 shows that the magnitude of the infinite order T = 0 ps response is larger than that of the first order T = 0 ps response. However, the infinite order DID approximation and the first order DID approximation are qualitatively equivalent. 175 0.004 ) 1/2 T=0 ps m 0.003 T=50 ps long T 5 ( OKE (excited, eq) Response function 0.002 0.001 0.000 0 2 4 6 8 10 12 14 t (ps) (a) 0.004 ) 1/2 m 0.003 5 ( Response function 0.002 0.001 0.000 0.0 0.2 0.4 0.6 0.8 1.0 t (ps) (b) Figure 4.13 Comparison of response functions using SHO-like approximation at T = 0 ps, T = 50 ps, long T (uncorrelated with exp(β ∆V )) with the equilibrium excited state OKE response function. Equation 4.88 is used toRt calculatefinite T responses, i.e. T = 0 ps, 50 ps. For the long T response, R = (1/m) 0 Cvv (t)dt h∇Π(0) · ∇Π(t)ie . The equilibrium OKE response is evaluated using R(t) = −β (d/dt) hΠ(0)Π(t)ie . All curves in this figure are calculated using the first order DID approximation. Upper panel (a) shows longer time scale and lower panel (b) shows only the first picosecond. 176 0.005 ) 1/2 R (T=0) full DID 0.004 m R (T=0) first order DID 5 0.003 ( Response function 0.002 0.001 0.000 0 2 4 6 8 10 12 14 t (ps) Figure 4.14 Comparison of full DID and first order DID response functions using SHO-like approximation at T = 0 ps. So SHO-like approximation is not good enough to study picosecond dynamics in liquids. The difference between SHO-like approximation and the INM approximation is the basis. In INM treatment, the INM basis is mutually orthogonal, whereas in the SHO- like treatment, the Cartesian basis is usually not orthogonal thus bringing in some error. 4.3.3 Non-INM numerical evaluation Although direct evaluation of the fundamental Poisson bracket is computing-challenging, it is possible to do so.197 The fundamental Poisson bracket has its own equation of motion. Define 3N × 3N Jacobian matrix JRP (t) = χ T (0,t), ∂ rkν (t2 ) χ j µ ,kν (t1,t2) = , (4.89) ∂ p j µ (t1 )  RP  ∂ r j µ (t) J (t) j µ ,kν = . (4.90) ∂ pkν (0) 177 Take partial derivative with respect to pkν (0) on the Newton’s equation, we have   ! ∂ ∂ 1 ∂ V r¨ j µ (t) = − (4.91) ∂ pkν (0) ∂ pkν (0) m ∂ r j µ R(t)  RP  ∂ r¨ j µ (t) 1 ∂ 2V ∂ rl ξ (t) ⇒ J¨ (t) j µ ,kν = =− ∑ · ∂ pkν (0) m l ξ ∂ r j µ (t)∂ rl ξ (t) ∂ pκν (0)   = − ∑ [D(t)] j µ ,l ξ JRP (t) l ξ ,kν . (4.92) lξ The equation of motion for the fundamental Poisson bracket is thus J¨ RP (t) = − D(t)JRP(t), (4.93) and the initial condition is 1 JRP (0) = 0, J˙ RP (0) = 1. (4.94) m One can solve the equation of motion for the fundamental Poisson bracket using a standard MD numerical integrator to propagate Newton’s equation, such as the Verlet algorithm: JRP (t + δ t) =2JRP (t) − JRP(t − δ t) − (δ t)2D(t)JRP(t), (4.95) 1 JRP (0) =0, JRP (−δ t) = − (δ t)1. (4.96) m Or the velocity Verlet algorithm: 1 JRP (t + δ t) =JRP (t) + δ t J˙ RP(t) + (δ t)2 [−D(t)JRP(t)], (4.97) 2 1 J˙ RP (t + δ t) =J˙ RP (t) + δ t[−D(t)JRP(t) − D(t + δ t)JRP (t + δ t)]. (4.98) 2 1 JRP (0) =0, J˙ RP (0) = 1. (4.99) m 178 I tested this numerical approach for our preferential solvation model and got noise. In this thesis, we are not going to pursue this numerical approach any further. 4.4 Short-time Expansion of Solute-Pump/Solvent-Probe Response To test if our hybrid INM/MD method should give us a correct short time behavior, we can perform a Taylor expansion of the response about instantaneous configuration R(T ). The two-dimensional solute-pump/solvent-probe response has the form220 R(0, T, T + t) = hA(0) {B(T ),C(T + t)}i , (4.100) where A, B,C are functions of configuration R only, .D E A(t) = A(R(t)) = eβ δ ∆V (t) eβ δ ∆V (0) , (4.101) e B(t) = B(R(t)) = Πxz (t), (4.102) C(t) = C(R(t)) = Πxz (t). (4.103) The purpose of this section is to show that the short time t expansion of the solute- pump/solvent-probe response is the same as that obtained using the INM approximation and SHO-like approximation. The response function using the INM approximation is   sin ωα t ∂ B ∂ C R(0, T, T + t) = A(0) ∑ . (4.104) α mωα ∂ qα T ∂ qα T +t and the response function using the SHO-like approximation is Z t 1 R(0, T, T + t) = Cvv (τ )dτ hA(0)∇R B(T ) · ∇RC(T + t)i . (4.105) m 0 179 4.4.1 Exact expansion The central part of the response function is the Poisson bracket {B(T ),C(T + t)}. Before showing the expansion of exact response function, we shall expand the configuration about R(T ), ˙ )+ t2 ¨ t 3 ... R(T + t) = R(T ) + t R(T R(T ) + R(T ) + · · · , (4.106) 2! 3! ∂ ∂ denoting ∇ ≡ ∇R = ∂R, ∇P = ∂P. ˙ ) = P(T ) , R(T (4.107) m ¨ ) = F(T ) , R(T (4.108) m   ... ∂ ∇V dR P(T ) R(T ) = − · = −∇∇V (T ) · . (4.109) ∂R m dt R(T ) m2 So, t t2 t3 R(T + t) = R(T ) + P(T ) + F(T ) − ∇∇V · P(T ) + · · · . (4.110) m 2m 3!m2 Then one can expand the Poisson bracket {B(T ),C(T + t)} about R(T ), {B(R(T )),C(R(T + t))} = ∇R B|R(T ) · ∇PC|R(T +t) (4.111)   ∂ R(T + t) = ∇R B|R(T ) · ∇R ,C|R(T +t) · , ∂ P(T ) ∇RC|R(T +t) = ∇RC|R(T ) + ∇R ∇RC|R(T ) · (R(T + t) − R(T )) 1 + ∇R ∇R ∇RC|R(T ) : (R(T + t) − R(T ))2 + · · · . (4.112) 2 Substitute the expansion of R(T + t) and the following fundamental Poisson bracket into 180 the above equation, ∂ R(T + t) t t3 = {R(T ), R(T + t)} = 1 − ∇∇V (T ) + · · · . (4.113) ∂ P(T ) m 3!m2 Finally, we get the exact short-time expansion of the Poisson bracket, t {B(T ),C(T + t)} = [∇B · ∇C]T m t2 + [∇B · ∇∇C · P]T m2 1 t3 − [∇B · ∇∇V · ∇C]T 3! m2 1 t3 + [∇B · ∇∇C · F]T 2 m2 1 t3 + [∇B · ∇∇∇C : PP]T + O(t 4 ), (4.114) 2 m3 and the exact short-time expansion of the solute-pump/solvent-probe response, t hA(0) {B(T ),C(T + t)}i = hA(0)∇B(T ) · ∇C(T )i m t2 + hA(0)∇B(T ) · ∇∇C(T ) · P(T )i m2 1 t3 − hA(0)∇B(T ) · ∇∇V (T ) · ∇C(T )i 3! m2 1 t3 + hA(0)∇B(T ) · ∇∇C(T ) · F(T )i 2 m2 1 t3 + 3 hA(0)∇B(T ) · ∇∇∇C(T ) : P(T )P(T )i + O(t 4). 2m (4.115) 181 4.4.2 INM approximation By means of diagonalizing the dynamical matrix D = ∇R ∇RV |R(0) , (in this section, we use the non-mass-weighted INM, the same as in chapter 1). Λ = UDUT =diag{mωα2 }, (4.116) ∂ qα Uα , j µ = , (4.117) ∂ r j µ R(0) ∂ r j µ UTjµ ,α = , (4.118) ∂ qα R(0) one obtains the independent harmonic INMs (α = 1, . . ., 3N) q =U(R(0)) · [R(t) − R(0)], (4.119) fα vα (0) qα (t) = 2 (1 − cos ωα t) + sin ωα t (4.120) mωα ωα fα 2 =vα (0)t + t +··· , (4.121) 2m pα (t) =m q˙α (t) = m vα (t). (4.122) Now we can expand the Poisson bracket in the INM basis, ∂ B(R(t1)) ∂ C(R(t2 )) {B(R(t1)),C(R(t2))} = ∑ jµ ∂ r j µ (t1 ) ∂ p j µ (t1) ∂ B(R) ∂ C(R) ∂ rkν (t2) =∑∑ , j µ kν ∂ r j µ t1 ∂ rkν t2 ∂ p j µ (t1) 182 where the fundamental Poisson bracket is ∂ rkν (t2) ∂ rkν (t2) ∂ qα (t2 ) ∂ pβ (t1 ) =∑ · · ∂ p j µ (t1) α ,β ∂ qα (t2 ) ∂ pβ (t1 ) ∂ p j µ (t1) ∂ qα (t2) ∂ r j µ (t1) = ∑ UTkν ,α (t2) · · α ,β ∂ pβ (t1) ∂ qβ (t1) ∂ qα (t2) T = ∑ UTkν ,α (t2) · ·U (t1 ) α ,β ∂ pβ (t1) j µ ,β sin ωα (t2 − t1 ) = ∑ UTkν ,α (t2) · δαβ · UTjµ ,β (t1) α ,β mωα sin ωα (t2 − t1 ) T = ∑ UTkν ,α (t2) · · U j µ ,α (t1 ). (4.123) α mωα In above derivation, the second equal sign used the “direct conditions” of the canonical transformation288 (r j µ , p j µ ) → (qα , pα ), ∂ pα ∂ r j µ ∂ pα ∂ p jµ = , =− , (4.124) ∂ p j µ ∂ qα ∂ r jµ ∂ qα ∂ p j µ ∂ qα ∂ r jµ ∂ qα = , =− . (4.125) ∂ pα ∂ r j µ ∂ pα ∂ p jµ ∂ qα (t2 ) sin ωα (t2 −t1 ) fα and the fourth equal sign ∂ pα (t1 ) = mωα is a consequence of qα (t) = mωα2 (1 − pα (0) cos ωα t) + mωα sin ωα t. Then the derivative with respect to the INM is ∂ B ∂ B(R) T ∂ qα t1 ∑ = U (t1), (4.126) jµ ∂ r j µ t1 j µ ,α ∂ C ∂ C(R) T ∂ qα t2 ∑ = U (t2 ). (4.127) kν ∂ rkν t2 kν ,α 183 The INM fundamental Poisson bracket is ∂ qα (t2) sin ωα (t2 − t1 ) χα ,β (t1 ,t2) ≡{qα (t1 ), qβ (t2 )} = = δαβ ∂ pβ (t1) mωα =χα (t2 − t1 )δαβ , (4.128) sin ωα t t 1 t3 χα (t) = = − ωα2 + · · · . (4.129) mωα m 3! m With a full INM treatment without assuming the observable only depends linearly on the coordinates, the Poisson bracket can be expanded into ∂ B ∂ C {B(T ),C(T + t)} = ∑ {qα (T ), qβ (T + t)} (4.130) α ,β ∂ qα T ∂ qβ T +t ∂ B sin ωα t ∂ C =∑ . (4.131) α ∂ qα T mωα ∂ qα T +t Within a linear INM theory (which assumes the observables only depend linearly on ∂ the INM coordinates), we find (using Equation 4.129) with notation ∇q = ∂q, ∂ B sin ωα t ∂ C t ∂ B ∂ C ∑ ∂ qα T mωα ∂ qα T m ∑ = α ∂ qα T ∂ qα T α 1 t3 ∂ B 2 ∂C − 2 ∑ 3! m α ∂ qα T mωα ∂ qα T +··· (4.132) t 1 t3 = [∇q B · ∇qC]T − [∇q B · Λ · ∇qC]T + · · · . (4.133) m 3! m2 For the full INM theory, the nonlinear INM terms result from the following higher-order terms ∂ C ∂ C ∂ 2C 1 ∂ 3C − = ∑ ∂ qα T +t ∂ qα T β ∂ qα ∂ qβ T q β (t) + ∑ q (t)qγ (t) + · · · , 2! β ,γ ∂ qα ∂ qβ ∂ qγ T β (4.134) 184 thus the higher-order terms contributing to the Poisson bracket are (using Equation 4.121)   ∂ B sin ωα t ∂ C ∂ C ∑ ∂ qα mωα ∂ qα − ∂ qα α T T +t T t2 ∂ B ∂ 2C = ∑ v (T ) m α ,β ∂ qα T ∂ qα ∂ qβ T β t3 ∂ B ∂ 3C + ∑ v (T )vγ (T ) 2m α ,β ,γ ∂ qα T ∂ qα ∂ qβ ∂ qγ T β t3 ∂ B ∂ 2C + 2∑ f (T ) + O(t 4) (4.135) 2m α ,β ∂ qα T ∂ qα ∂ qβ T β t2 t3 = [∇q B · ∇q ∇qC · v]T + [∇q B · (∇q ∇q ∇qC : vv)]T m 2m t3 + 2 [∇q B · ∇q ∇qC · f]T + O(t 4 ). (4.136) 2m So the full INM approximation to the Poisson bracket exactly matches the exact answer, Equation 4.114, through O(t 4). To prove this statement, we simply use the invariance of Λ ·∇qC = ∇B·∇∇V · these tensor products of gradients to a basis switch, for example, ∇q B·Λ ∇C. However, the linear INM approximation omits terms involving multiple gradients of observables, such as the second, the fourth and the fifth terms on the right-hand side of the following expansion. t {B(T ),C(T + t)} = [∇q B · ∇qC]T m t2 + [∇q B · ∇q ∇qC · v]T m 1 t3 − [∇q B · Λ · ∇qC]T 3! m2 t3 + [∇q B · ∇q∇qC · f]T 2m2 t3 + [∇q B · (∇q ∇q ∇qC : vv)]T + O(t 4). (4.137) 2m 185 4.4.3 SHO-like approximation The velocity autocorrelation function has a short-time expansion289 hv(0) · v(t)i t 2 h˙v2 (0)i Cvv (t) = = 1 − +··· . (4.138) hv2 (0)i 2 hv2 (0)i t 3 hF2 (0)i Z t So, dτ Cvv (τ ) = t − +··· . (4.139) 0 6 3NmkB T Insert the expansion of ∇C(T + t) (Equation 4.112) into the response   1 t 3 hF2 (0)i R(0, T, T + t) = t− + · · · hA(0)∇B(T ) · ∇C(T + t)i (4.140) m 6 3NmkB T t = hA(0)∇B(T ) · ∇C(T )i m t2 + 2 hA(0)∇B(T ) · ∇∇C(T ) · P(T )i m 1 t3 hF2 (0)i − hA(0)∇B(T ) · ∇C(T )i · 3! m2 3NkB T 1 t3 + hA(0)∇B(T ) · ∇∇C(T ) · F(T )i 2 m2 1 t3 + hA(0)∇B(T ) · ∇∇∇C(T ) : P(T )P(T )i + O(t 4). (4.141) 2 m3 Only the third term in the above equation appears different from the exact expansion 3 − 3!1 mt 2 hA(0)∇B(T ) · ∇∇V (T ) · ∇C(T )i. Although they look similar considering the follow- ing relation Z Z   1 −β V 1 h∇∇V i = dRe ∇∇V = − dR −β e−β V ∇V ∇V = β hFFi, (4.142) Q Q they are fundamentally different because the ensemble average of a product is not equal to a product of ensemble averages. This term probably introduces a relatively large error. Thus the SHO-like approximation only gives the short-time expansion exact through O(t 3). 186 4.5 2D Solute-Pump/Solvent-Probe Spectra for Our Model In the previous section, it was demonstrated that the hybrid INM/MD method is able to handle the ultrafast dynamics in our model at the long waiting time limit (T → ∞) and for finite waiting times, it was shown that the short-time expansion for the nonequilib- rium solute-pump/solvent-probe response is exact through O(t 4) (at least within a full nonlinear INM treatment). In this section, we are going to compute the full 2D solute- pump/solvent-probe spectra for our preferential solvation model. In order to have a good averaging, the simulated molecular system is chosen to be a single solute and 107 solvent atoms. For computational efficiency, the calculations are limited to first-order DID approximated polarizabilities. Figure 4.15 and Figure 4.16 display the full two-dimensional solute-pump/solvent-probe spectra for our 10% S and 50% S preferential solvation model calculated using the hybrid INM/MD method. For each waiting time T , the ultrafast intermolecular vibrational dynamics is evaluated in the basis of the INMs at that time, with the imaginary frequencies plotted on the negative frequency axis. However, from previous INM analysis of OKE spectra at the large waiting time limit, we know that the INM method is likely to overestimate the response for frequencies ω < 25 cm−1 , so the trustworthy region of the INM predictions to the 2D spectra is likely to be frequencies higher than 25 cm−1 . For the 10% S solvent mixture, we can clearly see the progression of the ultrafast vibrational spectra with increasing waiting time T . By contrast, for the 50% S solvent mixture, the progression of the ultrafast vibrational spectra is not so easy to observe. We can find some answer to why it is hard to observe the progression of the ultrafast vibrational spectra in 50% S mixture from Figure 4.17, which reveals the error bars for finite-T spectra for both 10% S and 50% S solvent mixtures and their comparison with the e–g differences. The statistical error at a certain frequency and time delay T is calculated in the following way. For every single configuration, we treat the response function in each histogram bin with associated range of frequency as a data point, and the sample size is 187 the number of configurations, Nconf . The statistical errors for response functions R(ω , T ), Rg (ω ) and ∆R(ω , T ) = R(ω , T ) − Rg (ω ) are defined as σR (ω , T ), σg (ω ) and σ∆R (ω , T ), respectively. The error bars plotted in Figure 4.17 are ±σ∆R (ω , T ) for fixed waiting time T. q 1 σR (ω , T ) = √ h[δ R(ω , T )]2i Nconf − 1 q 1 ≈√ hR2 (ω , T )i − hR(ω , T )i2 , (4.143) Nconf 1 q 2 and σg (ω ) = √ hRg (ω )i − hRg (ω )i2 , (4.144) Nconf q so σ∆R (ω , T ) = σR2 (ω , T ) + σg2 (ω ). (4.145) Accordingly an increase in the number of configurations leads to a smaller error. In the 10% S case, the e-g difference has very small statistical errors across the whole range of frequencies, whereas the T =15 ps spectrum has somewhat bigger error bars but still quite narrow. In the 10% S system there is a clean separation between the finite-T result and the e–g difference due to relatively small fluctuations. On the contrary, the 50% S system has much larger relative fluctuations in the finite-T spectra than the 10% S system, to the point that the 50% S error bars are greater than the difference between the finite-T result and the e–g difference (even though the spectra of the 50% S system are computed by averaging over four times more number of liquid configurations as in the 10% S system). The reason for higher signal-to-noise level in 10% S than in 50% S is unclear, except for the fact that 10% S has a larger relative change in strong solvent population than 50% S. Based on this observation, the solute-pump/solvent-probe spectroscopy may not be keen on picking up the puny change of dynamics in the 50% S system, thus this system may not be a good candidate to test our newly developed theory for the 2D spectra. We will focus on the 10% S system in the rest of the chapter. 188 -6 )) 4.0x10 -2 -6 /(m cm 3.0x10 -6 2.0x10 -6 1.0x10 4 0.0 ( -6 -1.0x10 ,T) -6 -2.0x10 -6 -3.0x10 R( -6 -4.0x10 -6 -5.0x10 10% S -6 -6.0x10 -50 N=108 0 0 5 10 50 20 /2 30 100 40 150 s) 50 c (p 10 (c 200 e- T 0 m g -1 ) Figure 4.15 Two-dimensional solute-pump/solvent-probe spectra for our preferential solvation system 10% S solvent mixtures. The 2D spectra were calculated by hybrid INM/MD methods, and the simulated system contains a single solute and 107 solvent atoms. The frequency ω corresponds to the ultrafast intermolecular vibrational dynamics (t) and the imaginary INM frequencies are plotted on the negative frequency axis. The spectra show how the ultrafast dynamics evolves with the waiting time T after the solute excitation. The large T limit (e–g) is the difference between excited- and ground-state INM OKE influence spectra. The results are averaged over 1.6 × 107 liquid configurations using first-order many-body polarizabilities. 189 -1 ) m (c c /2 T ( ps e-g 0 200 ) 10 50 150 40 100 30 20 50 10 0 0 -50 -6 5.0x10 )) -2 /(m cm 0.0 -6 -5.0x10 -5 4 -1.0x10 ( -5 -1.5x10 ,T) -5 50% S -2.0x10 R( N=108 -2.5x10 -5 Figure 4.16 Two-dimensional solute-pump/solvent-probe spectra for our preferential solvation system 50% S solvent mixtures. The 2D spectra were calculated by hybrid INM/MD methods, and the simulated system contains a single solute and 107 solvent atoms. The frequency ω corresponds to the ultrafast intermolecular vibrational dynamics and the imaginary INM frequencies are plotted on the negative frequency axis. The spectra show how the ultrafast dynamics evolves with the waiting time T after the solute excitation. The large T limit (e–g) is the difference between excited- and ground-state INM OKE influence spectra. The results are averaged over 1.6 × 107 liquid configurations using first- order many-body polarizabilities. 190 -6 5.0x10 e-g -6 /(m cm )) 4.0x10 T=15 -2 -6 3.0x10 -6 4 2.0x10 ( ,T) -6 1.0x10 R( 0.0 -6 -1.0x10 -50 0 50 100 150 200 -1 /2 c (cm ) (a) -6 5.0x10 0.0 /(m cm )) -2 -6 -5.0x10 -5 -1.0x10 50%S 4 -5 -1.5x10 ( N=108 -5 ,T) -2.0x10 e-g -5 -2.5x10 T=5ps R( -5 -3.0x10 -5 -3.5x10 -50 0 50 100 150 200 -1 /2 c (cm ) (b) Figure 4.17 Comparison of the fluctuations of solute-pump/solvent-probe responses at a finite waiting time T with that of the equilibrium e–g difference in INM OKE influence spectra for our preferential solvation system including a single solute and 107 solvent atoms. (a) 10% S solvent mixture, T =15 ps and e–g difference; (b) 50% S solvent mixture, T =5 ps and e–g difference. The results for 10% S and 50% S mixtures are calculated by averaging 1.6 × 107 and 8 × 107 liquid configurations, respectively, using first-order many- body polarizabilities. Note the two panels are different in scale. 191 4.6 Structural Information in Preferential Solvation 4.6.1 First-shell population dynamics The most important result of this chapter, Figure 4.18, depicts the real-frequency 2D solute- pump/solvent-probe spectra for the 10% S solvent mixture. (Note that the direction of the waiting time T axis is opposite of that in Figure 4.15.) The T -dependent ultrafast vibrational spectra evidently converge to the equilibrium e–g difference in the optical Kerr spectra within 50–100 ps. The time scale for this convergence is broadly consistent with the time scale for the solvation correlation function of the same 10% S solvent mixture (C(t) in Figure 3.3). However, to describe the relaxation observed by the 2D solute-pump/solvent- probe spectra quantitatively, one can define a normalized relaxation profile for the spectra at any chosen vibrational frequency ω : ∆R(ω ; T ) − ∆R(ω ; T = ∞) S(T ; ω ) = . (4.146) ∆R(ω ; T = 0) − ∆R(ω ; T = ∞) This profile S(T ; ω ) reflects the relaxation of the ultrafast dynamics that the spectroscopy is most sensitive to, at any given frequency. In Figure 4.19, the spectroscopic relaxation profiles S(T ) at fixed frequencies 50 cm−1 and 75 cm−1 are compared with the solute- solvent interaction energy relaxation C∆V (t), the first-shell population relaxation correla- tion function CnS (t) and their cross correlation C∆V,nS (t). These correlation functions are defined as follows. hδ ∆V (0)δ ∆V (t)ie C∆V (t) = , (4.147) h(δ ∆V )2 ie hδ nS (0)δ nS (t)ie CnS (t) = , (4.148) h(δ nS )2 ie hδ ∆V (0)δ nS (t)ie C∆V,nS (t) = , (4.149) h(δ ∆V )(δ nS )ie 192 where nS is the number of S solvents in the first solvation shell of the solute and its fluctuation is defined as δ nS ≡ nS − hnS ie. Thus, correlation function CnS (t) is a direct measure of the solvent structural relaxation. Figure 4.19 displays a remarkable agreement between this structural correlation func- tion and the spectroscopic relaxation profiles S(T ). The spectroscopic relaxation profiles at two frequencies (50 cm−1 and 75 cm−1 ) are virtually identical. In addition, the mixed energy/structure cross correlation function C∆V,nS (t) agrees quantitatively with the spectro- scopic profiles. However, the spectroscopic profiles are noticeably faster than the standard potential-energy solvation relaxation C∆V (t). In other words, the solute-pump/solvent- probe spectra seems to be able to single out the explicitly structural dynamics featured by the first-shell population relaxation rather than the energetic dynamics described by the potential-energy solvation relaxation. Therefore, the two-dimensional spectra offers a new molecular perspective on solvation dynamics by reporting on direct structural dynamics of the solvent. This structural dynamics has a different relaxation rate from the energetic dynamics measured by traditional time-dependent fluorescence. The error bars of the spectroscopic profiles S(T ) are essential — they are small enough to separate spectroscopic profiles from the potential-energy solvation relaxation. The determination of the error bars of S(T ) follows those formulas below. For any variable X = u/v, the standard deviation of X satisfies  2  2 ∂X ∂X σu2 σv2 u2 σX2 = σu2 + σv2 = + 4 ∂u u ∂v v v2 v  2 σX  σ 2 u  σ 2 v ⇒ = + . (4.150) u/v u v 193 Applying the above relation, we have the errors for S(T ), R(ω ; T ) − Re (ω ) S(T ; ω ) = (4.151) R(ω ; T = 0) − Re (ω )    2  2 σS(T ;ω ) 2 σR(ω ,T )−Re (ω ) σR(ω ,0)−Re (ω ) = + , (4.152) S(T ; ω ) R(ω , T ) − Re (ω ) R(ω , 0) − Re(ω ) where 2 2 2 σR( ω ,T )−Re (ω ) = σR (ω , T ) + σe (ω ), (4.153) 1   σe2 (ω ) = hR2e (ω )i − hRe(ω )i2 . (4.154) Nconf The error bars plotted in Figure 4.19 are ±σS(T ;ω ) . Figure 4.20 presents the evolution of strong solvent population distributions in the first shell and second shell for our preferential solvation systems 10% S and 50% S. The first- shell population varies more dramatically than the second-shell population in both systems and exhibits no significant visible change after T = 60 ps in 10% S and after T = 40 ps in 50% S. Figure 4.21 compares the structural correlation with the energetic correlation along with their cross correlation in the 50% S system. The difference between the first- shell population relaxation and the solute-solvent energy correlation in the 50% S system is not as discernible as that in the 10% S system. The indistinguishable structural dynamics and energetic dynamics in the 50% S solvent mixture also make it difficult to tell what microscopic information the 2D spectra pick up. It is important to note that the 2D nonequilibrium spectra are more sensitive to structural changes in some systems than in others. Besides noting that the 10% S has a better signal than 50% S system, one can switch the polarizabilities of S and W solvents (with αS = 0.73A˚ 3 and αW = 3.99A˚ 3 ) in the same preferential solvation system with 10% S solvent. The resulting 2D spectra are plotted in Figure 4.22 and look very different from the normal polarizability case plotted in Figure 4.18. Although the Hamiltonians, initial conditions and 194 hence the dynamics of the two 10% S systems are absolutely identical, now the primary structural change during the preferential solvation involves approximately a 130% increase in the spectroscopically dark S solvents in the first shell and only 12% decrease in the abundant spectroscopically bright background W solvents. It is not surprising to see almost no progression of the 2D spectra after the first 2 ps in Figure 4.22. So why could this spectroscopy serve as a good structural probe and why is the spectrum closely connected to the first-shell population dynamics? According to our previous T → ∞ analysis, the solvation shells beyond the first also contribute significantly to the difference spectrum. The key to answer these questions is related to how the solvation shell populations change. We will explore another piece of structural information in our preferential solvation model, the distributions of strong-solvent/solute/strong-solvent bond angles, in the following section. 4.6.2 Bond-angle distribution For the sake of understanding the solute-pump/solvent-probe spectroscopy, we look at the T -dependent bond-angle distribution of solvent pairs with respect to the solute, which re- flects the solvent orientational dynamics within the region close to the solute. Considering a solvent mixture with strong and weak kinds of solvent, we focus on strong solvent and denote θ as the bond angle between different solute-strong-solvent pairs in a given shell (or pair of shells), see Figure 4.23. We define the distribution of cos θ as follows. 1  ˆ j (T ) · Ω  ˆ k (T ) , ρ (cos θ )T = #SSpairs ∑δ cos θ − Ω (4.155) j