Computer Simulations of Viscosity Dependent Molecular Relaxation Processes

The dependence of excited-state lifetime, fluorescence quantum yield and isomerization 
rate of organic dye molecules on solvent viscosity has been a subject of numerous 
experimental and theoretical investigations. To explain the viscosity dependence of 
excited-state lifetime in this paper the temporal behavior of the excited state population 
is calculated for several models of the molecular relaxation process by a computer 
simulation incorporating both the stochastical motion of large molecular parts in the 
excited state and relaxation transitions. The described method is applicable to calculating 
the probability of changing the electronic state as a function of time and internal 
rotation coordinate.


INTRODUCTION
The dependence of the lifetime of excited electronic states of polyatomic molecules on the viscosity of the solvent has been a subject of numerous experimental and theoretical investigations for many [1][2][3][4] years.
Relations between fluorescence quantum yield b or excited state lifetime " and solvent viscosity r/have been found for various solutesolvent combinations.
Such dependences have been discussed on the basis of conformational changes of molecules by internal rotation which are hindered due to forces of the surrounding solvent on the solute.
Oster and Nishijima investigated the viscosity dependence of the fluorescence quantum yield of diphenyl-methane dyes. Their results could be approximated by the relation -I=A+BT/rl (A, B ll9 molecular parameter, T temperature). F/Srster and Hoffman were able to fit their experimental data for triphenyl-methane dyes by b Cr/, where the molecular parameter y was found to be approximately 0.67 for the dye crystal-violet.
As found from further experiments on the basis of ps-spectroscopy the excited state lifetime z can be related with the solvent viscosity r/ by z= ?r/v over a relatively broad viscosity range. -4 Herein the molecular parameter y is seen to vary between 0.2 and 1.4, dependent on the solute-solvent combination.
Whereas in all these experiments the lifetime z increases with increasing viscosity there has been reported only one case where a decrease of z with increasing r/ has been found. 7 It should be noted that in most cases the decay in time could be described with sufficient accuracy by a single exponential.
Several theoretical attempts have been tried in order to explain the experimental data. Oster and Nishijima discussed their data in terms of the free rotational diffusion model which yields y-1. From the F6rster-Hoffman theory (FH) a value of 2results for 3' and the decay in time can be described by exp {-kot-(tr2/rtz)t3}.
A different approach was followed by McCaskill and Gilbert 4 who treated the problem as a conformational change across an intramolecular barrier under consideration of a hydrodynamic drag and a random fluctuating force.
Using Kramers approximate solution for the appropriate unabbreviated Fokker-Planck equation (FP equation) the authors were able to calculate the parameter 3'. A generalization of the theory was given by Wilhelmi 9 and Martire and Gilbert, who took hopping processes between different electronic states into account. We present in our paper the application of the Langevin-simulation proposed by Martire and Gilbert in order to solve the appropriate FP equation which has been extended in its range of application by introducing a relaxation term describing transitions between different electronic states. We discuss various molecular models differing mainly in the relaxation term in the corresponding FP equation.

