MIXED QUANTUM/SEMICLASSICAL THEORY FOR SMALL-MOLECULE DYNAMICS AND SPECTROSCOPY IN LOW-TEMPERATURE SOLIDS by XIAOLU CHENG A DISSERTATION Presented to the Department of Physics and the Graduate School of the University of Oregon in partial fulfillment of the requirements for the degree of Doctor of Philosophy March 2013 DISSERTATION APPROVAL PAGE Student: Xiaolu Cheng Title: Mixed Quantum/Semiclassical Theory for Small-Molecule Dynamics and Spectroscopy in Low-Temperature Solids This dissertation has been accepted and approved in partial fulfillment of the requirements for the Doctor of Philosophy degree in the Department of Physics by: Dr. Miriam Deutsch Dr. Jeffrey A. Cina Dr. Jens U. No¨ckel Dr. Steven J. van Enk Dr. Marina G. Guenza Chair Advisor Member Member Outside Member and Kimberly Andrews Espy Vice President for Research and Innovation and Dean of the Graduate School Original approval signatures are on file with the University of Oregon Graduate School. Degree awarded March 2013 ii c©March 2013 Xiaolu Cheng iii DISSERTATION ABSTRACT Xiaolu Cheng Doctor of Philosophy Department of Physics March 2013 Title: Mixed Quantum/Semiclassical Theory for Small-Molecule Dynamics and Spectroscopy in Low-Temperature Solids A quantum/semiclassical theory for the internal nuclear dynamics of a small molecule and the induced small-amplitude coherent motion of a low-temperature host medium is developed, tested and applied to simulate and interpret ultrafast optical signals. Linear wave-packet interferometry and time-resolved coherent anti- Stokes Raman scattering signals for a model of molecular iodine in a 2D krypton lattice are calculated and used to study the vibrational decoherence and energy dissipation of iodine molecules in condensed media. The total wave function of the whole model is approximately obtained instead of a reduced system density matrix, and therefore the theory enables us to analyze the behavior and the role of the host matrix in quantum dynamics. This dissertation includes previously published co-authored material. iv CURRICULUM VITAE NAME OF AUTHOR: Xiaolu Cheng GRADUATE AND UNDERGRADUATE SCHOOLS ATTENDED: University of Oregon, Eugene, Oregon Nankai University, Tianjin, China DEGREES AWARDED: Doctor of Philosophy in Physics, 2013, University of Oregon Bachelor of Science in Physics, 2005, Nankai University AREAS OF SPECIAL INTEREST: Quantum Mechanics, Computational Physics Medical Imaging Numerical Simulation PROFESSIONAL EXPERIENCE: Graduate Research Assistant, University of Oregon, 2006 – 2013 Graduate Teaching Assistant, University of Oregon, 2005 – 2012 GRANTS, AWARDS AND HONORS: Meng Fellowship, University of Oregon, 2006 – 2007 Student Scholarship, Nankai University, 2001 – 2005 v PUBLICATIONS: Craig T. Chapman, Xiaolu Cheng, and Jeffrey A. Cina, “Numerical Tests of a Fixed Vibrational Basis/Gaussian Bath Theory for Small Molecule Dynamics in Low-Temperature Media” J. Phys Chem. A, 115, 3980, (2011). The first two authors contributed equally to this work. vi ACKNOWLEDGEMENTS I would like to thank my thesis advisor, Jeff Cina, who is my role model of a serious and happy researcher. vii To my parents and all my dear friends, who listened and talked to me during difficult times. viii TABLE OF CONTENTS Chapter Page I. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 II. NUMERICAL TESTS OF A FIXED VIBRATIONAL BASIS/GAUSSIAN BATH THEORY FOR SMALL MOLECULE DYNAMICS IN LOW-TEMPERATURE MEDIA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 II.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 II.2. Summary of FVB/GB and Setup for Numerical Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 II.2.1. FVB/GB Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 II.2.2. Coupled Harmonic Oscillators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 II.2.3. Are the Bath Wave Packets Gaussian? . . . . . . . . . . . . . . . . . . . . 17 II.2.4. Population Nonconservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 II.3. Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 II.3.1. J = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 II.3.2. J = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 II.3.3. Electronic Absorption Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . 29 II.4. Comments and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 III. VARIATIONAL FVB/GB AND MODEL OF MOLECULAR IODINE IN A 2D KRYPTON LATTICE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 III.1. Variational FVB/GB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 III.1.1. Time-Dependent Variational Principle in Quantum Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 III.1.2. Equations of Variational FVB/GB . . . . . . . . . . . . . . . . . . . . . . . . 42 III.1.3. Comparison of FVB/GB and Variational FVB/GB . . . . . . . . 46 III.2. The Model of Molecular Iodine in a 2D Kr Lattice . . . . . . . . . . . . . 49 ix Chapter Page III.2.1. Model System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 III.2.2. Equilibrium Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 III.2.3. System-Bath Partition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 III.2.4. Taylor Expansion of the Potential. . . . . . . . . . . . . . . . . . . . . . . . . 58 III.2.5. The Purity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 III.2.6. Cage Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 IV. LINEAR WAVE-PACKET INTERFEROMETRY AND TR-CARS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 IV.1. Linear Wave-Packet Interferometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 IV.1.1. Theoretical Expression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 IV.1.2. Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 IV.2. Tr-CARS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 IV.2.1. Theoretical Expression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 IV.2.2. Bath Wave Packets after Pulse-Driven Transitions. . . . . . . . . 90 IV.2.3. Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 V. DISCUSSION AND FUTURE DIRECTIONS . . . . . . . . . . . . . . . . . . . . . . . . . 97 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 A. FVB/GB EQUATIONS OF MOTION FOR A MULTIDIMENSIONAL BATH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 B. SOME INTEGRALS USED IN CALCULATING POTENTIAL ELEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 C. NUMERICAL STABILIZATION. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 REFERENCES CITED . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 x LIST OF FIGURES Figure Page II.1. Position and momentum expectation values for bath wave packets accompanying selected vibrational levels (J = 0.01) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 II.2. Real and imaginary parts of selected bath wave packet width parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 II.3. Populations of system vibrational levels ν = 0, 2, and 4 . . . . . . . . . . . . . . . 23 II.4. Coordinate expectation values and imaginary width parameters for the ν = 0, 2, and 4 Gaussian bath wave packets (J = 0.05) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 II.5. Selected level populations and the total population of all system vibrational states (J = 0.05) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 II.6. Fidelities of FVB/GB bath wave packets belonging to several system vibrational states (J = 0.05) . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 II.7. Expectation value of the position coordinate of the bath wave packet attached to the ν = 1 level for J = 0.05 . . . . . . . . . . . . . 27 II.8. Magnitude of the Gaussian portion and the Hermite polynomial prefactor of the exact bath wave packet (J = 0.05) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 II.9. Electronic absorption spectrum of a two-mode model . . . . . . . . . . . . . . . . . . 31 III.1. Total population of all system vibrational states in FVB/GB and in variational FVB/GB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 III.2. Fidelity of the overall wave function in FVB/GB and in variational FVB/GB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 III.3. 2D I2/Kr model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 III.4. Potential of Kr-Kr interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 III.5. Purity of the reduced density matrix for the 2D I2/Kr model . . . . . . . . . . 64 III.6. System potential in the excited (B) state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 III.7. System potential in the ground state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 IV.1. In-phase and in-quadrature interferograms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 IV.2. Absorptive portion of the linear susceptibility . . . . . . . . . . . . . . . . . . . . . . . . . 76 IV.3. Motions of two highest-frequency bath modes . . . . . . . . . . . . . . . . . . . . . . . . . 77 IV.4. Energy expectation value for the highest-frequency and second-to-highest-frequency bath modes . . . . . . . . . . . . . . . . . . . . . . . . . . 78 xi Figure Page IV.5. Expectation value of the highest-frequency and second-to-highest-frequency bath coordinates . . . . . . . . . . . . . . . . . . . . . . . . . 79 IV.6. Absolute value of tr-CARS dipole moment contribution M (0) 123(t32, t43) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 IV.7. Absolute value of tr-CARS dipole moment contribution N (0) 213(t32, t43) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 IV.8. Tr-CARS signal for the 2D I2/Kr model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 xii LIST OF TABLES Table Page III.1. Potential parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 xiii CHAPTER I INTRODUCTION Quantum mechanics has proven successful in describing atomic-scale systems and predicting relevant experimental results, such as the absorption spectrum of a hydrogen atom. For many decades, researchers have been intrigued by other phenomena that classical theory cannot fully explain but quantum theory can, especially those relating to wave characteristics. One of the central theoretical and practical questions is whether quantum mechanics can be applied to larger- scale systems, like molecules and biological structures, because indeed mother nature plays a lot of tricks in these complex systems that challenge a first- principles understanding. In chemical reactions and biological processes—on the picosecond and femtosecond timescales—transformations accessible to time- resolved experiments have been the subject of considerable interest. In addition to seeking to elucidate the dynamics underlying these processes, scientists have ambitions to control these reactions and transformations, by altering their spatial and temporal environments. Take a classical example: if we want to slow a car on the road (provided that we cannot control the driver), what we can do is to increase the vehicle’s coefficient of friction by scattering some gravel on the road—changing the car’s environment. In the quantum world, bath-mediated dynamics has been 1 an interesting topic investigated by many physicists and chemists. Poulin and Nelson monitored the photolysis of the triiodide anion I−3 in three very different pure organic molecular crystals (TBAT, TEAT, and TPPT) [1]. They found that the photofragments I−2 and I recombined subsequently in the constrained environments of TBAT and TEAT, but did not do so in the least constrained crystal, TPPT. Their results illustrate the intimate connection between lattice structure and reaction dynamics. It has also been demonstrated that in a pigment- protein complex Fenna-Matthews-Olson (FMO), dephasing noise from molecular and intermolecular vibrations may support electronic excitation energy transfer from the light-harvesting chlorosomes to the bacterial reaction center [2]. Along with these findings, a number of cage-mediated events and dephasing-assisted processes have inspired researchers to conduct more studies of the role of the environment in directing reactions and system dynamics. To study bath-mediated dynamics, we focus on vibrational coherence, which is in the picosecond regime, while electronic coherence is typically on the order of femto to pico seconds for condensed systems. The Apkarian group at UC-Irvine and the Schwentner group at Freie Universita¨t Berlin have performed several types of ultrafast optical experiments on dihalogen molecules embedded in rare gas solids to elucidate various aspects of coherent motions [3–14]. These are attractive systems for investigation, since they combine the simplicity of a diatomic target molecule with the complexity of its interaction with a well-structured surrounding medium. 2 Halogen diatomics are excellent target species, as abundant data exist on their electronic and vibrational states. For I2 in solid Kr, a vibrational dephasing time of hundreds of picoseconds is observed [7], and the dependence of the coherence decay rate on temperature and vibrational level is accurately measured and analyzed as well [9], which provide valuable information to establish a model of pure dephasing. A recent report of spectrally resolved, 4-wave mixing measurements in five resonant colors on the same system [15] necessitates deeper discussions of fundamental principles of quantum mechanics of molecular matter, among them: the distinction between quantum and classical coherent dynamics of a system entangled with the environment, event-driven decoherence, and environment selected coherent states. Rohrdanz and Cina demonstrated one effect of coherent bath dynamics on a system by showing the difference between two kinds of tr-CARS (time-resolved coherent anti-Stokes Raman scattering) signals from iodine molecules in a 1D Ar chain [16]: in one case there is a pre-pulse that excites a remote Ca atom in the chain, and in the other there is not. In their calculations, a multidimensional Gaussian wave packet was used to describe all solute and solvent nuclear degrees of freedom. Chapman and Cina first put forward a Fixed Vibrational Basis/Gaussian Bath (FVB/GB) theory to analyze small molecule dynamics in low-temperature media [17]. The improvement from Rohrdanz’s approximation to FVB/GB is that Gaussian wave packets are only used for the bath degrees of freedom, while the system is treated fully quantum mechanically. They make use of the fact that in 3 ultrafast experiments on dihalogens in noble gas matrices, a few high-frequency intramolecular degrees of freedom are driven to a large-amplitude motion, while motion in the low-frequency lattice modes is induced indirectly by the system vibration. In this dissertation, we further investigate FVB/GB theory and its applications. First we implement FVB/GB theory for a numerical test model of bilinearly coupled harmonic oscillators. The predictions of FVB/GB theory are compared with results obtained using an exact basis-set method. We find that for significant lengths of time FVB/GB does a good job of following the exact results. Meanwhile, several issues arise that require further attention, such as numerical instability and population nonconservation. In the light of Burghardt’s work [18] and earlier studies [19], we adopt a variational approach and derive equations of motion for the bath wave-packet parameters from Lagrange’s equations. The total population and energy are both shown to be conserved during propagation under this variational scheme. We then apply variational FVB/GB to a realistic model of molecular iodine in a 2D krypton lattice, evaluate the reduced system dynamics under nonstationary initial conditions, and calculate linear wave-packet interferometry and tr-CARS signals from this model sample. We aim to bridge the gap between experimental results and theoretical models of dynamics by taking advantage of simulations on a realistic model system. Besides ultrafast optical signals, we can calculate any physical quantities, e.g., Wigner 4 distribution function, bath energy, and coherence, which are significant in verifying or rejecting any proposed dynamical mechanism. Gu¨hr, Bargheer and Schwentner observed long lived coherent zone boundary phonons (ZBP) in ultrafast pump- probe spectra of I2 guest molecules in a Kr crystal [10]. These gave rise to a sharp peak single line in the Fourier transform of the pump-probe signal, which corresponds to the frequency at the boundary of the first Brillouin zone in a phonon dispersion curve of solid Kr along [100]. A mechanism of displacive excitation of coherent phonons (DECP) was proposed to explain the generation of ZBP: the electronic excitation of I2 changes the equilibrium position of some Kr atoms in the vicinity whose motion thereafter is parallel to the I2 vibration and therefore is decoupled from the intramolecular mode. Unlike propagating phonons whose amplitudes are damped relatively soon, zone boundary phonons have a very low group velocity, which makes them stay in the vicinity of the guest molecule for as long as 10 ps. Karavitis and Apkarian found a similar local mode from time- resolved coherent anti-Stokes Raman scattering [6]. Besides overtone beats of the ground state vibrational frequency, a beat at 41.5 cm−1 is discernable in the Fourier transform of their tr-CARS signal. They attempted an interpretation of this amplitude modulation as intimate coupling between the molecule and a local mode. Based on its linewidth, this local mode persists for as long as 100 ps. Different from the DECP scheme, Karavitis and Apkarian proposed a scenario in which the iodine molecule evolves on a dissociative but cage-bound potential and 5 the caging is described as a sudden process. After collision with the cage wall, iodine atoms lose most of their vibrational energy and the overdriven cage rebounds with a characteristic period. With the powerful tool of FVB/GB, we are well equipped to determine which dynamical mechanism is closer to reality. This dissertation is organized as following: in Chapter II, we provide a review of the FVB/GB theory and present numerical testing results on bilinearly coupled harmonic oscillators; the conditions of validity for the semiclassical approximations are also discussed, as well as the problem of population nonconservation; in Chapter III, we introduce the variational approach and the model of molecular iodine in a 2D krypton lattice; in Chapter IV, we describe the simulations of linear wave-packet interferometry and tr-CARS signals and exhibit the simulation results; in Chapter V, we summarize and the future directions of this work are addressed. Chapter II was published and co-authored with Craig T. Chapman and Jeffrey A. Cina. 6 CHAPTER II NUMERICAL TESTS OF A FIXED VIBRATIONAL BASIS/GAUSSIAN BATH THEORY FOR SMALL MOLECULE DYNAMICS IN LOW-TEMPERATURE MEDIA This work was published in the Journal of Physical Chemistry A, 115, 3980, (2011). It was initiated by Jeffrey A. Cina; Craig T. Chapman and Xiaolu Cheng performed the calculations; Jeffrey A. Cina was the principle investigator for this work. II.1. Introduction Because of the exponential scaling with system size that would be entailed in direct basis-set applications of quantum mechanics to spectroscopic and dynamical simulations of molecules in condensed media, we must often resort to reduced descriptions or to classical or semiclassical treatments of all or part of the system. In reduced dynamical descriptions, the simplifying assumption of an unperturbed near-equilibrium medium (bath) state is frequently invoked, which sacrifices the prospect of predicting or interpreting induced nonequilibrium coherent motion in the host liquid or solid. Classical or semiclassical approximations enable a full dynamical simulation of the environment of a target molecule, but even the 7 simplest quantum mechanical effects, such as zero-point motion, often cannot be rigorously handled and as a consequence such treatments must typically assume a high-temperature bath. There may be no magical short-term solution to these conundrums, but the limitations they impose accentuate the relevance of specific, chemically interesting situations in which they can be bypassed or avoided. In addition, recent advances [20, 21] in the basic description of the equilibration of a quantum mechanical system immersed in a larger quantum mechanical medium, which succeed in removing the longstanding subjective requirement of a lack of precise knowledge of the initial state of the combined system and environment, could be helpfully augmented with practical demonstration calculations on realistic molecular systems if these become possible. The behavior of small molecules in low-temperature solid crystals, including electronic and vibrational dynamics [3–9, 22, 23] as well as chemical reactions, [1, 24, 25] is affected in fundamental ways by interactions with the environment; these interactions can also generate coherent nonequilibrium environmental states that are not swamped by thermal agitation. [10, 11] Electronic and vibrational excitation-induced dynamics in systems of this kind are being studied by myriad time-resolved nonlinear optical techniques. An interpretation of the resultant spectroscopic signals in terms of the underlying quantum mechanical motion is of keen interest in its own right, for the possible demonstration of coherent control, [12, 8 26, 27] and for potential quantum-information applications. In addition, several key features of small-molecule chromophores embedded in cryogenic noble gas hosts may enable their comprehensive simulation with computationally efficient first-principles approaches that provide a rigorous quantum mechanical description of a small number of directly excited intramolecular degrees of freedom along with a well-defined, albeit approximate, treatment of the indirectly triggered semiclassical dynamics of the surrounding medium. Chapman and Cina have put forward two related quantum/semiclassical theories tailored to the specific properties of small molecules in rare-gas matrices. [17] Both theories take advantage of the timescale difference between high-frequency intramolecular vibrations and the lower frequencies of acoustic lattice phonons in a host crystal, along with the weak coupling between them. The strategy in both cases is to accept exponential scaling with respect to a small number of directly excited, mutually interacting, perhaps anharmonic intramolecular vibrations whose energy levels are widely spaced and hence few in number. For the numerous low- frequency, nearly (but not quite) harmonic, almost (but not strictly) independent lattice modes whose motion is partially (but perhaps not exclusively) indirectly excited through their anharmonic coupling to the high-frequency vibrations, recourse is made to a semiclassical description that scales according to a power law with respect to the size of a potentially macroscopic environment. [28–31] Chapman and Cina’s first scheme, called the fixed vibrational basis/Gaussian 9 bath (FVB/GB) theory, expands the time-dependent state of the system in a single basis of energy eigenstates determined at a frozen bath geometry. Their second, adiabatic vibrational basis/Gaussian bath (AVB/GB) approach makes use of a vibrational basis parametrized by the variable bath geometry; it takes advantage of a timescale separation between the system and bath, but does not make an adiabatic approximation per se. Both theories incorporate nonequilibrium bath dynamics with a semiclassical approximation, which is based on the short timescales of ultrafast spectroscopy rather than high temperature, by constraining the multidimensional bath wave packets associated with each system basis state to take the form of a time-dependent multidimensional Gaussian function. The equations of motion for the parameters specifying the various evolving bath wave packets, which must in practice be solved numerically, are simpler for FVB/GB than for AVB/GB, and it is the former version of the theory on which we focus here. The next section of this chapter provides a short review of the FVB/GB approach, details some aspects of its approximate handling of the bath wave packets and overall population conservation, and introduces the model Hamiltonian on which it will be tested. We then implement the theory numerically for a test case of coupled harmonic oscillators, using it to follow the time development of the state-defining bath parameters and to calculate an electronic absorption spectrum. The final section summarizes, compares the FVB/GB theory to related work, 10 and sets the stage for its pending application to simulations of the spectra and dynamics of dihalogen and other molecules embedded in multiatom rare-gas hosts. II.2. Summary of FVB/GB and Setup for Numerical Tests II.2.1. FVB/GB Equations A comprehensive exposition of the FVB/GB theory was given earlier; [17] here we shall summarize and then test the treatment of small-molecule dynamics in media by considering a pair of coupled vibrational degrees of freedom whose dynamics is governed by a well-defined nuclear Hamiltonian. In systems amenable to treatment by FVB/GB, there is a significant timescale separation between one or more high-frequency “system” oscillators and a lower-frequency, multimode “bath.” The mass-weighted normal coordinates of our test model’s single-mode system and bath are referred to as q and Q, respectively, and the corresponding momenta are referred to as p and P . The full nuclear Hamiltonian is h = hs + hb + u(q,Q), (II.1) where hs = (p 2/2) + us(q) and hb = (P 2/2) + ub(Q). 1 If q and Q are normal coordinates of the composite system, then a power series expansion of the system- 1In this chapter, all quantities are rendered “dimensionless” by expressing them as multiples of the appropriate combination of reference mass, distance, and time units—m0, d0, and t0, respectively—two of which may be chosen freely and the third of which is determined by the requirement that ~ is equal to m0d20t −1 0 ; the constant ~ takes a unit value in this system and does not appear explicitly. 11 bath interaction potential u(q,Q) begins with third-order terms proportional to q2Q and qQ2. The time-dependent nuclear ket of the sample can be expressed in general as a sum of tensor product states, |Ψ(t)〉 = ∑ ν |ν〉〈ν|Ψ(t)〉 = ∑ ν |ν〉|ψν(t)〉, (II.2) in which the system vibrational eigenstates obey hs|ν〉 = εν |ν〉 and the accompanying bath wave packets, |ψν〉 ≡ 〈ν|Ψ〉, carry all the time dependence. The equations of motion for the bath wave packets, i ∂ ∂t |ψν(t)〉 = (εν + hb + uνν(Qˆ))|ψν(t)〉+ ∑ ν¯ 6=ν uνν¯(Qˆ)|ψν¯(t)〉, (II.3) are obtained by substituting Eq. (II.2) into the time-dependent Schro¨dinger equation, i|Ψ˙(t)〉 = h|Ψ(t)〉, and taking the inner product with a specific vibrational eigenstate. Eq. (II.3) uses the notation uνν¯(Qˆ) = 〈ν|u(qˆ, Qˆ)|ν¯〉. Eq. (II.3) remains exact, but the molecular systems studied in the experiments of interest here consist of isolated chromophores that are externally driven to large- amplitude motion within a crystalline host. The bath is weakly disturbed from equilibrium either directly, by a laser-induced resonant transition under which some bath modes are Franck–Condon-active, or indirectly, through coupling to the laser-driven system vibration. Although some bath nuclei may experience large-amplitude displacements, the motion of any single nucleus would typically result from a superposition of smaller-amplitude motion in many spatially extended 12 modes. In this situation and in the absence of a thermal population at cryogenic temperatures, it is appropriate to make an approximate, Gaussian ansatz. In our 1D bath, the wave packets are therefore assigned the forms 〈Q|ψν〉 = exp [iαν(Q−Qν)2 + iPν(Q−Qν) + iγν ], (II.4) which are specified by time- and ν-dependent parameters Qν = 〈ψν |Qˆ|ψν〉/〈ψν |ψν〉, Pν = 〈ψν |Pˆ |ψν〉/〈ψν |ψν〉, αν = α′ν + iα ′′ ν , and γν = γ ′ ν + iγ ′′ ν . In the case of a multidimensional bath, [17] Qν and Pν become vectors, αν becomes a complex- symmetric matrix, and γν remains a scalar. For the Gaussian wave packets to retain their form while propagating in possibly anharmonic potentials ub(Q) + uνν(Q) and experiencing amplitude transfer from other levels as a result of uνν¯(Q), we enforce a locally quadratic approximation [32] under which the bath-coordinate-dependent quantities in Eq. (II.3) are written as second-order expansions about the centers Qν of the various bath wave packets. We must in particular express the summand in the position–space representation of Eq. (II.3) in terms of 〈Q|ψν(t)〉 by writing 〈Q|uνν¯(Qˆ)|ψν¯〉 = uνν¯(Q)〈Q|ψν¯〉〈Q|ψν〉〈Q|ψν〉, (II.5) and expanding uνν¯(Q)〈Q|ψν¯〉/〈Q|ψν〉 through second order in Q−Qν : uνν¯(Q) 〈Q|ψν¯〉 〈Q|ψν〉 ∼= eif(νν¯) [ uνν¯ + (Q−Qν) ( ∂uνν¯ ∂Q + ig1(νν¯)uνν¯ ) + (Q−Qν)2 ( 1 2 ∂2uνν¯ ∂Q2 + ig1(νν¯) ∂uνν¯ ∂Q + ig2(νν¯)uνν¯ )] . (II.6) 13 uνν¯ and its derivatives in Eq. (II.6) are to be evaluated at Q = Qν . The new quantities appearing there are f(νν¯) = (Qν −Qν¯)2αν¯ + (Qν −Qν¯)Pν¯ − γν + γν¯ , (II.7) g1(νν¯) = 2(Qν −Qν¯)αν¯ − Pν + Pν¯ , (II.8) and g2(νν¯) = i 2 g21(νν¯)− αν + αν¯ . (II.9) The approximation of Eq. (II.6) is justified by the expectation that none of the bath wave packets, which are only weakly shifted from the ground state, will differ greatly from the others. The ratio 〈Q|ψν¯〉/〈Q|ψν〉 is therefore nearly constant, and its weak spatial variation can be adequately captured by a quadratic expansion. In the illustrative case of a 1D bath, the parameter equations of motion are in this way found to be α˙ν = −2α2ν − 1 2 ∂2(ub + uνν) ∂Q2 − ∑ ν¯ 6=ν eif(νν¯) [ 1 2 ∂2uνν¯ ∂Q2 + ig1(νν¯) ∂uνν¯ ∂Q + ig2(νν¯)uνν¯ ] , (II.10) Q˙ν = Pν + 1 2α′′ν ∑ ν¯ 6=ν e−f ′′ (νν¯) [( ∂uνν¯ ∂Q − g′′1 (νν¯)uνν¯ ) sin f ′(νν¯) +g′1(νν¯)uνν¯ cos f ′(νν¯)] , (II.11) 14 P˙ν = −∂(ub + uνν) ∂Q + ∑ ν¯ 6=ν e−f ′′ (νν¯) [g′1(νν¯)uνν¯ sin f ′(νν¯) − ( ∂uνν¯ ∂Q − g′′1 (νν¯)uνν¯ ) cos f ′(νν¯) ] + α′ν α′′ν ∑ ν¯ 6=ν e−f ′′(νν¯) [( ∂uνν¯ ∂Q − g′′1 (νν¯)uνν¯ ) sin f ′(νν¯) + g′1(νν¯)uνν¯ cos f ′(νν¯)] , (II.12) γ˙′ν = −εν − α ′′ ν + P 2ν 2 − ub − uνν − ∑ ν¯ 6=ν e−f ′′ (νν¯)uνν¯ cos f ′(νν¯) + Pν 2α′′ν ∑ ν¯ 6=ν e−f ′′ (νν¯) [( ∂uνν¯ ∂Q − g′′1 (νν¯)uνν¯ ) sin f ′(νν¯) + g′1(νν¯)uνν¯ cos f ′(νν¯) ] , (II.13) and γ˙ ′′ ν = α ′ ν − ∑ ν¯ 6=ν e−f ′′ (νν¯)uνν¯ sin f ′(νν¯). (II.14) The equations of motion for the wave packet parameters of a multidimensional bath were published previously [17] but are compactly reproduced in Appendix A. II.2.2. Coupled Harmonic Oscillators To test the efficacy of the FVB/GB treatment of vibrational dynamics on the simplest possible two-mode system, we consider a pair of harmonic oscillators between which bilinear coupling is retained. Such a system is of course separable in normal modes, providing exact results with which to compare the FVB/GB predictions. In future applications to anharmonic multidimensional systems, the system coordinate will be identified with the highest-frequency normal mode, and as mentioned above, system–medium coupling will commence at higher order. By investigating a test case whose dynamics is readily soluble, either analytically or by 15 rudimentary basis-set methods, we are able to gauge both the ease of integration of the FVB/GB equations and the accuracy of their solutions. Our test Hamiltonian is of the form of Eq. (II.1) with us(q) = ω 2q2/2, ub(Q) = Ω2Q2/2, and u(q,Q) = JqQ. We set ω = 1.0 and Ω = 0.5 and compare the performance of FVB/GB with two different coupling strengths, J = 0.01 and 0.05. In each case, the initial state is a direct product of system- and bath-only ground states displaced along the system coordinate by δ = 1 = 21/2qrms |Ψ(0)〉 = D(δ)|ν = 0〉|n = 0〉 = ∞∑ ν=0 |ν〉 e−1/4 ( 1 2 )ν/2 1√ ν! |n = 0〉︸ ︷︷ ︸ |ψν(0)〉 , (II.15) with D(δ) = exp[−iδpˆ] and hb|n〉 = Ω(n + 1 2 )|n〉. The displacement generates nonzero amplitude in the bath wave packets accompanying higher-lying system vibrational levels while leaving their widths, spatial centers, and momenta identical to those of |n = 0〉. The initial condition in Eq. (II.15) is much like that following a short-pulse electronic transition in which the system vibration is Franck–Condon- active and the “medium” is not but the two modes experience nonvanishing coupling in the excited electronic state (section I.3.3). As noted above, the collection of bath wave packets |ψν(t)〉 defines the state (Eq. II.2) of the composite system. Under the Gaussian ansatz (Eq. II.4), the collection of parameters Qν , Pν , αν , and γν for the various vibrational indices ν specifies the state of the system plus the bath at a given time. To assess the validity of the FVB/GB treatment, we therefore compare these wave packet parameters with 16 the corresponding quantities Qν = 〈ψν |Qˆ|ψν〉 〈ψν |ψν〉 , (II.16) Pν = 〈ψν |Pˆ |ψν〉 〈ψν |ψν〉 , (II.17) α ′′ ν = 〈ψν |ψν〉 4〈ψν |(Qˆ−Qν)2|ψν〉 , (II.18) α′ν = α ′′ ν 〈ψν |[(Qˆ−Qν)(Pˆ − Pν) + (Pˆ − Pν)(Qˆ−Qν)]|ψν〉 〈ψν |ψν〉 , (II.19) γ′ν = −i ln [ 〈Q|ψν〉 |〈Q|ψν〉| ] − Pν(Q−Qν)− (Q−Qν)2α′ν , (II.20) γ ′′ ν = ln [ 1√〈ψν |ψν〉 ( pi 2α′′ν )1/4] , (II.21) calculated using the exact bath wave packets |ψν(t)〉 = 〈ν|Ψ(t)〉. II.2.3. Are the Bath Wave Packets Gaussian? A Gaussian wave packet remains Gaussian under evolution in an arbitrary quadratic potential. [33] Given its Gaussian initial form, the exact wave function 〈q,Q|Ψ(t)〉 with which we test our 2D model therefore remains Gaussian for all time, but this fact does not render the bath packets 〈Q|ψν(t)〉Gaussian nor imply that the FVB/GB description should be exact for the test model. Instead, a straightforward 17 analysis shows that if the exact 2D wave function is Ψ(q,Q; , t) = exp i(q − qt, Q−Qt) · A(t) ·  q − qt Q−Qt  + iBT (t) ·  q − qt Q−Qt + iC(t)  , (II.22) then the actual wave packets for the bath are given by 2 ψν(Q; t) = (piω)1/4√ 2νν!X ( X − ω X )ν/2 Hν ( Y √ ω 2 √ X(X − ω) ) exp ( Y 2 4X + Z ) , (II.23) where Hν(x) is a Hermite polynomial X = ω 2 − iAss, (II.24) Y = 2i(Q−Qt)Asb − 2iqtAss + iBs, (II.25) and Z = i(Q−Qt)2Abb − 2i(Q−Qt)qtAsb + i(Q−Qt)Bb +iq2tAss − iqtBs + iC. (II.26) Only the bath wave packet ψ0(Q; t) accompanying the vibrational ground state stays strictly Gaussian for all time, but those accompanying higher levels may remain approximately Gaussian. This occurs when the Hermite polynomial function in 2The integral ∫ ∞ −∞ dxe−ax 2+bxHn(x) = (pi a )1/2 [(a − 1)/a]n/2eb2/4aHn[b/(2(a(a − 1))1/2)] is helpful in this derivation. 18 Eq. (II.23) varies slowly in the range of Q values about which exp[Z + (Y 2/4X)] is peaked. II.2.4. Population Nonconservation The total population of the system-and-bath, 〈Ψ(t)|Ψ(t)〉 = ∑ ν 〈ψν(t)|ψν(t)〉, is the sum over the populations of the individual system vibrational levels. From the equations of motion (Eq. II.3) for the bath wave packets, it is easy to confirm that the rate of change of the total population d〈Ψ|Ψ〉 dt = ∑ ν ∑ ν¯ 6=ν rνν¯ , (II.27) vanishes because rνν¯ = −rν¯ν , where rνν¯ = 2Im〈ψν |〈ν|u|ν¯〉|ψν¯〉 is the rate of population transfer to ν from ν¯, resulting from the Hermiticity of the system–bath interaction potential. In the FVB/GB theory, this population-transfer balancing property is compromised because, as described in section I.2.1, uνν¯ and uν¯ν (and hence rνν¯ and rν¯ν) are subject to slightly different approximations. The explicit expression for the transfer rate under FVB/GB, rνν¯ ∼= √ 2pi α′′ν e−2γ ′′ ν Im { eif(νν¯) [uνν¯ + 1 4α′′ν ( 1 2 ∂2uνν¯ ∂Q2 + ig1(νν¯) ∂uνν¯ ∂Q + ig2(νν¯)uνν¯ )]} , (II.28) nonetheless makes it clear that population-transfer balancing should be approximately in force under the assumed conditions of weak system–bath interaction. In future applications of the FVB/GB theory to systems for which exact calculations are 19 impracticable, monitoring the extent of population conservation (as well as energy conservation [32]) can serve as an internal check of the accuracy of the approximate treatment. II.3. Numerical Results II.3.1. J = 0.01 We first consider a case of very weak system–bath coupling, J = 0.01, with a tensor product of the uncoupled system ground state displaced by δ = 1 (i.e., to the classical outer turning point) and the uncoupled bath ground state as the initial state (section I.2.2). The normal-mode frequencies ω± = ((ω2 + Ω2)/2 ± ((ω2 − Ω2)2/4 + J2)1/2)1/2 are barely distinguishable from ω and Ω for this value of the coupling constant. The FVB/GB parameters in Eq. (II.10) through (II.14) were propagated with a fourth-order Runge–Kutta algorithm using a fixed time step of δt = 1.6×10−6τs, where τs = 2pi/ω is the systems vibrational period. Eleven system vibrational states were included (ν = 0− 10). Plots of the spatial centers and average momenta of several bath wave packets for this very weak coupling case are shown in Figure II.1. As expected for a semiclassical treatment, the agreement between FVB/GB (red) and exact (blue) expectation values is excellent at short times but begins to deteriorate eventually, after hundreds of system vibrational periods. Notice that the bath- coordinate values move asymmetrically about Qν = 0 because of the asymmetrical 20 -0.08 -0.04 0 0.04 0.08 0 5 10 P 4 t/τs -0.03 -0.015 0 0.015 0.03 P 2 -0.02 -0.01 0 0.01 0.02 P 0 -0.08 -0.04 0 0.04 0.08 Q 4 -0.02 -0.01 0 0.01 0.02 Q 2 -0.04 -0.02 0 0.02 0.04 0.06 Q 0 700 705 710 t/τs Figure II.1. Time-dependent position and momentum expectation values for bath wave packets accompanying selected vibrational levels in the case of J = 0.01. Exact results are shown in blue, and FVB/GB values are plotted in red. Time is reckoned in multiples of the system vibrational period, τs = 2pi/ω. The mixed quantum/semiclassical theory behaves accurately at short times (left panels) but begins to fail—earlier for bath wave packets attached to higher vibrational states— after many hundreds of vibrational periods. initial condition. The real and imaginary parts of selected bath wave packet width parameters αν = α ′ ν + iα ′′ ν for the same case are plotted in Figure II.2. Exact values for these parameters are determined as explained in section I.2.2. Discrepancies between the exact and semiclassical values of the width parameters appear somewhat earlier than for the coordinate and momentum expectation values (Figure II.1.). In Figures II.1. and II.2. it can be seen that the breakdown 21 0.2497 0.25 0.2503 0.2506 0 5 10 α 4" t/τs 0.2498 0.25 0.2502 0.2504 α 2" 0.2498 0.2499 0.25 0.2501 0.2502 α 0" -0.0006 -0.0003 0 0.0003 0.0006 α 4’ -0.0002 0 0.0002 α 2’ -0.0002 -0.0001 0 0.0001 0.0002 α 0’ 300 305 310 t/τs Figure II.2. Same as for Figure II.1. for the real and imaginary parts of selected bath wave packet width parameters in the case of very weak system–bath coupling. Exact values of these parameters (blue) are compared with those predicted by the FVB/GB scheme (red). of the FVB/GB description occurs sooner for bath wave packets associated with higher system vibrational levels. The populations 〈ψν(t)|ψν(t)〉 = (pi/(2α′′ν))1/2 exp[−2γ ′′ ν ] (Eq. II.21) of a few system vibrational levels are graphed in Figure II.3. The fractional population changes are very small in this instance, but the mixed quantum/semiclassical approach correctly tracks them for several hundred vibrational periods. Not plotted here are the real parts γ′ν of the bath wave packet phase parameters; these match the exact values (given by Eq. (II.20) with Q = Qν) for about 1000 vibrational periods. 22 0.001578 0.00158 0.001582 0.001584 0.001586 0 5 10 po p 4 t/τs 0.07576 0.07578 0.0758 0.07582 0.07584 po p 2 0.6065 0.6066 0.6067 0.6068 po p 0 400 405 410 t/τs Figure II.3. Same as for Figure II.1. for the populations of system vibrational levels ν = 0, 2, and 4. In this very weak coupling case, the quantum/semiclassical theory (red) accurately portrays the exact transfer of population (blue) among different quantum levels for several hundred periods before starting to break down. The total population (the sum of 〈ψν(t)|ψν(t)〉 over all system levels ν) is very nearly conserved in this case but not exactly so (section I.2.4); it exceeds unity by 3×10−7 after about 800 periods. To gauge the overall accuracy of the FVB/GB treatment, one can evaluate the fidelity |〈Ψ(t)|ΨFVBGB(t)〉|/(〈ΨFVBGB(t)|ΨFVBGB(t)〉)1/2 of the combined time-dependent system–bath wave function given by Eq. (II.2); this quantity falls below unity by only 1× 10−7 after 400 vibrational periods. II.3.2. J = 0.05 Under more substantial but still fairly weak system–bath coupling, the approximate nature of the quantum/semiclassical treatment becomes apparent at shorter evolution times. This behavior is illustrated Figure II.4., which shows the time development of the position expectation value and the imaginary part of the width parameter for a few bath wave packets in a bilinearly coupled pair of harmonic oscillators with system frequency ω = 1, bath-mode frequency 23 0.21 0.25 0.29 0.33 0 5 10 15 20 α 4" t/τs 0.23 0.24 0.25 0.26 0.27 α 2" 0.246 0.248 0.25 0.252 0.254 α 0" -0.9 -0.6 -0.3 0 0.3 0.6 Q 4 -0.08 -0.04 0 0.04 0.08 Q 2 -0.1 0 0.1 0.2 0.3 Q 0 Figure II.4. Coordinate expectation values and imaginary width parameters for the ν = 0, 2, and 4 Gaussian bath wave packets in the J = 0.05 model (red) are compared to the corresponding exact values (blue). Ω = 0.5, coupling strength J = 0.05, and initial system displacement δ = 1. The normal-mode frequencies for this model are ω+ = 1.00166 and ω− = 0.49667. Numerical integration of the FVB/GB parameter equations of motion in this case was carried out using a fourth-order Runge–Kutta method with fixed time-step of δt = 1.6 × 10−7τs (the results were no different from those obtained with a longer time increment of δt = 1.6× 10−6τs). The truncated vibrational basis kept system states with ν from 0 to 10. A comparison of Figure II.4. to the evolution of the corresponding parameters in Figures II.1. and II.2. confirms that the same initial system displacement induces 24 0.9994 0.9996 0.9998 1 0 5 10 15 20 to ta l p op ul at io n t/τs 0.0014 0.0016 0.0018 po p 4 0.074 0.075 0.076 po p 2 0.606 0.609 0.612 0.615 po p 0 Figure II.5. Selected level populations and the total population of all system vibrational states under FVB/GB for the case of J = 0.05 (in red). Discernible differences from the exact state populations (blue) and significant departures from population conservation appear after about 10 vibrational periods. larger-amplitude motion in the bath in this more strongly coupled situation. Unsurprisingly, the parameter development predicted by the FVB/GB scheme (red) begins to disagree with the exact behavior (blue) at shorter times (on the order of tens of vibrational periods, again earlier for the parameters of wave packets attached to higher levels). The same pattern is followed in the populations of individual vibrational levels, the total population (Figure II.5.), and the fidelity fν = |〈ψν |ψFVBGBν 〉|/(〈ψν |ψν〉〈ψFVBGBν |ψFVBGBν 〉)1/2 of the bath wave packets (Figure II.6.). Because of the approximate nature of population 25 0.96 0.98 1 0 5 10 15 20 fid el ity 4 t/τs 0.998 0.9985 0.999 0.9995 1 fid el ity 2 0.9999992 0.9999996 1.0000000 fid el ity 0 Figure II.6. Fidelities of FVB/GB bath wave packets belonging to several system vibrational states with system–bath coupling of J = 0.05. conservation under FVB/GB, the unphysical deviations of the total population from unity seen in the bottom panel of Figure II.5. are not necessarily attributable to numerical integration error. The expanded vertical scale of that plot, along with the very small losses in fidelity of the individual bath wave packets shown in Figure II.6., suggest that the practical utility of the quantum/semiclassical approximation could persist beyond the initial appearance of discernible errors in the highest-lying bath wave packets. It is worth noting that, by truncating the vibrational basis at a particular level, the FVB/GB treatment as implemented makes a second approximation in addition to the semiclassical assumption of Gaussian bath wave packets. Although this secondary approximation seems minor because of the very small bath amplitude associated with the highest levels in the finite vibrational basis, there is reason to consider the possibility that it may contribute significantly to the breakdown of the FVB/GB theory. Figure II.7. addresses this question by comparing the 26 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 5 10 15 20 25 30 35 Q 1 t/τs FVB/GB, νm=10FVB/GB, νm=20 Exact Figure II.7. Expectation value of the position coordinate of the bath wave packet attached to the ν = 1 level for J = 0.05. The semiclassical FVB/GB prediction follows the exact value (blue) for a longer time when 21 states are included in the truncated vibrational basis (red) than when only 11 states are retained (green). expectation value of the position coordinate of the wave packet associated with ν = 1 calculated for bases of 11 and 21 vibrational states. The FVB/GB prediction follows the exact value (blue) for a longer time when 21 states are included in the truncated vibrational basis (red) than when only 11 states are retained (green). 3 The solutions of the parameter equations of motion Eqs. (II.10)–(II.14) appearing in Figure II.7. were obtained using fourth-order Runge–Kutta integration with a fixed time step of δt = 1.6 × 10−7τs and appear to be numerically converged by virtue of their agreement with solutions obtained using a 10-times-longer step. The greater longevity of FVB/GB with a larger basis seen in Figure II.7. is at least consistent with the apparent “breakdown propagation” from higher- to lower- level bath wave packets noted in the preceding Figures. Cutting off the vibrational basis at a certain νmax is equivalent to replacing the actual bilinear system–bath 3Exploratory calculations with an even larger vibrational basis of 41 states show still-longer- term agreement between FVB/GB and numerically exact bath wave-packet parameters. 27 coupling with a more complicated interaction potential that effectively turns off the transfer of amplitude between ψνmax(Q) and ψνmax+1(Q) by the ν¯ = νmax +1 term in Eq. (II.3). A putatively greater sensitivity on the part of FVB/GB to such a change in u(q,Q) than the numerically rigorous calculations based on diagonalization of the Hamiltonian in a finite basis could be due to the constrained Gaussian form of the bath wave packets under FVB/GB. Alternatively, the mixed quantum/semiclassical theory could simply be more sensitive to round-off error in its handling of the barely populated higher-lying levels. Because q extends over a larger range of values in the higher-lying system vibrational states, strengthening the interaction of higher vibrational states with the bath, it cannot be concluded, however, that the propagation of truncation errors alone determines the survival time of the quantum/semiclassical description (loosely defined as the approximate propagation time before which the bath wave packets of the FVB/GB theory begin to differ significantly from their exact counterparts or otherwise to exhibit nonphysical behavior). Despite the demonstration in previous section I.2.3 that the exact bath wave packets do not remain strictly Gaussian (except for ψ0(Q; t)), one might wonder whether, but for “edge effects” associated with round-off error or vibrational basis truncation, some of them would essentially retain this form for substantial periods of time. We do not attempt to analyze this question exhaustively at present, but Figure II.8. illustrates the fact that even the bath wave packets accompanying low 28 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -8 -6 -4 -2 0 2 4 6 8 a bs ol ut e va lu e, ν = 1 Q Hermite Polynomial Gaussian function wave packet Figure II.8. Magnitude of the Gaussian portion (blue) and the Hermite polynomial prefactor (red) of the exact bath wave packet 〈Q|ψ1(60τs)〉 (purple) for the case of J = 0.05. The FVB/GB ansatz cannot exactly describe this bath wave packet because it assumes a constant prefactor. vibrational levels can have significant non-Gaussian contributions. This Figure graphs the magnitude of the purely Gaussian portion (blue) and the Hermite- polynomial prefactor (red) of the exact bath wave packet associated with the first vibrational excited state (ψ1(Q; t); see Eq. II.23) at t = 60τs (purple) for the case of J = 0.05. The non-negligible variation in the magnitude of the H1 factor over the spatial width of the wave packet indicates that the Gaussian ansatz of FVB/GB can become discernibly approximate at times similar to the breakdown times noted in Figures II.4.–II.8. II.3.3. Electronic Absorption Spectrum One of the primary motivations for the development and testing of our mixed quantum/semiclassical theory is eventually to enable the rigorous calculation of time-resolved nonlinear optical signals from small molecules embedded in low- temperature crystalline matrices. As a simple initial test of this program, we 29 can use the FVB/GB framework to determine the linear absorption spectrum or, equivalently, the short-pulse linear wave packet interferometry signal [34, 35] for a solvated “molecule” with two electronic levels, g and e, a Hamiltonian H = |g〉hg〈g|+ |e〉he〈e| (II.29) and an electronic dipole moment operator µˆ = µ(|e〉〈g| + |g〉〈e|). The ground- state nuclear Hamiltonian hg is that for the uncoupled intramolecular (system) and medium (bath) degrees of freedom, given by Eq. (II.1) with us(q) = ω 2q2/2, ub(Q) = Ω 2Q2/2, and u(q,Q) = 0. The excited-state nuclear Hamiltonian he has the same form with us(q) = (ω 2/2)(q+ 31/2)2 + εe, ub(Q) = (Ω 2/2)(Q+ 31/2)2, and u(q,Q) = J(q + 31/2)(Q + 31/2) as a result of the Franck–Condon activity in both the molecule and the medium and a Duschinsky rotation between them. We set ω = 1, Ω = 1/3, and J = 1/20. The vibronic absorption spectrum is the stick spectrum in Figure II.9. obtained from I(ξ) = µ2 ∑ n,N |〈(n,N)e|(0, 0)g〉|2δ(ξ−E(n,N)e +E(0,0)g) using the exact vibronic states (|g〉|(n,N)g〉 and |e〉|(n,N)e〉) and the eigenenergies (E(n,N)g = ω(n+ 1/2) + Ω(N + 1/2) and E(n,N)e = ω+(n + 1/2) + ω−(N + 1/2) + εe with ω± = ((ω 2 + Ω2)/2± ((ω2−Ω2)2/4 + J2)1/2)1/2). The FVB/GB approximation to the spectrum was calculated by using the standard expression [36, 37] I(ξ) = µ2 pi Re ∫ ∞ 0 dteit[ξ+(ω/2)+(Ω/2)]〈(0, 0)g|e−i[he−(iΓ/2)]t|(0, 0)g〉, (II.30) and evaluating the requisite time-dependent wave packet overlap by means of our 30 0 0.05 0.1 0.15 0.2 0.25 0 1 2 3 4 5 6 7 in te ns ity frequency/ω FVB/GB Exact Figure II.9. Electronic absorption spectrum of a two-mode model calculated exactly (blue) and using FVB/GB wave packet propagation (red). The frequency axis is down-shifted by the bare electronic transition frequency εe. See the text for details. quantum/semiclassical propagation algorithm4 along with analytical expressions for the resulting integrals over the bath coordinate variable. In Eq. (II.30), we have introduced an ad hoc “radiative” decay rate of Γ = 0.15/τs to zero out the integrand by the time the FVB/GB approximation begins to fail; this artificial decay, or the limited valid duration of FVB/GB that requires its introduction, gives rise to the nonzero line widths seen in Figure II.9. Apart from this difference, there is excellent agreement between the exact and approximate absorption spectra. We emphasize that in future applications of FVB/GB to systems embedded in multimode media it is to be anticipated that the theory itself will properly account for optical dephasing due to a permanent loss of overlap between an original 4Numerical integration of the FVB/GB equations of motion was performed with a Runge– Kutta routine using a fixed time step, δt = 1.6 × 10−6τs, for a total integration time of 63.7τs . Time-dependent parameters were found for the bath wave packets accompanying 11 vibrational states, but only 10 (ν = 0 − 9) were included in the spectrum. The population of the excluded ν = 10 state is 3.5× 10−6. 31 ground-state wave packet and a copy of it propagated in an excited electronic state. This process will typically occur on a much shorter timescale (tens to hundreds of femtoseconds) than spontaneous emission and will also precede the breakdown of the semiclassical approximation. II.4. Comments and Conclusions The analysis and numerical tests of the mixed quantum/semiclassical FVB/GB method given here augur well for its successful application to larger, more complex molecular problems. We have found that in an appropriately chosen test model and for significant lengths of time it does a good job of tracking the bath-state parameters that define the overall state of the system and bath, according to Eq. (II.2). The calculation of section I.3.3 illustrates FVB/GB’s ability to simulate spectroscopic signals. Initial realistic applications can be envisaged with respect to the short-time dynamics of van der Waals complexes of dihalogen molecules in small rare-gas clusters (the longer process of cluster evaporation due to vibrational predissociation would not, however, be directly accessible to the theory). It is noteworthy that the separation between the system and bath would be effected by using the normal coordinates in a chosen electronic state. In that electronic manifold at least, system–bath coupling would begin at cubic order, which is one degree higher and 32 hence likely even weaker than the bilinear coupling considered in the present trial simulations. Halogen diatomics and other species in rare-gas hosts or clusters have been the subject of numerous previous theoretical studies. Of particular interest in comparison to the FVB/GB scheme is another, more general semiclassical framework, the forward–backward implementation of a semiclassical initial value representation (FB-IVR) [38, 39] that was recently used to calculate some dynamical properties of I2 molecules interacting with one or six Ar atoms. [40, 41] In that study, Tao and Miller used FB-IVR to calculate the probability density of iodine internuclear distances (expressed as a time-correlation function) at various times following the initial displacement of a vibrational wave packet. They found that the theory, which treated all degrees of freedom semiclassically, was able to generate a fringe structure in the probability density characteristic of quantum mechanical interference. The FVB/GB theory differs from FB-IVR in treating selected high-frequency degrees of freedom rigorously. In its present form, FVB/GB is more specialized in that it assumes an initial pure state that may be described as a tensor product of system and bath states and is targeted at a narrower class of molecular problems. The advantages of the present theory are its avoidance of a time-consuming phase- space average over the initial states and the fact that it generates a fully entangled state-ket for both the system and the medium in which it is embedded. This 33 approximate time-dependent wave function can in principle be used, together with analytical expressions for multidimensional Gaussian integrals, to evaluate arbitrary system or bath observables, correlation functions, or spectroscopic signals without recalculation. Our FVB/GB approach has both similar and complementary features to a local coherent-state approximation (LCSA) put forward by Martinazzo and co-workers. [42] That paper adopts an expansion in discrete position-variable states of a quantum mechanical subsystem and treats the bath as a product of coherent states (frozen Gaussians) instead of the variably squeezed state (arbitrary multidimensional Gaussian) used here. The LCSA ansatz successfully captured the subsystem dynamics of several system–bath combinations including a Morse- oscillator subsystem coupled to a multimode harmonic bath of various frequencies and coupling properties. Its use of a position-representation expansion of the subsystem’s time-dependent state suggests that the approximation may be best suited to problems in which the subsystem is of a similar or lower frequency than the surrounding medium. Mintert and Heller [21] recently outlined a general semiclassical approach to the simulation of quantum systems that can exchange both energy and particles with the surroundings. The theory sketched in that paper envisages a quantum mechanical system and a quantum mechanical bath occupying different spatial regions of an extended sample and describes the states of both using 34 multidimensional Gaussian wave packets. Among other interesting features, ref [21] gives consideration to the decomposition of a mixed-state system density matrix as either an incoherent sum or a coherent superposition of Gaussian states. In future theoretical work, it will be interesting to investigate whether the finite survival time of the semiclassical treatment based on Gaussian wave packets observed in section I.3 can be extended without sacrificing the advantage of power- law scaling with respect to the dimension of the bath. One possible avenue would be to generalize the bath wave packets to include components proportional to Gaussian functions times low-order Hermite polynomials. [43] Preliminary analysis indicates that such an enhancement would allow the theory to handle pure Duschinsky rotations, for which FVB/GB in its native form fails.5 Consideration could also be given to the use of bath-parameter equations of motion derived from a variational principle rather than a quadratic expansion of the nuclear Hamiltonian about the center of the wave packet. [44, 45] A prospective technique for intermolecular communication by coherent acoustic phonons [16] depends crucially on the generation of a nonequilibrium bath state following the vibrational or electronic excitation of one guest species (the transmitter) and the influence of the resulting propagating disturbance on another guest molecule (the receiver). Further refinement and testing of this process stand 5Consistently with FVB/GB’s failure under abrupt initiation of bilinear coupling between harmonic system and bath modes in the absence of vibrational displacement (i.e. pure Duschinsky rotation), calculations show that the theory in the form presented here breaks down rather quickly in the presence of both bilinear coupling and only a very small vibrational displacement. 35 to benefit from the use of FVB/GB to describe the launching and motion of lattice waves in the medium surrounding the initially excited chromophore, their subsequent impingement on the vibrational dynamics of the receiver species, and the detection of the latter, perhaps by time-resolved coherent anti-Stokes Raman spectroscopy. [5] Interesting examples of coherent nuclear motion in a surrounding medium initiated by the short-pulse electronic excitation of an embedded target molecule have been recorded in ultrafast pump–probe spectra of various dihalogens in solid matrices. [10, 13] In these experiments, sinusoidal oscillations in the pump–probe signal at frequencies characteristic of the zone-boundary phonons of the host lattice were visible long after vibrational-frequency oscillations had decayed away. Two general mechanisms of excitation of these local lattice modes may be considered: (a) coupling between nuclear degrees of freedom of the system and bath that differs in the ground and excited electronic states (analogous to the J term given below Eq. II.29) and (b) direct electron–phonon coupling, referred to as the displacive excitation of coherent phonons (DECP, akin to the inclusion of δQ 6= 0 in ub(Q) below Eq. II.29). Chapman [46] recently performed detailed comparisons of these alternative mechanisms using FVB/GB in calculations similar to those in section I.3.3. As the example in that section illustrates, these two mechanisms need not be mutually exclusive. Another purely nuclear excitation mechanism for zone-boundary phonons in these systems, along the lines of mechanism (a), 36 would not require a shift in vibration–phonon coupling upon electronic excitation of the guest molecule. It could simply be the case that the large-amplitude vibrational motion launched by short-pulse electronic absorption activates the higher-order terms comprising u(q,Q) (Eq. II.1), which may be unimportant unless the vibrational coordinate q takes large values. Detailed dynamical and electronic- structure calculations are needed to resolve fully the origin of zone-boundary phonon excitation upon small-guest electronic absorption in rare-gas matrices. Intriguing recent experiments from the Schwentner group [3, 12, 14] make use of wave packet interferometry to exert optical control over the interactions between the vibronic excitation of a bromine molecule and the lattice vibrations of the cryogenic krypton crystal in which it was embedded. By varying the optical fringe structure of the exciting phase-locked pulse pairs to concentrate spectral intensity on either the pure vibronic progression of bromine (or another dihalogen species) or the progression of accompanying phonon side bands, those experimenters were able to prepare and probe different types of nonstationary states in which the molecular vibration was or was not accompanied by a superposition of one-quantum states within a finite range of the acoustic-phonon bandwidth. Experiments of this kind are prime candidates for simulation by FVB/GB. The van Hulst research group has pioneered the application of wave packet interference methods to single molecules isolated in crystalline media. Most recently, that group reported the first successful fluorescence-detected single- 37 molecule nonlinear optical spectroscopy measurements.6 [47] The theoretical simulation of further measurements along this line, perhaps including multidimensional wave packet interferometry and state-reconstruction experiments [48, 49] on solvated chromophores, would evidently require methods possessing the twin capability of FVB/GB to provide both a rigorous quantum mechanical description of intramolecular degrees of freedom and a well-defined approximate treatment of the associated medium dynamics. 6A collaboration of groups at the universities of Wu¨rzburg, Bielefeld and Kaiserslautern recently performed 2D coherent electronic spectroscopy measurements of individual defects on a metal surface by detecting an action-signal based on spatially resolved photo-electron emission, rather than fluorescence (T. Brixner). 38 CHAPTER III VARIATIONAL FVB/GB AND MODEL OF MOLECULAR IODINE IN A 2D KRYPTON LATTICE One of the problems occurring in FVB/GB with the locally quadratic approximation, which is the total population being not conserved, can be solved by adopting a variational treatment. More generally, variational FVB/GB may have improved accuracy and numerical stability. In this chapter, we start from the general time-dependent variational principle and derive equations of motion for the bath wave-packet parameters. To apply variational FVB/GB to a realistic system, we establish a model of molecular iodine in a 2D krypton lattice. Computations performed on this 2D model consume about seven times less time than a 3D model would require, and yield relatively efficient and reasonable simulation results. III.1. Variational FVB/GB III.1.1. Time-Dependent Variational Principle in Quantum Mechanics The time-dependent variational principle is formulated by minimizing the action function given by S = ∫ t2 t1 Ldt, (III.1) 39 where the Lagrangian is taken as L = 〈Ψ|i ∂ˆ ∂t −H|Ψ〉. (III.2) (We set ~ = 1.) Here the time derivative operator includes two parts acting on the right or left side respectively, [18] i ∂ˆ ∂t = − i 2 (←− ∂ ∂t − −→ ∂ ∂t ) . (III.3) The variation of the action is δS = ∫ t2 t1 δLdt = ∫ t2 t1 ( i 2 〈δΨ|Ψ˙〉+ i 2 〈Ψ|δΨ˙〉 − i 2 〈Ψ˙|δΨ〉 − i 2 〈δΨ˙|Ψ〉 −〈δΨ|H|Ψ〉 − 〈Ψ|H|δΨ〉) dt, (III.4) with 〈Ψ|δΨ˙〉 = d dt 〈Ψ|δΨ〉 − 〈Ψ˙|δΨ〉, (III.5) and its complex conjugate 〈δΨ˙|Ψ〉 = d dt 〈δΨ|Ψ〉 − 〈δΨ|Ψ˙〉. (III.6) We can integrate Eq. (III.4) by parts to obtain δS = [ − i 2 〈δΨ|Ψ〉∣∣t2 t1 + ∫ t2 t1 ( i〈δΨ|Ψ˙〉 − 〈δΨ|H|Ψ〉 ) dt ] + c.c. (III.7) |δΨ〉 is small everywhere in the interval of time from t1 to t2 and vanishes at both t1 and t2. The integrated term in Eq. (III.7) is zero. The action is required to be 40 stationary with respect to |δΨ〉 between the fixed points t1 and t2, and therefore for an arbitrary |δΨ〉, we have Re [ i〈δΨ|Ψ˙〉 − 〈δΨ|H|Ψ〉 ] = 0. (III.8) Let the real and imaginary parts of 〈δΨ|x〉 be a(x) and b(x) respectively, and 〈x|i ∂ ∂t −H|Ψ〉 ≡ c(x) + id(x). Then Eq. (III.8) becomes∫ ∞ −∞ dxRe {[a(x) + ib(x)][c(x) + id(x)]} = ∫ ∞ −∞ dx [a(x)c(x)− b(x)d(x)] = 0. (III.9) If a(x) and b(x) are considered to be independent and completely arbitrary, then we have c(x) = d(x) = 0, and i|Ψ˙〉 = H|Ψ〉, (III.10) which is the Schro¨dinger equation. Other forms of the Lagrangian have been used too. Kramer and Saraceno do not assume the wave function is normalized and take the Lagrangian as [50] L = i 2 〈Ψ|Ψ˙〉 − 〈Ψ˙|Ψ〉 〈Ψ|Ψ〉 − 〈Ψ|H|Ψ〉 〈Ψ|Ψ〉 ; (III.11) Kerman and Koonin use the form of a complex Lagrangian [51] L = ∫ (dr)AΨ∗[i(∂/∂t)−H]Ψ; (III.12) Heller minimizes the functional I to formulate the variational principle [19] I = ∫ |HΨ− i∂Ψ/∂t|2dv, (III.13) 41 where dv implies integration over all coordinate space. These forms of variational principle all result in the Schro¨dinger equation, which is one of the fundamental postulates of quantum mechanics. III.1.2. Equations of Variational FVB/GB To derive the equations of motion for the bath wave-packet parameters, we expand the time-dependent nuclear ket as a sum of tensor product states, |Ψ(t)〉 = ∑ ν |ν〉〈ν|Ψ(t)〉 = ∑ ν |ν〉|ψν(t)〉, (III.14) where |ν〉 is the system eigen state, and make a Gaussian ansatz for the bath wave packets, as in Chapter II, 〈Q|ψν〉 = exp[i(Q−Qν)T · αν · (Q−Qν) + iP Tν · (Q−Qν) + iγν ]. (III.15) Then we have the Lagrangian in terms of wave-packet parameters L = ∑ ν 〈ψν |ψν〉 ( −1 4 Tr[α˙′ν · (α′′ν)−1] + P Tν · Q˙ν − γ˙′ν ) − 〈Ψ|H|Ψ〉, (III.16) where the expectation value of energy is 〈Ψ|H|Ψ〉 = ∑ ν 〈ψν |ψν〉 ( ν + 1 2 P Tν · Pν + 1 2 Tr[α′2ν · (α′′ν)−1 + α′′ν ] ) + ∑ ν 〈ψν |ub + uνν |ψν〉+ ∑ ν,ν¯ 6=ν 〈ψν |uνν¯ |ψν¯〉, (III.17) and the population of the νth vibrational level is 〈ψν |ψν〉 = e−2γ′′ν √ piN 2Ndet[α′′ν ] . (III.18) 42 The variational principle leads to the Lagrange’s equations of motion, one of which is d dt ∂L ∂γ˙′ν = ∂L ∂γ′ν , (III.19) and it gives the time derivative of the νth level population d dt 〈ψν |ψν〉 = − ∑ ν¯ 6=ν (i〈ψν |uνν¯ |ψν¯〉+ c.c) . (III.20) Along with this equation, the Lagrange’s equations for γ′′ν , Qν , Pν , α ′ ν , and α ′′ ν lead to the equations of motion for all the wave-packet parameters: P˙ν = Pν 〈ψν |ψν〉 ∑ ν¯ 6=ν (i〈ψν |uνν¯ |ψν¯〉+ c.c.) − 1〈ψν |ψν〉 ∂ ∂QTν [ 〈ψν |ub + uνν |ψν〉+ ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) ] ,(III.21) Q˙ν = Pν + 1 〈ψν |ψν〉 ∂ ∂P Tν ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) , (III.22) α˙′ν = 2(α ′′ ν) 2 − 2(α′ν)2 + 2α ′′ ν 〈ψν |ψν〉 [ 〈ψν |ub + uνν |ψν〉+ 1 2 ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) ] + 4α ′′ ν 〈ψν |ψν〉 · ∂ ∂α′′ν [ 〈ψν |ub + uνν |ψν〉+ ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) ] · α′′ν , (III.23) α˙ ′′ ν = −2α ′′ ν · α′ν − 2α′ν · α ′′ ν − α ′′ ν 〈ψν |ψν〉 ∑ ν¯ 6=ν (i〈ψν |uνν¯ |ψν¯〉+ c.c.) − 4〈ψν |ψν〉α ′′ ν · [ ∂ ∂α′ν ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) ] · α′′ν , (III.24) 43 γ˙′ν = −εν + 1 2 P Tν · Pν − Tr[α ′′ ν ] − Nb + 2 2〈ψν |ψν〉 [ 〈ψν |ub + uνν |ψν〉+ 1 2 ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) ] − 1〈ψν |ψν〉Tr [ α ′′ ν · ∂ ∂α′′ν ( 〈ψν |ub + uνν |ψν〉+ ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) )] + 1 〈ψν |ψν〉P T ν · ∂ ∂P Tν ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) , (III.25) γ˙ ′′ ν = Tr[α ′ ν ] + Nb + 2 4〈ψν |ψν〉 ∑ ν¯ 6=ν (i〈ψν |uνν¯ |ψν¯〉+ c.c.) + 1 〈ψν |ψν〉Tr [ α ′′ ν · ∂ ∂α′ν ∑ ν¯ 6=ν (〈ψν |uνν¯ |ψν¯〉+ c.c.) ] . (III.26) Nb is the number of bath modes. Q and P are vectors of Nb elements and α is an Nb by Nb matrix. It is worth mentioning that the total population is conserved in variational FVB/GB, because d dt ∑ ν 〈ψν |ψν〉 = ∑ ν,ν¯ 6=ν (−i〈ψν |uνν¯ |ψν¯〉+ i〈ψν¯ |uν¯ν |ψν〉) = ∑ ν,ν¯ 6=ν −i〈ψν |uνν¯ |ψν¯〉+ ∑ ν¯,ν 6=ν¯ i〈ψν |uνν¯ |ψν¯〉 = 0. (III.27) We switched indices ν and ν¯ for the second term. It also can be proved that the total energy is conserved. Let T ≡ 〈Ψ|i ∂ˆ ∂t |Ψ〉 = ∑ ν 〈ψν |ψν〉 ( −1 4 Tr[α˙′ν · (α′′ν)−1] + P Tν · Q˙ν − γ˙′ν ) , (III.28) and E ≡ 〈Ψ|H|Ψ〉. Then L = T −E. T is a linear function of α˙′ν , Q˙ν and γ˙′ν ; it also involves α ′′ ν , Pν and γ ′′ ν . For convenience of notation, we separate the parameters 44 into two groups by defining: χ1 = α ′ ν , Qν , γ ′ ν , χ2 = α ′′ ν , Pν , γ ′′ ν . (III.29) T is a function of χ˙1 and χ2, but not χ1 or χ˙2, with ∂T ∂χ1 = 0 and ∂T ∂χ˙2 = 0. The linearity of T with χ˙1 is crucial in this derivation, and enables us to simplify the following summation ∑ χ1 ∂T ∂χ˙1 χ˙1 = T. (III.30) The Lagrange’s equations lead to d dt ∂T ∂χ˙1 = − ∂E ∂χ1 , (III.31) and 0 = ∂T ∂χ2 − ∂E ∂χ2 . (III.32) We have ∂E ∂χ1 = − d dt ∂T ∂χ˙1 = − ∑ χ2 χ˙2 ∂2T ∂χ2∂χ˙1 . (III.33) The chain rule was used in the last step. The time derivative of the energy can be 45 shown to vanish by making use of the linearity of T with χ˙1: d dt E = ∑ χ1 χ˙1 ∂E ∂χ1 + ∑ χ2 χ˙2 ∂E ∂χ2 = − ∑ χ1χ2 χ˙1χ˙2 ∂2T ∂χ2∂χ˙1 + ∑ χ2 χ˙2 ∂T ∂χ2 = − ∑ χ2 χ˙2 ∂T ∂χ2 + ∑ χ2 χ˙2 ∂T ∂χ2 = 0. (III.34) III.1.3. Comparison of FVB/GB and Variational FVB/GB In the FVB/GB treatment, we started from the time-dependent Schro¨dinger equation (see Chapter II). To derive the equations for the time derivatives of the wave-packet parameters, we made locally quadratic approximations for both the bath and interaction potentials, and also expanded uνν¯(Q)〈Q|ψν¯〉/〈Q|ψν〉 through second order in (Q−Qν). In the variational FVB/GB approach, on the other hand, we use Lagrange’s equations to determine the behavior of the Gaussian wave-packet parameters. The equations of motion for the wave-packet parameters obtained by applying the variational principle are different from the previous equations of FVB/GB with a locally quadratic approximation. In the old version of FVB/GB, we evaluated the system-bath interaction and its derivatives at the wave-packet center position Qν , uνν¯(Qν), ∇uνν¯(Qν), ∇2uνν¯(Qν). The symbol ∇ represents a vector of derivatives with respect to the bath-coordinate 46 variables, ∇ = ∂ ∂Q . In the variational treatment, we evaluate the elements of the interaction operator matrix and its derivatives with respect to the center position Qν , momentum Pν , and the width parameter αν , 〈ψν |uνν¯ |ψν¯〉, ∂ ∂QTν 〈ψν |uνν¯ |ψν¯〉, ∂ ∂P Tν 〈ψν |uνν¯ |ψν¯〉, ∂ ∂α′ν 〈ψν |uνν¯ |ψν¯〉, ∂ ∂α′′ν 〈ψν |uνν¯ |ψν¯〉. The major improvement of the variational method is exhibited in the high- order coupling cases. For instance, for two coupled harmonic oscillators, when the interaction is of the form u = JqQ3 for 1D bath, and we have uνν¯ = JQ 3 1√ 2ω (√ νδν¯,ν−1 + √ ν + 1δν¯,ν+1 ) , (III.35) ∇uνν¯ = 3JQ2 1√ 2ω (√ νδν¯,ν−1 + √ ν + 1δν¯,ν+1 ) , (III.36) ∇2uνν¯ = 6JQ 1√ 2ω (√ νδν¯,ν−1 + √ ν + 1δν¯,ν+1 ) . (III.37) If we start with bath wave packets centered at zero, Qν = 0, then the right-hand side of each of the Eqs. (III.35)– (III.37) vanishes initially and the coupling terms in the FVB/GB equations of motion are all zero. With this anharmonic system-bath interaction, all the parameters keep their initial values and the wave packets do not evolve at all under the old FVB/GB scheme. Given this initial condition, FVB/GB with a locally quadratic approximation fails to propagate the bath wave packets. However, with the same interaction potential and the same initial condition, the variational treatment generates bath wave-packet propagation and can provide the 47 time-dependent total wave function. With Qν = 0, Pν = 0, and α ′ ν = 0 initially for all vibrational levels, the only nonzero coupling term is ∂ ∂QTν 〈ψν |uνν¯ |ψν¯〉 = 3α ′′ ν 2(α′′ν + α ′′¯ ν)2 √ pi α′′ν + α ′′¯ ν e−γ ′′ ν−γ′′¯ν × J√ 2ω (√ νδν¯,ν−1 + √ ν + 1δν¯,ν+1 ) . (III.38) Substituting this expression into the equation of motion for P˙ν , we can have the wave packets evolving with time. From this example, it is obvious that variational FVB/GB is more capable to deal with high order coupling cases than the old version of FVB/GB. We test the variational FVB/GB for an illustrative case of two coupled harmonic oscillators with the Hamiltonian H = p2 2 + 1 2 ω2q2 + P 2 2 + 1 2 Ω2Q2 + J2qQ 2, (III.39) where the system and bath frequencies are ω = 1 and Ω = 0.5, and the coupling coefficient is J2 = 0.0159 (such that J2qrmsQrms = 1/(14τs)). Initially the system coordinate is displaced by 1/ √ ω. FVB/GB and variational FVB/GB yield fairly similar results at short times. But under variational FVB/GB, the total population is conserved; under FVB/GB, it exceeds one at longer time (Figure III.1.). Variational FVB/GB also shows better performance by following the exact results longer, which can be seen from the fidelity (Figure III.2.). As defined in section II.3.1, the fidelity is fidelity = |〈Ψ(t)|ΨFVBGB(t)〉|/(〈ΨFVBGB(t)|ΨFVBGB(t)〉)1/2, (III.40) 48 1 1.005 1.01 1.015 1.02 1.025 0 1 2 3 4 5 to ta l p op ul at io n t/τs FVB/GB Variational Figure III.1. The total population of all system vibrational states in FVB/GB and in variational FVB/GB with Hamiltonian expressed in Eq. (III.39). 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 0 2 4 6 8 10 Fi de lity t/τs FVB/GB Variational Figure III.2. Fidelity of the overall wave function in FVB/GB and in variational FVB/GB with Hamiltonian expressed in Eq. (III.39). where |Ψ(t)〉 refers to the exact state ket calculated using a basis-set method. III.2. The Model of Molecular Iodine in a 2D Kr Lattice III.2.1. Model System The model for our calculation is a 2D lattice with 25 atoms: one iodine molecule and 23 krypton atoms. The reason we choose this number of atoms will be explained in detail in the following paragraphs. The total potential is taken to be a sum of 49 Parameter Values De 12547.194 cm −1 re 2.666 A˚ a0 4.99203 a1 17.6987 a2 52.1267 a3 −5.93197 a4 8.58628 A 8.57701× 107 cm−1 B 3.26046 A˚ −1 C 4.32949× 105 cm−1A˚6 D 1.43582× 107 cm−1A˚8 σ 3.74 A˚ ε 162.3 cm−1 Table III.1. Potential parameters pair-wise atom-atom interactions. The I-I interaction is represented by a modified Hulbert-Hirschfelder potential (X state),1 [16] VI−I(r) = De[(1 + a1ξ3 + a2ξ4)e−2a0ξ − 2(1 + a3ξ3 + a4ξ4)e−a0ξ], (III.41) with ξ = (r/re−1). For Kr-Kr, we use a Buckingham-type function of the form [52] VKr−Kr(r) = Ae−Br − C/r6 −D/r8. (III.42) The I-Kr interaction is approximated by the Xe-Kr Lennard-Jones potential, [53] VI−Kr(r) = 4ε ( σ12 r12 − σ 6 r6 ) . (III.43) Parameters for these potentials are listed in Table III.1. 1The I2 potential of the B state is also used in the simulations of ultrafast optical signals. We make no change in the I-Kr potential upon electronic excitation of I2. 50 People may wonder why we use 25 atoms. On the one hand, we want to include all the interactions between atoms in the bulk solid that are not negligible; on the other hand, in order to reduce the amount of computing time, we will wish to use the smallest possible number of atoms. Periodic boundary conditions are essential in order to mimic the presence of the infinite bulk surrounding our model system. A square box containing N atoms is treated as a unit of an infinite periodic lattice of identical boxes. By wisely choosing the number of atoms N and the size of the box, we only need count the interactions among atoms inside of the box. It is assumed that for r greater than a certain cut-off distance rc, any pair-potential and its first and second derivatives are negligible. The size of the box L should be larger than twice rc. For a particular atom i in the box, if its distance from an atom j (outside of the box) rij is smaller than rc, then the distance between atom i and the image of atom j, which is inside of the box, will be L−rij > rc. We switch the interaction on atom i from atom j’s image with that from j, calculate the interactions of two atoms with a distance smaller than rc, and neglect those with a distance larger than rc. Therefore, we always only need to count the atom j or its image, never both. The value of this cut-off distance rc is chosen to be the distance at which the Kr-Kr interaction energy has magnitude equal to one-tenth of the vibrational quantum ~ωKr: |VKr−Kr(rc)| = 0.1~ωKr = 0.1~ √ (∂2VKr−Kr/∂r2)eq mKr/2 = 2.32 cm−1, (III.44) from which we find rc = 8.095A˚. At this cutoff distance, both I-I and I-Kr 51 Figure III.3. The 2D I2/Kr model. The red dots represent iodine atoms and the blue are krypton atoms. Scale bar = 1 A˚. interactions have magnitudes smaller than one tenth of the dissociation energy. The length of the box side is an integer times the square root of the area per atom, L = n √ area per atom > 2rc. (III.45) The primitive cell for a 2D krypton lattice is a parallelogram and effectively contains one atom. Its area is 13.776 A˚2, which was obtained by carrying out a cooling simulation on an open 2D cluster of 100 krypton atoms. Our method for simulating the cooling process will be described in the next section. The smallest integer that satisfies the equation above is five. Since the sample is two dimensional, the total number of atoms is 25. The analysis presented here is based on krypton atoms. We replace two krypton atoms by one iodine molecule, and repeat cooling to obtain an equilibrated structure, while maintaining a constant sample area equal to 344.40 A˚2. (Figure III.3.) An abrupt truncation of the pair-wise interactions (simply setting them to be 52 zero for r greater than rc) would produce a small discontinuity leading to an infinite force. To avoid this discontinuity, potentials are conventionally shifted upwards by a small amount, which however leads to a discontinuity in the first derivative and an infinite second derivative. We wish to avoid this kind of shift because we will calculate the transition between electronic ground and excited states, in which the size of the shift may differ. The transition energy between these electronic potential energy is important and should not be obscured at the beginning of modeling the system. The approach we adopt is to smooth the potentials in a small range near rc by making the pairwise potential and its first and second derivatives continuous. The range we choose is [x1, x2], with x1 = rc, x2 = L/2− 0.1A˚. (III.46) For distances smaller than x1, the potential being used is still the general function as in Eqs. (III.41), (III.42) and (III.43); for distances between x1 and x2, we use an interpolation function g(r) given below; for distances larger than x2, the potential is zero. For Kr-Kr interaction, we have VKr−Kr(r) =  Ae−Br − C/r6 −D/r8, (r < x1) g(r), (x1 ≤ r ≤ x2) 0, (r > x2) (III.47) (Figure III.4.). The interpolation function is defined as g(r) = g0 + g1(r − x1) + g2(r − x1)2 + g3(r − x1)3 + g4(r − x1)4 + g5(r − x1)5. (III.48) 53 −8e−23 −6e−23 −4e−23 −2e−23 0 7.8 8 8.2 8.4 8.6 8.8 9 9.2 V K r− Kr /J r/(Å) VKr−Kr g(r) zero Figure III.4. Potential of Kr-Kr interaction. The red curve presents the potential in Eq. (III.42); the blue presents the interpolation function between x1 and x2 in Eq. (III.48). Here x1 = 8.095 A˚ and x2 = 9.179 A˚. Coefficients can be calculated by fulfilling the conditions of continuity at x1 and x2: g(x1) = g0 = VKr−Kr(x1) g′(x1) = g1 = V ′Kr−Kr(x1) g′′(x1) = 2g2 = V ′′Kr−Kr(x1) g(x2) = g ′(x2) = g′′(x2) = 0. Similar procedures are implemented for the iodine-iodine and iodine-krypton potentials. III.2.2. Equilibrium Configuration With explicit potential functions, we need to find the equilibrium configuration for these 25 atoms. If the total number of atoms is less than four, it is easy to solve this geometry problem by just setting the distance between each pair of atoms equal to the equilibrium distance. For example, the equilibrium value for the iodine- iodine internuclear distance is r12 = √ (x1 − x2)2 + (y1 − y2)2 = re = 2.666 A˚. For 54 a cluster of N atoms, there are N(N − 1)/2 bonds, and 2N − 3 vibrational degrees of freedom (periodic boundary conditions are not applied). For N = 2 or N = 3, we have N(N − 1)/2 = 2N − 3; since the number of coordinate equations is equal to the number of unknown variables, the set of equations can be solved. However, if the number of atoms N is bigger than or equal to four, the condition cannot be satisfied that every pair is at its equilibrium configuration, and therefore we have to do something else to obtain the minimum potential configuration. The method suitable for our case is analogous to a cooling process: start with a room temperature sample, cool it to nearly absolute zero, and obtain a crystal structure that minimizes the total energy, e.g., face-centered cubic for a pure krypton crystal. If we decrease the temperature too fast, the atoms might become stuck in a meta- stable state; but this situation can be avoided by cooling slowly and controlled for by repeating the cooling process several times with different initial conditions. We start with a uniform lattice, five atoms in a row and five in a column, with spatial intervals ∆x = ∆y = 3.7116 A˚. The density, i.e. the number of atoms per unit area, is the same as the density we found for 100 krypton atoms. The velocities of these atoms are randomly selected from a normal distribution with zero mean value and standard deviation equal to √ kBT/mKr. kB is Boltzmann constant and mKr is the atomic mass of krypton. The initial temperature is T = 300 K. The velocity of the center of mass is set to zero by adding the necessary identical increment to all atomic velocities. Then we run a classical molecular dynamics simulation on 55 these atoms. Every time step, the net force on each atom exerted by all others is calculated and therefore new velocities are updated as well as new positions. For ith atom, the force in the x-direction is F ix = − ∑ j ∂Vij ∂rij ∂rij ∂xi , (III.49) with rij = √ (xi − xj)2 + (yi − yj)2. Vij is the potential between the ith and jth atoms, as described in section III.2.1. In this molecular dynamics simulation, we adopt a velocity-Verlet algorithm [54]. The updated position and velocity in the x-direction after one time step ∆t are x(t+ ∆t) = x(t) + vx(t)∆t+ Fx(t) 2m ∆t2, (III.50) vx(t+ ∆t) = vx(t) + Fx(t+ ∆t) + Fx(t) 2m ∆t. (III.51) The same procedures are implemented for the y-direction. The time step we use is one femtosecond, which is short enough compared with the ground vibrational period of the I-I potential (∼ 150 fs), the vibrational period of the I-Kr potential (∼ 1.6 ps) and the period of the Kr-Kr potential (∼ 1.4 ps). We extract energy gradually by decreasing all the velocities of atoms by three percent every 10,000 time steps. The total time for cooling is 2 × 108 fs and therefore the final temperature is 300 K × (0.97)2×2×104 , which is less than 10−300 K. At this point, the lattice is cold enough and we get the configuration corresponding to the minimum total potential energy. Other cooling rates were tested, e.g., decreasing the velocities by 56 two percent every 10,000 steps, or three percent every 1000 steps, and they resulted in the same configuration. III.2.3. System-Bath Partition There is no universal definition of “system” versus “bath” as people partition them in different ways for different purposes [55]. We here take advantage of the timescale difference between the intramolecular vibrations and the crystal phonons, and separate them unambiguously by treating the highest-frequency normal coordinate as the system and the remaining as the bath. First, we need to determine the normal coordinates for our 2D model system. For a 1D harmonic oscillator, we know the restoring force is related to the displacement by Fx = −kx = mx¨ = −mω2x. It’s straightforward to generalize the force constant k from a scalar to a matrix, whose elements are the various second derivatives of the potential at the equilibrium configuration. For the model under study, it is a 50 by 50 matrix and the elements are ( ∂2Vtot/∂xixj ) eq , (i, j = 1, 2, · · · , 50). By diagonalizing the force constant matrix, we obtain 48 normal-mode frequencies, along with two zero frequencies corresponding to free translations of the center of mass. The normal coordinate of the highest-frequency mode (41 rad · THz) is taken to be the system coordinate, while the remaining 47 internal normal coordinates are taken to constitute the bath coordinates. We also obtain a unitary matrix U that can be used to transform Cartesian to normal coordinates, which enables us to write the total potential energy function in terms of the system coordinate q and 57 a 47-element bath-coordinate vector Q, Vtot(q,Q). It is then easy to decompose the potential function into system, bath, and interaction potentials: Vtot(q,Q) = Vs(q) + ub(Q) + u(q,Q), (III.52) with Vs(q) = Vtot(q, 0), ub(Q) = Vtot(0,Q)− Vtot(0, 0), u(q,Q) = Vtot(q,Q)− Vs(q)− ub(Q). (III.53) The form (Eq. III.52) is just that needed to implement the variational FVB/GB theory. It is important to emphasize that we are not making a normal mode approximation, but rather only making use of the definitions of the normal coordinates to specify the system and bath degrees of freedom. III.2.4. Taylor Expansion of the Potential In variational FVB/GB, we need to calculate terms of the form 〈ψν |ub|ψν¯〉, 〈ψν |uνν¯ |ψν¯〉, and their derivatives with respect to the bath wave-packet parameters. Recall that |ν〉 is an eigen state of system Hamiltonian and |ψν〉 is the corresponding bath wave packet. Since the total potential function is a sum of pair-wise atom- atom interactions, it takes the form of a rather complicated function of the normal coordinates. If the total potential function can be expressed as a power series in the various normal coordinates, however, use can be made of the multidimensional Gaussian nature of the bath wave packets to analytically calculate 〈ψν |ub|ψν¯〉 and 58 〈ψν |uνν¯ |ψν¯〉. In the light of this idea, we make a Taylor expansion of the bath and system-bath interaction potentials, and only keep terms up to fourth order. Coefficients can, in principle, be obtained by using a finite difference method. As a practical alternative, we simply fit the multidimensional bath (and system-bath interaction) potentials to quartic polynomials in the bath (or system and bath) coordinates. This approach is efficient because we need not fit the whole potential function, but only separate parts involving different combinations of up to four different coordinates. For example, a part of the bath potential only involves Q1, the lowest- frequency bath coordinate, Vfit1(Q1) = Vtot(q = 0,Q = [Q1, 0, · · · , 0]T )− Vtot(q = 0,Q = 0). This function is then fit to a polynomial of the form Vfit1(Q1) ≈ c1Q21 + c2Q31 + c3Q41. (III.54) Of course, since there are 47 bath coordinates, there are 47 parts of the bath potential involving only one bath coordinate: Vfit1(Q1), Vfit2(Q2), · · · , Vfit47(Q47). The part of potential involving both Q1 and Q2 is Vfit12(Q1, Q2) = Vtot(q = 0,Q = [Q1, Q2, 0, · · · , 0]T ) −Vtot(q = 0,Q = [Q1, 0, · · · , 0]T ) −Vtot(q = 0,Q = [0, Q2, 0, · · · , 0]T ) +Vtot(q = 0,Q = 0), (III.55) 59 which is approximated as Vfit12(Q1, Q2) ≈ c4Q21Q2 + c5Q1Q22 + c6Q21Q22 + c7Q31Q2 + c8Q1Q32. (III.56) There are 47×46/2 potential parts, like that above, involving two bath coordinates. Note that there is no bilinear term because Q1 and Q2 are normal coordinates. As we keep quartic polynomials, there are also parts of the bath potential involving three or four bath coordinates, e.g., Vfit123(Q1, Q2, Q3) and Vfit1234(Q1, Q2, Q3, Q4). The system-bath interaction is split into parts for fitting in the same way. It turns out we have about 300,000 potential functions to fit. The choice of fitting range is crucial. If the range is too small, the potential cannot show any curvature and it will be hard to accurately determine the coefficients for the fourth-order terms; if the range is too large, higher order terms like fifth or higher may play a significant role and the Taylor expansion approximation up to fourth order would break down. We therefore tested the range by comparing two different kinds of fits: one is what we are actually going to use—keeping terms up to fourth order; the other includes the fifth and sixth-order terms. There are some ranges where these two fits result in similar coefficients, and both fits promise small standard deviations from the full potential, which means that the fifth and sixth-order terms are not significant. We finally chose the largest range that satisfied this condition, which turned out to be approximately 60 [−0.2Qrms, 0.2Qrms], where, for instance, the root mean square value of lowest- frequency bath coordinate is Qrms1 = √ 1 2Ω1 . In order to carry out the fitting procedure, we adopt a Monte-Carlo-like approach. Every time step, small random changes are made to the coefficients and only accepted if they result in a smaller standard deviation from the true potential. The small changes to the coefficients are random numbers between 0 and 1 generated by the computer, adjusted by subtracting 0.5, and then multiplying by a prefactor. We also monitor the acceptance ratio during successive intervals of 1000 iterations. If the absolute values of the random numbers are quite small relative to difference between the current and future ultimate value of a given coefficient, then the acceptance ratio will roughly be one half because about half of the numbers generated are positive and half are negative, with one of these halves being in the right direction. On the other hand, if the random changes are most often too big to bring the potential parameter closer to its true value, then almost all of them will be rejected, and the acceptance ratio will be close to zero. We decided to set the acceptance percentage threshold to be ten percent imposed by this prediction. Let us say if the ratio is less than one in ten, most of those changes do not improve the fit, so we explore smaller changes, by making the prefactor ten times smaller; if the acceptance ratio is greater than one in ten, we increase the prefactor by a factor of ten. To be safe, the initial prefactor is assigned a fairly large, physically motivated value; for example, the initial prefactor associated 61 with c4 (in Eq. III.56) is 1 2 (~Ωrms1 + ~Ωrms2 )/(Qrms1 )2Qrms2 . We run the Monte Carlo calculation 1000 times, count the steps accepted, adjust the prefactor accordingly, and then repeat this procedure several times. Roughly speaking, the magnitude of the prefactor determines the precision of the coefficient. A random number between 0 and 1 is most likely to be in order of 0.1. The smaller the prefactor is, the more precise the coefficient is. Therefore, we can be sure that single precision accuracy of the coefficient is achieved and the fitting is finished when the prefactor has eventually been decreased to 10−8 times the coefficient itself. There are as many as 447,910 coefficients and terms in the Taylor expansion. Note that the approximate expansions are only made for bath and system-bath potentials. For the system, we keep the full potential. With the Taylor expansion coefficients of the potentials ready, we can evaluate 〈ψν |ub + uνν |ψν〉 and 〈ψν |uνν¯ |ψν¯〉 analytically. For example, for coupling terms like u = cnijq nQiQj, 〈ψν |uνν¯ |ψν¯〉 = cnij〈ν|qn|ν¯〉 ∫ ∞ −∞ dQ QiQj × exp[−i(Q−Qν)T · α∗ν · (Q−Qν)− iPTν · (Q−Qν)− iγ∗ν ] × exp[i(Q−Qν¯)T · αν¯ · (Q−Qν¯) + iPTν¯ · (Q−Qν¯) + iγν¯ ]. (III.57) Qi is the i th bath coordinate variable (47 in total) and Qν is the spatial center vector of the wave packet accompanying νth vibrational level. Solutions to such 62 integrals are shown in Appendix B. These potential elements can be substituted in the bath wave-packet parameter equations of motion, and then we can apply variational FVB/GB to this 2D lattice model in order to perform calculations of the overall wave function. III.2.5. The Purity As a preliminary application of our theory, we can evaluate the reduced system dynamics under nonstationary initial conditions. The overall state ket is presumed to be a sum of products of system eigen states and bath wave packets, |Ψ(t)〉 = ∑ ν |ν〉|ψν(t)〉. (III.58) The reduced system density matrix is the trace of the total density matrix over the bath degrees of freedom, σs = Trb[ρ] = Trb|Ψ〉〈Ψ| = ∫ ∞ −∞ dQ 〈Q|Ψ〉〈Ψ|Q〉 = ∑ νν¯ (∫ ∞ −∞ 〈Q|ψν〉〈ψν¯ |Q〉 ) |ν〉〈ν¯| = ∑ νν¯ 〈ψν¯ |ψν〉|ν〉〈ν¯|. (III.59) We define the purity as the trace of the square of the reduced density matrix of the system. It is one for pure states and is less than one for mixed states. Purity ≡ Trsσ2s = ∑ νν¯ |〈ψν¯ |ψν〉|2 (III.60) 63 0.97 0.975 0.98 0.985 0.99 0.995 1 1.005 0 5 10 15 20 25 30 35 40 Pu rit y t/τs dt=0.0016τsdt=0.0008τs 0.985 0.99 0.995 1 0 2 4 6 8 Figure III.5. The purity of the reduced density matrix as defined in Eq. (III.60) for the 2D I2/Kr model under nonstationary initial conditions. To calculate 〈ψν |ψν¯〉, we substitute A = iα∗ν − iαν¯ , B = 2iα∗ν ·Qν − iPν − 2iαν¯ ·Qν¯ + iPν¯ , C = −iQTν · α∗ν ·Qν + iP Tν ·Qν − iγ∗ν + iQTν¯ · αν¯ ·Qν¯ − iP Tν¯ ·Qν¯ + iγν¯ (III.61) into 〈ψν |ψν¯〉 = √ piN det[A] exp ( 1 4 BT · A−1 ·B + C ) . (III.62) The system coordinate is initially displaced by √ ~/ω (ω = 41rad · THz is the system frequency), and we include six populated vibrational states ν = 0 − 5 for the calculation. The behavior of the purity is presented in Figure III.5. Calculations 64 were carried out with a fourth-order Runge-Kutta algorithm using fixed time steps of dt = 0.0016τs and dt = 0.0008τs (τs = 2pi/ω = 153 fs). Results obtained with these two different time steps agree until 30τs. At first the purity oscillates with a period of about one system vibration τs, and then the oscillation frequency doubles after five system periods. It decreases with a damping rate of about 0.001 τ−1s (0.007 ps−1) for the first 10 τs and then have a steady oscillation around 0.988 for another 20 τs. Doublets starting to show up after first five periods can be attributed to overtone beats. At first, only adjacent vibrational states are relatively strongly coupled; as time goes by, 0–2, 1–3, 2–4, and other similar beats start to contribute more effectively to the purity. III.2.6. Cage Effect There are some peculiarities associated with our model system, one of which is that it is dissociation free. Because the iodine molecule is caged in the krypton matrix, it cannot dissociate as it does in the gas phase. When the iodine atoms separate from each other very widely, their interactions with nearest in-line krypton atoms are fairly strong. These krypton atoms will push the iodine atoms back together and hence iodine dissociation is effectively avoided. In the gas phase, the equilibrium internuclear distance of iodine in the electronic ground state is 2.666 A˚, while in a 2D krypton crystal, the equilibrium configuration, i.e. when the system and bath coordinates are all zero, corresponds to an iodine internuclear distance of 2.6647 A˚. If the iodine molecule is pumped to the electronic excited (B) 65 −5 0 5 10 15 20 25 −15000 −10000 −5000 0 X: 20.76 Y: −1.105e+004 System potential, excited state q V e /c m − 1 Figure III.6. The system potential Vs(q) in the excited (B) state. The system coordinate q vibrates in the indicated range. state, then it vibrates with the system coordinate q varying from zero to 20.76 units, corresponding to an iodine-iodine distance from 2.6647 A˚ to 3.690 A˚ (Figure III.6.). The unit of q is √ ~/ω, where ω = 41rad · THz is the system frequency. From the equilibrium configuration, we know that the nearest in-line krypton atoms are separated from each other by 10.765 A˚, greater than 3.690 A˚. The amplitude of krypton’s motions is quite small compared with iodine. In this case, the iodine atoms will not pass the nearest krypton atoms and the motion is confined to be intramolecular. We can compare the potential of iodine-iodine interaction in the gas phase and the system potential in our 2D crystal, both in the electronic ground state (Figure III.7.). These two potentials very nearly coincide at short distances. 66 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 −1.5 −1 −0.5 0 0.5 1 x 104 r/angstrom V I I/c m − 1 −10 −5 0 5 10 15 20 25 30 −2.5 −2 −1.5 −1 −0.5 0 x 104 q V s /c m − 1 Figure III.7. The system potential Vs(q) in the ground state (blue), compared with I-I interaction potential (red) as in Eq. (III.41). At longer distances, they deviate from each other: the free iodine molecular potential flattens toward a constant value determined by the dissociation energy; the system potential in the 2D crystal goes higher and steeper with larger q. Thus the system potential is bound and the iodine molecule will not dissociate as it does in the gas phase. It is worth noting that the system coordinate q mostly corresponds to the iodine vibration. Recall that q is the highest-frequency normal coordinate and that it is a linear combination of 50 mass-weighted Cartesian coordinates of 25 atoms (23 krypton and two iodine): q = c1x1 + c2y1 + c3x2 + c4y2 + · · ·+ c49x25 + c50y25. (III.63) The first four Cartesian coordinates belong to the two iodine atoms. The coefficients 67 are normalized: c21 + c 2 2 + · · ·+ c250 = 1, (III.64) and are very small numbers except the first four, c21 + c 2 2 + c 2 3 + c 2 4 = 0.999907, (III.65) which means that iodine intramolecular motion participates in 99.99% of the behavior of the system coordinate. This explains the high degree of similarity between iodine-iodine interaction potential and the system potential at short distances. 68 CHAPTER IV LINEAR WAVE-PACKET INTERFEROMETRY AND TR-CARS With the model of 2D I2/Kr established, we are now in a position to implement FVB/GB and perform simulations of ultrafast spectroscopies. Various kinds of ultrafast optical experiments have been performed on samples of dihalogen molecules embedded in rare gas solids, including pump-probe [3, 10, 11, 13, 53, 56], spectrally resolved transient grating and, most recently, resonant 5-color 4-wave mixing [15]. We focus on two of them: linear wave-packet interferometry [14] and time-resolved coherent anti-Stokes Raman scattering [6, 7, 9]. The former provides information on wave-packet dynamics in an excited electronic state and the latter is primarily sensitive to nuclear motion in the electronic ground state. Simulations of these two experiments are performed without detailed information on higher-lying states, on which limited data may exist. First we carry out simulations with frozen Gaussian wave packets for the bath, meaning the width parameter αν is not allowed to change with time. The “frozen” FVB/GB is a simplified theory, with which the results of the full, “thawed” FVB/GB theory can be compared. We have developed and employed procedures to stabilize the numerical calculations, which are described in Appendix C. 69 IV.1. Linear Wave-Packet Interferometry IV.1.1. Theoretical Expression In a linear wave-packet interferometry experiment, the B ← X transition of an iodine molecule is driven by two resonant phase-locked ultrashort pulses, and an induced two-field-dependent fluorescence emission is detected. The signal provides information on the dynamics of the molecule on the excited state electronic potential energy curve. It has been demonstrated that the in-phase and in- quadrature phase-locked interferograms can be combined to determine the linear susceptibility [35]. In particular, the imaginary part of the linear susceptibility can give us a linear absorption spectrum within the pulse frequency range. To formulate a theoretical description of the linear wave-packet interferometry signal, we start from a Hamiltonian that includes the two-state Hamiltonian for the sample and the interaction with the laser pulses, H(t) = H + V (t), (IV.1) H = |g〉Hg〈g|+ |e〉He〈e|, (IV.2) where g and e denote ground and excited electronic states respectively. The interaction V (t) = V1(t) + V2(t), (IV.3) 70 is given in the dipole approximation by Vj(t) = −µˆEj(t), (IV.4) in which µˆ = µ(|e〉〈g|+ |g〉〈e|) (IV.5) is the dipole moment operator of iodine. The electric fields of the two resonant pulses have the form, E1(t) = Ee −(t−t1)2/2σ2 cos[Ω(t− t1)], (IV.6) and E2(t) = Ee −(t−t2)2/2σ2 cos[Ω(t− t2) + ΩLtd + φ], (IV.7) with a delay time td = t2 − t1. (IV.8) Ω is the center frequency and ΩL is the locking frequency. For the special case Ω = ΩL, we would have Ω(t− t2) + ΩLtd + φ = Ω(t− t1) + φ, (IV.9) and the phase difference between the two pulses would be φ. To further illustrate the concept of locking frequency, we take the Fourier transforms of the electric 71 fields: E˜1(ω) = ∫ ∞ −∞ dt eiωtE1(t) = E 2 eiωt1 √ 2piσ2 [ e−(ω+Ω) 2σ2/2 + e−(ω−Ω) 2σ2/2 ] , (IV.10) E˜2(ω) = ∫ ∞ −∞ dt eiωtE2(t) = E 2 √ 2piσ2 [ e−(ω+Ω) 2σ2/2+i(ωt2+ΩLtd+φ) + e−(ω−Ω) 2σ2/2+i(ωt2−ΩLtd−φ) ] . (IV.11) In Eqs. (IV.10) and (IV.11), the first term is negligible compared with the second one. For the Fourier components at the locking frequency, we have E˜1(ΩL) ∼= E 2 √ 2piσ2e−(ΩL−Ω) 2σ2/2+iΩLt1 , (IV.12) E˜2(ΩL) ∼= E 2 √ 2piσ2e−(ΩL−Ω) 2σ2/2+i(ΩLt2−ΩLtd−φ) = e−iφE˜1(ΩL). (IV.13) The difference between E˜1(ΩL) and E˜2(ΩL) is just a phase factor e −iφ, and φ is defined as the phase locking angle. For other frequencies, E˜2(ω) ∼= ei(ω−ΩL)td−iφE˜1(ω), (IV.14) The phase difference of frequency components of the two fields is dependent on the particular frequency, the locking frequency, and the delay time. First order time-dependent perturbation theory leads to |Ψ(t)〉 = e−iH(t−t0)|Ψ(t0)〉 − i ∫ t t0 dτ eiH(τ−t)V (τ)e−iH(τ−t0)|Ψ(t0)〉. (IV.15) 72 The initial time t0 is long before the first pulse. The combined system and bath are initially taken to be in the overall ground state, |Ψ(t0)〉 = |g〉|0g〉|ψ0g〉. (IV.16) We define a reduced pulse propagator p(eg) ≡ 1 σ ∫ ∞ −∞ dτ eiHeτe−iHgτe−iΩτ−τ 2/2σ2 . (IV.17) Because fluorescence is a slow process (typically on the order of hundreds of nanoseconds), we can set the upper limit of integration to infinity. We have also made a rotating-wave approximation. The pulse propagators are defined by P1 = i 2 µEσ [ p(eg)|e〉〈g|+ p(ge)|g〉〈e|] , (IV.18) P2 = i 2 µEσ [ p(eg)e−iΩLtd−iφ|e〉〈g|+ p(ge)eiΩLtd+iφ|g〉〈e|] . (IV.19) The total state ket can be written in terms of pulse propagators as |Ψ(t)〉 = e−iH(t−t0)|Ψ(t0)〉+ e−iH(t−t1)P1e−iH(t1−t0)|Ψ(t0)〉 +e−iH(t−t2)P2e−iH(t2−t0)|Ψ(t0)〉. (IV.20) Population in the excited electronic state is 〈e|Ψ(t)〉〈Ψ(t)|e〉, with 〈e|Ψ(t)〉 = i 2 µEσe−iHe(t−t1)p(eg)e−iHg(t1−t0)|0g〉|ψ0g〉 + i 2 µEσe−iHe(t−t2)p(eg)e−iΩLtd−iφe−iHg(t2−t0)|0g〉|ψ0g〉. (IV.21) 73 The contribution to the excited state population due to interference between the e-state amplitudes generated by pulse 1 and 2 is Wφ(td) = 1 4 µ2E2σ2ei(ΩLtd+E0g td+φ)〈ψ0g |〈0g|p(ge)e−iHetdp(eg)|0g〉|ψ0g〉+ c.c. (IV.22) For the in-phase case, W0(td) = 1 2 µ2E2σ2Re [ ei(ΩLtd+E0g td)〈ψ0g |〈0g|p(ge)e−iHetdp(eg)|0g〉|ψ0g〉 ] ; (IV.23) for the in-quadrature case, Wpi/2(td) = −1 2 µ2E2σ2Im [ ei(ΩLtd+E0g td)〈ψ0g |〈0g|p(ge)e−iHetdp(eg)|0g〉|ψ0g〉 ] . (IV.24) The cosine and sine transforms of the population due to interference are given by Re [Wφ(ω)] = ∫ ∞ 0 dt cos(ωt)Wφ(t) (IV.25) and Im [Wφ(ω)] = ∫ ∞ 0 dt sin(ωt)Wφ(t). (IV.26) The absorptive linear susceptibility is expressed as [35] χ′′(ΩL ± ω) = 1 4pi2|E˜1(ΩL ± ω)|2 { Re [W0(ω)]± Im [ Wpi/2(ω) ]} . (IV.27) The dispersive linear susceptibility is approximately χ′(ΩL ± ω) ∼= 1 4pi|E˜1(ΩL ± ω)|2 { Re [ Wpi/2(ω) ]∓ Im [W0(ω)]} . (IV.28) 74 -0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 1 2 3 4 5 6 In te rfe re nc e Si gn al (a rb. un its ) t/τs φ=0 φ=pi/2 Figure IV.1. The in-phase and in-quadrature interferograms calculated with frozen Gaussians. The center wavelength of both pulses is 580 nm and the temporal width of the pulses is 20 fs. IV.1.2. Simulation Results First we simulate linear wave-packet interferometry signals on the model of 2D I2/Kr (described in Chapter III) with two phase-locked pulses both centered at 580 nm and having a FWHM of 20 fs. The locking frequency is the same as the center frequency. The in-phase (Eq. IV.23) and in-quadrature (Eq. IV.24) interferograms calculated with frozen-Gaussian variational FVB/GB are shown in Figure IV.1. We include 34 vibrational levels of the excited state and keep couplings between any level ν and its adjacent levels ν¯ = ν±4 (see more in Appendix C). The unit of time τs is one ground state vibrational period (153 fs), while the excited state vibrational period is 251 fs for νe = 0 and 314 fs for νe = 26. The initial peak has a width about one third of τs, corresponding to the time for loss of overlap between the wave packet prepared by the first pulse and that prepared by the second. There is a small peak around 2τs, corresponding to the round-trip time of the wave packet 75 0 0.0002 0.0004 0.0006 0.0008 0.001 -4 -3 -2 -1 0 1 a bs or pt io n (ar b. un its ) (ω-ΩL)/ωs 20fs 10fs Figure IV.2. The absorptive portion of the linear susceptibility χ′′ calculated using frozen-Gaussian FVB/GB. The temporal width of both pulses is 20 fs (or 10 fs in another case) and the center wavelength of both pulses is 580 nm. generated by the first pulse. Compared with gas phase cases [35], interferograms of our 2D I2/Kr model do not exhibit subsequent peaks at the following excited state vibrational periods. We suggest an interpretation of these missing bumps as resulting from a cage effect—short dephasing time prevents the wave packet evolving on the excited state for more than one cycle from having sizable overlap with the one just being excited from the ground state. The absorptive part of the linear susceptibility is calculated according to Eq. (IV.27) and presented in Figure IV.2. In a small range near the locking frequency, peaks show up corresponding to excitations to different vibrational levels νe of the excited state. The wavelength 580 nm corresponds to a transition of νg = 0 in the ground state to νe = 13 in the excited state. The Franck–Condon overlap 〈νe|0g〉 has a maximum at νe = 26, transition to which would require the energy of a photon with wavelength of 535.8 nm. Absorption grows with frequency 76 Figure IV.3. Motions of two highest-frequency bath modes. The left is the highest-frequency bath mode and the right is the second-to-highest-frequency bath mode. higher than the locking frequency. The absorption spectrum obtained using 10-fs pulses is similar to that with 20-fs pulses, except that the linewidths of some peaks are a little wider. We do not observe phonon side bands in the absorption spectrum. The linewidths of the absorption peaks are fairly wide, making it difficult to discern phonon side bands between subsequent zero-phonon peaks. However, we can investigate the generation of zone boundary phonons (ZBP) by watching the bath response after excitation by the first laser pulse. ZBP have zero group velocity and the highest frequency, which correspond to the highest-frequency mode (47th) for our model of 2D I2/Kr. The second-to-highest bath frequency is close to the highest. We assign these two modes as ZBP, with wave numbers of 59.2 cm−1 and 58.3 cm−1, respectively. Figure IV.3. shows the motions of these two bath modes, which are mostly parallel to the vibration of iodine molecule. Gu¨hr, Bargheer 77 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 1 2 3 4 5 6 7 8 En er gy /(- h ω s) t/τs 46th bath mode 47th bath mode Figure IV.4. The energy expectation value for the highest-frequency and second- to-highest-frequency bath modes after the first pulse. and Schwentner reported ZBP of 51.5 cm−1 in pump-probe spectra of I2 in a Kr crystal [10]; Karavitis and Apkarian observed a local mode of 41.5 cm−1 in tr-CARS for the same sample [6]. We calculate the energy involving only the 47th bath coordinate and momentum. E47 = 〈Ψ|P 247/2 + V (Q47)|Ψ〉/〈Ψ|Ψ〉 = ∑ ν 〈ψν |P 247/2 + V (Q47)|ψν〉/ ∑ ν 〈ψν |ψν〉. (IV.29) A similar calculation is performed for E46. The energies of these two bath modes (Figure IV.4.) start with values 0.4% greater than the zero-point energies because the normal coordinates involve iodine coordinates, though a very small percentage, and the excitation of iodine makes these bath modes weakly Franck– Condon active. Interestingly, the energies begin to increase markedly after one excited-state vibrational period, and then present an oscillation as well as a continual increase. We can also watch the coordinate trajectory of these two 78 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 8 Q /   √- h /ω s t/τs Q46Q47 Figure IV.5. The expectation value of the highest-frequency and second-to-highest-frequency bath coordinates. 〈Ψ|Q47|Ψ〉/〈Ψ|Ψ〉 =∑ ν 〈ψν |Q47|ψν〉/ ∑ ν 〈ψν |ψν〉. Similar for Q46. bath modes (Figure IV.5.). They start to move noticeably after one excite-state period, then oscillate at the system frequency rather than their characteristic mode frequency for one cycle, and then the oscillation slows down a little. These behaviors suggest a non-negligible coupling between ZBP and the iodine molecular vibration. Using the root mean square value of Q and the equilibrium value of q in the excite state, we find that the interaction energy proportional to qQ is about 0.0021~ωs for the 46th bath mode and 0.0011~ωs for the 47th mode. Among 47 bath modes, these two modes are not the most strongly coupled to the intramolecular motion. The interaction energy (proportional to qQ) is 0.017~ωs for the 28th bath mode, which is the most strongly coupled to the iodine vibration and involves more distinct motions of the nearest in-line Kr atoms. However, from the behaviors of energies and coordinate trajectories of the 46th and 47th bath 79 modes, we cannot agree that ZBP are decoupled from the iodine intramolecular vibration. IV.2. Tr-CARS Time-resolved coherent anti-Stokes Raman scattering (tr-CARS) is another ultrafast optical experiment whose signal we wish to simulate and interpret. For simplicity, we imagine the molecular system starting from the ground vibrational level of its ground electronic state. The first pulse pumps the molecule to its excited electronic state. The second pulse is the Stokes pulse, which has a longer center wavelength than the pump. It dumps the molecule back to and generates coherent motion on the ground electronic potential energy surface. After a variable delay, the third pulse probes the system. These three resonant laser pulses induce a third-order molecular dipole moment. A nonlinear optical signal is detected in the wave-vector matched direction as a function of the delay time between the second and the third pulses. During this delay time t32, several vibrational levels of the ground state are coupled by the bath and undergo vibrational decoherence, which is the process we are interested in. IV.2.1. Theoretical Expression We derive a theoretical expression of the tr-CARS signal from a sample of iodine molecules embedded in a krypton crystal. The Hamiltonian consists of a 80 time-independent portion and a time-dependent interaction with three laser pulses: H(t) = H + V (t), (IV.30) where V (t) = V1(t) + V2(t) + V3(t), (IV.31) with the laser pulses being numbered 1, 2 and 3. The time-independent portion of Hamiltonian for our system is H = |g〉Hg〈g|+ |e〉He〈e|, (IV.32) in which Hg and He are the nuclear Hamiltonian in the ground and the excited states, respectively. The matter-field interaction for the j-th pulse is written in the dipole approximation, Vj(t) = −µˆEj(t), (IV.33) where the dipole moment operator is µˆ = µ(|g〉〈e|+ |e〉〈g|), (IV.34) and the transition dipole µ is assumed to be independent of nuclear coordinates. The electric field for the j-th pulse is taken to have a Gaussian envelope, Ej(t) = Ej exp [ − 1 2σ2j (t− tj)2 ] cos [Ωj(t− tj) + φj] , (IV.35) where σj is the temporal width, Ωj is the center frequency, φj is the optical phase, and tj is the pulse arrival time. The center frequencies of the pulses are chosen to 81 drive the B ← X electronic transition of I2. Note that tj is location dependent, with tj = t ′ j + kj · r/Ωj, where r points from some origin to a point within the sample, t′j is the pulse arrival time at the origin, and kj is the wave vector. We treat the interactions with the laser pulses as perturbations. In the interaction picture, the state ket is |Ψ˜〉 = eiHt|Ψ(t)〉, (IV.36) (note that ~ ≡ 1) and the Schro¨dinger equation becomes i ∂ ∂t |Ψ˜(t)〉 = V˜ (t)|Ψ˜(t)〉, (IV.37) where V˜ (t) = eiHtV (t)e−iHt. Therefore, |Ψ˜(t)〉 = |Ψ˜(t0)〉 − i ∫ t t0 dτ V˜ (τ)|Ψ˜(τ)〉. (IV.38) Third order time-dependent perturbation theory leads to |Ψ˜(t) ∼= |Ψ˜(t0)〉 − i ∫ t t0 dτ V˜ (τ)|Ψ˜(t0)〉 − ∫ t t0 dτ2 V˜ (τ2) ∫ τ2 t0 dτ1 V˜ (τ1)|Ψ˜(t0)〉 +i ∫ t t0 dτ3 V˜ (τ3) ∫ τ3 t0 dτ2 V˜ (τ2) ∫ τ2 t0 dτ1 V˜ (τ1)|Ψ˜(t0)〉. (IV.39) Reverting to the Schro¨dinger picture, we obtain an expression for the state after 82 the third pulse: |Ψ(t)〉 ∼= e−iH(t−t0)|Ψ(t0)〉 − i ∫ t t0 dτ eiH(τ−t)V (τ)e−iH(τ−t0)|Ψ(t0)〉 − ∫ t t0 dτ2 e iH(τ2−t)V (τ2)e−iHτ2 ∫ τ2 t0 dτ1 e iHτ1V (τ1)e −iH(τ1−t0)|Ψ(t0)〉 +i ∫ t t0 dτ3 e iH(τ3−t)V (τ3)e−iHτ3 ∫ τ3 t0 dτ2 e iHτ2V (τ2)e −iHτ2 × ∫ τ2 t0 dτ1 e iHτ1V (τ1)e −iH(τ1−t0)|Ψ(t0)〉. (IV.40) Next, we define pulse propagators Pj(t; τ) ≡ −i ∫ t t0 dτ eiH(τ−tj)Vj(τ)e−iH(τ−tj). (IV.41) Let the initial time t0 be long before the arrival of the first pulse. The lower limit of integration can be replaced by −∞. We make a rotating-wave approximation and obtain Pj(t; τ) ∼= iµEj 2 [∫ t −∞ dτ eiHe(τ−tj)|e〉〈g|e− 1 2σ2 j (τ−tj)2 e−iΩj(τ−tj)−iφje−iHg(τ−tj) +H.c.] . (IV.42) Then at time t, the state of the system takes the form |Ψ(t)〉 ∼= e−iH(t−t0)|Ψ(t0)〉+ ∑ j e−iH(t−tj)Pj(t; τ)e−iH(tj−t0)|Ψ(t0)〉 + ∑ jk e−iH(t−tk)Pk(t; τ2)e−iH(tk−tj)Pj(τ2; τ1)e−iH(tj−t0)|Ψ(t0)〉 + ∑ jkl e−iH(t−tl)Pl(t; τ3)e−iH(tl−tk)Pk(τ3; τ2)e−iH(tk−tj)Pj(τ2; τ1) ×e−iH(tj−t0)|Ψ(t0)〉. (IV.43) 83 Note the nesting of pulse propagators occurring in this expression: in the second order terms, for example, the first argument of Pj is the second argument of Pk. We may further define reduced pulse propagators p (eg) j (t; τ) ≡ 1 σj ∫ t −∞ dτ eiHe(τ−tj)e − 1 2σ2 j (τ−tj)2 e−iΩj(τ−tj)e−iHg(τ−tj), (IV.44) which are dimensionless operators acting on the vibrational state. Then Pj(t; τ) = iµEjσj 2 [ |e〉〈g|e−iφjp(eg)j (t; τ) +H.c. ] . (IV.45) Note also that ( p (eg) j (t; τ) )† = p (ge) j (t; τ). (IV.46) We take |Ψ(t0)〉 to be the state ket for the ground vibrational level in the ground electronic state, |Ψ(t0)〉 = |g〉|0g〉|ψ0g〉. (IV.47) Then, |Ψ(t)〉 = e−iHg(t−t0)|g〉|0g〉|ψ0g〉 +|e〉 ∑ j e−iHe(t−tj) iµEjσj 2 e−iφjp(eg)j (t; τ)e −iHg(tj−t0)|0g〉|ψ0g〉 −|g〉 ∑ jk e−iHg(t−tk) µ2EjEkσjσk 4 eiφk−iφjp(ge)k (t; τ2)e −iHe(tk−tj) ×p(eg)j (τ2; τ1)e−iHg(tj−t0)|0g〉|ψ0g〉 −|e〉 ∑ jkl e−iHe(t−tl) iµ3ElEkEjσlσkσj 8 e−iφl+iφk−iφjp(eg)l (t; τ3)e −iHg(tl−tk) ×p(ge)k (τ3; τ2)e−iHe(tk−tj)p(eg)j (τ2; τ1)e−iHg(tj−t0)|0g〉|ψ0g〉. (IV.48) 84 This third-order expansion takes full account of any temporal pulse overlap that may occur. The time-dependent molecular dipole moment is responsible for the coherent signal emission. The emission field from an individual chromophore is [57] e(R, t) = (µ¨(t4)× n)× n c2|R− r| , (IV.49) where r is a location within the sample, and the location of a detector R is far from the sample. n ≡ n4 is a unit vector pointing from r to R (the detector is positioned in the same direction of the emitted light). Since the detector is far from the sample, we can approximate n as pointing along R from an origin within the sample. t4 = t − |R − r|/c is the time the emitted field is generated at the chromophore. We define t′4 as the emission time from the origin of the sample and write t4 = t ′ 4 + n4 · r/c. (IV.50) Because µ(t) ∼ e−iΩ4t, we have µ¨(t) ∼ −Ω24µ(t). As the optical line widths are small compared to the electronic transition frequency, we use µ¨(t) ∼= −Ω¯2µ(t) for all the chromophores, and ~Ω¯ is the electronic transition energy. The direction of the induced dipole is approximately perpendicular to the propagation direction of the emitted light. Therefore (µ× n)× n can be simplified with (−µ). Approximating |R− r| ∼= R in the denominator, we obtain the net field emitted from the sample E(R, t) ∼= ρΩ¯ 2 Rc2 ∫ V dr µ(t4). (IV.51) 85 where ρ is the density—the number of iodine molecules per volume. The expectation value of the induced dipole moment is 〈Ψ(t)|µˆ|Ψ(t)〉. Contributing to tr-CARS are the terms involving all three pulses, µtrCARS(t4) = ∑ jkl (Mjkl(t4) +Njkl(t4)) + c.c., (IV.52) where Mjkl(t4) = − i 8 µ4EjEkElσjσkσle −iφl+iφk−iφj〈ψ0g |〈0g|eiHg(t4−t0)e−iHe(t4−tl) ×p(eg)l (t4; τ3)e−iHg(tl−tk)p(ge)k (τ3; τ2)e−iHe(tk−tj) ×p(eg)j (τ2; τ1)e−iHg(tj−t0)|0g〉|ψ0g〉, (IV.53) and Njkl(t4) = i 8 µ4EjEkElσjσkσle iφl+iφk−iφj〈ψ0g |〈0g|eiHg(tl−t0)p(ge)l (t4; τ)eiHe(t4−tl) ×e−iHg(t4−tk)p(ge)k (t4; τ2)e−iHe(tk−tj)p(eg)j (τ2; τ1)e−iHg(tj−t0)|0g〉|ψ0g〉. (IV.54) Note that the reduced pulse propagators retain location dependence because they involve the position-dependent pulse arrival times tj = t ′ j + kj · r. If we require that the angle between two laser pulses be small enough so that (nk − nj) · r/c is negligible compared with one vibrational period (nj is a unit vector pointing the same direction of kj, with nj/c = kj/Ωj), then 2pirθjk| cos θ|/c 150 fs, (IV.55) 86 where r is taken to be the laser spot size about 35 µm [16], θjk is the angle between nk and nj, and θ is the angle between r and (nk − nj) (which can be from 0 to pi). The inequation (IV.55) leads to θjk  11.7◦. The variation across the sample of the time delays between pulses (geometrical distortion) is therefore negligible on the timescale of vibrations if the laser beams are at most a few degrees from collinearity. With this requirement, we are safe to approximate (tj − tk) as (t′j − t′k) in the reduced pulse propagators and evaluate them at the sample origin. We also approximate e−iHe(nk−nj)·r/c as e−iεe(nk−nj)·r/c (and so for other similar factors). Since geometrical distortion can be neglected on a vibrational timescale, this approximation is also guaranteed by the requirement of small angle between pulses. We can write e−i(εe−εg)kj ·r/Ωj as e−ikj ·r by making use of the fact that these pulses are resonant with the electronic transition (k4 is the signal-field wave vector). Now we have Mjkl(t4) = − i 8 µ4EjEkElσjσkσle −iφl+iφk−iφje−i(k4−kl+kk−kj)·rM (0)jkl (t ′ 4), (IV.56) where M (0) jkl (t ′ 4) = 〈ψ0g |〈0g|eiHg(t ′ 4−t0)e−iHe(t ′ 4−t′l)p(eg)l (t ′ 4; τ3)e −iHg(t′l−t′k) ×p(ge)k (τ3; τ2)e−iHe(t ′ k−t′j)p(eg)j (τ2; τ1)e −iHg(t′j−t0)|0g〉|ψ0g〉, (IV.57) is not position dependent. M (0) jkl (t ′ 4) is the overlap of two wave packets: one is the initial wave packet staying in the ground state for a time t′4−t0; the other is a three- field-dependent wave packet: the initial wave packet stays in the ground state for 87 a time t′j − t0, then is excited by pulse j, evolves in the excited state for a time t′k − t′j, then is dumped back to the ground state by pulse k, evolves for another time t′l− t′k, then is excited again by pulse l, and finally evolves in the excited state for t′4 − t′l. Integrating Mjkl over the sample volume, we obtain a sharply peaked function centered at the emitted light wave vector k4 equal to kl−kk+kj. Tr-CARS collects signals in the direction k1 − k2 + k3. Therefore, M123 and M321 contribute to tr-CARS dipole moment. We can see from the expression for M321 that the third pulse acts on the initial wave packet first, then the second, and the first pulse acts last. Only when the three pulses are nearly overlapped is M321 nonnegligible. Another group of third order dipole moments is of the form Njkl(t4) = i 8 µ4EjEkElσjσkσle iφl+iφk−iφje−i(−k4+kl+kk−kj)·rN (0)jkl (t ′ 4), (IV.58) where N (0) jkl (t ′ 4) = 〈ψ0g |〈0g|eiHg(t ′ l−t0)p(ge)l (t ′ 4; τ)e iHe(t′4−t′l)e−iHg(t ′ 4−t′k) ×p(ge)k (t′4; τ2)e−iHe(t ′ k−t′j)p(eg)j (τ2; τ1)e −iHg(t′j−t0)|0g〉|ψ0g〉. (IV.59) Unlike M (0) jkl (t ′ 4), N (0) jkl (t ′ 4) is the overlap of a one-pulse wave packet and a two-pulse wave packet. The one-pulse wave packet is the initial one excited by pulse l after t′l − t0 and evolving in the excited state for t′4 − t′l. The two-pulse wave packet is the initial one excited by pulse j, dumped by pulse k, and then evolving in the ground state for t′4 − t′k. The vector determining the phase factor of Njkl(t4) is k4 ∼= kl +kk−kj. Hence N213 and N231 contribute to the tr-CARS dipole moment. 88 In N213, the second pulse acts before the first. Since the first two pulses are nearly overlapped in tr-CARS experiments, N213 can be nonnegligible. N (0) 231 is the overlap of two wave packets: one is pumped by pulse 1 and then evolves in the excited state; the other is pumped by pulse 2 and dumped by pulse 3. Because the interval between the first two pulses is short, these two wave packets are similar before t′3. After t′3, one wave packet evolves in the excited state, while the other does so in the ground state. Since these wave packets move very differently, their overlap N (0) 231 is typically small, and N231 makes only a small contribution to the tr-CARS signal. Overall, the dipole moment generating tr-CARS is µtrCARS(t4) = [M123(t4) +M321(t4) +N213(t4) +N231(t4)] + c.c. (IV.60) The net emitted field is E(R, t) ∼= ρΩ¯ 2 Rc2 ∫ V dr µ(t4) = ρΩ¯2 Rc2 ∫ V dr [(M123(t4) +M321(t4) +N213(t4) +N231(t4)) + c.c.] = ρΩ¯2 Rc2 ∫ V dr { − i 8 µ4E1E2E3σ1σ2σ3e −iφ3+iφ2−iφ1e−i(k4−k3+k2−k1)·r × ( M (0) 123(t ′ 4) +M (0) 321(t ′ 4) +N (0)∗ 213 (t ′ 4) +N (0)∗ 231 (t ′ 4) ) + c.c. } = ρV Ω¯2 8Rc2 µ4E1E2E3σ1σ2σ3f(k4 − k3 + k2 − k1) {−ie−iφ3+iφ2−iφ1 × ( M (0) 123(t ′ 4) +M (0) 321(t ′ 4) +N (0)∗ 213 (t ′ 4) +N (0)∗ 231 (t ′ 4) ) + c.c. } . (IV.61) If the integration were carried out over the whole space, we would obtain a delta function of (k4−k3 +k2−k1). However, we integrate over the sample volume, and therefore obtain a sharply peaked function f centered at k4 equal to k1 − k2 + k3. 89 k4 is the wave vector of the emitted light. The direction of R, pointing from an origin within the sample to the detector, is the same as the direction of k4 and should be along with k1 − k2 + k3 to obtain maximum tr-CARS signals. The tr-CARS signal is the time-integrated intensity of the emitted field, S ∝ ∫ ∞ t′3−3σ3 dt′4 E 2(R, t). (IV.62) It is worth noting that there can be emitting field even before the arrival time of the third pulse. Hence we set the lower limit of integration to be t′3−3σ3 in order to include all portions of the signal. The square of Eq. (IV.61) includes φ-dependent terms, which will average to zero over many laser shots since the experiments do not use phase-locked pulses. Then we obtain an expression for the signal S ∝ ∫ ∞ t′3−3σ3 dt′4 ∣∣∣M (0)123(t′4) +M (0)321(t′4) +N (0)∗213 (t′4) +N (0)∗231 (t′4)∣∣∣2 . (IV.63) IV.2.2. Bath Wave Packets after Pulse-Driven Transitions For one electronic state, wave packet accompanying each vibrational level is approximated as a Gaussian function. After a laser pulse, electronic transition occurs between one electronic state and another, i.e. |g〉 → |e〉 or vice versa. Then the wave packet accompanying each vibrational level in the new electronic state is not a single Gaussian function, but a weighted sum of Gaussian wave packets of the old electronic state. To implement FVB/GB and propagate wave packets in the new electronic state, we write each wave packet as a single Gaussian function by finding the population of each vibrational level, the spatial center of the wave packet, the 90 momentum, the variance, the position-momentum correlation, and the phase of the wave packet. Let’s start with a vibrational coherent state of the electronic ground state, |Ψ〉 = |g〉 ∑ νg |νg〉|ψνg〉. (IV.64) After one pulse, the state ket is P (τ2; τ1)|Ψ〉 = iµEσ 2 (p(eg)e−iφ|e〉〈g|+ p(ge)eiφ|g〉〈e|)|Ψ〉 = |e〉iµEσe −iφ 2 ∑ νg p(eg)|νg〉|ψνg〉. (IV.65) P (τ2; τ1) is the pulse propagator introduced in previous sections. Using the excited state vibrational ket as basis, we express the state ket as P (τ2; τ1)|Ψ〉 = |e〉 ∑ νe |νe〉|ψνe〉, (IV.66) and find |ψνe〉 = iµEσe−iφ 2 ∑ νg 〈νe|p(eg)|νg〉|ψνg〉. (IV.67) Recall that the reduced pulse propagator is p(eg) = 1 σ ∫ τ2 −∞ dτ1 e iHeτ1e−iHgτ1e−iΩτ1−τ 2 1 /2σ 2 . (IV.68) Then the bath wave packet for each vibrational level νe is 〈Q|ψνe〉 = ∑ νg f(νe, νg)〈νe|νg〉〈Q|ψνg〉, (IV.69) 91 where f(νe, νg) ≡ iµEe −iφ 2 ∫ τ2 −∞ dτ1 e −τ21 /2σ2+i(Eνe−Eνg−Ω)τ1 . (IV.70) Since 〈Q|ψνg〉 is approximated as a Gaussian function, 〈Q|ψνe〉 is a sum of νmaxg + 1 Gaussian functions. We approximate it as one single Gaussian 〈Q|ψνe〉 = exp[i(Q−Qνe)T ·ανe · (Q−Qνe) + iP Tνe · (Q−Qνe) + iγνe ] and find the wave-packet parameters: 〈ψνe|ψνe〉 = ∑ νg ,ν¯g f ∗(νe, νg)f(νe, ν¯g)〈νg|νe〉〈νe|ν¯g〉〈ψνg |ψν¯g〉, (IV.71) 〈ψνe|Q|ψνe〉 = ∑ νg ,ν¯g f ∗(νe, νg)f(νe, ν¯g)〈νg|νe〉〈νe|ν¯g〉〈ψνg |Q|ψν¯g〉, (IV.72) Qνe = 〈ψνe|Q|ψνe〉 〈ψνe |ψνe〉 , (IV.73) 〈ψνe|P |ψνe〉 = ∑ νg ,ν¯g f ∗(νe, νg)f(νe, ν¯g)〈νg|νe〉〈νe|ν¯g〉〈ψνg |P |ψν¯g〉, (IV.74) Pνe = 〈ψνe |P |ψνe〉 〈ψνe |ψνe〉 , (IV.75) 〈ψνe|(Q−Qνe)(Q−Qνe)T |ψνe〉 = ∑ νg ,ν¯g f ∗(νe, νg)f(νe, ν¯g)〈νg|νe〉〈νe|ν¯g〉 ×〈ψνg |(Q−Qνe)(Q−Qνe)T |ψν¯g〉,(IV.76) α ′′ νe = ( 4〈ψνe|(Q−Qνe)(Q−Qνe)T |ψνe〉 〈ψνe|ψνe〉 )−1 , (IV.77) 92 〈ψνe |(Q−Qνe)(P − Pνe)T |ψνe〉 = ∑ νg ,ν¯g f ∗(νe, νg)f(νe, ν¯g)〈νg|νe〉〈νe|ν¯g〉 ×〈ψνg |(Q−Qνe)(P − Pνe)T |ψν¯g〉, (IV.78) α ′ νe = α ′′ νe〈ψνe|(Q−Qνe)(P − Pνe)T |ψνe〉+ 〈ψνe |(P − Pνe)(Q−Qνe)T |ψνe〉α ′′ νe 〈ψνe|ψνe〉 , (IV.79) γ ′′ νe = − 1 2 ln ( 〈ψνe|ψνe〉 √ 2Nbdet[α′′νe ] piNb ) , (IV.80) γ ′ νe = Im ln ∑ νg f(νe, νg)〈νe|νg〉〈Qνe|ψνg〉  . (IV.81) IV.2.3. Simulation Results We simulate tr-CARS signals on the model of 2D I2/Kr under conditions similar to those adopted in prior experiments [7]. The center wavelength is 580 nm for pump and probe pulses, and 600 nm for the Stokes pulse. The temporal width is 80 fs for the pump and probe, and 40 fs for the Stokes pulse. The arrival time of the Stoke pulse is the same as that of the pump pulse (t′2 = t ′ 1). The tr- CARS signal obtained using frozen-Gaussian variational FVB/GB only involves the biggest dipole moment contribution M (0) 123. Contributions from other portions of the tr-CARS dipole moment are relatively small. Absolute values of the wave-packet overlaps M (0) 123 (Eq. IV.57) and N (0) 213 (Eq. IV.59) at selected delay time t32 are shown in Figures IV.6. and IV.7. For N (0) 213, we have not included the emission before the arrival time of the third pulse. The upper limit of the third pulse propagator t′ is 93 0 0.01 0.02 0.03 0.04 -1 -0.5 0 0.5 1 1.5 2 |M 12 3(0 ) | t43/τs t32=0 t32=2.25τs t32=5.5τs t32=10.75τs Figure IV.6. The absolute value of tr-CARS dipole moment contribution M (0) 123(t32, t43), calculated according to Eq. IV.57. 0 0.0001 0.0002 0.0003 0.0004 0.0005 0 0.5 1 1.5 2 |N 2 13 (0) | t43/τs t32=0 t32=2.25τs t32=5.5τs t32=10.75τs Figure IV.7. The absolute value of tr-CARS dipole moment contribution N (0) 213(t32, t43), calculated according to Eq. IV.59. set to be infinity. The non-sequential dipole moment N213 is on the order of one percent of the magnitude of M123. Therefore we can neglect the contribution from N (0) 213. Contributions from the other two dipole moments M321 and N231 are small too (section IV.2.1), and we have the tr-CARS signal written as S ∝ ∫ ∞ t′3−3σ3 dt′4 ∣∣∣M (0)123(t′4)∣∣∣2 . (IV.82) We include 82 vibrational levels of the excited state and 8 levels of the ground 94 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0 5 10 15 20 25 30 35 40 45 trC AR S sig na l (a rb. un its ) t/τs Figure IV.8. Tr-CARS signal for the 2D I2/Kr model with frozen Gaussians, calculated according to Eq. (IV.82). state, and only keep couplings between any level ν and its adjacent levels ν¯ = ν ± 4 (see more in Appendix C). Figure IV.8. shows the tr-CARS signal calculated according to Eq. (IV.82). It presents an oscillation at a regular period of one ground state vibrational period (153 fs). For the first 15 periods, the signal decreases with a vibrational dephasing rate of about 0.01ps−1. Then the signal starts to grow a small amount after 20 periods. Using the same pulses, Kiviniemi and coworkers find a decoherence rate of 0.01ps−1 from their tr-CARS signal of hundreds of picoseconds (about 1000 ground-state periods) [7]. It is worth mentioning that the temperature for their experiment is 2.6 K. For our model of 2D I2/Kr, the temperature is 0 K, which a temperature can be approximated as if it is lower than the Debye temperature of a krypton crystal (71.9 K). Alternatively, the temperature with which the energy kBT is lower than one quantum of the highest-frequency bath- mode energy ~Ω47 can be considered to suit for our model of 0 K. It turns out that this critical temperature is 85.2 K. Longer time calculations of tr-CARS signals are 95 to be performed in order to watch the dynamics of wave packets in the electronic ground state. 96 CHAPTER V DISCUSSION AND FUTURE DIRECTIONS A mixed quantum/semiclassical theory for small-molecule dynamics in low- temperature solids has been framed, tested and applied to a realistic 2D model. It uses the system energy eigenkets as a basis and approximates the bath wave packets to have a Gaussian form, and therefore is termed the fixed vibrational basis/Gaussian bath approach (FVB/GB). We first tested FVB/GB on a model of bilinearly coupled harmonic oscillators and it exhibits promising excellence by predicting the exact results well for a certain amount of time. A variational treatment is proven to have some advantages over the locally quadratic approximation. In a 1D bath case, the variational method has been shown capable handling the high order interactions (with coupling terms proportional to qQ2) more accurately. In addition, the total population and energy are conserved in the variational FVB/GB theory. A 2D I2/Kr model is established, including 23 krypton atoms and one iodine molecule with periodic boundary conditions. The bath potential and system- bath interaction are written as a quartic function of the normal coordinates, and coefficients are found by carrying out a Monte–Carlo–like simulation. We implement FVB/GB numerically for this model, evolve the wave packets on the ground and 97 excited electronic potential energy surfaces, and obtain an approximated total wave function. The purity of the reduced system density matrix is calculated, and shows a long dephasing time of about hundreds of picoseconds. We also simulate the linear wave-packet interferometry and tr-CARS signals with frozen Gaussians. The behaviors of energies and coordinate trajectories of the 46th and 47th bath modes (assigned as zone boundary phonons) after excitation by the first laser pulse suggest a non-negligible coupling between ZBP and the iodine molecular vibration. The future direction of this work is to make good use of the total wave function and further investigate the role of the bath in quantum dynamics. 98 APPENDIX A FVB/GB EQUATIONS OF MOTION FOR A MULTIDIMENSIONAL BATH It is straightforward to derive the parameter equations of motion for the FVB/GB parameters in the case of a multidimensional bath via the procedure outlined in ref [17]. Using the auxiliary functions f(νν¯) = (Qν −Qν¯)T · αν¯ · (Qν −Qν¯) + P Tν¯ · (Qν −Qν¯) + γν¯ − γν , (A.1) g1(νν¯) = 2αν¯ · (Qν −Qν¯) + Pν¯ − Pν , (A.2) and g2(νν¯) = αν¯ − αν + i 2 g1(νν¯)g T 1 (νν¯), (A.3) these equations are found to be α˙ν = −2α2ν − 1 2 ∇∇T (ub + uνν)− ∑ ν¯ 6=ν eif(νν¯) { 1 2 ∇∇Tuνν¯ + i 2 [(∇uνν¯)gT1 (νν¯) + g1(νν¯)(∇Tuνν¯)] + iuνν¯g2(νν¯) } , (A.4) Q˙ν = Pν + (2α ′′ ν) −1 · ∑ ν¯ 6=ν e−f ′′ (νν¯) {uνν¯g′1(νν¯) cos f ′(νν¯) + [ ∇uνν¯ − uνν¯g′′1 (νν¯) ] sin f ′(νν¯) } , (A.5) 99 P˙ν = −∇(ub + uνν) + α′ν · (α ′′ ν) −1 · ∑ ν¯ 6=ν e−f ′′(νν¯) {uνν¯g′1(νν¯) cos f ′(νν¯) +[∇uνν¯ − uνν¯g′′1 (νν¯)] sin f ′(νν¯) } − ∑ ν¯ 6=ν e−f ′′(νν¯) { [∇uνν¯ − uνν¯g′′1 (νν¯)] cos f ′(νν¯) −uνν¯g′1(νν¯) sin f ′(νν¯)} , (A.6) γ˙′ν = 1 2 P Tν · Pν − εν − ub − uνν − Trα ′′ ν − ∑ ν¯ 6=ν e−f ′′ (νν¯)uνν¯ cos f ′(νν¯) +P Tν · (2α ′′ ν) −1 · ∑ ν¯ 6=ν e−f ′′ (νν¯) {uνν¯g′1(νν¯) cos f ′(νν¯) +[∇uνν¯ − uνν¯g′′1 (νν¯)] sin f ′(νν¯) } , (A.7) and γ˙ ′′ ν = Trα ′ ν − ∑ ν¯ 6=ν e−f ′′ (νν¯)uνν¯ sin f ′(νν¯). (A.8) The symbol ∇ represents a vector of derivatives with respect to the bath-coordinate variables, to be evaluated at values of Qν . Eqs. (A.4)–(A.8) reduce to Eqs. (II.10)– (II.14) in the main text for the case of a 1D bath. 100 APPENDIX B SOME INTEGRALS USED IN CALCULATING POTENTIAL ELEMENTS The integral of a Gaussian function times a polynomial can be calculated analytically. We have I = ∫ ∞ −∞ dQ exp[−QT · A ·Q+BT ·Q+ C] = √ piN det[A] exp [ 1 4 BT · A−1 ·B + C ] , (B.1) I2 = ∫ ∞ −∞ dQ exp[−QT · A ·Q+BT ·Q+ C]Qi = 1 2 (A−1B)iI, (B.2) I3 = ∫ ∞ −∞ dQ exp[−QT · A ·Q+BT ·Q+ C]QiQj = ( 1 2 A−1 + 1 4 A−1BBTA−1 ) ij I, (B.3) I4 = ∫ ∞ −∞ dQ exp[−QT · A ·Q+BT ·Q+ C]QiQjQk = ( 1 4 A−1ik (A −1B)j + 1 4 A−1jk (A −1B)i + 1 4 A−1ij (A −1B)k + 1 8 (A−1B)i(A−1B)j(A−1B)k ) I, (B.4) 101 and I5 = ∫ ∞ −∞ dQ exp[−QT · A ·Q+BT ·Q+ C]QiQjQkQl = ( 1 4 A−1ik A −1 jl + 1 4 A−1jk A −1 il + 1 4 A−1ij A −1 kl + 1 8 A−1il (A −1B)j(A−1B)k + 1 8 A−1jl (A −1B)i(A−1B)k + 1 8 A−1kl (A −1B)i(A−1B)j + 1 8 A−1ik (A −1B)j(A−1B)l + 1 8 A−1jk (A −1B)i(A−1B)l + 1 8 A−1ij (A −1B)k(A−1B)l + 1 16 (A−1B)i(A−1B)j(A−1B)k(A−1B)l ) I. (B.5) 102 APPENDIX C NUMERICAL STABILIZATION Because computers can only store limited digits for a real number, sometimes numerical calculations are not stable and we have to invoke some procedures in order to carry out accurate and stable simulations. In the equations of motion for the bath wave-packet parameters (Eqs. III.21–III.26), the coupling term 〈ψν |uνν¯ |ψν¯〉/〈ψν |ψν〉 is proportional to the ratio of amplitudes of two vibrational levels cν¯/cν . If one vibrational state has zero amplitude, cν = 0, then the coupling term goes to infinity and we can not perform the propagation. In practice, there are states with very small amplitudes that the computer can not determine accurately. For example, if the most populated state ν¯ has an amplitude of 0.8, and there is a state ν with an amplitude of 10−50, then this number 10−50 is not accurate because only 16 digits are kept during calculations for double-precision numbers. To avoid an unreasonably large coupling, which is proportional to cν¯/cν = 0.8× 1050 in this example, we set any amplitude smaller than 10−8 times the largest amplitude to be 10−8 × cmax so that the least population is 10−16 times the most population. In this example, the amplitude of the state ν is changed from 10−50 to 0.8× 10−8. For two states ν and ν¯ are far away from each other, e.g., ν = 0 and ν¯ = 20, their amplitudes can be very different and the ratio cν¯/cν can be a very large number. 103 Ideally, this may not cause a problem because the coupling is also proportional to uνν¯ , which is nearly zero for states apart from each other. Then, we have uνν¯ (about zero) times a big number cν¯/cν , and obtain a reasonable coupling. However, an element of the matrix uνν¯ (ν, ν¯ = 1, 2, · · · , 47) is not accurate when its absolute value is smaller than 10−16 times the maximum absolute value in this matrix. The actual number stored in a computer may not be small enough to eliminate the large coupling term. Accordingly, we decide to include only the couplings between near states, i.e. for state ν, we keep its couplings with ν¯ = ν ± 4 while couplings with ν¯ > ν + 4 and ν¯ < ν − 4 are forced to be zero. Lastly, we regularize the imaginary width parameter matrix α′′ν . Numerical difficulties in the propagation of Gaussian wave packets have been reported to arise in connection with the inversion of matrices [58]. When the time step in integration is not small enough, the determinant of α′′ν can be zero during propagation and therefore the inverse can not be calculated. To avoid this happening, we set a low limit for an eigen value of the imaginary width parameter matrix α′′ν , which is one twenty-fifth of its initial value. If one of the eigen values ever becomes smaller than its threshold, it is changed to be the smallest value allowed automatically. 104 REFERENCES CITED [1] P. R. Poulin and K. A. Nelson, Science 313, 1756 (2006). [2] M. B. Plenio and S. F. Huelga, New J. Phys. 10, 113019 (2008). [3] M. Fushitani, M. Bargheer, M. Gu¨hr, and N. Schwentner, Phys. Chem. Chem. Phys. 7, 3143 (2005). [4] Z. Bihary, M. Karavitis, and V. A. Apkarian, J. Chem. Phys. 120, 8144 (2004). [5] M. Karavitis, R. Zadoyan, and V. A. Apkarian, J. Chem. Phys. 114, 4131 (2001). [6] M. Karavitis and V. A. Apkarian, J. Chem. Phys. 120, 292 (2004). [7] T. Kiviniemi, J. Aumanen, P. Myllyperkio¨, V. A. Apkarian, and M. Pettersson, J. Chem. Phys. 123, 64509 (2005). [8] M. Gu¨hr and N. Schwentner, J. Chem. Phys. 123, 244506 (2005). [9] M. Karavitis, T. Kumada, I. U. Goldschleger, and V. A. Apkarian, Phys. Chem. Chem. Phys. 7, 791 (2005). [10] M. Gu¨hr, M. Bargheer, and N. Schwentner, Phys. Rev. Lett. 91, 085504 (2003). [11] M. Fushitani, N. Schwentner, M. Schro¨der, and O. Ku¨hn, J. Chem. Phys. 124, 024505 (2006). [12] H. Ibrahim, M. He´jjas, and N. Schwentner, Phys. Rev. Lett. 102, 088301 (2009). [13] M. Gu¨hr and N. Schwentner, Phys. Chem. Chem. Phys. 7, 760 (2005). [14] M. Fushitani, M. Bargheer, M. Gu¨hr, H. Ibrahim, and N. Schwentner, J. Phys. B: At. Mol. Opt. Phys. 41, 074013 (2008). [15] D. Segale and V. A. Apkarian, J. Chem. Phys. 135, 024203 (2011). [16] M. A. Rohrdanz and J. A. Cina, Mol. Phys. 104, 1161 (2006). 105 [17] C. T. Chapman and J. A. Cina, J. Chem. Phys. 127, 114502 (2007). [18] D. V. Shalashilin and I. Burghardt, J. Chem. Phys. 129, 084104 (2008). [19] E. J. Heller, J. Chem. Phys. 64, 63 (1976). [20] N. Linden, S. Popescu, A. J. Short, and A. Winter, Phys. Rev. E 79, 061103 (2009). [21] F. Mintert and E. J. Heller, Eur. Phys. Lett. 86, 50006 (2009). [22] V. Senekerimyan, I. Goldschleger, and V. A. Apkarian, J. Chem. Phys. 127, 214511 (2007). [23] I. U. Goldschleger, V. Senekerimyan, M. S. Krage, H. Seferyan, K. C. Janda, and V. A. Apkarian, J. Chem. Phys. 124, 204507 (2006). [24] T. Momose, M. Fushitani, and H. Hoshina, Int. Rev. Phys. Chem. 24, 533 (2005). [25] A. Borowski and O. Ku¨hn, J. Photochem. Photobiol. A 190, 169 (2007). [26] C. S. Guiang and R. E. Wyatt, J. Chem. Phys. 112, 3580 (2000). [27] K. Ohmori, Proc. Jpn. Acad., Ser. B 84, 167 (2008). [28] P. Huo, S. Bonella, L. Chen, and D. F. Coker, Chem. Phys. 370, 87 (2010). [29] Q. Shi and E. Geva, J. Chem. Phys. 129, 124505 (2008). [30] O. V. Prezhdo, Theor. Chem. Acc. 116, 206 (2006). [31] J. M. Riga, E. Fredj, and C. C. Martens, J. Chem. Phys. 124, 064506 (2006). [32] E. J. Heller, J. Chem. Phys. 62, 1544 (1975). [33] D. J. Tannor and E. J. Heller, J. Chem. Phys. 77, 202 (1982). [34] N. F. Scherer, R. J. Carlson, A. Matro, M. Du, A. J. Ruggiero, V. Romero- Rochin, J. A. Cina, G. R. Fleming, and S. A. Rice, J. Chem. Phys. 95, 1487 (1991). [35] N. F. Scherer, A. Matro, L. D. Ziegler, M. Du, R. J. Carlson, J. A. Cina, and G. R. Fleming, J. Chem. Phys. 96, 4180 (1992). [36] E. J. Heller, J. Chem. Phys. 68, 2066 (1978). [37] R. G. Gordon, Adv. Magn. Reson. 3, 1 (1968). 106 [38] K. Thompson and N. Makri, Phys. Rev. E 59, R4729 (1999). [39] X. Sun and W. H. Miller, J. Chem. Phys. 110, 6635 (1999). [40] G. Tao and W. H. Miller, J. Chem. Phys. 130, 184108 (2009). [41] G. Tao and W. H. Miller, J. Chem. Phys. 131, 224107 (2009). [42] R. Martinazzo, M. Nest, P. Saalfrank, and G. F. Tantardini, J. Chem. Phys. 125, 194102 (2006). [43] G. D. Billing, J. Chem. Phys. 110, 5526 (1999). [44] R. D. Coalson and M. Karplus, J. Chem. Phys. 93, 3919 (1990). [45] B. Lasorne, M. A. Robb, and G. A. Worth, Phys. Chem. Chem. Phys. 9, 3210 (2007). [46] C. T. Chapman, Ph.D. dissertation (University of Oregon, 2010). [47] D. Brinks, F. D. Stefani, F. Kulzer, R. Hildner, T. H. Taminiau, Y. Avlasevich, K. Mu¨llen, and N. F. van Hulst, Nature 465, 905 (2010). [48] T. S. Humble and J. A. Cina, J. Phys. Chem. B 110, 18879 (2006). [49] J. A. Cina, Annu. Rev. Phys. Chem. 59, 319 (2008). [50] P. Kramer and M. Saraceno, Lecture Notes in Physics: Geometry of the Time- Dependent Variational Principle in Quantum Mechanics (Springer-Verlag, Berlin, 1981). [51] A. K. Kerman and S. E. Koonin, Ann. Phys. 100, 332 (1976). [52] T. Kiljunen, M. Bargheer, M. Gu¨hr, and N. Schwentner, Phys. Chem. Chem. Phys. 6, 2185 (2004). [53] R. Zadoyan, Z. Li, C. C. Martens, and V. A. Apkarian, J. Chem. Phys. 101, 6648 (1994). [54] D. Frenkel and B. Smit, Understanding Molecular Simulation (Acamemic Press, San Diego, 1996). [55] M. Bonfanti, G. F. Tantardini, K. H. Hughes, R. Martinazzo, and I. Burghardt, J. Phys. Chem. 116, 11406 (2012). [56] Z. Bihary, R. Zadoyan, M. Karavitis, and V. A. Apkarian, J. Chem. Phys. 120, 7576 (2004). 107 [57] J. D. Jackson, Classical Electrodynamics (New York, Wiley, 1962). [58] I. Burghardt, M. Nest, and G. A. Worth, J. Chem. Phys. 119, 5364 (2003). 108