Citation |

- Permanent Link:
- http://ufdc.ufl.edu/AA00047454/00001
## Material Information- Title:
- Photodynamics of extended polyatomic systems
- Creator:
- Srivastava, Deepak, 1957-
- Publication Date:
- 1988
- Language:
- English
- Physical Description:
- vii, 164 leaves : ill. ; 28 cm.
## Subjects- Subjects / Keywords:
- Approximation ( jstor )
Coordinate systems ( jstor ) Electric fields ( jstor ) Electronics ( jstor ) Fourier transformations ( jstor ) Lighting ( jstor ) Molecular interactions ( jstor ) Photons ( jstor ) Potential energy ( jstor ) Vibrational states ( jstor ) Dissertations, Academic -- Physics -- UF Lasers ( lcsh ) Molecular dynamics ( lcsh ) Physics thesis Ph. D Quantum theory ( lcsh ) - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt )
## Notes- Thesis:
- Thesis (Ph. D.)--University of Florida, 1988.
- Bibliography:
- Includes bibliographical references.
- General Note:
- Typescript.
- General Note:
- Vita.
- Statement of Responsibility:
- by Deepak Srivastava.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
- Resource Identifier:
- 025018456 ( ALEPH )
20233917 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

PHOTODYNAMICS OF EXTENDED POLYATOMIC SYSTEMS By DEEPAK SRIVASTAVA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1988 TO MY WIFE ACKNOWLEDGEMENTS I would like to express my thanks and appreciation to some of the many people who have encouraged and supported me in this work. The foremost amongst them is my advisor Dr. David A. Micha, who constantly provided the guiding thought throughout the course of this work and also provided a part of the funding for my financial support. I am also thankful to many graduate students, postdoctoral associates, and faculty members at the QTP who have contributed in some way or the other towards increasing my knowledge in the field of theoretical and computational chemical physics. Finally, I would like to express my appreciation for all the help I have received from the secretaries and the staff at the QTP and at the Physics Department that has made my stay at the University of Florida a pleasant experience. iii TABLE OF CONTENTS page ACKNOWLEDGMENTS............................... iii ABSTRACT...................................... . vi CHAPTERS 1 INTRODUCTION. .................................... 1 2 COLLISIONAL TIME-CORRELATION FUNCTION APPROACH TO LIGHT-MOLECULE INTERACTIONS.... 8 2.1 Transition Ratesfrom Time-Correlation Functions .......... 9 2.2 Radiation Field TCFs for Thermal Light............... ...13 2.1 Radiation Field TCFs for Laser Light ................ 17 3 MOLECULAR TCFs IN THE ELECTRONIC ADIABATIC REPRESENTATION........ ...22 3.1 Electronic Adiabatic Representation for Molecular TCFs . . .. 24 3.2 Complex Time TCFs and their Fourier Transforms........ ... 28 3.3 Properties of the Complex Time TCFs and their Fourier Transforms .................. . 34 4 FACTORIZATION OF TCFs INTO PRIMARY AND SECONDARY REGION TCFs......... ... 39 4.1 Sufficient Conditions for the Factorization of the Amplitude . . 40 4.2 Norm Conserving Time-Dependent Variational Principle(TDVP) 43 4.3 Time-Dependent Self-Consistent Field (TDSCF) Approximation 45 4.4 Primary and Secondary Region Factorization of TCFs ...... 50 5 TIME-DEPENDENT EFFECTIVE FREQUENCY OPERATORS 52 5.1 Intramolecular Processes: Strongly Coupled Primary and Secondary Regions ...... ...53 5.2 Intermolecular Processes: Weak Coupling Limit ......... . 62 5.3 Relation to the Generalized Langevin Dynamics Approach. . 64 5.4 Treatment of Many Degrees of Freedom in the Primary Region Dynamics ..................... 69 iv 6 PATH-INTEGRAL APPROACH TO THE PRIMARY REGION DYNAMICS ............ ...73 6.1 Numerical Matrix Multiplication (NMM) Method and its Shortcomings. ............................. 75 6.2 The Complex Time Propagator for the Primary Region Dynamics ................ . 78 6.3 A More Efficient Transformation from Real to Complex Time . 87 7 APPLICATIONS AND RESULTS. ........................ 91 7.1 Electronic Absorption Spectra of a Model Two-Potential Probelem ................... 92 7.2 Photodissociation of CH3I from its Ground Vibrational State. . 100 7.3 Photodissociation of CH3I from Initially Excited Vibrational States. ............... 120 8 DISCUSSION AND CONCLUSIONS................. 137 APPENDICES A SECOND ORDER TRANSITION RATES. ................. 152 B NUMERICAL ANALYTIC CONTINUATION ............... 155 C SCF EIGENVALUES AND EIGENFUNCTIONS........... ..157 BIBILIOGRAPHY................................... 159 BIOGRAPHICAL SKETCH. ................................ 164 V Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy PHOTODYNAMICS OF EXTENDED POLYATOMIC SYSTEMS by DEEPAK SRIVASTAVA December 1988 Chairman: Prof. David A. Micha Major Department: Physics A general time-dependent quantum mechanical approach to the interaction of visible and UV light with extended polyatomic systems is presented in this work. The interaction of light with a large polyatomic system is treated as a two-step process: a photon absorption excites electronic transitions in the target system, and this is followed by a dynamical evolution of the system on the excited potential energy surface. Transition rates for photon-molecule interaction processes are described in terms of time-correlation functions (TCF) of transition operators. These are written as a product of TCFs of the target and of the light source. The time evolutions in these two TCFs are carried out independently, and the statistical properties of the light source are explicitly incorporated in the formalism. The procedure is developed for resonance absorption-emission and for scattering of thermal, chaotic and coherent light. The time evolution in a large polyatomic system is treated within a molecular TCF approach. For a general two-surface electronic excitation problem, a transformation of these molecular TCFs from the real arguments to complex time is vi described so that their computation by a numerical path-integral method become feasible. The complex time dependence in these TCFs is described by the dynamical evolution of an amplitude under the influence of a frequency operator. A time-dependent self-consistent field (TDSCF) approximation is used for large polyatomic system to factor the molecular TCFs into primary and secondary region TCFs. The primary and secondary regions are modeled by considering general primary motion coupled to many harmonic degrees of freedom in the secondary region. The strong coupling between the primary and the secondary regions has a general dependence on the variables of the primary region, whereas it has linear and bilinear terms in the variables of the secondary region. The weak coupling limit is considered by dropping the bilinear terms in the coupling. The complex time propagators for the secondary region dynamics for both of these cases are constructed analytically, while an efficient path-integral method is developed for the dynamics of the primary region. As applications, the path-integral method for the primary region dynamics is used for computing the electronic absorption spectra of a model two-potential problem, and the TDSCF approach to the dynamics of large systems is used for studying the photodissociation dynamics of CH3I from ground and excited vibrational states of the ground potential energy surface. The cross-sections for both of these cases are compared with other available works with good results. Extensions of the method to the other extended systems, with many degrees of freedom, and involving tunneling and barrier crossing in the primary region dynamics, are briefly discussed. Vii CHAPTER 1 INTRODUCTION Theories of light absorption-emission, elastic and Raman scattering are described by time-dependent perturbation expressions in terms of transition dipole moments, and polarizability tensors [Kramers and Heisenberg, 25], respectively. The interaction of atoms and molecules with quantized radiation field was first introduced by Dirac [Dirac, 27]. Since then many developments have been made along the same lines and have included higher order corrections in these theories to increase accuracies [Herzberg, 50]. The basic approach of these developments is to concentrate on small atoms, accurately compute all of their eigenvalues and eigenstates, and use them in a summation according to the Fermi Golden Rule [Messiah, 62] for computing the total cross-section for the required photointeraction process. In case of polyatomics, using the Born-Oppenheimer approximation [Balint-Kurti, 75], the electronic part of the Schrodinger equation is first solved for fixed nuclear positions to obtain electronic potential energy surfaces; these are then used for the nuclear part of the problem. Given electronic potential energy surfaces for a small polyatomic system, one can still obtain all the low-lying vibrational and rotational eigenstates and eigenvalues to calculate cross-section. However, cross-sections processes which involve rearrangement in the molecule due to the absorption of light are not easily described by the method of summing over all the accessible quantum states of the system. This is mainly because one needs exact scattering wavefunctions of many channels in the excited state of the system, and the number of such available channels can be very large even for small polyatomics. Still, the time1 2 independent quantum mechanical method, based on summation over all the state to state information, and known as the stationary coupled channel (CC) approach, has been used successfully to produce total cross-sections for light-molecule interaction processes involving rearrangement in the target [Balint-Kurti and Shapiro, 81; Heather and Light, 83]. However, as the size of the system starts to increase, methods such as these involve a larger and larger number of available quantum states of the system, and therefore soon become too cumbersome to be of any practical use. Alternatively, efforts have been made to change the approach from the timeindependent quantum mechanical theories, which involve summations over all of the stationary states of the system, to the time-dependent quantum mechanical descriptions which are based on the dynamical evolution of the system from given initial conditions to the required final conditions [Lee and Heller, 79; Berens and Wilson, 81; Thirumalai and Berne, 83; Mukamel et al., 85; Srivastava and Micha, 87; Coalson, 87]. These methods do not require summations over all the states of the system, and therefore are better suited for dealing with large molecular systems. The basic feature of these approaches is to start from the expression for the total cross-section in the energy domain, use the Fourier transform definition of the energy conserving delta function to make a transformation to the time domain, and use the time-dependent Schrodinger equation for the dynamical evolution of the system between given initial and required final conditions. In recent years, many advances have been made in solving the time-dependent Schrodinger equation for nuclear vibrational and rotational motion with given initial conditions for an arbitrary general potential in many dimensions [Kosloff, 88 for a recent review]. They can be catagorized into the following broad groups: (a) semi-classical methods [Marcus, 72; Miller, 74; Mukamel and Jortner, 76; Micha, 83; Billing and Jolicard, 3 84; Gray and Child, 84; Henriksen, 85] are generally based on running classical trajectories for the nuclear motions, whereas the electronic transitions are treated quantum mechanically; (b) approximate quantum mechanical approaches based on wavepacket dynamics [Heller, 75, 81; Agrawal and Raff, 82; Mowrey and Kouri 85; Sawada and Metiu, 86], where the mean position and the spread of the wavepacket follow a canonical set of equations which are derived from a timedependent variational principle [Heller, 81]; (c) straightforward integration of the time-dependent Schrodinger equation by discretizing the phase space and using a Fast Fourier Transform (FFT) algorithm [Feit, Fleck and Steiger, 82; Kosloff and Kosloff, 83; Hellsing, Nitzan and Metiu, 86; Kosloff, 88]; and (d) path integral approaches [Thirumalai and Beme, 83; Miller, Schwartz and Tromp, 83; Behrman et al., 83; Coalson, Freeman and Doll, 85] which provide complex time propagators between given initial and final positions in the coordinate space. The present developments in the time-dependent descriptions still have the following problems which must be addressed. 1. By restricting the radiation field part of the interaction Hamiltonian to simple monochromatic circular wave forms, these approaches do not explicitly incorporate the statistical and dynamical properties of the light source in the formalism. 2. Semiclassical methods for solving the time-dependent Schrodinger equation are based on classical trajectories which do not sample the non-classical regions of the potential energy surfaces. 3. Approximate quantum mechanical techniques based on wavepacket dynamics require local quadratic forms for potentials, hence they exclude extensive delocalization of wavepackets. 4 4. The numerical path-integral methods, for a potential of arbitrary functional form, are not useful in the calculation of the required real time dynamics because of the oscillatory nature of the real time propagator. 5. Straightforward integration of the time-dependent Schrodinger equation by discretizing the phase space and using FFT for propagation are practical only for problems with few spatial dimensions. The time-correlation function approach to thermodynamic and electromagnetic response properties of extended molecular systems is well known in statistical mechanics [Forster, 75; McQuarrie, 76; Lovesey, 80]. Recently, there have also been efforts to describe molecular processes occuring in a target following its collision with a projectile, using a collisional time-correlation function (TCF) approach [Micha, 77, 79, 81, 86]. Since, as before, this approach does not require summations over all the states of the target, it is useful in dealing with many degrees of freedom in extended molecular targets [Vilallonga and Micha, 87; Micha and Parra, 88]. In this work we present a formalism and applications of a timedependent quantum mechanical approach to light-molecule interactions in large systems with the following desirable features. 1. We treat the radiation field and the molecular target of the interaction Hamiltonian on equal footing so that, at any stage in the formalism, explicit statistical and dynamical properties of the light source are easily incorporated. 2. Within a mean field approximation, we separate the degrees of freedom of an extended molecular system into a primary region ( containing a few important degrees of freedom along which all significant dynamical events take place) and a secondary region (with all of the remaining degrees of freedom). 3. For a general potential of arbitrary functional form in the primary region, we use 5 a complex time path-integral method which can incorporate important quantum effects such as tunneling and barrier crossing in this region. 4. The computational efforts for the method do not increase quadratically with an increase in the range and the dimensionality of the potential. In chapter 2 we describe transition rates for general light-molecule interaction processes in terms of the TCF of the Lippman-Schwinger t operator. In the dipole approximation to the light molecule interaction, the t operator is expanded into terms containing factors for the radiation field and the molecular target. The transition rates for all first order processes such as resonance absorption and emission are written in terms of the molecular dipole TCF of the target, and the electric field TCF of the light source. For all second order processes, such as multiphoton absorption-emission, elastic and Raman scattering, the transition rates are written in terms of the TCFs of the molecular polarizability tensors. The statistical properties of the light source are incorporated in the formalism, and specific cases of thermal, chaotic and coherent light are discussed. In chapter 3 we introduce electronic potential energy surfaces by using the Born-Oppenheimer approximation to separate electronic and nuclear motions in a molecular target, and write a general molecular TCF in the electronic adiabatic representation. We describe a transformation from real to complex time of these TCFs, and discuss some of their general properties. The complex time dependence in the TCF, for a general light-molecule interaction process involving electronic transitions in the target, is shown to be equivalent to the dynamical evolution of an amplitude (equal to the initial eigenstate times the transition dipole) induced by a frequency operator (an operator depends upon the energy difference between the upper and the lower potential energy surfaces divided by h). 6 In chapter 4 we describe the time evolution of the amplitude introduced in the previous chapter through a norm conserving time-dependent variational principle (TDVP). Considering fast and localized electronic tranisitions in a large molecular target, we use a time-dependent self-consistent field (TDSCF) approximation to factor the molecular dipole TCF of the target into self-consistent primary and secondary region TCFs. In chapter 5 we describe time-dependent effective frequency operators for many intra- and intermolecular photointeraction processes by considering an arbitrary general motion in the primary region dynamics coupled linearly and bilinearly to many degrees of freedom in a secondary region dynamics undergoing harmonic motions. We show a direct relationship between our TDSCF approach and the phenomenological Generalised Langevin Equation (GLE) [Mori, 65] approach to the stochastic dynamics in large systems. We explicitly prove a FluctuationDissipation Relationship (FDR) between the fluctuation force and the memory kernel in our TDSCF approach for linearly coupled primary and secondary region dynamics. We also develop the formalism for many degrees of freedom in the primary region dynamics coupled linearly and bilinearly to many harmonic degrees of freedom in the secondary region dynamics. In chapter 6, we describe details of an efficient path-integral method for the primary region dynamics. Using a discretized coordinate space, the computational method does not depend upon any specific functional form of the primary region potential and also avoids problems of other methods which scale quadratically with respect to the range and the dimensionality of the potential. In chapter 7 we describe applications and results of the TCF approach to the interaction of visible and UV light with extended molecular systems. We first de- 7 scribe an application to the electronic absorption spectra of a model one-dimensional problem for two-potential energy surfaces. The complex time computations for the TCF in a single spatial dimension, using the complex time numerical path-integral method, are performed so as to provide a check of the accuracy of the numerical method for the primary region dynamics. We study the photodissociation dynamics of a realistic polyatomic system (CH3I) in an electronically excited state within a TDSCF approximation by factoring the TCF for the process into self-consistent primary and secondary region TCFs. The resulting total TCF has been used to compute the total cross-sections for the photodissociation of CH3I, which are then compared with results from other methods with good results. We have also extended the method to the computations of total cross-sections for many vibrationally excited states of CH3I on the ground potential energy surface. Finally, in chapter 8 we discuss advantages and distinctive features of our approach and mention possibilities for further improvements and future extensions of our work. CHAPTER 2 COLLISIONAL TIME-CORRELATION FUNCTION APPROACH TO LIGHT-MOLECULE INTERACTIONS We treat the interaction of light with a molecular system as a photon scattering process through a collisional TCF approach [Srivastava and Micha, 87]. Starting from the transition rates for a collision process in terms of the Lippmann-Schwinger transition operator t [Messiah, 62] in section 2.1, we derive the transition rates for the interaction of light with a molecular system in terms of the TCF of these T operators. In the dipole approximation to the light-molecule interaction, the T operator is written as an expansion with terms containing radiation field and molecular system factors. For first order processes such as resonance absorption and emission, transition rates are derived in terms of molecular dipole TCF of the target and the electric field dipole TCF of the light source. For second order processes such as elastic and inelastic (Raman) scattering due to polarization fluctuations in the target, transition rates are written in terms of the polarizability tensor TCFs of the molecular system and the polarizability tensor TCFs of the light source. In section 2.2 we treat the interaction of molecules with thermal light, which is given in the photon number representation by a Bose-Einstein distribution of the initial number of photons in a particular mode. Section 2.3 considers the interaction of a molecule with a single mode laser. The laser light source is treated by considering the radiation field to be in the coherent state representation. The first order transition rates for laser-molecule interaction processes are derived when a) the laser is operating below the threshold (chaotic light with a poisson distribution 8 9 for the mean number of photons), and b) the laser is operating above the threshold (coherent light with a statistical distribution for the phase and the mean number of photons given by the experimental conditions). 2.1. Transition Rates from Time-Correlation Functions We consider the interaction of visible and UV light with a molecular system. The treatment is general so that the molecular system can be as large as a protein molecule, a polymer chain, an adsorbate on a surface, or a solid surface. We denote with fSm, HR, and H1 the Hamiltonians of the molecular system, the radiation field, and their interaction, respectively. The total Hamiltonian is expressed as f = HM + HR + HI- (2.1) Eigenstates j) and I r), corresponding to the noninteracting molecular system Hamiltonian Hm and the radiation field Hamiltonian HR, respectively, are given by HMj) Eu1j) (2.2) HRjr)= Err). We treat the interaction of light with a molecular system as photon scattering, during which the radiation field goes from an initial state I r) to a final state I r/). Later on we shall describe various processes like absorption or Raman scattering by specifying I r) and I r/). Transition operators t+ describing the scattering satisfy the Lippmann-Schwinger equation [Messiah, 62] 11 + H^ I T+, (2.3) E+i-HI-HR where E -- 0 is a limit to be taken after matrix elements of these operators have been calculated. Rates of transition R, averaged over the statistical distribution w1 10 of the initial target states, wr of the initial radiation states, and summed over final states of the target and the radiation, Ij/) and Ir/), respectively, are expressed as R = + Euj [ w,|(r'j'| t+|r j)|26(E - E') (2.4) i rr' where E = Er + E. is the total energy, and the Dirac b function imposes energy conservation. By using the Fourier transform of the 6- function 6(E - E') - 2 ] e(E-E')t dt, (2.5) 6(E- E) =27rh _-0 (-5 and the completeness relation Er)('r'I = 1, (2.6) it r Eq.(2.3) can be written as [Micha, 81] R = dt ((t+(t)tt+(O)))M,R (2.7) where the time-dependence in the T+ is given by t+ exp _-Hot) t+exp(+ ot) (2.8) and Ho = HM + HR, is the total Hamiltonian for the molecular system and the radiation field without interaction. The double bracket < ... >M,R indicates the quantal average over the initial state IJ) of the molecular system with a statistical distribution wj, and over the initial state Ir) of the radiation field with a statistical distribution wr. So far we have not specified the interaction or the coupling Hamiltonian H1. Treating the interaction Hamiltonian in the electric dipole approximation, H1 = -E.D (.9 (2.9) 11 where f is the dipole vector of the molecular system, and f is the electric field vector of the light source. The transition operator T of Eq.(2.2) is written up to the second order, with an implicit + sign, as = Z b~+ E~, +^ID 77 -HM- H (2.10) 7( E + i H- m - HR where ((, r, () = (x, y, z), and the time-dependence of the t operator is given by Eq.(2.8). Since HAM and HR commute with each other, it is straightforward to take matrix elements of the first term in Eq.(2. 10) and to separate them into radiation field and molecular system factors. In the second term we use the following expression 1 1 +O a b = - ( dA 2 A2 (2.11) a + b 7r _-00 a +2 b2 + A2 valid also for complex values of a and b to write 1 ih +o MM MR E H1 - +00 du .- u2 - (2.12) E +ie-fHM -H H 1 -o 12 _i h2U2 M12 - h2U2 where MM = EM + Zie - fM (2.13a) MR = E R + ZEM - HR (2.13b) E = EM + ER (2.13c) +E (2.13d) and we have replaced A = zhu. Using these expressions in the second term of Eq.(2.10), to factor it into a radiation field part and a molecular system part, and the result in Eq.(2.6), we write the transition rates R up to the second order as 12 R = Ed ((b ,(t)b))) ((',(O ,()M))R - + du j dt ((b (t) ()))A(u.14) -oo -oo (ET ,(t)i(u, 0)))R - 1 +oo j+ u' d - du' dt (( ',,(u', t)bC(o)))m .4 7r ~ ,O 7 (2.14a) {{ (U', t) kC(0))) R h2 00 d +0 du' +0dt 7r-o -K1 ~ oo f -o u1(',t) '?((U,0))) M ( ( ',, , t ) ,7((u, O)) R where the polarizability tensor operator for the molecular system is defined by AlmU 2,g(u, 0) = D, Jb fhU bc (2.15a) the polarizability tensor operator for the radiation field is defined by ,, (U, 0)= MR E2C ,(2.15b) and the time-dependence in these operators are given for example by ic (u, t) = eC-HM iTiC (u, 0) e+IMI (2.16) The expression for R has four terms, each of them corresponding to different photointeraction processes. The first term describes first order processes such as resonance dipolar absorption and emission of light. The second and the third terms are cross-terms the contribution of which will depend upon the statistical properties 13 of the light source used in the experiment. For thermal light in the photon number representation, they become zero and do not contribute to the overall transition rates. In the later sections, we shall, however, see that the contribution of these terms cannot be ignored if the light source used is a coherent light. The fourth term describes second order processes such as elastic and Raman scattering mediated by polarization fluctuations in the target. The transition rates R expressed in Eq.(2.14), as opposed to the standard expressions [Dirac, 47; Heitler, 54], are specially useful for the treatment of lightmolecule interactions in large systems because they do not involve summations over excited states of the molecular system. Furthermore, Eq.(2.10), which is used to separate in the transition operator the radiation field and molecular system variables, can also be used to derive a similar separation of the transition operator to order higher than the second. Thus, a similar formulation of the transition rates for higher order photointeraction processes can be systematically derived. Lastly, the main purpose of separating the TCF of the t operator into the TCF of the radiation field and the TCF of the molecular system is to uncouple the time evolution of the molecular system from that of the radiation field. This has a twofold advantage: firstly, one can use time-dependent quantum mechanical techniques to describe the time evolution of the molecular system without referring to the time evolution of the radiation field, and secondly, we can easily incorporate the statistical properties of the light source into the dynamics. 2.2. Radiation Field TCFs For Thermal Light When the light source is maintained at a temperature T, it is described by the Bose-Einstein statistical distribution function W,(w8, T), of the number of photons 14 in a mode s with energy hL,. Then it is easier to work in the photon number representation of the radiation field. The electric field operator in Eq.(2.14) is decomposed as Eg (1) = E+ ) (t) + E- (0 (2.17) where +( t)2 eg &ge'i (2.18) k and E t)= + t (2.19) where k is the photon wavevector, eg is the ( component of the polarization vector, and cii the annihiliation operator of a photon of energy hwo = hkc and polarization (. In the photon number representation, a general radiation field state with several modes s is written as I{(k,,)N-}) where N, is the number of photons in mode s; then E +) ( , 2oV Ceia't,j {... (ks .) N.-1 ...) (2.20) Using Eqs.(2.17-2.20) in Eq.(2.14), one can write the transition rates for first and second order processes. In the photon number representation, the second and third terms of Eq.(2.14) do not contribute, because they do not conserve the number of photons in the initial and the final radiation state. Therefore the transition rate is written as 15 R = R11 + R22 (2.21) where R11 is the transition rate for the first order processes, in which a single electric field operator at an initial time t = 0, is correlated to its adjoint at a later time t. Similarly, R22 is the transition rate for second order processes, where the time-dependent polarization tensor operator, gC(u, 0), which is second order in the electric field operators of Eq.(2.17), is correlated to its adjoint at a later time t. In the same notation R12 and R21 are the cross-terms, the contributions of which are zero in the photon number representation. The transition rates R11 for first order processes are given by +00 R11= dt M(t R ,(t) ,(2.22) where fA (t) = KK, (t)t b (0))) (2.23) is the molecular dipole TCF of the extended system, and f, (t)= E) (t) (0))R + E (t) (0) is the electric field TCF of the light source. The first term in Eq.(2.24) corresponds to resonance absorption, whereas the second term corresponds to spontaneous and stimulated emissions. Other cross-terms in R11 do not contribute because they do not conserve the number of photons. Using Eqs.(2.18-2.20) in the above equation we write 16 fj, (t) 2 o 6 WV [ Nsw e".r C, kV + (N, + 1)wT e-SW't] (2.25) and the transition rate RI, is given by 2RV W VE8 I dt [N, cos(wkt) (2.26) 2Re((bD (t)t b (0)))M + Reiwisi((D( (t)t D (0)))A M}] where the statistical probability distribution function for a thermal light WV' is the Bose-Einstein distribution for N, photons with wavevector T, and frequency wj, and the summation in Eq.(2.26) over the number of photons can be carried out analytically. We have also used in Eq.(2.26); ((D (-t) D(0)))m = ((D ( t)D (O)))*, to ensure that the transition rate R1 remains real. For the second order processes, we write the transition rates R22 as h2 +00 +00 +00 R22= 72 duf du Jdt 17C ''VQ-' -00 -0 f , (U' U, t ) ff7 17 (U' U, t) (2.27) where f"7, (u', u,) , KKf (U t)t t (u, 0) (2.28) is the molecular polarizability tensor TCF of the target, and fA ,1( u', u, t) , (u', t)t e(u, 0) ))R (2.29) 17 is the electric field polarizability tensor TCF of the light source. Using Eqs.(2.172.20) in Eq.(2.29), the electric field polarizability tensor TCF of the light source can also be further decomposed into terms representing elastic and Raman scattering [Appendix A]. These terms when substituted back in Eq.(2.27), will give transition rates for second order processes when thermal light is interacting with a large molecular system. Equations (2.22) and (2.27) for the transition rates R11 and R22 are computationally useful, because they do not involve the summations over the excited states of the target, and the corresponding molecular TCF can be evaluated by time-dependent techniques. 2.3. Radiation Field TCFs For Laser Light When the light source is a laser, the statistical distribution of its electric field factorizes into a distribution for the electric field amplitude, times a distribution for its phase [Louisell, 73; Sargent, Scully and Lamb, 74; Loudon, 83]. If the laser is operated below the threshold for coherent emission, the amplitude fluctuations of a given mode, which are proportional to the mean number of photons in that mode, have an exponentially decaying distribution about the most probable amplitude (which occurs at zero), whereas fluctuations in the phase, caused by spontaneous emissions, are quite large and therefore the phase can have any value between 0 and 27r. If the laser is operated above the threshold, both the amplitude and the phase have sharply peaked distribution about their most probable values, which are given by experimental conditions. For such a light source it is necessary to include in the description the statistical distributions of both the amplitude and the phase of the electric field, and it is advantageous to work in the coherent state representation 18 [Klauder and Sudarshan, 68; Loudon, 83]. A general radiation state in the coherent state representation is written as I {a,}), with ak, {a,})= a I {a, }) where the eigenvalue a, is a complex number whose square modulus gives the mean number of photons in the s mode of the laser. Using equations (2.17-2.20) for the creation and annihiliation terms of the electric field operator, and using a single mode coherent states I as), we write 2 +) (t) as) = i E( 6 se-w aslas) (2.30) The radiation field TCF of Eq.(2.22), in the coherent state representation, can be written as f () t(t) t j ())R (Etj) (t) Et - (0))) Etl +t) 0)R + )R + ET (t) E (0))) + E (t) (0) 2.31) Unlike f R (t) in the photon number representation, cross-terms in Eq.(2.31) are not zero, because we do not need to conserve the number of photons in the initial and the final radiation states I r) and Ir'), respectively. As we shall see very shortly, however, it is still possible to have a laser in which the statistical distribution of the phase is such that it again can make the contributions of these cross-terms equal to zero. Using equations (2.18), (2.19) and (2.30) in Eq.(2.31), we write f R (t) in a single mode coherent state representation as f/1(t) = Ja f d2a,P (a,) x [-a2 e + (ja~l + 1)eCWri + a e +ius - c4 eaC+iwt] (2.32) 19 where P(a,), a real statistical weight for the complex amplitude a, is given by experimental conditions. The two dimensional integral f d2a, [...] extends over the entire complex plane of a. Light sources of experimental significance have different statistical weights P(a8). In the first and the fourth terms of Eq.(2.32), we 2, a2 have checked that the statistical averages of a , and .f are real for several weight factors P(a,) given in the literature [Schubert and Wilhelmi, 86]. Therefore, the radiation field TCF fJ,(t) in Eq.(2.32) satisfies the relation f ,(t) f (-t), and hence, the resulting transition rates of experimental significance, are real from Eq.(2.32). As examples of laser sources, with given statistical weights, we consider the radiation field TCFs for chaotic and coherent light. To compare Eq.(2.32) with the equivalant expression in the photon number representation in Eq.(2.25), we have first considered the case of chaotic light, which is the same as that of a laser operating below the threshold. In this case the phase of the light source is randomly distributed between 0 and 27r, therefore the weight factor P(a,) does not depend on the phase of a, and hence, the contributions from the first and the fourth terms of the integrand in Eq.(2.32) are zero. Thus, below the threshold fj,(t) takes a simpler form fg (t) = 2-" ' + ((N,) + 1)e i-'] (2.33) which agrees with Eq.(2.25) for a single mode, insofar as the average is (N,) = IfdIa'I |al P (IaI) (2.34) 0 20 Substituting Eqs.(2.32) and (2.33) in Eq.(2.21), we can write the first order transition rate R11 for the interaction of a single mode laser with a molecular system, when the laser is operating below the threshold, as 00 Richa*tc = , e E2 dt[(N,) cos(Lq~t) x 2Re(((bt (t)t bt (0)))M ) (2.35) Re(e-i"';.i(Ab (t 4g (0)))M)] If the laser is operating above the threshold, the rate of transition involves instead the statistical distribution P(a,), which depends on the amplitude I a, as well on the phase , = Arg(a,), and is given by +00 Coherent = 2hfV 2 J dtf d2a, p (a,) x [2 las12 {cos(W t)} - cos(wkt - 2p,) x ((A (W) bt (0))) M (2.36) + e-iw ((b (t)t b (0))) M] This can be rewritten as Rcoheent = RChaotic 2- E Jdt C 0 x [A, cos(wk t)Re{((bD (t)t bE (0)))M} (2.37a) + iB, sin(Lq t)Im{((b (t)t b (0)))M}] 21 with A, = d2a, P (a,) Ia,12 cos (2,p) (2.37b) B, = d2a, P (a,)la,12 sin (2o,) (2.37c) showing an additional rate due to coherence. The term containing B, can be shown to be zero for statistical distributions of experimental interests [Schubert and Wilhelmi, 86]. CHAPTER 3 MOLECULAR TCFs IN THE ELECTRONIC ADIABATIC REPRESENTATION In the previous chapter, transition rates for all first and second order photointeraction processes are expressed as integrals over time of the product of the molecular TCFs of the target and the radiation field TCFs of the light source. The time-dependence in the radiation field TCFs is expressed for many physically interesting light sources in terms of exponentials, oscillating with the frequency of the light source used. Therefore, most of the transition rates are finally written as the Fourier transforms of the molecular TCFs of the target. In this chapter, we describe general features of these molecular TCFs and their Fourier transforms. Indeed, the task of a general formalism should normally include describing a light source within a frequency range that extends from the far infrared region on one side, to the ultravoilet region on the other. A complete discussion of such an extended frequency range, however, is beyond the scope of this work, and we will describe mainly the interaction of visible or UV light with large molecular systems. As shown in table 3-1, the light source in this frequency range excites nuclear vibrations, accompanied by electronic transitions between different states of the system [McQuarrie, 83]. Thus, our work involves mainly the electronic excitations of the molecular system, followed by a time-dependent evolution according to the upper potential energy surface, which in some cases can also cause the dissociation of the system. We have also considered nuclear motions in the ground electronic state, which occur generally when the molecule is in a thermal bath. 22 23 Table 3-1 The Electromagnetic Spectrum And The Corresponding Molecular Processes Region Wave # (/cm) Process Microwave 0.033 - 3.300 Rotations of Polyatomic Molecules Far Infrared 3.30 -330 Rotations of small molecules Infrared 330 - 3300 Vibrations of flexible bonds Visible and Ultravoilet 3300 - 330000 Electronic transitions with vibrations of flexible bonds [From, "Quantum Chemistry" by D. A. McQuarrie; 83, pp.438] In section 3.1, we introduce electronic potential energy surfaces for the motion of the nuclei, and write a general molecular TCF of the target in the adiabatic representation. Due to the difficulty involved in computing a real time propagator with realistic potential energy surfaces, we describe a transformation to complex time TCFs in section 3.2. Section 3.3 discusses general features of these complex time TCFs, and writes them in a form which is used in the later part of this work. 24 3.1. Electronic Adiabatic Representation for Molecular TCFs We consider the evaluation of a general positive time TCF for two different operators Z and b. In case of the first order processes, they are the molecular transition dipole operators, and in case of the second order processes such as elastic and Raman scattering, they become the molecular polarizability tensor operators. A general positive time molecular TCF for operators A and b is written as f* (t) = K (t) t (0))) (3.1) where ((..)M is the quantal average over the probability distribution of the initial state of the molecular system, and the time-dependence in the operator A is given by A et) J*f'mA(O) e+LfHmt (3.2) Here Hm the total Hamiltonian for the extended molecular system involves the nuclear as well as the electronic motions in the system, and is written as 3N-6 h2 a2 h g2 a2 HSl = - T2Mn 0R2 + V(r, R) (3.3) n=1 A i=1 2 with r = (rj, r2, r3,...) is a set of s generalized coordinates for the motion of electrons, and R = (R1, R2, R3, ...) is a set of 3N -6 generalized coordinates for the vibrational motions of the nuclei. We use the Born-Oppenheimer approximation to separate the electronic and nuclear motions in the target. The description is suitable for the ground and first few electronically excited states of the target, where such a separation in the electronic and nuclear motion generally yields reasonably good results. In the Born-Oppenheimer approximation [Balint-Kurti, 75] the electronic 25 wave functions are the solutions of the following equations Sh2 a2 + V(r, R)] 4j(r, R) = Ej(R) Oj(r, R) (3.4) This equation is parametric in the nuclear variables R, and {Ej(R)} is a set of electronic potential energy surfaces for the motion of the nuclei. The total initial wave function of the target is then written as I'Ij(R)) = Vj(R) koi(R)) (3.5) where 10i(R)), the initial electronic state, depends parametrically on the nuclear variables R, I,) is the nuclear vibrational wave function on the ith electronic potential energy surface, and the explicit dependence of the electronic wave functions 0i (R)) over the electronic coordinates r is suppressed, because they are integrated out in the next step. Then, if WI is the probability of finding the nuclear vibrational wave function in the vibrational state IV,,) of the I'h electronic potential energy surface, the TCF fAp(t) as defined in Eq.(3.1) can be written as fAA (t) = W, (0i (R) I ( I e+ Ht (A ) (3.6) e I'h t I)| I)O4 (R)) Introducing a complete basis set in the coordinate representation E IR) (R I j)I j(R))(0j(R)I(Vj IR) (R I = 1, (3.7) and using the adiabatic approximation for electronic potential energy surfaces [Olson and Micha, 82] expressed as HfMl d(R) ~ q4;(R))Hi(,PR) (3.8) we write Eq.(3.6) as fAB(t) = 5[ WV1J dR V'f(R)e+1H(AtR)t (A(R))ij I I39) e- hi'f) (b(R))ji Vj(R) 26 where -i(, PR) is the Hamiltonian for the nuclear motion on the ith potential energy surface, and PR is the set of the momentum operators conjugate to the operators in the set R. In the time-dependent formulation, we interpret the TCF fA(t) in Eq.(3.9) as follows: at an initial time, t'= 0, the amplitude of the nuclear motion on the th electronic surface is given by (b),, 1 01) . At a later time, t' > 0, this amplitude propagates on the excited electronic surface with the nuclear Hamiltonian f-j = j(fOR). Finally, at the required time t' = t, the time evolved amplitude on the jih electronic surface is overlapped with the complex conjugate of (A)je(i/hjt) I VI). In case the system has initially been prepared in a pure vibrational state I VI) of the ith potential energy surface, the initial probability distribution function WI becomes unity, and the TCF in Eq.(3.9) describes a general light-molecule interaction process involving two different electronic states i and j. If the system, however, is not initially in a pure state I Vb), then the finite temperature formalism, in which the initial state of the system has been prepared in a thermal equilibrium, can be introduced simply by selecting a suitable probability distribution function AII. Such a probability distribution function is usually described by the Boltzmann distribution expressed as Wi = , where 0 is the inverse temperature, kBT, at which the initial state of the system has been prepared, and Ei, the stationary eigen state on the ith potential energy surface, is given by Hiloi) = (3.10) For a two surface problem using Eq.(3.10) in Eq.(3.9), the TCF fAi(t) is expressed in the following equivalent forms: 27 fA (t) = ' Tr[e h (At+i#h) e-(A) (b)i3] (3.11a) 1 1 dR dR' (Rj e+ J(t+i~h) R') A(R') (R'j e- I /R) B(R) (3.11 b) W~i e /iL(A) e- (B); I'i) (3.11c) where '-j = '-j(R, PR) is the nuclear Hamiltonian operator on the tIh potential energy surface, and the Tr[...] in Eq.(3.1Ia) is to be evaluated over the initial nuclear vibrational state OIi) of the system. Equation (3.1 1b) rewrites this trace and all the operators in the coordinate representation, which is useful especially if the initial nuclear vibrational state of the system is not known, and the propagators on the it and the jth electronic surfaces are not difficult to obtain. The form in Eq.(3.1 Ic) is used when the initial nuclear vibrational state 101y) and the corresponding stationary state eigenenergy Ei are not difficult to obtain by time-independent techniques, while the important effects, such as dissociation and rearrangement, are occuring due to the time-dependent dynamics on the excited surface. Equations (3.11a), (3.11b), and (3.11c) present a general two-operator TCF in the adiabatic representation, where (A),j and (B),, are the transition moment operators between two different electronic states i and j. If the frequency range of the light source is in the infrared region, then it excites nuclear vibrations and rotations only on the ground electronic surface of the system. In our work, this corresponds to the case j = i, with i refering to the ground electronic state of the system. Transition moments (A);j and (B),, in this case become the net dipole, or the net polarization of the system in its ground electronic state, and the TCF as 28 expressed in Eq.(3.1 la) acquires a symmetric form which can also be evaluated in our time-dependent formulation. 3.2. Complex Time TCFs and their Fourier Transforms The usual route for computing the transition rates for various different lightmolecule interaction processes that emerges from the discussion till now, is to first compute the time dependence in the TCF as defined in Eq.(3.11), and then take its Fourier transform to get the required transition rate. The route as suggested, however, is not very easy to follow for any physically interesting problem involving realistic potential energy surfaces. This is primarily due to the presence of a real time excited potential propagator in the expression for the TCF in Eqs.(3.1la, b, c). The form of this propagator is more apparent in Eq.(3.1 1b), where it is expressed in the coordinate representation. Computations of the matrix elements of a pure real time propagator, involving realistic potentials, is a very difficult problem in itself, so that to our knowledge, till now, there is no known numerical or analytical method that even tries to address this problem for a completely general potential. This is also the main reason that path-integral methods [Feynmann and Hibbs, 65; Schulman, 81] for evaluating a propagator in real time dynamics, though very simple to understand and physically appealing to the intuition, have remained basically a formal device, and have not been used in problems involving realistic potentials. Recently, there has been a resurgence of interests in these methods [Miller, Schwartz and Tromp, 83; Behrman, Jongeward and Wolynes, 83; Thirumalai and Berne, 83; Doll, 84], and significent advances are made towards developing numerical techniques for the evaluations of rapidly oscillating integrals in many dimensions [Chang and Miller, 87; Doll, Coalson and Freeman, 87]. 29 Alternatively, there have been developments such as time-independent quantum mechanical coupled channel methods [Shapiro and Bersohn, 80; Heather and Light, 83], the time-dependent quantum mechanical wave-packet methods [Heller, 81; Mowrey and Kouri, 85; Sawada and Metiu, 86], and a host of other semiclassical methods [Marcus, 72; Miller, 74; Mukamel and Jortner; 76, Micha 83; Billing and Jolicard, 84; Gray and Child, 84; Henriksen, 85] which are generally based on running classical trajectories for the nuclear motions, while the electronic transitions are treated quantum mechanically. Each of these methods have their own advantages and disadvantages, and have been applied to the electronic absorption and photodissociation problems with varying degrees of success. Finite element grid methods which are based on the path-integral approach to computing a propagator in the coordinate or the phase spaces [Thirumalai, Bruskin and Berne, 83; Coalson, 85; Jackson and Metiu, 85], and are not restricted to any specific form for the potential, have not been applied to any of these light-molecule interaction processes involving realistic potential energy surfaces in more than one dimension. This can be explained very briefly by considering the coordinate space matrix elements of a real time propagator after a time interval 2e as +00 K (X, X', 2E) = dX" (X le- L|X") (3.12) -00 ( X"|e~ h IX) where ft is a general one dimensional time-independent Hamiltonian along the variable X.For E -- 0, a typical expression for a short time approximation for the matrix elements appearing in the integrand of Eq.(3.12) [Schulman, 81], is written as 30 (X~e- 1|2X') = (2_ exp{- (X X')2 (3.13) E _+- 0i (V (X) + V (X'))} The first term in the exponent, which corresponds to the short time approximation for the kinetic energy term in the propagator, increases as the square of the difference (X - X'). The oscillating nature of the exponentials therefore increases rapidly as one considers larger and larger values of (X -X'), whereas the magnitude of these rapidly increasing oscillations remain constant over the entire coordinate space. All of the numerical techniques for computing the integrals in Eq.(3.12), which are based on approximating the limits of the infinite integrals in Eq.(3.12) to be within some large finite region of the coordinate space, therefore do not converge to any physically meaningful result. One approach to avoid this problem is to use a transformation which shifts the problem of computing the real time TCF to another problem, in which instead one computes the same TCF for complex times, and in the end, after computing the Fourier transform of this complex time TCF, one uses some suitably defined inverse transformation to extract the required physical quantity of interest. The real to the complex time transformation of the computations in Eq. (3.12), the real short time E -+ oo is replaced by an equivalent complex short time K -- oo, where K has a fixed non-zero negative imaginary part. As a result, the rapidly oscillating exponentials in Eq.(3.12), which have a uniform maginutde in the entire coordinate space, are replaced by the rapidly oscillating exponentials the magnitude of which also decay rapidly, as one considers larger and larger values (X - X'). Therefore, the numerical evaluation of the infinite integrals, within sufficiently large finite limits in the coordinate space in Eq.(3.12), does not cause any problem. 31 Next, we describe the transformation from real to the complex time for a molecular TCF fAi(t) as written in Eq.(3.11), and derive the required inverse transformation that gives us the required Fourier transform of the original TCF. If initially the molecular system is in a pure nuclear vibrational state I,,) of the 7th potential energy surface, fAb(t) of Eq.(3.llc) can then be written as fAe (t) = E eJ-)i( (0I(A)ij )(/,j I(t)jiV01) (3.14) I,J where 1j) is the final nuclear vibrational state on the Jth potential energy surface, and wjj = (ej - E1)/h. The Fourier transform of Eq.(3.14), defined as +00 fgh (w) = 1 J dt et''fAi (t) (3.15) gives the familiar Golden rule result fg (w) = Z(iI(A)J)(jl(B)iIlki) 6(W - wj) (3.16) I,J where the Dirac 6 function conserves the energy, by requiring that the energy hwji of the electronic transition between the initial and the final nuclear vibrational states IV,,) and 10j) is equal to the energy hw given by the light source. The Fourier transform fAb(w), as defined in equations (3.15), requires that the computations are done along the real time axis in the complex time plane. The transform to an axis in the complex time plane, which is parallel to the real time axis, but with a non-zero constant imaginary part, is achieved simply by replacing in Eq.(3.14) the real time t by a complex time r = t - ihu, i.e we use fgb(t) = fAb(r) with 32 fAj (t) = iwe-i''' he iJ (v |(A)i)j ,jv ) (3.17) 1,J where u is a parameter which always remains real. Using the definition in Eq.(3.15) for the Fourier transform, we can see that the Fourier transform of the complex time TCF f.A (t) is given by LS (W) e-wjuh (4j (A),jI'J) (3.18) I, J (0J A()ji 10) 6 (w - Wij) From equations (3.16), (3.18), and using the property of the Dirac 6 function requiring w7j = o, it follows that fAj(w), which defines the transition rates for different light-molecule interaction processes in the previous chapter, is given by fA (L) = wuhfA (w) (3.19) where fAA(w) is the Fourier transform of the complex time TCF fA'(t) as defined in Eq.(3.17). For a constant u this transform is taken along an axis, which is parallel to the real time axis but has a non-zero constant imaginary part. In case the initial state of the molecular system is not a pure vibrational state, but is at thermal equilibrium at the inverse temperature kB!, we use Eq.(3.llb) of the previous section and steps similar to Eqs.(3.14-3.19), to show that if the complex time TCF for a finite temperature two-surface photoexcitation problem is given by 33 f 3 (t) = +Tr [C*(ih ) (A C-*th)i (3.20) with f = - u, and u as the artificially introduced parameter, then fAf(w)= e"A'j3(w), where fA.6(w) is the Fourier transform of Eq.(3.20) as defined in Eq.(3.15). Some of the other works along these lines treat the artificially introduced parameter u as temperature; the computations of a complex time propagators are then suitably referred to as real time finite temperature computations. In case of an actual finite temperature computation, due to a thermal distribution for the initial state of the system, the requirement of a positive 3 in the argument of the lower surface propagator in Eq.(3.20) imposes therefore an upper limit on the choice of the parameter u. In our work, however, we do not use time-dependent computations for the propagator on the lower potential surface; therefore, we have not needed any constraint on the size of the parameter u imposed by the size of the natural 3 for the problem. We continue to call u simply an artificially introduced parameter, the size of which depends only upon the following competing factors. Increasing the value of u increases the damping of the upper surface oscillations and therefore increases the numerical stability of the dynamical computations. Increasing the value of u, on the other hand, also decreases the size of the numbers we are dealing with in the computations; therefore, beyond a certain value, they start to affect the overall accuracy of the computations. For each individual problem, one has to choose a compromise between these two competing factors and select a suitable range, such that u is sufficiently large for an effective damping for the upper surface oscillations, and at the same time sufficiently small so that the numbers in the computations are not too small to affect the overall accuracy. In the end, the accuracy of a particular 34 computation is measured by checking the stability of a computation with respect to the variations in the value of u. The transformation of a computation along the real time axis to another computation along a parallel axis, but with a non-zero imaginary component u, is not unique. There are other transformations, such as rotating the computational axis in the complex time plane, that can achieve similar results. In fact, some of them may later prove to be more efficient than the transformation used in this work. We discuss this further in the later chapters. Lastly, the property of the 6 function for the conservation of energy in Eqs.(3.16) and (3.18), which has been used in Eqs.(3.17) and (3.20) to get the inverse transformation between f4A(w) and Ij g(w), holds true only for a sharp spectra. Other cases, where one has to incorporate also a phenomenological line shape in the problem, such as molecular line-broadening effects, we can use modifications of the similar transformation [Coalson, 85] to achieve similar results. 3.3. Properties of the Complex Time TCFs and their Fourier transforms. Some properties of the complex time TCF fA.,(t), which we use later in this work, follow from its explicit expression which we write once again as f~b = M 'I() + ( (fl h (3.21) The initial and the excited surface Hamiltonians fi- and 7ij, respectively, are the operators over the nuclear variables {, PRI of an extended molecular system. Sometimes these Hamiltonians may contain constant terms like hAI and hA, respectively, such that the difference hAij= h(Aj - A,) may correspond simply to 35 a vertical displacement of the upper surface with respect to the lower one. These constant terms affect only the position of the line shape function in the frequency space. Therefore, they can be taken care of, very easily, in the definition of the Fourier transform fA.A(w). If the upper and the lower surface Hamiltonians are written as )Rj= k, - hAj and fi= H, - hA,, respectively, 'A in Eq.(3.22) can then be written as fAb (t) = e-"F AA (t) (3.22) with its Fourier transform given by fAb (w) = eA (-A ) . (3.23) The shifted Fourier transform F 6(w - Aj,) is given by +00 FAP (w') dt eI''ifA (t) (3.24) - 00 with o' = w - Aji, and the shifted TCF FAA(t) is written in a compact form as FAA(t) = Z W '(A) e- I) r (b) (3.25) 1 Here we have used H, I j) = EI 1), and the frequency operator nj(E) which describes the time-evolution is given by ( t- E1) A, (E1) = (3.26) 36 where fIj is an operator over only the nuclear variables {R, PR} for the j A electronic surface. For two different operators A and b, we define a symmetric form for the TCF as A 2 = (#f (t) + FbA (t)) (3.26) which is invariant under the interchange of operators A and . Replacing t by -t in Eq.(3.25) FAs(-t)= W1 K ( e+ib,(E1)r* (B) i1) (3.28) we see that FA"(-t) =FA(t)*. Now considering the Fourier transform F-A,(w) of Eq.(3.25) as 0 +00 {()dt'eiwi FA&(t) + dt'e'LI PA (t) (3.29) replacing in the first term t' by -t', and using the property FA (-t') = we get +00 FPjj (w) dt' {e-i''F,1 (t)* + e~IFj (t)} (3.30) 0 For B = A this gives +00 F,jj()= ] dt'Re teTiw"FAA (t)j (3.31) 0 37 and for f $ A, we use the symmetrized form, which is invariant to the interchange of the operators A and 5, to write Fjg (w) = Jdt'Re {e+iWI'F (t)} (3.31) 0 Equations (3.31) and (3.32) for the Fourier transforms FAA(w) and k. ( respectively, ensure that the physical properties represented by these Fourier transforms remain real. Furthermore, for B # A, we can show that the symmetrized TCF FA ,(t)can also be written as Fjb (t) = {F2 () - FA (t) - F (t) (3.33) where Fag(t), FA(t) and Fbj(t) are the autocorrelation functions corresponding to the operators C, A and b, respectively, and C = A + b. Therefore, for all physically interesting light-molecule interaction processes, from now on, we consider the evaluation and uses of the complex time autocorrelation function FAA(t), which we again write in a compact form as FAA(t) = W1 (Tj (0)1T;(r)) (3.34) I where the time-dependent amplitude I Tj (t)), which is evolving according to the frequency operator Q,(Ej) is given by lTj;(t)) = exp(-njt)|Tjg(O)) with IT j,(O) = (A),,IiI'). The overlap (T ;(O) I Tj;(r)), completely describe the timedependence in the TCF FAA(t) for a complex time r = t - iuh. To compute a vibrational infrared or a rotational microwave spectra, we merely need to replace the frequency operator n, by a different frequency operator n2,= (HI, - Ej)/h, 38 and redefine the transition dipole moment operator (A)3, as (A)jj; the total dipole operator of the system in its i"h electronic state. The physical interpretation of Eq.(3.34) is as follows: at an initial real time t = 0 and the complex time 7 = -inh, the amplitude I T 7-iu'h)) is almost identical in shape, but of smaller magnitude in comparison to I T~j(O)). For later times t > 0, the propagation of the amplitude occurs according to the frequency operator n2,, and a fall off in the value of the overlap (Tjg(O) I Tjg(r)) is observed, as the amplitude I Ti(r)) moves away from its initial position and begins to change its shape. The short time behavior of the overlap primarily decides the position and the shape of the overall line shape in the high frequency region, whereas the long time behavior of the overlap decides its shape in the low frequency region. If the overlap decays to zero at intermediate times and shows some recurrences in the long time behavior, one finds fine structures overlaying the otherwise smooth spectra. CHAPTER 4 FACTORIZATION OF TCFs INTO PRIMARY AND SECONDARY REGION TCFs We consider the evaluation of a general complex time dipole-dipole TCF Fb'b(t) for light-molecule interaction processes, involving electronic transitions between two potential energy surfaces. In particular, we are interested in application to the photoexcitation and dissociation dynamics of extended molecular systems. The usual problem one encounters in dealing with extended molecular systems is the large number of degrees of freedom involved. If we were to apply any quantum molecular dynamics method for computing a molecular dipole TCF of an extended system with realistic potential energy surfaces, the first and the foremost consideration we would have to think of is the computational efficacy of the method for a multivariable problem. It has been observed that the computational efforts for most of the methods, where there is more than one degree of freedom in the problem, increase as the square of the number of degrees of freedom involved. Therefore, all the methods that try to compute the dynamics, simultaneously along all the coupled degrees of freedom, become more and more impractical as the number of degrees of freedom start to increase from more than one or two. To avoid this problem, we treat the dynamics of an extended molecular system in the spirit of the well-known mean field approximation, which factorizes the full dynamics of an extended molecular system into many single variable problems, with time-dependent self-consistent couplings between the motions along these uncoupled variables. As a result of such factorization the computational efforts 39 40 in a multivariable dynamical problem do not scale as the square of the number of the variables involved. In section 4.1, we describe the sufficient conditions for, a) defining fast and localized photoexcitation processes in a large molecular system, interacting with visible or UV light, and b) factorizing the molecular dipole TCF for these cases into self-consistent primary and secondary region TCFs. Section 4.2 describes a timedependent variational principle (TDVP) for the amplitudes I TI(t)), introduced in the previous chapter. Section 4.3 derives a norm conserving time-dependent self-consistent field (TDSCF) approximation for factoring these amplitudes into self-consistent primary and secondary region factors, and finally in section 4.4, we use the results of sections 4.1, 4.2, and 4.3 to factor a dipole-dipole TCF Ftt(t) into self-consistent primary and secondary region TCFs. 4.1. Sufficient Conditions for the Factorization of the Amplitudes In general, light-molecule interaction processes in extended molecular systems are classified into two broad groups. First, the intramolecular process, in which a single molecule involving several vibrational and rotational degrees of freedom absorbs light in a given frequency range, depending upon which it dynamically evolves on the ground, or excited electronic surface. Second, the intermolecular process, in which an extended molecular system absorbs light in a given frequency range, followed by a dynamical evolution of the system in a given electronic state. For both cases, if the absorption of light in the extended molecular system is a fast process, we can define a sufficient condition to describe a localized region in the large molecular system, such that this region plays a more dominant role in the light-molecule interaction as compared to the rest of the system. If the frequency 41 of the light source is in the visible or UV region, the sufficient condition that defines such a localized region, is that the transition dipole operator (b(R))ij of the system, is an operator over the nuclear variables of the localized region only. In other words, the transition dipole operator, for a localized electronic excitation process in a large molecular system, is written as (b(R))ij = (f(X))ij (4.1) where R = {X, Q} is a set of generalized nuclear variables for the entire molecular system, with X = {X1, X2, X3, ...} a subset that defines the variables of the the localized region within which the important dynamics takes place, and Q = {Qi, Q2, Q3, ...} the subset that contains all the remaining variables of the system. We use "primary" region for the region covered by the elements of the subset X, and "secondary" region for the region covered by the elements of the subset Q. With the definitions of the primary and the secondary regions of an extended molecular system at an initial time t = 0, in place, the sufficient condition for the factorization of the amplitude I TI(0)) into factors corresponding to the primary and secondary region dynamics for a later time t > 0, is described by writing the initial nuclear vibrational state 0 t1) on the i'h potential surface as 01) = I 101) (4.2) where p corresponds to the primary region dynamics and s to the secondary region dynamics. The stationary state eigenvalue E1, corresponding to the sufficient condition in Eq.(4.2), is then written as E, = E7+Ej + EP', where the last term is due to the coupling between the primary and the secondary regions on the ith potential energy surface. Furthermore, it is not difficult to see that the sufficient condition in Eq.(4.2), for the factorization of the amplitude into self-consistent 42 primary and secondary region factors, is the validity of the time-independent Hartree (TH) approximation for the stationary state vibrational wavefunction on the ith potential surface [Bowman, 86]. Given the sufficient condition in Eq.(4.1) for the definition of a localized region in an extended molecular system is satisfied, it is not difficult to find a time-independent Hartree approximation to the stationary state eigenfunction 01) on the ith potential energy surface, such that IT, (0)) = 1k' (0)) Irn' (0)) (4.3) where I ('(0)) = (b(X))ij I Of) is the initial value of the amplitude for the primary region dynamics, and I n'(O)) =1 ;bl) is the initial value of the amplitude for the secondary region dynamics. As we shall see later, equations (4.1) and (4.2) completely define the sufficient conditions for the factorization of the timedependent amplitude into self-consistent motions corresponding to the primary and secondary region dynamics. Therefore, these are also the sufficient conditions for factoring the molecular dipole TCF for an extended system into primary and secondary region TCF factors. If the frequency of the light source is in the infrared region, the sufficient condition for describing a localized process in an extended system in Eq.(4.1) can be generalized by using a series expansion for the total dipole moment of the system in a single electronic state to write it as (D(R));i = Z(b'(X))ji PI(Q) (4.4) where (b'(X)u is the net dipole moment of the primary region when the molecular system in the zih electronic state, and P1(Q) is a polynomial in the variables of the secondary region. The description in the later chapters of this work, with slight 43 modification, can also incorporate this generalization of the sufficient conditions in Eqs.(4.1) and Eq.(4.2), to include the interaction of infrared light with extended molecular systems. 4.2. Norm Conserving Time-dependent Variational Principle (TDVP) The time-dependent amplitude IT(t)) which defines the time-dependence of the total TCF in Eq.(3.34) is a solution of the following equation of motion (Zi - A2) IT (t)) = 0 (4.5) where the initial condition |To) = IT(0)), defines the normalization N at an initial time t = 0, as N = (To I To). For a two surface problem, we have dropped in Eq.(4.5) specific references to the subscript j in the frequency operator n, and to the subscript I in the amplitude I TI). When the exact solution IT) = IT(t)) of equation (4.5) is complicated, it is sometime desirable to replace it by an approximate solution I) = jI(t)) with the initial conditions |To) = IT(0)) at t = 0. We construct an approximate solution to Eq.(4.5) by using a time-dependent variational principle (TDVP), the general form of which is similar to those given by Dirac [Dirac, 26, 30] and Frenkel [Frenkel, 34] in the earlier days, and by Mclachlan [McLachlan, 64] in later modifications, see for example [Langhhoff, Epstein and Karplus, 72] for a review. In our work, however, we also impose the conservation of the total norm N by using the well known Lagrange undetermined multiplier technique. The variational functional is constructed in steps, by first introducing the Lagrangian 44 L [T, T*] = ( KT ) + (T IT)) (4.6) where i = i - 9, and the form of the Lagrangian in Eq.(4.6) is invariant under time reversal. From this, the action between the initial and the final times follows as S [T, T*] = dt' L [T, T*] (4.7) ti and finally the norm conserving functional for Eq.(4.5) is written as S' [TA, T-] = S[T, T-]+ I dt' A (t') (TAIT) (4.8) ti where in our notation the subscript A in the approximate solution ITA) refers to the undetermined Lagrange multiplier A, which is solved by using the constraint that the norm N defined as N = (TA(t) I TA(t)) is a conserved quantity. We solve equation 6S'[T*, TA] = 0 with fixed boundary conditions, and treat the variations 6T* and 6T\ independently of each other. We choose the variation 6TA = 0 and allow for a T*, to get (6TAII)TA + AT,\) = 0 (4.9) which is equivalent to VTA + ATA = , (4.10) The subscript A in the approximate solution TA refers to the additional constraint, which needs to be satisfied solving for A. Using a transformation to remove A from Eq.(4.10), we write the trial solution T\ as 45 TA = ) exp[i J dt' A (t')] (4.11) iI where 41 satisfies V4' = 0 (4.12) and the normalization condition N = (TAI TA) gives tN exp[-2] dt' {Im( A (t'))} ] = (4.13) tl Since Re(A) does not appear in the normalization constraint in Eq.(4.13), we are free to choose it equal to zero. Then the norm conserving solution to Eq.(4.5) is given by TA = N' ( ) A (4.14) where 4A(t) satisfy Eq.(4.12), and the normalization constant N is constructed from the initial conditions. 4.3. Time-dependent Self-consistent Field (TDSCF) Approximation Expressing the norm conserving time-dependent variational principle (TDVP), of the previous section, in a compact form as (STAIDTA + ATA) = 0 (4.9) The frequency operator Q, which governs the time evolution of the exact amplitude IT(t)), is of the form 46 Q = QP + Q-1 + QP8 (4.16) where the superscript p refers to the variables of the primary region, s to the variables of the secondary region, and ps to the coupling between the primary and the secondary regions. We write the time-dependent self-consistent field (TDSCF) approximation to Eq.(4.15) as I TXA() [t)W() (4.17) where the factors (t) and j(t) correspond to the self-consistent dynamics of the primary and the secondary regions, respectively. Treating the variations 6 (t)* and 6b(t)* independently, and choosing either 6b(t)* = 0, or 6b(t)*= 0, we get the following coupled equations i~l)+ iRQS) -( 7A ~) + A (I-~) = 0 (4.18a) i{(|$;+ i()i )- (( ()i+ ((|) = 0 (4.18b) Using the transformations (=('e I dt' (4.19a) W - ( I ) 1 q =' exp i ~~ di' (4.19b) in equations (4.18a) and (4.18b), respectively, we write them for (' and 17' in a form which is more easy to implement. Further using other transformations, 47 S= exp dt'A(t') (4.20a) 'i dt'A (t') (4.20b) we remove the normalization constraint A from the differential equations corresponding to '(t) and q'(t), respectively. We find the differential equations corresponding to the self-consistent primary and secondary region amplitudes (t) and 7(t), respectively, in their final form as - ( - ( P1) = 0 (4.21a) i(I - P - 77 = 0 (4.21b) Imposing the normalization constraint N = (TA(t) I TA(t)), and using the transformation equations (4.19a,b) and (4.20a,b), we write N = (77I1 r) exp[- { f RO + 07-7) (4.22) tIt + ( ~ }Odt' - 4 dt'Im { A (t') Subtracting complex conjugate of Eq.(4.18a) from Eq.(4.18a), and using the result in Eq.(4.22), we express the normalization constraint as 48 (423 exp -Idt'In { A (t')} ] (4.23) where again as in the previous case we are free to choose Re(A) = 0, since the normalization constraint does not involve Re(A). Finally, using the transformations in equations ( 4.19a,b) and (4.20a,b) in the TDSCF approximation in Eq.(4.17), we get N = r exp[+i { (4.24) ti + + 21m (A (t'))} dt'] where again we use Eq.(4.18a) to write + (4.24) i(MI) (QII _ P -11 A which together with the normalization constraint in Eq.(4.23), and the choice Re(A) = 0, gives the norm conserving TDSCF approximation in its final form as T A = N 2 exp i dt' K (7 ) (4.26) ( 1() (77171)2 ( 77 7) The amplitudes (t) and 7(t), corresponding to the self-consistent primary and secondary region dynamics, respectively, satisfy the following coupled equations: 49 -- (t) = 0 (4.27a) ii) - n'.1 (t) rj = 0 (4.27b) where the time dependence in the effective frequency operators f (t) and "f f (t) are given by jP = np + (4.28a) Eff (r7 1r/) nff = n + (f ) (4.28b) Similar to the time-dependent Hartree (TDH) factorization of the total wave function [Schatz, 77; Gerber, Buch and Ratner, 82; Kugar, Meyer and Cederbaum, 87; Makri and Miller, 87] for a coupled many-mode problem, we also achieve a formal separation of the total amplitude, T(t), into self-consistent primary and secondary region factors (t) and 7(t), respectively. The time-dependence in each of these factors is governed by an effective frequency operator, which includes an average of the interaction frequency operator, over the other factor. The essential feature of the TDSCF approximation is that, while it provides the formal factorization of the dynamics of a multivariable system into primary and secondary region dynamics, it still permits the self-consistent energy exchange between these two regions through the effective couplings in the frequency operators. An essential requirement in our work is, that we must have a norm conserving form for the amplitude T(r) for a complex time r > 0, and a correct phase factor in Eq.(4.26). This is because we are interested in computing the overlap of the complex time amplitude with its initial 50 value at r = 0; therefore the phase factor an integral of the effective coupling from 0 to 7, does not get canceled by its own value at 7 = 0. The norm conserving form for the amplitude in Eq.(4.26) is essential because the effective frequency operators ,f () and niff(T) are non-hermitian for any complex time r > 0. 4.4. Primary and Secondary Region Factorization of TCFs Using the TDSCF factorization of the amplitude T(r) as shown in Eqs.(4.2628), in Eq.(3.34) of the previous chapter, we write the complex time molecular TCF as Ftb (t) = T W1N1f1 (r) g, (r) e i"&() (4.29a) where we have replaced the operator A by a molecular dipole operator 1, restored the subscript I corresponding to the sum over the initial nuclear vibrational state IV,,), and the complex time 7 is given by r = t - ihu. We define the "primary region TCF" fI(7r), as a dimensionless overlap fj (r) = (0) 1' (7)) (4.29b) and the "secondary region TCF" gl(r), also as a dimensionless overlap gj (r) =(T' (0)I7 (-)) (4.29c) The dimensionless primary region TCF fj(r) and the secondary region TCF gj(r), together with an overall phase factor al(7), which is an integral of the selfconsistent coupling from 0 to r given by 51 aj (r) = dt' (4.29) 0 completely describe the time-dependence in the total TCF Fbb(t). The normalization constant N, in Eq.(4.29a), which is computed from the given initial conditions 6i(0) and rI(O) as according to N1 = (6, (0) 77, (0) (' (0) 7,' (0)) (4.29d) provides correct dimensions to the overall TCF. CHAPTER 5 TIME-DEPENDENT EFFECTIVE FREQUENCY OPERATORS In the present chapter, we consider a general nuclear vibrational Hamiltonian in a given electronic state, and solve the equations of motion for the secondary region dynamics. In section 5.1, we describe the time-dependent effective frequency operators by considering many degrees of freedoms in the secondary region interacting strongly with a single degree of freedom in the primary region. We model the strongly coupled primary and secondary regions couplings which are of arbritrary analytic form in the variables of the primary region, and have a linear and bilinear dependence in the variables of the secondary region. Section 5.2 describes the effective frequency operators for a weak coupling limit, and ignores the terms in the coupling which are bilinear in the variables of the secondary region. The solutions to the equations of motion for the secondary region operators, in both of these cases, are reduced to analytic forms, which are then used to construct the propagators for the secondary region dynamics. The primary region dynamics will be described in the next chapter by using a numerical path-integral method. In section 5.3 we derive an operator equation of motion for the primary region dynamics, show its relationship to the well known phenomenological Generalized Langevin equation (GLE) [Mori, 65], by constructing the fluctuation force term and the dissipative memory kernel, and prove the Fluctuation-Dissipation Relationship (FDR) [Van Hove, 54] explicitly for a general intermolecular process. Section 5.4 considers many degrees of freedom in the primary region dynamics interacting with a set of coupled harmonic oscillators in the secondary region dynamics. 52 53 5.1. Intramolecular Processes: Strongly Coupled Primary and Secondary Regions In a general intramolecular photointeraction process, following the absorption of light in the visible and UV region, we consider electronic transitions between fixed electronic states i and j. Therefore, for notational brevity throughout this chapter, we avoid writing i and j as the specific references to the initial and the final electronic states, respectively. Consider a general nuclear vibrational Hamiltonian H in the given electronically excited state as consisting of the following terms: 0 = +V" (k (5.1a) 2mx is the Hamiltonian for the primary region dynamics with an arbitrary potential V'( ) given either in a parametrized analytic form based on experimental information, or on a grid in coordinate space from ab initio or semi-empirical electronic structure calculations; p2 =t + -mkwkQk (5.1b) k 2mk is the Hamiltonian for the secondary region expressed as a sum of k(k = 1, 3...N) harmonic oscillators, such that the mass mk and the frequency Wk of the kh harmonic oscillator correspond to the mass and the frequency of the kih normal mode vibration in the molecule. The coupling between the primary and the secondary regions is expressed as Hp = { PkAkXPX + Vk ( ) Qk + E QkVk'k (k) Qk (5.1c) k k' I 54 where the arbitrary functions Vk(Z) and Vkk'(Z) are operator dependent forces and potential curvatures, respectively, and the bilinear momentum coupling parameter AkX has dimensions of inverse mass. The arbitrary functions Vk(Z) and Vkk'(X) are also given either in an experimentally parametrized analytic form, or on a grid in coordinate space as computed in an electronic structure calculation. Following the absorption of light in the visible or UV region, one finds that the primary region dynamics, describing for example rearrangement in the system due to the formation or breaking of a bond, is strongly coupled to the rest of the molecule. Therefore, we have added to the coupling the terms linear and bilinear in the operators of the secondary region. We also note that the formulation does not depend upon the strength of the various terms in the coupling, and in the next section we consider the weak coupling limit by ignoring the first and the last terms in Eq.(5.1c). For a single degree of freedom in the primary region dynamics coupled to many degrees of freedom in the secondary region, we write the norm conserving TDSCF approximation for the amplitude TI(t) as T' (t) = N2 e+i(i) ( (t) k (t)(5.2) I 7k= kI 2kI)} Here the secondary region amplitude I(t) in Eq.(4.26) of the previous chapter is written as a product over k, and the amplitude r/kI(t) corresponds to the kh normal mode motion in the secondary region dynamics. The time-dependent effective frequency operator f2fJ f(t) for the primary region dynamics is given by 55 nff (t) = - + Pk (t, [71]) ALx-x k + Vk(X) qk (t, [ii])+ I Vkk( ) q' (t, [1]) (5.3) 2k qkIk(t,[i])Vk'k(-C)qk (t, [7]) - Epk1] 2kjk' where we have omitted specific indices within the functional dependence for notational brevity, and on the ith potential energy surface we have written the stationary state eigenvalue EI as EI = EP + E- + EP", with E-= k E, and Ej8 = Ek Eyk. The time-dependent coefficients, which couple the primary region dynamics to the normal mode motions in the secondary region, are the normalized expectation values of the corresponding operators within the amplitudes coresponding to the secondary region dynamics. For example, the coefficient q (t, [1]) is of the form kI4 kI) q(t,[]) =(7kI17) (5.4) The time-dependent effective frequency operator ,'gff(t) for the motion in the secondary region dynamics is given by n (= ne f (t) (5.5) k where the effective frequency operator for the kih normal mode motion in the secondary region is written as 56 (t) [ + MkD'(t, [)Q PkAkXPX (t,[i]) +) k (t,[, ) Qk + ek (E [k _] E Ejk] (5.6a) where i means any 27k, with k' 5 k. The time-dependent effective frequency Wk(t, [6j), the effective force Vk(t [(, i]), and the effective energy gk(t' [ ,]) between the kIh normal mode in the harmonic bath, the primary region, and the rest of secondary region, are written as Wk (t, k + vkk }(5.6b) T)- k t k(,[] + tkk t[6)+ Vikk (t, [6])) qk, (t, [ND k'$k (5.6c) gk (t, [, q]) = {Pk' (t, [7]) Ak'XPX (t, []) k'$k + Vk' (t, [6]) qk' (t, [N]) + !v'' (t, [)) , (t, [i]) (5.6d) + qk" t,9]Vk"k' k k2 k"54k'54k 57 The time-dependent coefficients that couple the kik mode of the secondary region to the primary region dynamics are defined, for example, by k'k W I vk'k le) Vk'k (t, = (5.6c) The time-dependent effective phase a(t, [77, s]), defined in Eq.( 4.29d) and needed to compute the total TCF Ftt(t) in Eq.(4.29a), is given by a(t, [ , 7]) = dt' Z{pk (t', [7]) AkXPX (t', [a]) 0 k + Vk (t' qk (t', [711 + 1 tkk ('[ ]) qI (t, [,q]) (5.7) + qk ', [ ) k ' , [ q)ki(, ) 21 k'ok The form of the effective frequency operator ff(t) in Eq.(5.6) for the kth mode of the secondary region dynamics, is that of a linearly forced parametric oscillator. Operators Qk(t) and Pk(t), corresponding to the kih mode motion in the secondary region dynamics, evolve according to the following equations of motion Vk = i (t), Qk] (5.8a) kn= i (t), 1k (5.8b) where [...]c are quantum commutators. Using 2f(t) from Eq.(5.6a), it is easy to show that Qk(t) satisfies 58 k + ,k (t) Qk = AkXnX (t) _ k (t) (5.9) rnk where we have dropped the functional dependence on and i for notational convenience. Equation (5.9) is the equation of motion of a linearly forced parametric oscillator, therefore we solve it by separating it into a homogeneous equation and a nonlinear forcing term. We are using this approach, because it is suitable for application to computations in the complex time plane, and also because it is useful in deriving a parallel between the TDSCF dynamics and the GLE approach to stochastic dynamics. We introduce uk(t) and vk(t), two linearly independent solutions of the homogeneous equation +k()= 0 (5.10) Tt2+ CkI{t)] k t)} = with the initial conditions Uk(0) = 0, and Vk(O) = 1. The full solution to the operator Eq.(5.9) is written as Pk (0) Qk (t) = Qk (0) vk (t) + Uk Uk (t) + Yk (t) (5.11a) where Qk(0) and Pk(O) are, respectively, the initial values of the operators Qk(t) and Pk(t) at the initial time I = 0, and the the last term yk(t) is an integral of the form Yk (t) dt' A {k xtx _ i' k - )I Uk (t -t') (5.11b) MkWk 0 59 Using equations (5.11a) and (5.11b) with the definitions in Eq.(5.4) for the timedependent coffecients of the coupling between the primary and the secondary region dynamics, we express the coefficients needed to compute the effective frequency operator nff (t) for the primary region dynamics as qk (t) = qk (0) vk (t) + Pk (0)Uk (t) + Yk (t) (5.12a) MkWk Pk (t) Pk (0)ik (t) + mkqk (0) i'k (t) + mkyk (t) - AkXmkPX (1) (5.12b) q2(t) = q (0) v2 (t) + U (t) + y (t) qk k k2 2Uk )+kk mkwk + 2 qk (0) Vk (t) + Pk (0) U k (t) yk (t) (5.12c) MkWk where the initial values qk(O), Pk(O), and q2(0) are defined as according to (0)= (r/kI (0) IQ21r/k (0)) qk ( ()= kI (0) IrkI (0)) (5.12d) At any time t the time-dependent effective frequency operator nPf f(t) depends upon the previous history of the primary and the secondary region dynamics through the term yk(t), which involves integral of the nonlinear forcing term from 0 to t, therefore, in our TDSCF approach we refer to yk(t) as a "memory" term. Furthermore, in a later section of this chapter, we show a direct relationship between the TDSCF approach and the GLE approach to stochastic dynamics, by deriving a 60 fluctuation-dissipation relationship between the memory kernel and the fluctuation force. It appears from Eq.(5.llb), that the computation of the memory term Yk(t) at any time t, involves knowing the integrand at all previous times t' < t, as well as at the present time t' = t. It is, however, not difficult to see that the choice of the initial conditions, uk(O) = 0 and vk(O) = 1, for the solutions of the homogeneous Eq.(5.10), make the integrand in the memory term yk(t) equal to zero at t' = t. Therefore, this choice of the initial conditions is significant in the sense that it enables us to perform the TDSCF computation, for a nonlinearly coupled primary and the secondary regions dynamics, without going through iterations at each time. In other words, these initial conditions have made it possible to perform our computations in a step by step progression along a desired axis in the complex time plane, while at the same time maintaining the self-consistency in the dynamics. We also note that the linearly independent solutions uk(t) and VOk(t) which satisfy the homogenous Eq.(5.10) for the kth normal mode motions in the secondary region, are completely independent of the simultaneous solutions uk'(t) and vk'(t) for all k' 3 k. Therefore, for any secondary region made up of N coupled harmonic oscillators, we have to simultaneously solve only 2N uncoupled homogeneous equations of motion to achieve the self-consistent solution of the full dynamics. The initial values q,(O), and pk(O) defined as according to Eq.(5.12d), are the normalized expectation values of the corresponding operators within the amplitude 77k(O). They can also be thought of as the mean position and momentum corresponding to the kih normal mode motion in the secondary region dynamics. If the secondary region is initially prepared at thermal equilibrium in the presence of the primary region, then we can use a suitable thermal distribution for these ini- 61 tial values to introduce a temperature dependence in the dynamics of the primary region operators. We end this section by using equations (5.5-5.6) for the definition of the effective frequency operator ngff(t) for the kih mode motions in the secondary region dynamics, and equations (5.11-5.12) as the solutions for the operator Eq.(5.9), to construct the propagator for the kth mode motion in the secondary region dynamics [Yamashita and Miller, 85] as Kk ^ ^ MkWk iMk ' = ihukp(t) 2huk (t) {(Qk + Qk)V (t) - 2Qk - Tfkw jfk (t') Uk (t') dt' MkLkJ 0 2 t -M2k J fk (ti) Uk (t - t') dt' (5.13a) 0 t t' -2W mw Af (ti) Uk (t - t') dt']f fk (t") Uk (t"l) dt" I 0 0 + (Elk+ Ejk) t] where the effective force fk(t), in the integrand, is given by fk (t) = mkAkXpx (t) - k k (t) (5.13b) and the time evolution of the amplitude corresponding to the k'h normal mode motion in the secondary region dynamics is described by 77 (Qk,t) = dQ'k (QkIKk (OkQ't)IQ') 7k1 (Q',o) (5.14) -010 62 This is about as far as one can proceed in general, i.e, without numerically solving the homogeneous Eq.(5.10) for a given process. 5.2. Intermolecular Processes: Weak Coupling Limit Intermolecular photointeraction processes such as overtone line broadening, dephasing, and energy transfer in an extended molecular system, following a localized absorption of light in the visible and UV region, are often described by the time evolution of a localized primary region, interacting weakly with a collection of coupled harmonic normal mode motions in the secondary region. In the formulation of the previous section, we drop the bilinear coupling terms, and keep the term which is linear in the operator dependence of the secondary region dynamics and has an arbitrary analytic form in the operator dependence of the primary region. In other words, the effective frequency operator Of f(t) for the primary region dynamics simplifies to f(- E + V qk (t, [ri]) - (5.15) k where the effective coupling of the primary to the secondary region is also defined by the example in Eq.(5.4). The time-dependent effective frequency operator for the kih mode of the secondary region simplifies to (mkwkQi + - E0(5.6)pk ef (t 1 L+ + +M k k~ t [ ) k - EI ) (5.16) 63 and the time-dependent effective phase needed to compute the total TCF is written as a (t, [ , r]) = dt' IV (t, [ )qk (t, [r)- E, (5.17) The form of the operator in Eq.(5.16) is that of a linearly forced harmonic oscillator. It is not very difficult to see that the homogeneous Eq.(5.10) for the functions Uk(t), and t'(t) reduces to the equation of motion for a harmonic oscillator with frequency wk and mass mk. Therefore, the functions uk(t) and vk(t) with the initial conditions, uk(0) = 0 and vk(0) = 1, simply reduce to 1k(t) = sin(wkt) and Vk(t)= coS(wkt), respectively. The time-dependent effective coupling qk(t, [7]), which couples the primary region dynamics to the kth mode motions in the secondary region dynamics is given by qk (t) = qk (0) COs Wk t + Pk(0) sin wkt + Yk (t) (5.18) where the "memory" term yA(t) for the weak coupling case simplifies to t Yk(t) 1 fIdt' Vk (t') sin wk (t - t') (5.19) Mkok 0 Here again, we notice that the computations for the memory term in Eq.(5.19) involve the time dependence in the effective potential V (t) only for the previous times t' < t. Therefore, the TDSCF computations for the primary region dynamics can be performed in a step by step progression along any chosen time axis. Finally, the propagator for the secondary region dynamics, in the weak coupling case, is expressed in an analytic form as [Schulman, 81] 64 Kk ( A' ( mkwk Z imk j( 2 ) , 2Qp Q'S~ t~') 27wih sin Wkt/eP 2h, sin Wkt I k + qk)CS - 2QQ' - 2Q' dt'v k (t') sin Wkt' - 2Qk Mkwk I Mkwk 0 x dt'k (') sin y (t -t) - m 2 (5.20) 00 t t' x Idt' Vk (') sin Wk (t - t') IdtitV k (t") sin wkt"} 0 0 + (E k + E pk) ] 5.2. Relationship to the Generalized Langevin Equation Approach The method of factoring the dynamics of a large molecular system into a "primary region" consisting of a few important degrees of freedom, interacting with the rest of the system labeled as the "secondary region", is well known in statistical mechanics, where one often needs to deal with the dynamics of a "system" in thermal equilibrium with a large heat reservoir called "bath". Among many different approaches to deal with the factorization of a large system into coupled primary and secondary region dynamics, the approach based on the Generalised Langevin Equation (GLE) [Mori, 1965], solving it under different approximations for the primary region dynamics, is well covered in the literature. The terms which affect the dynamics of primary region variable, due to the coupling with a large heat reservoir, are the "fluctuation force" term and the "memory" term. The form of the fluctuation force term is decided by the dynamical and the statistical properties of the secondary region, and depends largely upon the way the secondary region is brought to a thermal equilibrium in the presence of the primary region. 65 If the primary region is coupled to an infinite heat reservoir for a sufficiently long time, the primary region also attains a steady state due to a relationship between the memory kernel and the flutuation force, this relationship is often referred to as the Fluctuation-Dissipation Theorem [Van Hove, 54; Berne, 71]. The GLE is a stochastic differential equation with additive fluctuation force and dissipative memory terms. Starting with the time-dependent effective frequency operator n ~f(t) in Eq.(5.15), we derive an operator equation for the primary region dynamics, show its relationship with the GLE by constructing the fluctuation force and the memory kernel terms, and finally derive the Fluctuation-Dissipation Relationship(FDR) for our TDSCF dynamics. The time evolution of the operators and Px for the primary region dynamics is described by the following operator equations = i [f f(t), (5.21a) X= i [f f(t), Px (5.21b) Using Eq.(5.15) for the effective frequency operator, Eq.(5.21b) is written as x = [V" (k) ,fx] + [i/ () ,Px qk(t) (5.22) k where we avoid the functional dependence in qk(t) for notational simplicity. From Eq.(5.18) the effective coupling qk(t) is written as 66 qk(= q (0) + 1 Skk (0) COS Wkt + pk( ) Sin C14t Mk~Jk I rkwk V k (t _ 1 mkk J 1 dt' (t') coswk (t - t') (5.23) MkWk MkWk 0 0 where we have used integration by parts in the memory term in Eq.(5.19), and substituted the result in Eq.(5.18). Using Eq.(5.23) in Eq.(5.22) and rearranging the result we get S_ [Mod () ] [Vk ( X c k , 2 COS Wk (t - t') 0 2MkmXwk + Rk (t) } (5.24) where the commutator involving the modified frequency operator UVod(t) given by fAod pX, + k _m V (t)Vk -E (5.25a) describes the time-dependence in the uncoupled primary region dynamics in the presence of a stationary secondary region, and the fluctuation force Rk(t) is given by Rk (t) = ]c ( 1 (,q (0) + Vk (0)) coswkt + Pk (0) sinwktI. hc MkWk / mkwk ) (5.25b) 67 The integrand in Eq.(5.23) has been simplified to the form in Eq.(5.24) by using the equation of motion of the operator V1k() to write j'k = Vk (5.25c) / 2mlLJ with (0..) =0) . (5.25d) ( (0) (0)) The form of Eq.(5.24), a general operator equation for the primary region dynamics in the presence of a weakly coupled secondary region, is an analogue of the GLE in our TDSCF approach. We have considered a coupling which is linear in the variables of the secondary region, and has a general analytic form Vk(X) in the variables of the primary region. Therefore, Eq.(5.24) for the primary region dynamics in our TDSCF approach is more general than the usual GLE given in the literature. As an example, we consider a simple case of Vk(Z) = FjX, with Ik a free force constant, and use this in Eq.(5.24) to write it as PX -- i [ 0 ,x P, -Idt'M (t - t') K= R(t) 0 (5.26) where M (t -t') = Mk (t - t') (5.27a) k and 68 R (t - t') Z Rk (t - t') (5.27b) k with AMk (t - t') = Wkc0s o (t - t') (5.28a) MkWk and Rk (t) = k (qk (0) + k 2 x (0) cos (Wkt) + Pk (0) sin (wkt) (5.28b) \ Mkwk / mkwk ) respectively. Equation (5.26) is the analogue of the usual GLE for a bilinearly coupled system-bath problem, where the fluctuation force term R(t) and the memory kernel M(t - t') are given in Eqs(5.27a) and (5.27b), respectively. It is not difficult to see that the primary region dynamics depends upon the initial values of the secondary region operators Qk(0) and Pk(O) only through the fluctuation force term R(t), which contains the normalized expectation values qk(0) and pk(O), given for example by q7 (0) 1 k (0) 1k (0)) (5.29) q (0) (77k (0) 177 (0)) Therefore, the fluctuation force term R(t) can be given a stochastic interpretation by observing that these expectation values are described by a distribution function which depends upon the initial temperature of the system. If the initial state of the secondary region is an equilibrium state, then the distribution of all the bath operators is the same and consequently R(t) is Gaussian (Central Limit Theorem). The mean value of the resulting Gaussian depends upon the form of the distribution function. If initially the secondary region has been brought to the equilibrium in the 69 presence of the primary region, a suitable choice of the distribution function, which give the mean value of the fluctuation force equal to zero, is a displaced Boltzman distribution function [Zwanzig, 73; Lindenberg and Seshadri, 82] expressed as W(q(O), p(o)), 3) = exp(- J[ + ) mkw{qk(o) + (0)}) k5 2mk 2 mkwk (5.30) The thermal averaged correlation function of the force term R(t) then is written as (R(t)R(t'))en = S 5/ dq() dp(O) W(q(O), p(O)), f) (5.31) k (531 x Rk (qk(0), Pk(0),1t) Rk'(qk'(0), N, 0) where Rk(qk(O),pk(O),t) and Rk,(qk'(0),pk,(0),t') are as defined in Eq.(5.28b). Using the distribution function W(q(O), p(O)) in Eq.(5.31), we obtain the correlation function of the fluctuation force (R(t)R(t')), = kBT mk CoSwk (t -t') (5.32) k M~ which is the same as (R (t) R (t')}e = kBT M (t - t') . (5.33) This shows a FDR between the memory kernel M(t - t') and the fluctuation force R(t). For the more general operator equation in (5.24) for the primary region dynamics, we get a similar relationship between the memory kernel and the fluctuations in the force. 70 5.4. Treatment of Many Degrees of Freedom in the Primary Region Dynamics To treat many degrees of freedom I = ltoL in the primary region dynamics, we generalize the Hamiltonian operators of section 5.1, and write L ^ 2 = m S+ V0(i) , (5.34a) fts = -+ 2mkw, (5.33b) k=1 ftP =S PkAkIPI + Vk(X)Qk - QVk'k(x)Qk } (5.33c) k I k'#k where i= {X1, 2, ... L } are the operators of position variables of the primary region, and k runs over all the harmonic modes in the secondary region dynamics. Using the TDSCF approximation TI = N ___ IkI(Qk) exp(ia(t)) (5.34) ( Ig k' 12kIkg the time-dependent frequency operator for the primary region dynamics is given by ff t) - fI - EP+ f (4 D f( = [E-- E+{(Pk (t) AkIA k I + qk(t)Vk(f) + lq(t)Vkk(f) (5.35) qk' (t)Vk'k(X)qk(t) - Ek} k'#k and the time-dependent effective frequency operator for the kih mode of the secondary region dynamics is given by eff t) [ + Mk '(t) Qk+ PkAkIPl (t) + f) (t) Qk +- k(t)-E -EP (5.36a) 71 with Wk = 2 + rk 2 (5.36b) i7k (=k (t) + + k'k (t) + Vkk' (t) qk (t) (5.36c) k'#k k Pk' (t) A PI () + V k' (t) qk' (t) k'$k I 1 :qk" (t) Vk"k' (t) qk' (t) } (5.36d) k"#k'34k where the time-dependent coefficients of the coupling between the primary and the secondary region dynamics are defined as in sections 5.1 and 5.2. The effective frequency operator for the kih mode of the secondary region dynamics is still similar to that of an operator for a linearly forced parametric oscillator. Therefore, the solution to the equations of motions for the secondary region dynamics, obtained as described in section 5.1 and substituted in Eq.(5.36), give the multi-variable time-dependent frequency operator ff(t). One can use Eq. (5.35) to deal with the dynamics of a multi-variable primary region coupled to a harmonic secondary region. If, however, the computations of a complex time propagator in two or three dimensions are prohibitive, a further simplification is achieved by using an SCF assumption for the primary region, so that 7 = r (5.37) 72 in Eq.(5.34) to write the TDSCF approximation as T I~ rr ______ kI( Qk) T ' 2 (CII I rkI)' exp('a(t)) (5.38) (gI ~ k k2I ~ The time-dependent effective frequency operator for the primary region dynamics is then written as ff (t) = fnf (t) (5.39a) with eff(t)= + o ( 1,t) + pk(t)AkPJ h 21771 k + k ( ,t) qk (t) + ,1 kk (kt) q2 (t) 2'$k + 1 Pk (t) 'kI'PL' (t) - Ek } - E ] (5.39b) where v(XI, t), Vk(XZ t) and vk'k(X, t) are defined for example by v(-1)I(1+1)I (-1)I(+1) (5.39c) The effective frequency operator for the motion along the 11h degree of freedom of the primary region can be used to compute its complex time propagator, and the computational efforts for the full primary region propagator scales linearly with the number of degrees of freedom in the primary region dynamics, instead of quadratically, as would happen if no SCF assumption was made for the primary region dynamics. CHAPTER 6 PATH-INTEGRAL APPROACH TO THE PRIMARY REGION DYNAMICS The path-integral approach to the primary region dynamics is described mainly by considering the time-evolution of the primary region amplitude I (t)) according to i |(t)) = n'p(t)| (t)) ,(6.1) where the time-dependent effective frequency operator nY(t) has been described in the previous chapter. We write the integral solution of Eq.(6.1) in coordinate space as +00 (X ()) = dXo KP (X, Xo, t) (Xo I(0)) (6.2) -00 where KP(X, X0, t), the real time propagator for the primary region dynamics, is written as +00 +00 KP(X,Xo,)= dXN-1... dX1(X e-6_,t XN-1 -00 -00 .. X1 e-aX0 (6.3) In Eq.(6.3), we have divided the total time t into N short time steps of length At =7 where N -+ oo. The time-dependence in the frequency operator ni during the ith step in time for i going from 0 to N - 1, is described in Eqs(5.3) and (5.15) of the previous chapter. In the Lagrangian formulation of the propagator the 73 74 multidimensional integral in Eq.(6.3) is also written as a functional integral, which is equivalent to a summation over all possible paths between the initial and the final positions in coordinate space. The number of all such possible paths is infinite; therefore, since their inception in the early sixties, the path-integral (PI) methods [Feynmann and Hibbs, 65], though very appealing to the physical intuition, have mostly been used as formal devices. In recent years, these methods have also been shown to be a convenient numerical tool in the studies of many-body quantum mechanical systems, and have been used in the computations of imaginary as well as complex time propagators. All of the recently developed numerical approaches to the multidimensional integral in Eq.(6.3), or equivalently to a functional integral in the Lagrangian formalism, are in general separated into two broad groups. Firstly, there are the discretized PI formulations, in the cartesian coordinate space [Thirumalai, Bruskin and Berne, 83], and in the discretized phase space by using the Fast Fourier Transform (FFT) alogorithms [Kosloff and Kosloff, 83; Hellsing, Nitzan and Metiu, 86]. Secondly, there are the Monte Carlo based approaches to the multidimensional integrals; by Fourier Coefficient path-integration (FPI) methods [Coalson, Freeman and Doll, 86], by random walk methods [Miller, Schwartz and Tromp; 83], and by using Metropolis Monte Carlo schemes [Behrman, Jongeward and Wolynes, 83]. In section 6.1, we briefly review a discretized PI formulation in the cartesian coordinate space known as the Numerical Matrix Multiplication (NMM) method [Thirumalai, Bruskin and Berne, 83], and discuss its shortcomings due to which it becomes impractical for use in physical problems involving more than one degree of freedom. Section 6.2 describes the evaluation of the complex time memory term y(r) introduced in the previous chapter, and introduces a restricted matrix 75 multiplication (RMM) method to efficiently compute a complex time propagator for the primary region dynamics. By deforming the integration contour in the complex time plane, we describe in section 6.3 an efficient transformation from real to complex time, and show that the efficiency of the RMM method can further be increased by using this new transformation. 6.1. Numerical Matrix Multiplication (NMM) Method and its Shortcomings For simplicity, we consider only a Cartesian Hamiltonians in one dimension, i.e, the form of the general Hamiltonian is expressed as = +V (k) (6.4) 2mx where mx is the reduced mass for the degree of freedom along which k and PX are the associated position and momentum operators. The solution to the timedependent Schrodinger equation for the time-independent Hamiltonian of Eq.(6.4) is obtained by writing the coordinate space propagator as K (X, X', fl) = (X e X') (6.5) where [according to Thirumalai, Bruskin and Berne; 83, and also for simplicity] we are considering the propagator along the imaginary time axis such that t = -ioh. For a time-independent Hamiltonian, using exp(-2KH) = exp(-KH).exp(--k), it is easy to show that +00 K (X, X')K) I dX" K (X, X"; K) K (X", X'; tc) (6.6) -00 76 Starting with K(X, X'; K), one can iterate Eq.(6.6) N times to get K(X, X'; 2AK) in terms of K(X, X'; K). If K is chosen such that K = 4, then N iterations yield the desired propagator in Eq.(6.5). If N is sufficiently large then K is small enough to write the short time approximation for K(X, X'; K) as K (X,X'; K) = exp (X - X') - {V(X) + V (XI)} 2r h K 2h'K 2 (6.7) The infinite integral in Eq.(6.6) is first replaced by a finite integral; this can be done firstly, because the diagonal terms in Eq.(6.7), for which X'= X, are exponentially damped by the contributions from the potential energy factor exp{--KV(X)} for large values of X under the presumption that one is dealing with binding potentials, and secondly, because the far off-diagonal terms are damped by the kinetic energy factor exp{-(mx/2h2K)(X-X')2}. A grid is constructed of M+1 equally spaced points over this finite range. If A is the spacing between these equally spaced grid points, and using the trapezoidal rule for the finite integral in the coordinate space, each integration in Eq.(6.6) then becomes a simple numerical matrix multiplication expressed as M+1 K (iA, jA; 2K) = A K (iA, lW; K) K (lA, jA; K) (6.8) 1=1 where the real positive nature of K ensures the numerical stability of the procedure. The advantages of using the NMM for the integral in Eq.(6.6) are its ability to deal with any binding potentials, and the fast convergence towards the final imaginary time Oh because during each iteration the initial imaginary time K increases by a factor of 2. 77 The main disadvantage of the NMM method is that the number of grid points needed to span the appropriate region of a one dimensional potential is typically of the order of 102. In one dimension, each NMM step for a 100 x 100 matrix requiring about 106 mathematical operations, is very well within the limits of present day computers. However, if one were simply to extend this procedure to a two dimensional problem, the matrix size of 104 x 104, with each NMM step requiring roughly 1012 operations, starts to reach the limits of present day computing capabilities. Therefore, a straightforward extension of the NMM method to physical problems, involving more than one or two degrees of freedom, is impractical. One suggested way of improving upon this is to use some higher order quadrature instead of the simple trapezoidal rule. Given the appropriate weight factors for the selected grid points in the coordinate space, one uses a far smaller number of the total grid points in each NMM step. This improvement, however, has never been implemented in any numerical computations till now, and its usefulness is more in doubt if one is using the NMM method for computing the propagator for a time-dependent Hamiltonian. Thirumalai and Berne recommend the use of the NMM method only for attractive or binding potentials, because otherwise the number of grid points to effectively sample a repulsive potential in the coordinate space becomes too large. So that even in one dimension this increase in the number of grid points in the NMM procedure is reflected in a quadratic decrease in the computational efficacy of the overall method. Lastly, the fast convergence to the final imaginary time 0 for the problem is achieved only for time-independent Hamiltonians; the method as described above does not work for a time-dependent Hamiltonian. 78 6.2. The Complex Time Propagator for the Primary Region Dynamics The NMM method for a one dimensional time-independent Hamiltonian, can be easily extended to computations in the complex time plane by replacing in the above description the imaginary time -i0h by a complex time r = t - iuh with u as a constant real parameter. Therefore, the short imaginary time K = is replaced by the short complex time e =g. The complex time propagator matrix for one dimensional time-independent Hamiltonian is then written as K(r) = AN-1.K(E)N-1 (6.9) where if N is sufficiently large and e sufficiently small, the short time approximation for the propagator; Kij(e) = K(iA, jA, e), is written as MX (m 2mA2( *2 + S)= I2 i fe) D (i - j) exp 2 - _ 2 (Vi + V) (6.10) where A is the spacing between any two successive grid points in a (M + 1) x (M + 1) matrix over a range R given by R = MA, and V and V are the potentials at the ith and the jA grid points respectively. The contributions from the decay factor D(i - j), given by D (I - exp 2 N XA (i J) ] , (6.11) are equal to one for the on-diagonal terms for which i = j, and decay rapidly as exp(-C(i - j)2) for i , j, where the decay constant C depends upon the choice of the parameters M and N. The contributions from the off-diagonal terms in the 79 short time complex matrix K(e) can be increased or decreased by decreasing or increasing the value of the decay constant C. In the TDSCF approach the effective frequency operators, which govern the time-evolution in the primary region dynamics, are time-dependent. Therefore, the NMM method in which the initial complex time e increases by a factor of 2 during each iteration, does not work. However, the complex time-evolution in the primary region amplitude (X, 7-) can still be written in a discretized pathintegral formulation by discretizing a finite region in the coordinate space within which the primary region amplitude (X, r) evolves on the upper potential energy surface. For a primary region dynamics, evolving with a time-dependent effective frequency operator, the discretized PI formulation is then written as a time-ordered matrix product K(r) = I(N-1) K(N-1)(E).K(N-2)(E)...K(0)(E) (6.12) where e =g, and K(')(E) is a complex time matrix at the 1ih step in time with 1 going from 0 to N - 1. The progression from the 1th step to the (I + 1)1h step is simply a finite complex element matrix multiplication. We note that the complex time TCF Fbb(-r) and its Fourier transforms Fb,b(w) are computed on a contour along which the the parameter u, which defines the imaginary component of the complex time r = t - iuh, remains a constant. As shown in Fig.1, the complex time TCF and its Fourier transform are computed along a path C along which the parameter u remains a constant and the real time t' varies from 0 to t -- oo. To effectively compute the Fourier transform along the path C, one needs complex time TCF for many values of t' between 0 and t. The complex time 7'= t' - iuh corresponding to each of these values of t' is therefore divided by an integer N and the time-ordered matrix multiplication for 80 the primary region dynamics is performed at each instant along the path C'. To complete the discussion about the computations of a complex time propagator for the primary region dynamics, it is natural that we now consider the evaluations of the complex time effective frequency operators along the path C'. The time-dependence in the effective frequency operator n P(7r') along a path C' in the complex time plane, as described in equations (5.3) and (5.15) of the previous chapter, is through the coupling between the primary and secondary region dynamics. For simplicity, we consider intermolecular processes in which the effective coupling at the 11h step along C', is written as q (71) = q (0) cos (wr1) + P (0) sin (wrl) + y (r) (6.13a) rnw where the memory term y(7-1) is an integral in the complex time plane TI y (r1) - dc v (t,) sin w (r - t,) (6.13b) 0 and at the 11h step 7- = le' with E' = . Provided the time-dependent effective force v(tc) is an analytic function in the complex time plane, the value of the complex integral in the memory term is independent of the path C'. Among many possible approaches to the evaluation of the integral in the complex time plane we have successfully used two alternate ways to compute the memory term in Eq.(6.13b). 81 N -uli IMT I ReT C, C2 1 T''t'-iU? T-t-iuE C Figure 1. Computations of the primary region propagator in the complex time plane. 82 The first approach uses the analytic property of the memory term in the complex time plane; according to which, the memory term is first computed along the negative imaginary time t, = -iu'h for many different values of 0 < u' < i, and then uses these values at u' = u, to analytically continuing the memory term to the complex time plane at T- = t; - ituh by using a numerical Pa-de' approximant (Schlessinger's Point Method [Appendix B]) scheme. Numerical analytic continuation for the complex time computations has also been used for computing the Flux-Flux TCFs [Miller et al., 83] for reaction rates , and the dipoledipole TCFs [Thirumalai and Berne; 83] for spectroscopic properties. However, the numerical analytic continuation methods though accurate for short real time processes, with t < uh, are not very accurate for physical processes involving long real times t >> uh. Therefore, in our second approach to the computations of the integral in Eq.(6.13b), we avoid the numerical analytic continuation and write the integration of the memory term in the complex time plane as a sum of the integrals along the paths C' and the path C, i.e (...)dt =(...)dc+ (...)dte (6.14) 0 CC2 where as shown in Fig. 1, the path C' at the real time t = 0, along the negative imaginary axis involves a real integral from 0 -- -iulh, and is followed by the path C2 along which the imaginary time -iujh remains constant while the real time varies from 0 -+ tj. Substituting in the first term tc = -iu'h, and in the second term tc = t'- iulh, the memory term y(r) in Eq.(6.13b) at the complex time r = tj - iulh is written as 83 U' y(rj) - du' v (-iu'h) sinhu (uj - u') h 0 dt' v (t' - iuih) sinw (t1 - t') (6.15) 0 and the integrals in the above equation are evaluated numerically. Using the same procedure in the complex integrals in Eq.(5.20), we can also simplify the propagator for the secondary region dynamics as K"S (Q, Q' r') = eXp[ M (27rih sin Lor' 2h sin Wr' x { (Q2 + Q'2) cos wr' - 2QQ' - 2Q1 (7') 2Q' 2 MW 2 (r) - ;22 13 (7')} + (El - Ej7) r'] (6.16a) where I1 (r') = - du' v (-iu'h) sinh (wu'h) 0 + dti" v (t" - iuh) sin (t" - wuh) (6.16b) 0 I2 (r') = - du' v (-iu'h) sinhw (u - u') h 0 + J dt" v (t" - 1uh) sinw (t - t") (6.16c) 0 84 and 13 (r') = Jdu' (-iu'h) sinhw (u - u') h du"v (-iu"h) sinhwu"h 0 0 U - dt"v (t" - iuh) sinw (t' - t" ) du'v (-iu'h) sinhwu'h 0 0 + Jdt"v (t" - iuh) sinw (t' - t") 0 x dt"'v (t"' - iu'th) sino (t"' - iuh) (6.16d) 0 The complex time computations of the primary region dynamics for a dipole-dipole TCF and its Fourier transform are illustrated in Fig. 2. The dipole-dipole TCF Fb)(r') is computed along C for many points r' - iuh. The path-integral method, through an efficient restricted matrix multiplication (RMM) method (to be described in the next paragraph) for the primary region dynamics, is used along the path C' joining 0 and r'. The memory terms needed to perform the RMM procedure along C' are computed numerically according to Eq.(6.15) by using the force term v(7) at the grid points shown in Fig.2. The grid points needed to compute the memory term at any complex time 71 are the points for which the corresponding complex times are less than rj; therefore, the force terms at these grid points are already known from the previous computations. 85 IMT t' ReT C1 Figure 2. Complex time computations for the TCF and its Fourier transform; discretization of the complex time plane. 86 Lastly, we explain the efficient Restricted Matrix Multiplication (RMM) method mentioned in the preceding paragraph. As mentioned before, even in a one dimensional problem, the straightforward numerical matrix multiplication (NMM) [Thirumalai and Berne, 83] is suitable mostly for binding or attractive potentials where the range of the matrix sampling the potential is fixed by the range of the potential. Therefore, for the nonbinding and repulsive potentials, the method becomes unsuitable even in one dimensional problems. We use two properties of the short time approximation to the complex matrix K(e) in Eq.(6.10) to efficiently incorporate a linear increase with the range of the potential. Firstly, the short time approximation for the complex time matrix element Kij(e) as written in Eq.(6.10) is invariant to the interchange of i and j, i.e Kij(c) = KT(E). Therefore, by using the property that the product of two symmetric matrices is a symmetric matrix, we reduce the comptational efforts in any NMM iteration by a factor of two. More importantly, the second property comes from separating, in the short time approximation for the complex time matrix element Kjj(e), a completly real decay factor D(i - j). For diagonal elements, i.e for j = i, the decay factor D(0) is one. For j : i the decay factor, D(i - j) cx exp(-C(i - j)2), starts to cause a rapid decay in the magnitudes of the off-diagonal matrix elements as one goes away from the diagonal. Since suitable grid spacings, are generally achieved by keeping the decay constant C within a range 0.1 < C < 1.0 [Thirumalai, Bruskin and Berne, 83; Coalson, 85], the corresponding decay factor for far offdiagonal matrix elements, i.e (i - J) ;> 15 is typically of the order of - 10-40. This means that the contributions from the far off-diagonal matrix elements rapidly approach zero as one goes further away from the diagonal. Instead of using the 87 full matrix in an NMM step, we use the rapidly decaying property of the decay factor D(j - i), and suggest a natural cut-off in the contributions from the far off-diagonal matrix elements. The name Restricted Matrix Multiplication (RMM) hence refer to the NMM procedure where the contributions from the far off-diagonal elements beyond a suitable cut-off are chosen as zero. If M is the suitable cut off beyond which we are ignoring the contributions from the off-diagonal elements in a (Al + 1) x (M + 1) matrix, the number of mathematical operations needed to perform a restricted matrix multiplication (RMM) step are typically of the order of ~ (M + 1) x M2 as opposed to the operations ~ (M + 1)3 needed in an equivalent NMM step. Therefore, as the size of the matrix starts to increase from (MI+1) x (M+1) to (M+L+1) x (M+L+1), the RMM step for the computations of the same accuracy needs only - L x VI extra operations, whereas for the similar increase in the matrix size the NMM step will need L extra operations. In typical applications both M and L are generally of the order of ~ 100 whereas the cut-off parameter M is generally of the order of - 10. Therefore, specially for unbound potentials, the efficacy of the RMM procedure is manyfold as compared to that of the NMM procedure. 6.3. A More Efficient Transformation from Real to Complex Time As mentioned before the transformation from real to complex time along a contour for which the imaginary part -iuh remains a constant, is not unique. We also note that having separate computation paths, for the complex time Fourier transform along C and for the RMM steps along C', is not very efficient. This is because the computations being done along two separate contours require information on many more points on the complex time grid, as compared to a situation where C 88 and C' coincide, i.e where the Fourier transform computation is done along the RMM computation path. By deforming the path C for the complex time Fourier transform, we show that there exists another transformation, for which the Fourier transform contour overlaps with the RMM contour, making the computations for this new transformation more efficient. The line shape function for a general two surface photointeraction process is related to the complex time Fourier transform tihw; +00 iI Fbb (w) = dt' e"oF ) (t' -- iuh) (6.17) -00 of the complex time dipole-dipole TCF F(r') with ' = - iuh. We can rewrite Eq.(6.17) as FDD (w) = lim Re dt' e-~"hPbL, (t' - ilh) (6.18) 0 t - 00 where we have used F(-t' - iuh) =F(t' - iuh)*, and for constant u we have taken the exponential factor containing u inside the integral. As shown in Fig.(3), deforming the integration path C we rewrite the integral as dz'eIwz'Fj (z') = dz'esz'fgb (z') + Jdz'eswz'fbp (z') (6.19) C Ci Cc where the complex integration variable z' = t' - iuh. 89 IUT t ReT C re-ict -ui re C Figure 3. Deformation of the Fourier transform contour. 90 In the first term we use z' = -iu'h to show that /dz'eoz'F b (z') = ih fdu' ewno'Fb (-iu'h) (6.20) Ci 0 where Fb p(-iu') is real; therefore, the contribution to the line shape function from the first term in the deformed path, Ci, is equal to zero. In the second term, we use r = Vt2 +u2 and a =arcsin u with r -+ oo for constant small angle a, corresponding to t - oo for constant u, and we rewrite the second term as e-0a dz'ewz' fg, (z') = J dz' esoz' f (z') (6.21) Cc 0 Replacing the integration variable z' = r' exp(-ia) for constant a, and rearranging the integrand we get Fbb (w) = lim e-C J dr' e-Wr'-e' F b, (r'e-i-) (6.22) 0 r -+ oo where the Fourier transform integration and the complex time dipole-dipole TCF F(r'exp(-ia)) for a constant positive angle a, are evaluated along the same path Cc. We also note that, instead of first translating the integration range in the lower half of the complex time plane by a constant amount -iuh, and then deforming the resultant path as C = C1+Cc to achieve the new integration path in Eq.(6.22), we could have directly obtained the desired transformation of the integrals merely by rotating the real time axis by a constant negative angle a = - arcsin where u is a real parameter. CHAPTER 7 APPLICATIONS AND RESULTS In this chapter we present applications of the TCF approach to the interaction of light with a polyatomic system. Section 7.1 describes application to the electronic absorption spectra for a model one-dimensional problem with two electronic potential surfaces. The upper and the lower potential energy surfaces are modeled by displaced Morse potentials [Coalson, 85] in one dimension. Computations of the complex time TCF and its Fourier transform for this model are performed to provide a check of the convergence of the numerical path-integral method. In section 7.2, we use the dipole-dipole TCF of an extended molecular system within a TDSCF approximation to study the photodissociation dynamics of methyl iodide (CH3I) on the potential energy surface of an electronically excited state. Realistic potentials energy surfaces for this well known example are taken from a Gaussian Wave-packet Dynamics (GWD) study of the same problem [Lee and Heller, 82]. Using the TDSCF approximation in chapters 4 and 5, the photodissociation dynamics of CH3I has been factored into a primary and a secondary region dynamics. The complex-time evolution of the primary region amplitude is performed through an efficient discretized path-integral method, whereas the complex time propagator for the secondary region amplitude is constructed analytically according to Eqs.(6.16a,b,c,d) in the previous chapter. The total cross-section for the photodissociation of CH3I from its ground vibrational state, using our TDSCF approach, is compared to that of the GWD approach of Lee and Heller. The TDSCF approach to photointeraction processes in many dimension, coupled with a path91 92 integral method for the primary region dynamics, is independent of any particular functional form of the initial amplitudes for the primary and the secondary region dynamics, and of any particular functional form for the primary region potential. An additional application to the photodissociation dynamics of vibrationally excited CH3I, for which the initial values of the primary and the secondary region amplitudes are not restricted to simple Gaussian functions, are described for many vibrationally excited states in section 7.3. 7.1. Electronic Absorption Spectra of a Model Two-Potential Problem For a two surface electronic absorption problem the complex time dipole-dipole TCF fb b(t) of Eq.(3.1 lb) at the temperature (kB/)- is written in the coordinate representation as fbb (t) = dX dX' (X e+Ki/3+ X D (X') (X' e- X) D (X) (7.1) where ft and l, are the Hamiltonians for the lower and the upper potential energy surfaces, the complex time r = t - iuh, and Z# is the partition function at the initial inverse temperature kBI. Electronic excitations in a single nuclear variable problem, such as in diatomics, are described by one dimensional potential energy curves. We model the lower and the upper potential energy curves as two displaced one dimensional Morse oscillator potentials [Reimers, Wilson, and Heller, 83; Thirumalai and Berne, 85; Coalson; 85], expressed as V (X) = V ( 1 - e-"X , (7.2a) 93 Upper and Lower Potential Energy Surfaces 3 -4 0 X (a.u 8 12 Figure 4. Model one dimensional Morse potentials for the ground _ , and the excited - - -- electronic potentials for a two state electronic absorption problem (see text). 4- 140 120 100 80 60 40 20 0 -20 U W.. |

Full Text |

PAGE 1 PHOTODYNAMICS OF EXTENDED POLYATOMIC SYSTEMS By DEEPAK SRIVASTAVA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1988 PAGE 2 TO MY WIFE PAGE 3 ACKNOWLEDGEMENTS I would like to express my thanks and appreciation to some of the many people who have encouraged and supported me in this work. The foremost amongst them is my advisor Dr. David A. Micha, who constantly provided the guiding thought throughout the course of this work and also provided a part of the funding for my financial support. I am also thankful to many graduate students, postdoctoral associates, and faculty members at the QTP who have contributed in some way or the other towards increasing my knowledge in the field of theoretical and computational chemical physics. Finally, I would like to express my appreciation for all the help I have received from the secretaries and the staff at the QTP and at the Physics Department that has made my stay at the University of Florida a pleasant experience. PAGE 4 TABLE OF CONTENTS page ACKNOWLEDGMENTS iii ABSTRACT vi CHAPTERS 1 INTRODUCTION 1 2 COLLISIONAL TIME-CORRELATION FUNCTION APPROACH TO LIGHT-MOLECULE INTERACTIONS. ... 8 2.1 Transition Ratesfrom Time-Correlation Functions 9 2.2 Radiation Field TCFs for Thermal Light 13 2.1 Radiation Field TCFs for Laser Light 17 3 MOLECULAR TCFs IN THE ELECTRONIC ADIABATIC REPRESENTATION 22 3. 1 Electronic Adiabatic Representation for Molecular TCFs .... 24 3.2 Complex Time TCFs and their Fourier Transforms 28 3.3 Properties of the Complex Time TCFs and their Fourier Transforms 34 4 FACTORIZATION OF TCFs INTO PRIMARY AND SECONDARY REGION TCFs 39 4.1 Sufficient Conditions for the Factorization of the Amplitude . . 40 4.2 Norm Conserving Time-Dependent Variational Principle(TDVP) 43 4.3 Time-Dependent Self-Consistent Field (TDSCF) Approximation 45 4.4 Primary and Secondary Region Factorization of TCFs 50 5 TIME-DEPENDENT EFFECTIVE FREQUENCY OPERATORS 52 5.1 Intramolecular Processes: Strongly Coupled Primary and Secondary Regions 53 5.2 Intermolecular Processes: Weak Coupling Limit 62 5.3 Relation to the Generalized Langevin Dynamics Approach. . . 64 5.4 Treatment of Many Degrees of Freedom in the Primary Region Dynamics 69 IV PAGE 5 6 PATH-INTEGRAL APPROACH TO THE PRIMARY REGION DYNAMICS 73 6.1 Numerical Matrix Multiplication (NMM) Method and its Shortcomings 75 6.2 The Complex Time Propagator for the Primary Region Dynamics 78 6.3 A More Efficient Transformation from Real to Complex Time . 87 7 APPLICATIONS AND RESULTS 91 7.1 Electronic Absorption Spectra of a Model Two-Potential Probelem 92 7.2 Photodissociation of CH3I from its Ground Vibrational State. . 100 7.3 Photodissociation of CH3I from Initially Excited Vibrational States 120 8 DISCUSSION AND CONCLUSIONS 137 APPENDICES A SECOND ORDER TRANSITION RATES 152 B NUMERICAL ANALYTIC CONTINUATION 155 C SCF EIGENVALUES AND EIGENFUNCTIONS 157 BIBILIOGRAPHY 159 BIOGRAPHICAL SKETCH 164 PAGE 6 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy PHOTODYNAMICS OF EXTENDED POLYATOMIC SYSTEMS by DEEPAK SRIVASTAVA December 1988 Chairman: Prof. David A. Micha Major Department: Physics A general time-dependent quantum mechanical approach to the interaction of visible and UV light with extended polyatomic systems is presented in this work. The interaction of light with a large polyatomic system is treated as a two-step process: a photon absorption excites electronic transitions in the target system, and this is followed by a dynamical evolution of the system on the excited potential energy surface. Transition rates for photon-molecule interaction processes are described in terms of time-correlation functions (TCF) of transition operators. These are written as a product of TCFs of the target and of the light source. The time evolutions in these two TCFs are carried out independently, and the statistical properties of the light source are explicitly incorporated in the formalism. The procedure is developed for resonance absorption-emission and for scattering of thermal, chaotic and coherent light. The time evolution in a large polyatomic system is treated within a molecular TCF approach. For a general two-surface electronic excitation problem, a transformation of these molecular TCFs from the real arguments to complex time is VI PAGE 7 described so that their computation by a numerical path-integral method become feasible. The complex time dependence in these TCFs is described by the dynamical evolution of an amplitude under the influence of a frequency operator. A time-dependent self-consistent field (TDSCF) approximation is used for large polyatomic system to factor the molecular TCFs into primary and secondary region TCFs. The primary and secondary regions are modeled by considering general primary motion coupled to many harmonic degrees of freedom in the secondary region. The strong coupling between the primary and the secondary regions has a general dependence on the variables of the primary region, whereas it has linear and bilinear terms in the variables of the secondary region. The weak coupling limit is considered by dropping the bilinear terms in the coupling. The complex time propagators for the secondary region dynamics for both of these cases are constructed analytically, while an efficient path-integral method is developed for the dynamics of the primary region. As applications, the path-integral method for the primary region dynamics is used for computing the electronic absorption spectra of a model two-potential problem, and the TDSCF approach to the dynamics of large systems is used for studying the photodissociation dynamics of CH3I from ground and excited vibrational states of the ground potential energy surface. The cross-sections for both of these cases are compared with other available works with good results. Extensions of the method to the other extended systems, with many degrees of freedom, and involving tunneling and barrier crossing in the primary region dynamics, are briefly discussed. vii PAGE 8 CHAPTER 1 INTRODUCTION Theories of light absorption-emission, elastic and Raman scattering are described by time-dependent perturbation expressions in terms of transition dipole moments, and polarizability tensors [Kramers and Heisenberg, 25], respectively. The interaction of atoms and molecules with quantized radiation field was first introduced by Dirac [Dirac, 27]. Since then many developments have been made along the same lines and have included higher order corrections in these theories to increase accuracies [Herzberg, 50]. The basic approach of these developments is to concentrate on small atoms, accurately compute all of their eigenvalues and eigenstates, and use them in a summation according to the Fermi Golden Rule [Messiah, 62] for computing the total cross-section for the required photointeraction process. In case of polyatomics, using the Bom-Oppenheimer approximation [Balint-Kurti, 75], the electronic part of the Schrodinger equation is first solved for fixed nuclear positions to obtain electronic potential energy surfaces; these are then used for the nuclear part of the problem. Given electronic potential energy surfaces for a small polyatomic system, one can still obtain all the low-lying vibrational and rotational eigenstates and eigenvalues to calculate cross-section. However, cross-sections processes which involve rearrangement in the molecule due to the absorption of light are not easily described by the method of summing over all the accessible quantum states of the system. This is mainly because one needs exact scattering wavefunctions of many channels in the excited state of the system, and the number of such available channels can be very large even for small polyatomics. Still, the time1 PAGE 9 2 independent quantum mechanical method, based on summation over all the state to state information, and known as the stationary coupled channel (CC) approach, has been used successfully to produce total cross-sections for light-molecule interaction processes involving rearrangement in the target [Balint-Kurti and Shapiro, 81; Heather and Light, 83]. However, as the size of the system starts to increase, methods such as these involve a larger and larger number of available quantum states of the system, and therefore soon become too cumbersome to be of any practical use. Alternatively, efforts have been made to change the approach from the timeindependent quantum mechanical theories, which involve summations over all of the stationary states of the system, to the time-dependent quantum mechanical descriptions which are based on the dynamical evolution of the system from given initial conditions to the required final conditions [Lee and Heller, 79; Berens and Wilson, 81; Thirumalai and Berne, 83; Mukamel et al., 85; Srivastava and Micha, 87; Coal son, 87]. These methods do not require summations over all the states of the system, and therefore are better suited for dealing with large molecular systems. The basic feature of these approaches is to start from the expression for the total cross-section in the energy domain, use the Fourier transform definition of the energy conserving delta function to make a transformation to the time domain, and use the time-dependent Schrodinger equation for the dynamical evolution of the system between given initial and required final conditions. In recent years, many advances have been made in solving the time-dependent Schrodinger equation for nuclear vibrational and rotational motion with given initial conditions for an arbitrary general potential in many dimensions [Kosloff, 88 for a recent review]. They can be catagorized into the following broad groups: (a) semi-classical methods [Marcus, 72; Miller, 74; Mukamel and Jortner, 76; Micha, 83; Billing and Jolicard, PAGE 10 3 84; Gray and Child, 84; Henriksen, 85] are generally based on running classical trajectories for the nuclear motions, whereas the electronic transitions are treated quantum mechanically; (b) approximate quantum mechanical approaches based on wavepacket dynamics [Heller, 75, 81; Agrawal and Raff, 82; Mowrey and Kouri 85; Sawada and Menu, 86], where the mean position and the spread of the wavepacket follow a canonical set of equations which are derived from a timedependent variational principle [Heller, 81]; (c) straightforward integration of the time-dependent Schrodinger equation by discretizing the phase space and using a Fast Fourier Transform (FFT) algorithm [Feit, Fleck and Steiger, 82; Kosloff and Kosloff, 83; Hellsing, Nitzan and Metiu, 86; Kosloff, 88]; and (d) path integral approaches [Thirumalai and Berne, 83; Miller, Schwartz and Tromp, 83; Behrman et alÂ„ 83; Coalson, Freeman and Doll, 85] which provide complex time propagators between given initial and final positions in the coordinate space. The present developments in the time-dependent descriptions still have the following problems which must be addressed. By restricting the radiation field part of the interaction Hamiltonian to simple monochromatic circular wave forms, these approaches do not explicitly incorporate the statistical and dynamical properties of the light source in the formalism. 2. Semiclassical methods for solving the time-dependent Schrodinger equation are based on classical trajectories which do not sample the non-classical regions of the potential energy surfaces. 3. Approximate quantum mechanical techniques based on wavepacket dynamics require local quadratic forms for potentials, hence they exclude extensive delocalization of wavepackets. PAGE 11 4 4. The numerical path-integral methods, for a potential of arbitrary functional form, are not useful in the calculation of the required real time dynamics because of the oscillatory nature of the real time propagator. 5. Straightforward integration of the time-dependent Schrodinger equation by discretizing the phase space and using FFT for propagation are practical only for problems with few spatial dimensions. The time-correlation function approach to thermodynamic and electromagnetic response properties of extended molecular systems is well known in statistical mechanics [Forster, 75; McQuarrie, 76; Lovesey, 80]. Recently, there have also been efforts to describe molecular processes occuring in a target following its collision with a projectile, using a collisional time-correlation function (TCF) approach [Micha, 77, 79, 81, 86]. Since, as before, this approach does not require summations over all the states of the target, it is useful in dealing with many degrees of freedom in extended molecular targets [Vilallonga and Micha, 87; Micha and Parra, 88], In this work we present a formalism and applications of a timedependent quantum mechanical approach to light-molecule interactions in large systems with the following desirable features. 1. We treat the radiation field and the molecular target of the interaction Hamiltonian on equal footing so that, at any stage in the formalism, explicit statistical and dynamical properties of the light source are easily incorporated. 2. Within a mean field approximation, we separate the degrees of freedom of an extended molecular system into a primary region ( containing a few important degrees of freedom along which all significant dynamical events take place) and a secondary region (with all of the remaining degrees of freedom). 3. For a general potential of arbitrary functional form in the primary region, we use PAGE 12 5 a complex time path-integral method which can incorporate important quantum effects such as tunneling and barrier crossing in this region. 4 . The computational efforts for the method do not increase quadratically with an increase in the range and the dimensionality of the potential. In chapter 2 we describe transition rates for general light-molecule interaction processes in terms of the TCF of the LippmanSchwinger T operator. In the dipole approximation to the light molecule interaction, the T operator is expanded into terms containing factors for the radiation field and the molecular target. The transition rates for all first order processes such as resonance absorption and emission are written in terms of the molecular dipole TCF of the target, and the electric field TCF of the light source. For all second order processes, such as multiphoton absorption-emission, elastic and Raman scattering, the transition rates are written in terms of the TCFs of the molecular polarizability tensors. The statistical properties of the light source are incorporated in the formalism, and specific cases of thermal, chaotic and coherent light are discussed. In chapter 3 we introduce electronic potential energy surfaces by using the Bom-Oppenheimer approximation to separate electronic and nuclear motions in a molecular target, and write a general molecular TCF in the electronic adiabatic representation. We describe a transformation from real to complex time of these TCFs, and discuss some of their general properties. The complex time dependence in the TCF, for a general light-molecule interaction process involving electronic transitions in the target, is shown to be equivalent to the dynamical evolution of an amplitude (equal to the initial eigenstate times the transition dipole) induced by a frequency operator (an operator depends upon the energy difference between the upper and the lower potential energy surfaces divided by h). PAGE 13 6 In chapter 4 we describe the time evolution of the amplitude introduced in the previous chapter through a norm conserving time-dependent variational principle (TDVP). Considering fast and localized electronic tranisitions in a large molecular target, we use a time-dependent self-consistent field (TDSCF) approximation to factor the molecular dipole TCF of the target into self-consistent primary and secondary region TCFs. In chapter 5 we describe time-dependent effective frequency operators for many intraand intermolecular photointeraction processes by considering an arbitrary general motion in the primary region dynamics coupled linearly and bilinearly to many degrees of freedom in a secondary region dynamics undergoing harmonic motions. We show a direct relationship between our TDSCF approach and the phenomenological Generalised Langevin Equation (GLE) [Mori, 65] approach to the stochastic dynamics in large systems. We explicitly prove a FluctuationDissipation Relationship (FDR) between the fluctuation force and the memory kernel in our TDSCF approach for linearly coupled primary and secondary region dynamics. We also develop the formalism for many degrees of freedom in the primary region dynamics coupled linearly and bilinearly to many harmonic degrees of freedom in the secondary region dynamics. In chapter 6, we describe details of an efficient path-integral method for the primary region dynamics. Using a discretized coordinate space, the computational method does not depend upon any specific functional form of the primary region potential and also avoids problems of other methods which scale quadratically with respect to the range and the dimensionality of the potential. In chapter 7 we describe applications and results of the TCF approach to the interaction of visible and UV light with extended molecular systems. We first de- PAGE 14 7 scribe an application to the electronic absorption spectra of a model one-dimensional problem for two-potential energy surfaces. The complex time computations for the TCF in a single spatial dimension, using the complex time numerical path-integral method, are performed so as to provide a check of the accuracy of the numerical method for the primary region dynamics. We study the photodissociation dynamics of a realistic polyatomic system (CH3I) in an electronically excited state within a TDSCF approximation by factoring the TCF for the process into self-consistent primary and secondary region TCFs. The resulting total TCF has been used to compute the total cross-sections for the photodissociation of CH3I, which are then compared with results from other methods with good results. We have also extended the method to the computations of total cross-sections for many vibrationally excited states of CH3I on the ground potential energy surface. Finally, in chapter 8 we discuss advantages and distinctive features of our approach and mention possibilities for further improvements and future extensions of our work. PAGE 15 CHAPTER 2 COLLISIONAL TIME-CORRELATION FUNCTION APPROACH TO LIGHT-MOLECULE INTERACTIONS We treat the interaction of light with a molecular system as a photon scattering process through a collisional TCF approach [Srivastava and Micha, 87]. Starting from the transition rates for a collision process in terms of the Lippmann-Schwinger transition operator T [Messiah, 62] in section 2.1, we derive the transition rates for the interaction of light with a molecular system in terms of the TCF of these T operators. In the dipole approximation to the light-molecule interaction, the T operator is written as an expansion with terms containing radiation field and molecular system factors. For first order processes such as resonance absorption and emission, transition rates are derived in terms of molecular dipole TCF of the target and the electric field dipole TCF of the light source. For second order processes such as elastic and inelastic (Raman) scattering due to polarization fluctuations in the target, transition rates are written in terms of the polarizability tensor TCFs of the molecular system and the polarizability tensor TCFs of the light source. In section 2.2 we treat the interaction of molecules with thermal light, which is given in the photon number representation by a Bose-Einstein distribution of the initial number of photons in a particular mode. Section 2.3 considers the interaction of a molecule with a single mode laser. The laser light source is treated by considering the radiation field to be in the coherent state representation. The first order transition rates for laser-molecule interaction processes are derived when a) the laser is operating below the threshold (chaotic light with a poisson distribution 8 PAGE 16 9 for the mean number of photons), and b) the laser is operating above the threshold (coherent light with a statistical distribution for the phase and the mean number of photons given by the experimental conditions). 2.1. Transition Rates from Time-Correlation Functions We consider the interaction of visible and UV light with a molecular system. The treatment is general so that the molecular system can be as large as a protein molecule, a polymer chain, an adsorbate on a surface, or a solid surface. We denote with Hm, Hr, and H] the Hamiltonians of the molecular system, the radiation field, and their interaction, respectively. The total Hamiltonian is expressed as H = H m + Hr + Hj. (2.1) Eigenstates | j) and | r), corresponding to the noninteracting molecular system Hamiltonian Hm and the radiation field Hamiltonian Hr, respectively, are given by Hu\i) = Ej\j) ( 2 . 2 ) A*|r) = E, |r). We treat the interaction of light with a molecular system as photon scattering, during which the radiation field goes from an initial state | r) to a final state | r/). Later on we shall describe various processes like absorption or Raman scattering by specifying | r) and | r/). Transition operators T+ describing the scattering satisfy the Lippmann-Schwinger equation [Messiah, 62] T+ = Hj + Hj K .-T+, E + ieH m Hr (2.3) where e Â— > 0 is a limit to be taken after matrix elements of these operators have been calculated. Rates of transition R, averaged over the statistical distribution wj PAGE 17 10 of the initial target states, w r of the initial radiation states, and summed over final states of the target and the radiation, | jr) and |r/>, respectively, are expressed as (2.4) Â»' rr' where E = E T + Ej is the total energy, and the Dirac 8 function imposes energy conservation. By using the Fourier transform of the 8Â— function 6{E -E') = -^-r e^ E E ^ dt, (2.5) 2tt h J_ oo and the completeness relation j' r ')U' r '\ = ( 2 6 ) )' r> Eq.(2.3) can be written as [Micha, 81] 1 f + Â°Â° = dt ((T + (t)'T + (0)))M,R (2.7) where the time-dependence in the T + is given by f + = exp(-^of) f + exp(+^ 0 f) (2.8) ri h and H 0 = Hm + Hr, is the total Hamiltonian for the molecular system and the radiation field without interaction. The double bracket PAGE 18 11 where D is the dipole vector of the molecular system, and E is the electric field vector of the light source. The transition operator T of Eq.(2.2) is written up to the second order, with an implicit + sign, as T = ~ XI EriDv 1 E + ie Â— Hm Â— Hr -rEcD, < U C ( 2 . 10 ) Â£ vC where (Â£, 77 , Q = (x, y, z), and the time-dependence of the f operator is given by Eq.(2.8). Since Hm and Hr commute with each other, it is straightforward to take matrix elements of the first term in Eq.(2.10) and to separate them into radiation field and molecular system factors. In the second term we use the following expression 1 1 r +o Â° Â— -7 = / d\ a + b 7 r J_ 00 a 2 + A 2 b 2 + A 2 Â’ valid also for complex values of a and b to write ( 2 . 11 ) E + ie Â— Hm Â— Hr = Â•Â» r du K J00 Mm Mr M 2 m Â— fi 2 u 2 M 2 r Â— h 2 u 2 ( 2 . 12 ) where M m = E M + ie M H m (2.13 a) Mr = E R + ie M Hr (2.13 b) E = E m + E r (2.13 c) e = e M + e R , (2.13 d) and we have replaced A = ihu. Using these expressions in the second term of Eq.(2.10), to factor it into a radiation field part and a molecular system part, and the result in Eq.(2.6), we write the transition rates R up to the second order as PAGE 19 12 p + OO R = A I dt {( D \,( t ) D ( ( 0 )))M ({ E \,( t ) E t ( 0 ))) R J-0 O X) -tEZT f' *?C J ~Â°Â° p+oo / +oo -oo du I dt (( D ^,( t ) x ^( u , 0)))m (( 4 (Â«)V(Â«.o)Â»j! Â• + /*+oo r-f-oo ^ 2^2 du> I dt (( i \, c ( U ', t ) Et ( 0 ))) R V vC n 2 ' + Â°Â° + ^EE/ du f + Â°Â° du ' f + Â°Â° dt 2 i y_ 00 y.oo 7T vC ^C '7 00 (2.14) ((X'J'C'( u, ^)Xi?c( u Â’ 0 )))a/ ((^'C'( u,,< )^c( u Â’Â°)))Â« where the polarizability tensor operator for the molecular system is defined by XÂ„ c (u, 0 )-OÂ„ a ^_" v K c the polarizability tensor operator for the radiation field is defined by (2.15a) Mr *, c (u,0)-Â£,^_ fiV < Â» (2.15fc) and the time-dependence in these operators are given for example by X, c ( U ,() = e-^"Â‘Â«(x.O) e+S*"' (2.16) The expression for i? has four terms, each of them corresponding to different photointeraction processes. The first term describes first order processes such as resonance dipolar absorption and emission of light. The second and the third terms are cross-terms the contribution of which will depend upon the statistical properties PAGE 20 13 of the light source used in the experiment. For thermal light in the photon number representation, they become zero and do not contribute to the overall transition rates. In the later sections, we shall, however, see that the contribution of these terms cannot be ignored if the light source used is a coherent light. The fourth term describes second order processes such as elastic and Raman scattering mediated by polarization fluctuations in the target. The transition rates R expressed in Eq.(2.14), as opposed to the standard expressions [Dirac, 47; Heitler, 54], are specially useful for the treatment of lightmolecule interactions in large systems because they do not involve summations over excited states of the molecular system. Furthermore, Eq.(2.10), which is used to separate in the transition operator the radiation field and molecular system variables, can also be used to derive a similar separation of the transition operator to order higher than the second. Thus, a similar formulation of the transition rates for higher order photointeraction processes can be systematically derived. Lastly, the main purpose of separating the TCF of the T operator into the TCF of the radiation field and the TCF of the molecular system is to uncouple the time evolution of the molecular system from that of the radiation field. This has a twofold advantage: firstly, one can use time-dependent quantum mechanical techniques to describe the time evolution of the molecular system without referring to the time evolution of the radiation field, and secondly, we can easily incorporate the statistical properties of the light source into the dynamics. 2.2. Radiation Field TCFs For Thermal Light When the light source is maintained at a temperature T, it is described by the Bose-Einstein statistical distribution function W 8 {u} a , T), of the number of photons PAGE 21 14 in a mode s with energy ha> s . Then it is easier to work in the photon number representation of the radiation field. The electric field operator in Eq.(2.14) is decomposed as Et(t) = EÂ± +) (i) + iÂ£ \t) (2.17) where and (2.18) Â£<-Â»(()= (^ +l (())* (2.19) where k is the photon wavevector, e ^ is the Â£ component of the polarization vector, and a ^ the annihilation operator of a photon of energy = hkc and polarization Â£. In the photon number representation, a general radiation field state with several modes s is written as |{(Â£ J 6) Ar Â‘}) where N 3 is the number of photons in mode s; then % +) (*) lUU,)*'}) = i E (^) ! % ( 2 . 20 ) Using Eqs.(2. 17-2.20) in Eq.(2.14), one can write the transition rates for first and second order processes. In the photon number representation, the second and third terms of Eq.(2.14) do not contribute, because they do not conserve the number of photons in the initial and the final radiation state. Therefore the transition rate is written as PAGE 22 15 R = R u + R 2 2 ( 2 . 21 ) where Rn is the transition rate for the first order processes, in which a single electric field operator at an initial time t = 0, is correlated to its adjoint at a later time t. Similarly, R 22 is the transition rate for second order processes, where the time -dependent polarization tensor operator, <^(u, 0), which is second order in the electric field operators of Eq.(2.17), is correlated to its adjoint at a later time t. In the same notation R\ 2 and R 2 \ are the cross-terms, the contributions of which are zero in the photon number representation. The transition rates Rn for first order processes are given by where +oo i?ll = ^ f dt (t) , tt'-OO ( 2 . 22 ) (2.23) is the molecular dipole TCF of the extended system, and / { J |,(<) = ((^r ) (<)Â£< +) ( Â«)))* + ((4Â° W^Â’w))* (2.24) is the electric field TCF of the light source. The first term in Eq.(2.24) corresponds to resonance absorption, whereas the second term corresponds to spontaneous and stimulated emissions. Other cross-terms in Rn do not contribute because they do not conserve the number of photons. Using Eqs.(2. 18-2.20) in the above equation we write PAGE 23 16 26nV' + (N a + 1W e-'^ 4 ] (2.25) and the transition rate R\\ is given by ^ = ^EÂ»-EE / * w Â«Â»ko ( 2 26 ) s it, Â£ o 2i? e Â«Â£> Â£ (()' Â£> Â£ (0 )))m + jee{e-'"Â‘.Â‘((i> Â£ (<)Â’ i>< ((>)))Â«}] where the statistical probability distribution function for a thermal light W a is the Bose-Einstein distribution for N s photons with wavevector k s and frequency u^, and the summation in Eq.(2.26) over the number of photons can be carried out analytically. We have also used in Eq.(2.26); ((D^ (-t)D^(O ))) m = (0^(t)D^(O)))*, to ensure that the transition rate Rn remains real. For the second order processes, we write the transition rates R 22 as +00 +00 +00 R 22 = ^EE / du ' J du J dt VC V'C'Â—oo Â—00 Â—00 < 2 27 ) where /Â«, Â„ c (u',u,t) = ((*Â«. (Â«'.<)' X,<(Â“.0))) M (2.28) is the molecular polarizability tensor TCF of the target, and KcvC ( u> Â’ u Â’*) = ((^'C' ( 2 29 ) PAGE 24 17 is the electric field polarizability tensor TCF of the light source. Using Eqs.(2.172.20) in Eq.(2.29), the electric field polarizability tensor TCF of the light source can also be further decomposed into terms representing elastic and Raman scattering [Appendix A]. These terms when substituted back in Eq.(2.27), will give transition rates for second order processes when thermal light is interacting with a large molecular system. Equations (2.22) and (2.27) for the transition rates R u and R 22 are computationally useful, because they do not involve the summations over the excited states of the target, and the corresponding molecular TCF can be evaluated by time-dependent techniques. 2.3. Radiation Field TCFs For Laser Light When the light source is a laser, the statistical distribution of its electric field factorizes into a distribution for the electric field amplitude, times a distribution for its phase [Louisell, 73; Sargent, Scully and Lamb, 74; Loudon, 83]. If the laser is operated below the threshold for coherent emission, the amplitude fluctuations of a given mode, which are proportional to the mean number of photons in that mode, have an exponentially decaying distribution about the most probable amplitude (which occurs at zero), whereas fluctuations in the phase, caused by spontaneous emissions, are quite large and therefore the phase can have any value between 0 and 27 t. If the laser is operated above the threshold, both the amplitude and the phase have sharply peaked distribution about their most probable values, which are given by experimental conditions. For such a light source it is necessary to include in the description the statistical distributions of both the amplitude and the phase of the electric field, and it is advantageous to work in the coherent state representation PAGE 25 18 [Klauder and Sudarshan, 68; Loudon, 83]. A general radiation state in the coherent state representation is written as | {a s }), with a ^ | {o 3 })= a s | {a,})%, where the eigenvalue a s is a complex number whose square modulus gives the mean number of photons in the s mode of the laser. Using equations (2.17-2.20) for the creation and annihilation terms of the electric field operator, and using a single mode coherent states | a s ), we write The radiation field TCF of Eq.(2.22), in the coherent state representation, can be written as Unlike f^(t) in the photon number representation, cross-terms in Eq.(2.31) are not zero, because we do not need to conserve the number of photons in the initial and the final radiation states | r) and |r'), respectively. As we shall see very shortly, however, it is still possible to have a laser in which the statistical distribution of the phase is such that it again can make the contributions of these cross-terms equal to zero. Using equations (2.18), (2.19) and (2.30) in Eq.(2.31), we wnte in a single mode coherent state representation as (2.30) /& m = (( 4 _> Â« 4 +) (Â°)))Â„ + (( 4 Â° Â« 4 "Â’ (Â°>))Â„ (2.31) (2.32) PAGE 26 19 where P(a s ), a real statistical weight for the complex amplitude a s , is given by experimental conditions. The two dimensional integral f d 2 a*[...] extends over the entire complex plane of a s . Light sources of experimental significance have different statistical weights P(a s ). In the first and the fourth terms of Eq.(2.32), we factors P(a s ) given in the literature [Schubert and Wilhelmi, 86]. Therefore, the radiation field TCF in Eq.(2.32) satisfies the relation and hence, the resulting transition rates of experimental significance, are real from Eq.(2.32). As examples of laser sources, with given statistical weights, we consider the radiation field TCFs for chaotic and coherent light. To compare Eq.(2.32) with the equivalant expression in the photon number representation in Eq.(2.25), we have first considered the case of chaotic light, which is the same as that of a laser operating below the threshold. In this case the phase of the light source is randomly distributed between 0 and 27r, therefore the weight factor P(a 3 ) does not depend on the phase of a 3 , and hence, the contributions from the first and the fourth terms of the integrand in Eq.(2.32) are zero. Thus, below the threshold takes a simpler form which agrees with Eq.(2.25) for a single mode, insofar as the average is have checked that the statistical averages of a 2 , and a* 2 are real for several weight OO 0 (2.34) PAGE 27 20 Substituting Eqs.(2.32) and (2.33) in Eq.(2.21), we can write the first order transition rate R n for the interaction of a single mode laser with a molecular system, when the laser is operating below the threshold, as OO R c nÂ°Â“ c = Â£^7 E f Â£ / <Â“[(*Â•> Â“*("Â£.0 ( 0 x2Re(((D l (t)'D ( (0))) M ) (2.35) (()*Â£Â£( 0)))m)) If the laser is operating above the threshold, the rate of transition involves instead the statistical distribution P(o s ), which depends on the amplitude | a s |, as well on the phase = Arg(a s ), and is given by TtCoherent K u U k. +oo J dt J d 2 a s P(a s ) ^ Â— OO i2 2he 0 V x [2 |Q a |^ {cos(uj-t)} cos(<^< 2 Â£t)Im{((DÂ£ (f) f (0)))ji/}]
M) V* *(A-)g* (t, hi) <}] ^ k^k 1 where we have omitted specific indices within the functional dependence for notational brevity, and on the i th potential energy surface we have written the stationary state eigenvalue Ej as Ej = E p + Ej + E ps , with Ej= ^2 k Ej, and E ps = Yhk The time-dependent coefficients, which couple the primary region dynamics to the normal mode motions in the secondary region, are the normalized expectation values of the corresponding operators within the amplitudes coresponding to the secondary region dynamics. For example, the coefficient ql(t, [?/]) is of the form Â«i(*.fo]) h^lQlh* 7 ) (5.4) The time-dependent effective frequency operator Cl s e ^(t) for the motion in the secondary region dynamics is given by 6!//Â« = EÂ“Â«//W < 5 5 ) k where the effective frequency operator for the k th normal mode motion in the secondary region is written as |