MODEL OF MOTION AND CHECK COMPUTATIONS
The internal rotation of large molecular parts can be described within the framework of classical statistics by solving the appropriate FP where Pl(q, qS, t) de dq5 is the probability of finding a molecule in the excited state $1 at time at the configurational coordinate of internal motion and with angular velocity ft. U1() represents the energy of the molecule in the excited state $1, which depends on and might be considered at least in principle to be known from calculations in the Born-Oppenheimer approximation. I is the reduced moment of inertia of the molecule with respect to the internal rotation and a represents a drag coefficient which is taken to be proportional to the viscosity r/. The transition rate klo(q, qS), which might depend on and qS, takes into account that the molecule can leave the excited state. In our treatment we assume that only the motion along one configurational coordinate is decisive for the radiationless deactivation and furthermore that the solute-solvent interaction can merely be described by a hydrodynamic drag--represented by a friction term--and a random fluctuating force F(t) in the equation of motion where the latter is due to buffeting by solvent molecules. According to the derivation of the FP equation (1) this force has to be -correlated: f oo F(t)F(t + s) dt= 2kTIa8(s). (2) Our treatment is therefore limited to such cases, where the characteristic time of a configurational change is considerably larger than the correlation time of the motion of the solvent molecules.
In this case an approximate solution of the FP equation (1) a where each actual realization of F(t) has been generated on the computer using one of the standard techniques. 4 '11 In dealing with the full FP equation containing relaxation transitions one is often not interested in the probability distribution P(, cb, t) but in the total probability N(t) of finding a molecule in the excited state at time t: This occupation probability of the excited state can on the other hand be obtained by recording the distribution of the times at which the individual hopping processes to the ground state occur. In order to consider transitions to the ground state and to obtain an approximate solution for N(t) the Langevin standard techniques #'l have to be modified in such a way that besides the simulation of the random fluctuating force a simulation of the stochastic hopping process has got to be incorporated. For the additional simulation there exist two possible methods: 1. During each time step (t, +/xt) a homogeneously distributed random number y (0, 1) is generated and compared with the instantaneous transition probability ko((t), b(t))" at. If during one particular time step y exceeds koAt, the hopping process is assumed to take place actually and the time is recorded. 2. The second method, which allows some considerable savings in computational time starts from a subensemble of molecules, which all have the same or nearly the same trajectories. For this subensemble an accumulated transition probability 'o ko((t'), (o(t')) dt' is compared with a random number-lny (where y is homogeneously distributed in (0, 1)), representing a value for the accumulated transition probability and thus for a given trajectory a random actual lifetime, which will be recorded. The obvious advantage of this method is that only one random number has to be generated to simulate the hopping process.
In this way the decay function of the excited state population can easily be obtained by recording the actual lifetime of each trajectory in the excited state which yields N(t) (compare Equation (4)).
For the one dimensional damped harmonic oscillator on which a random fluctuating force is acting, Langevin equation and the corresponding FP equation have been solved. 8 At this special case the quality of our method of simulation could be checked. in the computations are typical for molecules in which a group of the size of about one benezene ring is taking part in the internal rotation.
Since we treat the internal rotation classically we have to restrict ourselves to relative soft torsional modes (up to O0 cm -l at room temperature). The numerical solution of the FP equation (1) by simulation is a computer time consuming procedure. To avoid the computer time from growing too excessively we had to restrict most of the calculations to comparatively small excited state lifetimes (<20 ps) which means that the parameter ko(, if) was in most simulations chosen to be larger compared with typical molecular situations. Treating in the following section special molecular models we discuss how one has to proceed in order to apply the method to systems with longer excited state lifetimes, having in mind a possible fitting to experimental data.

FORSTER-HOFFMAN MODEL
First we apply the method described above to the well-known F6rster-Hottman model (FH), 6 in which displaced harmonic vibrations in both the ground state and the excited state are assumed (Figure 2). The relaxation process is thought to consist of an IC-process (ko) independent on a f-dependent radiationless process, the transition rate of which can be described by: This assumption seems to be quite reasonable, as in the range of very large r/ the internal rotation will freeze and the relaxation to the ground state will be completely determined by ko.
For the depopulation of the excited state FH gave the following ad hoc equation: dN(t)= -(ko + kl(q tpo)2)Nl(t) dt (7) In our opinion Eq. (7) seems to be questionable in as much a derivation of it from the appropriate FP equation seems impossible.
For this reason we treated the FH-model by solving the appropriate FP equation (1) Figure 3, where both a and k as parameters have been varied. It is clearly to be seen that after reaching a "thermal equilibrium" the decay functions can be described by a single exponential with sufficient accuracy. Only in the case of N(t) very large values of k a deviation from the exponential behavior within the first ps can be recognized. Such high values of k yield however very short excited state lifetimes (1,..., 2 ps). Typical values of the excited state lifetime cover the region of 10 ps up to 500 ps, which can be simulated by taking relative small values of k.
In these cases the temporal behavior of the decay can be accurately described by a single exponential. The exp (-t3) dependence of the decay in time calculated by FH can therefore not be confirmed by us.

MODEL OUT-OF-PLANE ROTATIONS
This model has very much similarity with the FH-model. However in the present model we consider a molecular system in which the equilibrium positions in the ground state So and the excited state S are identical whereas the transition probability is once again given by" k,o() ko + kl((P (Po) 2 (7) which is reasonable for the same reasons mentioned in the discussion of the FH model. In establishing this model we have in mind a number of plane molecules which exhibit a considerable decrease of fluorescence quantum yield in lower viscous media which can be explained by a deactivation through out-of-plane vibrations. The treatment with the method of simulation once again yields decay functions which can be described accurately by single exponentials. The dependence of the excited state lifetime, which has been derived from the slope of the logarithmic decay function In (N(t)/N(0)), on the viscosity is plotted in Figure 4. It is to be seen that over a viscosity range of one decade the functional dependence can be approximated by z Cr/v with y 0.6 in the particular case of the chosen parameters (ko, k,, a, I, to).

NON-ADIABATIC TRANSITIONS AT CROSSING-POINTS OF THE ELECTRONIC ENERGY
In connection with the investigation of the viscosity dependence of the isomerization rate of cyanine dyes non-adiabatic transitions at potential crossing-points are discussed in some papers. '3 COMPUTER SIMULATIONS OF VISCOSITY 127 FIGURE 4 Excited state lifetime (model out-of-plane rotation) versus viscosity. Parameters: o 5 ps-" I 5 10 -44 kgm 2" ko =0, kl 13.9 ps-1.
In this model it is assumed that the transition probability to the ground state (or lower excited state) is negligible for all co-ordinates o with exception of o +, where a crossing of the electronic energies of the two different states occur ( Figure 5). The transition probability P(5) has been calculated to "12 (8) where b is the angular velocity of the internal rotation at the crossing point, e12 is the energetic split-off of the potential energies in the crossing-point when perturbation theory has been taken into consideration, F is given by: dU dU2 (9) In our computation we have set C "-47r2e22/F the influence of this parameter is plotted in Figure 6. The decay functions plotted here show once again an exponential behavior. The dependence of the excited state lifetime on the viscosity is presented in Figure 7. In the functional dependence we can distinguish three regions" (i) Relative high r/: to/c << l: The curve can be approximated by r r/v with y 0.79 in the particular case.
(ii) Medium r/: to/al: Here the lifetime is almost independent of r/.
(iii) Small r/: to/a > Here r increases rapidly with decreasing r/, which can be understood within the concept of the thermal activation out of a potential well.
To prove the thermal character of the activation of the non-adiabatic  Parameters: as in Figure 6, C 0.6. The Eq. (10) is a strong indication for a thermal activation process out of a potential well. For the case of thermally activated adiabatic conformational changes across a potential barrier ' (11) where v and Voj are vibrational wavefunctions of the considered soft torsional mode in the excited state and the ground state respectively. E, Eoj are the total vibronic energies. Alo is factorized into the electronic part of the transition moment and into matrix elements of promoting modes. 5 It is assumed that Alo is a (, qS)-independent factor. The overlap integral is calculated within the semiclassical approximation. 16 The co-ordinate * is the one (or the ones) for which the angular momentum p(*)=p* is conserved p(i)=p(j). The FCfactor is then given by: where F IdU/dq dUo/dql*; e(i), eo(j) are the vibrational energies of the torsional vibration in the excited state and in the ground state respectively, i, j are vibration quantum numbers. Transitions take place at co-ordinates q* for all states i,j with Eli Eo.
Let us consider two cases.
COMPUTER SIMULATIONS OF VISCOSITY 131 (a) In the first one we assume that energy resonance can be reached for a large number of vibrational levels. We demand the energy resonance with vibrations of the electronic ground state to be fulfilled within an energy range of some kT, or in other words, for all thermically occupied vibrational levels of the excited state energy resonance might be reached. The electronic energy gap is overcome mainly by vibrational quants of "hard" modes, quasi resonance is finally reached due to a sufficient broadening of the vibrational levels which is supposed to result from interaction with a heat bath.
In our computations for this model, which have been carried out in the same manner as described above and taking expression (11) for klo, under the assumption of displaced harmonic potentials once again exponential decay function have been obtained. The decay time as a function of viscosity shows a decrease with increasing r/(resp, a) over a relative broad viscosity range.
(b) In the second case we assume a very sharp resonance. That means there exists only one "transition co-ordinate" R and only one "transition angular momentum" p p R (if the angular frequencies of the torsional vibrations in the excited state and the ground state are ditterent). Thus the transition rate ko has to be altered in this case to: kl={; nst frP=tpR'p=PRotherwise Because of the resemblance with the model of non-adiabatic transitions at potential crossing-points a similar viscosity dependence of the excited state lifetime might be expected. Figure 8 actually confirms this expectation. It should be emphasized that in the computations of this model excited state lifetimes of 20,..., 120 ps have been simulated, which lie within the range of typical decay times of a number of dyes. TM