Article pubs.acs.org/JPCA

Fourth-Order Perturbative Model for Photoinduced Internal Conversion Processes Brian P. Molesky and Andrew M. Moran* Department of Chemistry, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina 27599, United States S Supporting Information *

ABSTRACT: Essential to the functionality of numerous biological and synthetic molecular systems is the ability to rapidly convert electronic excitation energy into heat. Such internal conversion (IC) transitions often cannot be described by traditional second-order kinetic theories because of time-coincident electronic and nuclear relaxation processes. Here, we present a perturbative fourth-order phenomenological model for photoinduced IC that incorporates effects associated with finite laser bandwidths and nonequilibrium nuclear motions. Specialized knowledge of first-principles computational methods is not required, and many parameters can be obtained with standard spectroscopic measurements. The model is applied to the IC processes that precede electrocyclic ring-opening in α-terpinene. It is shown that the primary factor governing the shape of the population decay profile (Gaussian versus exponential) is the rate at which the wavepacket approaches the geometry corresponding to degeneracy between the excited states. Other parameters such as the displacement in the promoting mode and the thermal fluctuation amplitudes affect the sensitivity of the IC dynamics to motion of the wavepacket but do not alter the basic physical picture. Finally, we suggest a wavepacket representation of the IC process to visualize correlations between population-transfer dynamics and the amount of energy transferred from the system to the bath.

I. INTRODUCTION Internal conversion (IC) transitions abound in nature and technology. Molecules including DNA nucleobases, natural melanin pigments, provitamin D, and carotenoids rapidly convert electronic excitation energy into heat by way of IC.1−7 Such ultrafast relaxation protects many biological systems from photodamage in addition to facilitating photosynthesis. Studies of IC mechanisms have a long history in chemical physics.8−17 Nonetheless, innovative experimental approaches continue to uncover new insights into IC processes.18−27 The development of experimental techniques and growing consensus on the ubiquity of conical intersections in molecules have also motivated further theoretical work in this area.28−39 Modern theoretical methods provide valuable physical insights, particularly for molecules in the gas phase, where computational resources can be focused on the molecular degrees of freedom. However, fully capturing the complexity of systems in condensed phases remains a challenge despite steady progress in this field of research. In this paper, we present a phenomenological model for IC in which many of the parameters can be obtained by spectroscopic measurements. The model improves on (perturbative) second-order rate formulas by incorporating effects associated with finite laser bandwidths and nuclear relaxation processes that occur on the same time scale as IC. We are motivated by the ability to capture basic physical insights without specialized knowledge of first-principles methods. Phenomenological models developed with similar motivations have proven to be valuable tools for interpreting © 2013 American Chemical Society

experimental data acquired for molecules in condensed phases.40,41 A general aspect of such models is the use of a reduced description, where only a handful of “interesting” degrees of freedom are treated explicitly, and the remainder of the system is relegated to a heat bath composed of displaced harmonic oscillators.42,43 Written in this way, rate formulas for electron- and energy-transfer processes can be parametrized using a variety of spectroscopic measurements.41,44−49 Our development of a phenomenological model for IC is motivated by several recent investigations of DNA components and ring-opening systems.51−56 In this work, the model will be applied to one of these systems, alpha-terpinene (α-TP), which exhibits particularly unusual relaxation dynamics.53 To understand the basic relaxation scheme for α-TP, it is instructive to consider the potential energy surfaces of cyclohexadiene (CHD) shown in Figure 1.19,50 Photoexcitation near 270 nm populates a ππ* excited state (1B in CHD). Population then rapidly transfers into a lower-energy, optically forbidden state (2A in CHD) characterized by a doubly excited lowest unoccupied molecular orbital (LUMO).19,50 The branching space around the conical intersection involves ring deformation and CC stretching motions. In CHD, the two paths around the conical intersection are symmetrically equivalent. However, this symmetry is broken in derivatives of cyclohexadiene such as α-TP, thereby enabling the detection of separate wavepacket Received: August 7, 2013 Revised: November 28, 2013 Published: December 2, 2013 13954

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

specific aspects of the system (e.g., vibronic couplings, solvation) influence this particular IC mechanism. Below, we demonstrate the utility of the model by addressing a basic question about α-TP. Can a pair of trajectories in α-TP (i.e., parallel IC channels) give rise to dynamics resembling the population decay profile extracted from our measurements? This is an intriguing point because the time scales associated with the two components are so similar. The relevance of such time scales can be understood by considering that a set of integro-differential equations must generally be solved to describe the IC process43 d Pa(t ) = dt

∑ {∫

t

dt ′Wab(t − t ′)Pb(t ′) −

0

b≠a

∫0

t

dt ′Wba(t − t ′)Pa(t ′)}

(1)

Figure 1. Potential energy surfaces computed for CHD in ref 19. States 1B and 2A correspond to single and double HOMO-to-LUMO excitations, respectively. (b) Structure of α-TP. The electronic states of α-TP are thought to be similar to those found in CHD.19,50 However, the symmetry of the two paths around the conical intersection, which are indicated by the red arrows in panel (a), is broken by the alkyl substituents in α-TP.

Nonexponential kinetics may arise whenever populationtransfer dynamics are comparable to the time scale of electronic dephasing. The rate kernels Wab(t) and Wba(t) encode the (temporal) waveform of the coherence (i.e., dephasing dynamics) between the initial and final states. The interesting aspect of anisotropy decay in α-TP is that two components with fairly similar time constants are detected, but only one exhibits a pronounced nonexponential shape. The model developed here explores how this can occur when similar dephasing rates are associated with the two trajectories.

trajectories. We hereafter refer to the two lowest-energy excited states of α-TP as S1 and S2, which are the analogues of states 2A and 1B in CHD, respectively. It was recently suggested by Sension and co-workers, who have carried out extensive studies of CHD and its derivatives in solution,7 that a closely related system, provitamin D, undergoes a direct IC transition from the ππ* state to the ground state (i.e., relaxation may not proceed through an intermediate dark state as in CHD).57 Discussion regarding whether or not a direct S2 → S0 IC channel also contributes (either partly or fully) in α-TP is beyond the scope of this work. Pinpointing the dominant relaxation mechanism (S2 → S0 versus S2 → S1 → S0) is a challenge because clear spectroscopic signatures of sequential IC transitions in α-TP will be obtained only if a significant amount of population accumulates in the intermediate state S1. That is, in order for a straightforward assignment to be made, the first transition (S2 → S1) in the two-step sequential process (S2 → S1 → S0) must occur more slowly than the second transition (S1 → S0). In the present work, we focus only on the physics associated with IC between excited states (S2 → S1). The transient absorption anisotropy decay of α-TP can be decomposed into a dominant Gaussian-like component with a 100 fs HWHM and a less prominent 130 fs exponential component.53,57 The origin of these two components is not clear a priori. One possibility is that each component represents a different path around the conical intersection between excited states S2 and S1. Consideration of this mechanism takes inspiration from the parallel relaxation pathways detected in a closely related system, provitamin D, which is also distorted from C2 symmetry in both the ground and excited states.57 A second possibility, previously considered by us,53 is that the two components correspond to depopulation of states S2 and S1, respectively. Finally, the decay profile may represent a single trajectory in which nonexponential dynamics originate in timecoincident IC and nuclear relaxation processes. The first possibility now seems likely based on the agreement between the time scales of the anisotropy decay in the deep UV and the decay of an excited-state absorption nonlinearity associated with S2 in the visible spectral range.58 The reduced description presented here provides a framework for thinking about how

II. FORMULATION OF MODEL This section presents models that can be used to simulate the IC dynamics at both second-order and fourth-order in perturbation theory. We begin by partitioning the Hamiltonian into system and bath components. Time correlation function approaches are then used to generate rate formulas. The (conventional) second-order rate formula treats the special case in which the nuclear coordinates relax on a time scale that is fast compared to IC, whereas effects related to photoinduced nuclear motions are incorporated at fourth-order. Finally, we summarize key assumptions that are particular to the present model. IIA. Hamiltonian. The three components of the full Hamiltonian, H = Hsys + Hbath + Hsys−bath, are given by ∞

Hsys = |a⟩⟨a|

∑ |α⟩⟨α|[Ea + ℏαωp]

α=0 ∞

+ |b⟩⟨b|

∞

∑ |β⟩⟨β|[Eb + ℏβωp] + |c⟩⟨c| ∑ |ξ⟩⟨ξ|[Ec + ℏξωp] β=0

ξ=0

(2)

Hbath =

1 2

⎡ p2

⎤ + miωi2qi2 ⎥[|a⟩⟨a| + |b⟩⟨b| + |c⟩⟨c|] ⎢⎣ mi ⎥⎦

∑⎢ i

i

(3)

and Hsys − bath = |a⟩Q a⟨a| + |b⟩Q b⟨b| + |c⟩Q c⟨c|

(4)

As shown in Figure 2, the electronic structure of the model system resembles that of α-Tp in that a is the ground state, b an optically active state, and c an optically forbidden state with intermediate energy. Associated with each electronic state is a manifold of vibrational energy levels corresponding to the promoting mode: α is the number of vibrational quanta in state a; β and δ belong to state b; and ξ represents vibrational levels in state c. Notably, the vibronic indices associated with different 13955

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

population (β = δ) or vibronic coherence (β ≠ δ) within electronic state b. The final two interactions with the nonadiabatic coupling transfer population from state b to state c. The waveform of the coherence between states b and c (in t3) governs the kinetics in a manner similar to the dephasing dynamics that take place within the domain, t′, in eq 1. The interesting aspect of the fourth-order description is that the waveform in t3 depends on the time-evolving geometry of the nonequilibrium system in τ. IIB. Second-Order Rate Formula. Before considering the fourth-order process, it is useful to examine a traditional second-order rate formula obtained using the correlation function approach and notation particular to this work. The second-order rate constant is obtained by time-integrating a correlation function involving the nonadiabatic coupling operator V̂ K (2) = Figure 2. Model system considered in this work. Electronic states a, b, and c represent the ground state, an optically bright state, and an optically dark state, respectively. Greek symbols denote vibrational levels associated with each of the electronic states.

2 ℏ2

∫0

∞

dt1 ⟨V̂ (0)V̂ (t1)⟩

(5)

For illustration, we consider a transition between vibrational level β of electronic state b and vibrational level ξ of electronic state c. In this basis, the rate constant is given by Kc(2) ξ , bβ =

electronic states are not interchangeable because the equilibrium position of the (harmonic) promoting mode is displaced (see Figure 2). System−bath interactions are described with a Brownian oscillator model, which has been discussed in detail elsewhere.42,43,49,59 In the present model, it is assumed that fluctuations in the electronic energy levels are large compared to fluctuations of the vibrational energy levels. Thus, the collective coordinates in Hsys−bath depend only on the indices of the electronic states. These Brownian oscillator coordinates are related to the harmonic modes in the bath by Qa = ∑imiω2i daiqi, where dai is the displacement of mode i (i.e., coupling strength) with respect to the equilibrium position in the ground state. Dynamics are introduced by letting the system interact with an external radiation field, Ĥ rad−mat = −μ · E and the nonadiabatic coupling operator (see Appendix A). At secondorder, the IC transition between states b and c can be modeled without explicitly considering light absorption under the assumption that nuclear relaxation is fast compared to IC. In contrast, the fourth-order Feynman diagram shown in Figure 3 indicates that state b is populated after two interactions with the laser field; t1 is the time interval in which the system absorbs light. The system then evolves in t2 in either a vibronic

2 ℏ2

∑ ∑ Pbβ ∫ β

∞

0

ξ

dt1 ⟨b|⟨β|V̂ (0)|ξ⟩|c⟩⟨c|⟨ξ|V̂ (t1)|β⟩|b⟩ (6)

where Pbβ is a Boltzmann population. By assuming that the fluctuations in the electronic energy gap are large compared to fluctuations in the vibrational energy levels, the rate can be written as Kc(2) ξ , bβ = ⟨exp−[i

2 ℏ2

∫0

∑ ∑ Pbβ |Vbβ ,cξ|2 ∫ β t1

ξ

dt ′Ucb(t ′)]⟩

0

∞

dt1 exp[i(ωcb + ωξβ )t1]

(7)

where the solvation operator, Ucb, represents fluctuations of the energy gap with respect to the mean value, ℏωcb, at the equilibrium geometry of state b.42 The correlation function in eq 7 can be evaluated with a cumulant expansion by treating the bath that surrounds the system as a collection of displaced harmonic oscillators.42,59 For systems in solution, low-frequency, thermally driven motions of the solute and solvent are often handled using a single collective bath coordinate (see Appendix C). Application of the cumulant expansion yields42

Figure 3. (a) Double-sided Feynman diagram for the second-order rate constant. (b) Double-sided Feynman diagrams corresponding to the two fourth-order rate kernels. The time intervals ti separate perturbative interactions with an external electric field, Ĥ rad−mat, and the nonadiabatic coupling V̂ . In W1j(τ,t3) and W2j(τ,t3), the interval between the peak of the laser pulse and the first interaction with V̂ , τ, is closely related to t2; τ = t2 for a delta function laser pulse. 13956

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A 2 ℏ2

Kc(2) ξ , bβ =

∑ ∑ Pbβ |Vbβ ,cξ|2 ∫ β

0

ξ

Article

In the process of interest, two field−matter interactions occur before nonadiabatic coupling induces population flow between excited states. Correlation functions involving three time intervals must be evaluated at fourth-order to describe this sequence. The fourth-order rate constant for photoinduced IC can be written as

∞

dt1

* (t1) − g *(t1) exp[i(ωcb + ωξβ + λbb /ℏ + λcc /ℏ)t1−gbb cc + 2gbc*(t1)]

(8)

where gbc * = 0 under the assumption of uncorrelated thermal fluctuations in states b and c. The first-order terms in the cumulant expansion yield the reorganization energies in the argument of the exponential, λbb + λcc, because ωcb corresponds to the equilibrium geometry of state b. Using the definition of (2) the nonadiabatic coupling given in Appendix B, Kcζ,bβ becomes ℏωp

Kc(2) ξ , bβ =

μp

⟨c|

∫0

* (t1) dt1 exp[i(ωcb + ωξβ + λbb /ℏ + λcc /ℏ)t1 − gbb (9)

where it is assumed that the system possesses a single promoting mode, p, and ωξβ = (ξ − β)ωp. (2) It is useful to consider a simplified version of Kcξ,bβ for physical insight. We choose a form that parallels rate formulas for electron transfer because such models are well-known in the chemical physics community.44 First, it is assumed that the frequency of the promoting mode is large compared to kBT in order to restrict the summation over vibrational quanta β to β = 0. In addition, the line shape functions can be written as gbb(t) = kBTλbbt2/ℏ2 if the fluctuation amplitude in ωbc is large compared to the inverse of the bath’s relaxation time.42 Evaluation of the integral in eq 9 yields 2

Kc(2) ξ , b0 =

ωpℏ

|V bcel|2

2μp

π (λbb + λcc)kBT

∫0

∞

dt 2

W1j(τ , t3) = 2 ∑ ∑ ∑ ∑ Paα

∫0

∑ |⟨ξ|1⟩|2

⎡ (ℏω − ℏω − λ − λ )2 ⎤ bc ξ0 bb cc ⎥ exp⎢ − 4(λbb + λcc)kBT ⎢⎣ ⎥⎦

dt1

∫0

∞

dt3 ⟨Ĥ rad−mat(t1)

where Ĥ rad−mat represents the field−matter interaction and V̂ is the nonadiabatic coupling; tj are intervals between interactions with Ĥ rad−mat and V̂ (see Figure 3). Below, rate kernels associated with K(4) will be used to solve a set of integrodifferential equations of the type given in eq 1. Non-Markovian effects will be treated by integrating over only t1 and t2, thereby retaining the waveform associated with the coherence between the two states involved in the nonradiative transition in t3. Here, the term “non-Markovian” is used to indicate that the dephasing time scale is not short compared to the time scale of IC. The key mathematical issue is that the upper limit of the integral over t3 should not be infinite in the non-Markovian regime.43 Rate kernels associated with the two Feynman diagrams shown in Figure 3 are readily obtained by substituting the nonadiabatic coupling operator for the final two field−matter interactions in standard four-wave mixing response functions.42 The rate kernels are

β ⟨ξ|β − 1⟩|2

− gcc*(t1)]

∞

(11)

ξ ∞

∫0

+ ⟨Ĥ rad−mat(0)V̂ (t1 + t 2)V̂ (t1 + t 2 + t3)Ĥ rad−mat(t1)⟩

∂ |b⟩ ∂R p β + 1 ⟨ξ|β + 1⟩ +

2 Re ℏ4

V̂ (t1 + t 2)V̂ (t1 + t 2 + t3)Ĥ rad−mat(0)⟩

2

∑ ∑ Pbβ |− β

K (4) =

∞

dt1

∫0

α ∞

β

δ

μbβ , aα μaα , bδ Vbδ , cξ , jVcξ , bβ , j ℏ4

ξ

Re{

dt 2 exp[−i(ωbat1 + ωβαt1 + ωβδ t 2 + ωbc , jt3

+ ωβξt3)− f1j (t1, t 2 , t3)] E(τ − t 2 − t1) E*(τ − t 2)}

(12)

ξ

and W2j(τ , t3) = 2 ∑ ∑ ∑ ∑ Paα

(10)

Equation 10 shows that thermal fluctuations enter both electron-transfer and IC processes in similar ways. The key difference between IC and electron-transfer processes is that they are subject to different types of vibronic progressions. That is, under the present approximations, electron transfer and IC involve overlap integrals of the type ⟨ξ|0⟩ and ⟨ξ|1⟩, (2) respectively; it may be instructive to compare Kcξ,b0 to eqs (V.8) and (V.9) in ref 44. IIC. Fourth-Order Rate Formula. The second-order rate formula obtained in Section IIB generally provides an inadequate description of light-activated IC processes that (2) occur on the 100 fs time scale. First, Kcξ,b0 fails to incorporate photoinduced motions along Franck−Condon active coordinates. Second, evaluating the correlation function with an infinite integration limit tacitly assumes that the dephasing of the coherence between the initial and final states is fast compared to the transfer of population.43 The phenomenological model presented in this section aims to overcome these limitations without introducing an intractable number of empirical parameters.

∫0

∞

dt1

∫0

α ∞

β

δ

μaα , bδ μbβ , aα Vbδ , cξ , jVcξ , bβ , j

ξ

ℏ4

Re{

dt 2 exp[−i(ωabt1 + ωαδ t1 + ωβδ t 2 + ωbc , jt3

+ ωβξt3) − f2j (t1, t 2 , t3)] E*(τ − t 2 − t1) E(τ − t 2)}

(13)

where μaα,bδ are transition dipoles, and the line shape functions f1j(t1,t2,t3) and f 2j(t1,t2,t3) are discussed in Appendix C. We have included an index j in anticipation of treating a system that can relax by way of parallel IC channels (i.e., trajectories). In principle, each trajectory j can possess fully independent parameters. However, as discussed in the following section, differences in particular parameters are introduced only if needed to capture prominent aspects of the measured signals. The electric field, which is regarded as scalar here, is given by E(t ) = ε exp[−iωLt − w 2t 2/2ℏ2]

(14)

The above forms of W1j(τ,t3) and W2j(τ,t3) can be simplified to reduce computational expense and enhance physical insight. For brevity, we outline steps only for W1j(τ,t3). To begin, W1j(τ,t3) is written with explicit vibrational overlap integrals as 13957

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A 2 μp2

W1j(τ , t3) =

Article

and

∑ ∑ ∑ ∑ Paα |μba |2 |Vbcel ,j|2 ⟨β|α⟩⟨α|δ⟩ α

β

δ

ξ

∫0

⟨δ|∂/∂R p|ξ⟩⟨ξ|∂/∂R p|β⟩Re{

∞

dt1

∫0

d Pcj(τ ) = dτ

∞

dt 2

−

exp[−i(ωbat1 + ωβαt1 + ωβδt 2 + ωbc , jt3 + ωβξt3) − f1j (t1 , t 2 , t3)]E(τ − t 2 − t1) E*(τ − t 2)}

2 μp2

β

δ

∫0

∫0

∞

dt1

∞

dt 2exp[−i(ωbat1 + ωβαt1 + ωβδt 2 + ωbc , jt3 + ωβ t3)

− f1j (t1 , t 2 , t3)]E(τ − t 2 − t1) E*(τ − t 2)}

(16)

The index α is additionally eliminated by assuming that the frequency of the promoting mode is large compared to kBT/ℏ (Pa0 = 1). W1j(τ,t3) can then be written as W1j(τ , t3) =

2 μp2

∑ ∑ |μba |2 |V bcel , j|2 ⟨β|0⟩⟨0|δ⟩⟨δ|∂/∂R p exp(iHct3/ℏ) β δ ∞

∫0

∂/∂R p|β⟩Re{

dt1

∫0

∞

dt 2 exp[− i(ωbat1 + ωβ 0t1 + ωβδt 2

+ ωbc , jt3 + ωβ t3) − f1j (t1, t 2 , t3)]E(τ − t 2 − t1) E*(τ − t 2)}

(17)

In the final step, we substitute the auxiliary function, Φδβ(t3) (defined in Appendix B) into W1j(τ,t3) ωp W1j(τ , t3) = ∑ ∑ |μ |2 |Vbcel ,j|2 ⟨β|0⟩⟨0|δ⟩Φδβ (t3) ℏμp β δ ba

∫0

∞

dt1

Re{

∫0

∞

dt 2 exp[−i(ωbat1 + ωβ 0t1 + ωβδt 2

+ ωbc , jt3 + ωβ t3) − f1j (t1 , t 2 , t3)]E(τ − t 2 − t1) E*(τ − t 2)}

(18)

Analogous operations are carried out on W2j(τ,t3) to obtain ωp W2j(τ , t3) = ∑ ∑ |μ |2 |Vbcel ,j|2 ⟨β|0⟩⟨0|δ⟩Φβδ(t3) ℏμp β δ ba

∫0

Re{

∞

dt1

∫0

∞

dt 2 exp[−i(ωabt1 + ω0δt1 + ωβδt 2

+ ωbc , jt3 + ωβ t3) − f2j (t1 , t 2 , t3)]E*(τ − t 2 − t1) E(τ − t 2)}

(19)

With the rate kernels in hand, population dynamics involving states b and c are computed using the following set of integrodifferential equations, d Pbj(τ ) = − dτ +

∫0

∫0

τ

dt Wcb , j(τ , τ − t )Pbj(t )

τ

dt Wbc , j(τ , τ − t )Pcj(t )

τ

dt Wbc , j(τ , τ − t )Pcj(t )

(21)

(22)

The rate kernel Wcb,j(τ,t3) is made independent of the incident light intensity under the assumption of a “normalized excitation probability”, which is justified by the correspondence of second-order and fourth-order rate constants at large τ (see Supporting Information). This approximation makes it possible to simulate nonradiative population dynamics using eqs 20 and 21. Sensitivity to photoinduced nuclear motions originates in the dependence of Wcb,j(τ,t3) on the time τ. The rate kernels are continually “updated” to reflect the time-evolving geometry of the system. In the calculations presented below, it is assumed that Wbc,j and Wcb,j are related by detailed balance. This is an adequate approximation when dephasing in t3 (i.e., the domain of integration variable t′ in eqs 20 and 21) is fast enough to suppress oscillatory behavior. Finally, we comment on the decision not to assume that dephasing between states b and c is fast compared to the time scale of the nonradiative transition (i.e., Markovian dynamics are not assumed). It is safest not to make this assumption for αTP because the two time scales are separated by less than an order of magnitude. However, as shown in the Supporting Information, calculations of the population dynamics become much simpler if such a separation in time scales is invoked. When dephasing is fast compared to the IC transition, an analytical form of the rate constant can be obtained and population dynamics can be computed with ordinary differential equations instead of integro-differential equations. This may be the best approach in situations where IC dynamics are fast and qualitative insights are sought. IID. Key Approximations and Limitations of Model. This section further discusses two key assumptions that are particular to the present model. Notably, the reduced description employed in this work also has many inherent approximations, which have been discussed in detail elsewhere.42,43,49,59 The Feynman diagrams presented in Figure 3 assume a sequential process in which light absorption precedes the nonradiative transition between excited states. This assumption is valid when the laser pulse duration is short compared to the time scale of the IC transition. Models for conventional transient absorption spectroscopies invoke similar approximations in the regime of well-separated pump and probe pulses (e.g., doorway−window model).42 The 100 fs time scale of the fastest transition in α-TP is roughly 7 times larger than the laser pulse duration employed in ref 53. Thus, it is reasonable to limit our consideration to the sequential radiative and nonradiative processes indicated in Figure 3. A primary goal of the model is to describe the interplay between time-coincident nuclear relaxation and IC dynamics. The dependence of the rate (at which state c is populated) on the laser intensity is not of interest. Therefore, we have introduced a normalized excitation probability to remove dependence on the laser intensity. The condition is defined as

ξ

⟨δ|∂/∂R p|ξ⟩⟨ξ|exp(iHct3/ℏ)∂/∂R p|β⟩Re{

dt Wcb , j(τ , τ − t )Pbj(t )

Wcb , j(τ , t3) = W1j(τ , t3) + W2j(τ , t3)

∑ ∑ ∑ ∑ Paα |μba |2 |Vbcel ,j|2 ⟨β|α⟩⟨α|δ⟩ α

τ

where

(15)

The propagator for nuclear motion in the time interval, t3, is next written to the right of the bra, ⟨ξ|, to make clear how the sum over index, ξ, can be removed using the closure relation, ∑ξ|ξ⟩⟨ξ| = 1, W1j(τ , t3) =

∫0

∫0

(20) 13958

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A |ε|2 |μba |2 π [w 2(2λbbkBT + w 2)]1/2

Article

data. The parameters of the absorbance line shape were determined in ref 53 using a cumulant expansion; however, the symmetric CC stretching mode is treated explicitly as part of the system here for consistency with the model for IC presented above. That is, the intramolecular promoting mode is not incorporated with a cumulant expansion in the present work. Rather, the absorbance line shape is computed using

=1 (23)

where |ε|2 is the laser intensity, w the spectral width of the laser pulse, and λbb the reorganization energy for an optical transition to state b. In the Supporting Information, this approximation is motivated by correspondence between the second-order and fourth-order rate formulas at times τ that are long compared to the nuclear relaxation time scale Λ−1 bb,j. This correspondence is based on three key assumptions. First, it is assumed that the laser is tuned near the peak of the electronic absorbance spectrum. Second, the slow modulation regime of line broadening, (2λbbkBT)1/2 ≫ ℏΛbb,j, must be satisfied (i.e., Gaussian line shapes).42 Third, the IC transition rate must be slower than the time scale of dephasing between states b and c (i.e., a Markovian process). It is acknowledged that the third assumption is marginally satisfied in α-TP, where the dephasing and IC time scales differ by less than an order of magnitude. A more general way of eliminating the dependence on the laser intensity should be developed for systems that do not meet these three requirements.

σ(ω) = Re

∫0

∞

dt ⟨0|0(t )⟩ exp[i(ω − ωba)t − gbb(t )] (24)

where the form of the dynamic overlap integral ⟨0|0(t)⟩ is given in the Supporting Information.60,61 The absorbance spectrum constrains the line shape function (gbb(t)), the transition frequency (ωba), the frequency of the promoting mode (ωp), and the displacement in the promoting mode (d). Parameters are given in Table 1. Table 1. Parameters Used in Model Calculations parameter

value 35 400 cm−1 2.0 1450 cm−1 6.9 amu 1500 cm−1 −1 37 630 cm−1 425 cm−1 trajectory 2

ωba/2πc db ωp/2πcc μpc λbb/hcd ηd ωL/2πce w/hce trajectory 1 a

III. MODEL CALCULATIONS In this section, we apply the above model to a three-level system with properties resembling those of α-TP. The goal of these calculations is to establish general connections between IC dynamics and the model parameters. In particular, we will consider the time scale of overdamped wavepacket motion, displacement in the promoting mode, and the amount of correlation between energy levels of the excited states. IIIA. Establishing a Set of Parameters. Figure 4 overlays the measured absorbance spectrum of α-TP with the spectrum of the laser employed in our experimental work.53 Calculated versions of these spectra are shown along with the experimental

parameter

value

parameter

value

d Λ−1 bb,1 ωbc,1/2πca Velbc,1a

500 fs 3000 cm−1 225 nm−1

d Λ−1 bb,2 ωbc,2/2πca Velbc,2a

200 fs 1500 cm−1 180 nm−1

a c

Equations 18 and 19. bSections III and IV in Supporting Information. Appendix B. dAppendix C. eEquation 14.

The region of parameter space to be explored in our model calculations is further constrained by considering the experimentally determined population decay profile (i.e., not the raw signal). In Figure 4b, the population of state b, Pb1(τ) + Pb2(τ), is overlaid on the empirical decay profile F(τ ) = 0.5 exp[− (τ /120 fs)2 ] + 0.5 exp[−τ /130 fs] (25)

which was determined in ref 53. The profile F(τ) cannot be “fit” in the traditional sense because of computational expense. Therefore, parameters of the model were manually adjusted until good agreement was achieved. All parameters associated with each of the two trajectories j were initially set equal, and differences in particular parameters were introduced only if needed to capture the shape of the empirical decay profile. This procedure produces well-constrained parameters and is sufficient for investigating the basic physics of IC in α-TP. Effects of particular parameters on the population dynamics are discussed below. IIIB. Mechanics of Ultrafast Internal Conversion. We begin by investigating the interplay between energy detuning ωbc,j and the time scale of the bath Λ−1 bb,j. For context, it is useful to first consider the second-order rate constant presented in Section II, which makes clear how dynamics in the promoting mode and the bath coordinate are related to energy detuning (2) effects. From the perspective of the bath coordinate, Kcξ,b0

Figure 4. (a) Experimental and calculated laser spectra and linear absorbance line shapes for α-TP. (b) Population decay profile for the ππ* state (state b in Figure 2) determined experimentally is overlaid with calculated dynamics. It is assumed that the empirical decay profile represents the sum of two wavepacket trajectories, where 80% of the population relaxes by way of trajectory 1. Comparing these measured and calculated quantities defines the region of parameter space relevant to α-TP. Parameters are given in Table 1. 13959

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

Figure 5. Calculations reveal an interplay associated with the energy gap between levels b and c, ϕbc = ωbc,1/2πc = ωbc,2/2πc, and the time scale of the bath Λ−1 bb,j. Population of state b, Pbj(τ), computed for trajectories 1 (black) and 2 (red) are compared for energy gaps, which are (a) less than and (b) greater than the reorganization energy λbb/hc = 1500 cm−1. (c) Population differences Pb2(τ) − Pb1(τ) are plotted for various values of ϕbc. All other parameters are given in Table 1.

Figure 6. (a) Motion of the photoinduced wavepacket near the point of intersection in the bath coordinate has a strong influence on the population decay profile. The energy gap corresponding to the peak of the wavepacket jδEbβ,cξ(τ) is plotted for trajectories (b) j = 1 and (c) j = 2. The shaded area encompasses energy gaps ranging over ±1 standard deviations in the asymptotic width of the wavepacket. These plots represent the special case ωβξ = 0. The dynamic energy gap jδEbβ,cξ(τ) simply shifts in integer-multiples of ℏωp for other vibronic internal conversion channels.

with the position of the wavepacket because the two surfaces possess the same curvature. Overall, the system dissipates an amount of energy roughly equal to 4λbb,j during the relaxation process. To make the physics more transparent, we have obtained the following expression for the energy gap corresponding to the peak of the wavepacket

predicts inverted, activationless, and normal kinetic regimes similar to those found in electron-transfer processes.44 Whereas energy barriers must be surmounted in the bath dimension, the quantized promoting mode opens IC channels involving energy differences that are potentially much larger than the free energy gap between states b and c, thereby facilitating IC when ℏωbc ≫ kBT. The displacement in the promoting mode controls the amount of vibrational wave function overlap between the initial and final states. Notably, the vibronic progression in IC possesses a different shape than that associated with electron transfer because the perturbative part of the Hamiltonian is an operator in the space of the promoting mode. Figure 5 examines energy detuning effects from the perspective of the bath coordinate by keeping parameters of the promoting mode fixed. The difference in populations for trajectories 1 and 2, Pb2(τ) − Pb1(τ), is computed over a wide range of frequencies, ωbc,1 = ωbc,2. Of course, IC is generally slower for trajectory 2 because of its smaller nonadiabatic coupling (note that Velbc,1 > Velbc,2 in Table 1). However, the sensitivity of the population difference to the value of ωbc suggests that the time scale of the bath Λ−1 bb,j also plays an essential role. The calculations show that the population difference increases as ωbc decreases. Moreover, the population difference is most sensitive to the energy gap between states b and c when it is comparable to the reorganization energy λbb. The interplay between ωbc,j and Λ−1 bb,j can be interpreted by considering the free-energy surfaces associated with states b and c plotted with respect to the bath coordinate (Figure 6a). The wavepacket relaxes toward the equilibrium configuration in state b and may pass over the point of intersection between states b and c (depending on the values of ωbc,j, Λ−1 bb,j, and λbb,j). Notably, these diabatic surfaces intersect because the coupling is perturbative (i.e., this is not a conical intersection). With respect to the bath coordinate, the energy gap changes linearly

j

δEbβ , cξ(τ ) = ℏωbc , j + ℏωβξ + 4λbb[exp( −Λbb , jτ ) − 1] (26)

In the Supporting Information, δEbβ,cξ(τ) is derived by integrating the rate kernel in the “short time approximation”, a procedure that has been discussed previously in the context of four-wave mixing spectroscopies.42,62,63 Upon excitation, the wavepacket moves toward the region of coordinate space where the free-energy gap and reorganization energy are equal (i.e., the point of degeneracy in the space of the bath). To illustrate this point, Figure 6 plots jδEb0,c0(τ) for both trajectories; treatment of other vibronic relaxation channels, where ωβξ ≠ 0, will simply shift jδEbβ,cξ(τ) by integer multiples of ℏωp. The dynamic transition probability maximizes when the wavepacket traverses the point of intersection between the free energy surfaces of states b and c because the rate kernel Wcb,j(τ,t3) oscillates (in t3) with the lowest frequency in this region of space. The key distinction between the two trajectories is that a significant portion of the wavepacket overlaps with the point of intersection at τ < 500 fs for trajectory 1, whereas the wavepacket slides past the point of intersection before τ = 200 fs for trajectory 2. Departure from the point of intersection in trajectory 2 gives rise to an exponential-like “tail” in the population decay profile. The effect of such wavepacket dynamics on the population decay profile is fairly subtle in α-TP because the extremely large dimensionless displacement produces good overlap between j

13960

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

Figure 7. Calculations show that the population-transfer dynamics slow down as the amount of correlation between levels b and c increases. The population of state b, Pbj(τ), is computed for trajectories (a) 1 and (b) 2 over a range of the correlation parameter η. (c) Population differences Pb2(τ) − Pb1(τ) are plotted for various values of η. All other parameters are given in Table 1.

Figure 8. Calculations show that the population-transfer dynamics first speed up then slow down as the dimensionless displacement for the promoting mode, d, increases. Population of state b, Pbj(τ), computed for trajectories (a) 1 and (b) 2 over a range of displacements d. (c) Population differences Pb2(τ) − Pb1(τ) are plotted for various values of d. All other parameters are given in Table 1.

numerous pairs of initial and final vibronic levels (see Figure S1 in Supporting Information). The calculations presented in Figure 7 show that anticorrelation in levels b and c enhances the rate in addition to producing a Gaussian-like population decay profile. Population dynamics Pbj(τ) are generated with correlation parameters, η, ranging from η = −1 to η = 0.4. The onset of a Gaussian-like shape at early times becomes more pronounced as η decreases. The reason for this is that the transition becomes more sensitive to the position of the wavepacket in the bath coordinate in the anticorrelated regime (η < 0). To illustrate this point, eqs C3 and C4 show that amplitudes of the line shape functions that depend on both t2 and t3 increase as η decreases because they scale as (1 − η). In contrast, positive correlation (η > 0) suppresses sensitivity to overdamped wavepacket motion in τ and reduces tolerance for energy detuning between states b and c. Such energy detuning effects are understood by considering that the coherence between states b and c dephases more rapidly (in t3) when η decreases. Increasing the dephasing rate reduces the number of oscillations that occur within the rate kernel Wcb,j(τ,t3) before the amplitude fully decays. Thus, anticorrelation enhances the transition rate by inflating the amount of thermal noise; transitions between nondegenerate states are promoted by this additional bandwidth in the low-friction regime.64 Notably, this argument applies when the wavepacket is not located near the point of degeneracy in the bath coordinate (see Figure 6). At the geometry where states b and c are degenerate, population transfer speeds up when the dephasing rate decreases (i.e., the same is true for an activationless electron-transfer process at equilibrium).44,64 As discussed above, the extremely large dimensionless displacement in the promoting mode in α-TP (d = 2) enables rapid transitions between numerous pairs of vibronic levels by enhancing vibrational overlap intergrals. Beginning with an

extremely small value of d = 0.1, Figure 8 shows that the population-transfer dynamics speed up as d increases. For both trajectories, the rate of IC maximizes near d = 1.3 then slows down as the value of d increases further. This turnover in the rate of population flow represents a shift in the vibronic progression toward larger differences in the number of quanta associated with the initial and final vibronic levels. With respect to energy detuning, the most efficient relaxation channels available to trajectories 1 and 2 involve differences of a small number of vibrational quanta because the frequencies ωbc,1 and ωbc,2 correspond to roughly 2 and 1 vibrational quanta in the 1450 cm−1 promoting mode, respectively. Thus, extremely small displacements (near d = 0.1) induce slow population flow because transitions that annihilate single vibrational quanta dominate the IC dynamics. Increasing d produces appreciable wave function overlap between initial and final vibronic states that differ by more than one vibrational quantum number. However, matrix elements between pairs of vibronic levels involving differences of a few vibrational quanta begin to decrease at d = 1.3, where the vibronic progression favors transitions between states with energy differences that are much larger than ℏωp. To further illustrate this point, examples of vibronic progressions are computed at second-order in Figure S1 of the Supporting Information. In summary, the calculations presented in this section suggest three basic physical insights. First, the population decay profile transitions from a Gaussian-like shape to an exponential shape when the wavepacket passes through the region of space corresponding to the intersection between the free energy surfaces for states b and c (see Figure 6). Second, sensitivity to wavepacket motion requires either weakly correlated or anticorrelated fluctuations in the energies of states b and c; the best agreement with the empirical profile is found for fully anticorrelated fluctuations (see Figure 7). It is not presently clear how such anticorrelation should be interpreted in terms of 13961

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

Figure 9. The sum of wavepackets ξ(τ,qud,qod) is displayed with respect to the dimensionless coordinates qud and qod. The amplitudes of the wavepackets reflect the amount of population in state b, which corresponds to the ππ* excited state in α-TP. Trajectories 1 and 2 are associated with motion toward negative and positive values of the collective bath coordinate qod, respectively. The amount of energy transferred from the system to the bath is linear in qod, where unit displacement corresponds to 2λbb (3000 cm−1). Fairly small amplitude oscillations occur with respect to qud (i.e., the promoting mode) because the 23 fs period of the vibration is comparable to the 20 fs fwhm width of the electric-field envelope. ∞

the electron configurations that contribute to states b and c in α-TP. In general, the fluctuation amplitude in ωcb,j will increase (i.e., amount of correlation decreases) if states b and c involve dissimilar electron configurations. For related reasons, coherences between Frenkel excitons dephase more rapidly when they involve different molecular sites (i.e., suppression of exchange narrowing).45,65 Finally, the calculations in Figure 8 predict a turnover in the rate with respect to the displacement in the promoting mode. On the basis of this result, we conclude that the extremely large displacement in the promoting mode in α-TP (d = 2) suppresses sensitivity to ℏωcb,j by enabling transitions between a large number of initial and final vibronic states.

fud, j (t1 , t 2 , qud) =

(28)

The wavepacket corresponding to the overdamped collective bath coordinate is computed with fod, j (t1 , t 2 , qod) =

1 (4πλbbkBT )1/2

⎧ [q − u (t , t )]2 ⎫ ⎪ ⎪ od, j 1 2 exp⎨ − od − iωbat1 − gbb , j(t1)⎬ ⎪ ⎪ 4λbbkBT ⎭ ⎩

(29)

where uod,j(t1,t2) can be expressed in terms of quantities given above as

IV. WAVEPACKET REPRESENTATION OF DYNAMICS This section suggests a wavepacket representation that can be used to visualize the relationship between population flow and energy dissipation. The kinetic model outlined above treats one overdamped harmonic mode and one underdamped harmonic mode associated with each of the two trajectories. Expressions for both classes of photoinduced wavepackets are known within the approximations made in the present work.42 The underdamped wavepacket, which corresponds to the promoting mode, is given by42

uod, j(t1 , t 2) = −

i2kBTλbb exp( −Λbb , jt 2)[1 − exp(−Λbb , jt1)] ℏΛbb , j

− λbb + λbb exp( −Λbb , jt 2)[1 + exp(−Λbb , jt1)]

(30)

Effects related to the finite bandwidth of the laser pulse are incorporated by convoluting f ud,j and fod,j with the electric field envelopes, i.e.,

1 hud, j(t1 , t 2 , qud , pud ) = 2 2 1/2 ⟩⟨qud ⟩) 2π (⟨pud

ψud, j(τ , qud) = Re

∫0

∞

dt1

∫0

∞

dt 2 fud, j (t1 , t 2 , qud)

E(τ − t 2 − t1) E*(τ − t 2)

⎧ [qud − u ud(t1 , t 2)]2 [pud − z ud(t1 , t 2)]2 ⎪ − − exp⎨ 2 2 ⎪ ⟩ ⟩ 2⟨qud 2⟨pud ⎩ −iωbat1 − gbb , j(t1)}

∫−∞ hud,j(t1, t2 , qud , pud ) dpud

(31)

and ψod, j(τ , qod) = Re

(27)

where qud and pud are the position and momentum of the wavepacket, respectively, ⟨q2ud⟩ and ⟨p2ud⟩ are the variances; the time intervals t1 and t2 are defined in Figure 3. For brevity, the quantities uud(t1,t2) and zud(t1,t2) are given in the Supporting Information. We integrate over momentum for the purpose of visualizing the dynamics in the space of qud

∫0

∞

dt1

∫0

E(τ − t 2 − t1) E*(τ − t 2)

∞

dt 2 fod, j (t1 , t 2 , qod) (32)

The quantities ψud,j(τ,qud) and ψod,j(τ,qod) can be combined to produce two-dimensional wavepackets associated with each trajectory. The information contained in the model can be intuitively visualized using 13962

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

⎡4 ξ(τ , qud , qod) = ⎢ Pb1(τ ) ψod1(τ , −qod) ⎣5 ψud1(τ , qud) +

⎤ 1 Pb2(τ ) ψod2(τ , qod) ψud2(τ , qud)⎥ ⎦ 5

descriptions (i.e., at the level of Fermi’s golden rule). These effects may give rise to Gaussian-like population decay profiles under certain conditions. Our calculations suggest that the rate at which the wavepacket approaches (and departs from) the point of intersection in the bath coordinate (see Figure 6) is the primary factor governing the shape of the population decay profile. Other parameters such as the displacement in the promoting mode and correlation between excited states affect the sensitivity of IC dynamics to the wavepacket’s position but do not alter the basic physical picture. To convey clear physical insights, the representation of the wavepackets presented in Figure 9 reveals correlations between the IC dynamics and the amount of energy transferred from the system to the bath. The fourth-order model presented here will be useful for ultrafast spectroscopic investigations of IC dynamics in molecules, when traditional second-order rate theories cannot be applied because of time-coincident electronic and nuclear relaxation processes. Of course, as in any phenomenological model, choices made regarding particular parameters and/or the inclusion of multiple trajectories should be motivated by the system under consideration and the data available. Here, we have constrained parameters using a linear absorbance spectrum and transient grating signals. Additional spectroscopic measurements can be used to add detail to the model. In particular, resonance Raman experiments can provide Franck− Condon factors for intramolecular modes that are hidden underneath the broad absorbance line shape.41,66 Explicit treatment of additional modes will enable rigorous investigations of the dissipative processes that accompany IC.24 Vibrational modes that are not associated with large nonadiabatic couplings can be incorporated in the bath as underdamped Brownian oscillator coordinates. Parameterization of nonadiabatic couplings in systems with more than one dominant promoting mode will be a challenge. In these cases, ab initio calculations may be needed to evaluate matrix elements involving derivatives of the electronic wave functions.

(33)

where the coefficients 4/5 and 1/5 represent the relative weights of the two trajectories determined in Figure 4b. ξ(τ,qud,qod) is motivated by the display of correlations between the geometry (the wavepacket’s position) and the dynamics of population flow (the wavepacket’s amplitude). In principle, the dynamics of the wavepacket involve three dimensions; the promoting mode and collective bath coordinates associated with each of the two trajectories. However, in ξ(τ,qud,qod), the wavepackets are represented in a two-dimensional space to facilitate visualization; wavepackets associated with trajectories 1 and 2 move toward negative and positive values of the collective bath coordinate qod, respectively. Notably, this distinction in the direction of motion is not made for the promoting mode qud because this coordinate represents specific bond length and angle changes within the molecule (i.e., the direction is not arbitrary). Moreover, departure of the wavepacket in separate directions along qod is physically motivated by the possibility of parallel relaxation pathways in α-TP.57 Figure 9 presents wavepacket dynamics computed using eq 33 and the parameters in Table 1. Wavepackets associated with both trajectories are initially located near the Franck−Condon geometry at τ = 0 fs (qod = qud = 0). Departure of the wavepackets in separate directions manifests as an asymmetry (with respect to qod) near τ = 80 fs. Resolution in the peaks of the two wavepackets occurs near 160 fs, where a significant reduction in amplitude is also observed because of population flow from state b to state c. The wavepackets exhibit oscillations with small amplitudes in the underdamped dimension qud because the 20 fs laser pulse duration is comparable to the 23 fs period of the vibration (see movie in Supporting Information). The amount of energy transferred from the system to the bath is linear in qod, where unit displacement corresponds to 2λbb (3000 cm−1). The wavepacket moves slightly faster along trajectory 2 (Λ−1 bb,2 = 200 fs) and is able to transfer more than 3000 cm−1 to the bath before state b is fully depopulated. In contrast, slower motion in trajectory 1 (Λ−1 bb,1 = 500 fs) results in less energy dissipation. As indicated in Figure 6, the distance traversed by the wavepacket governs its “footprint” at the point of intersection between the diabatic surfaces in the space of the bath coordinate, where the dynamic transition probability is greatest. For example, with ωβξ = 0 (see eq 26), geometries corresponding to these energy degeneracies are reached when amounts of energy equal to 3000 and 1500 cm−1 have been transferred from the system to bath for trajectories 1 and 2, respectively.

■

APPENDIX A

Perturbative Non-Adiabatic Coupling

This appendix defines the perturbative interaction that induces non-radiative transitions between electronic states in the present model. The focus here is mainly on notation because background on IC transitions can be found elsewhere (see ref 24 and Chapter 9 of ref 66). We begin with a brief discussion of the perturbation’s origin to make clear the particular form of the matrix elements employed above. The electronic and nuclear wavefunctions, ψ(r;R) and χ(R), are separated based on their relative time scales of motion; ψ(r;R) is a function of the electronic coordinates r with parametric dependence on the nuclear coordinates R. Terms in the Schrödinger equation that originate in the nuclear kinetic energy operator prevent the attainment of product functions ψ(r;R) χ(R). For example, the action of the nuclear kinetic energy operator on such a product function yields66

V. CONCLUDING REMARKS Several recent experimental studies in our laboratory have examined ultrafast IC processes in DNA components and ringopening systems.51−56 In order to establish a framework for thinking about these observations, we have developed a fourthorder description of photoinduced IC processes and applied it to a model system resembling α-TP in this work. The model captures the interplay between the electronic populationtransfer dynamics and relaxation along Franck−Condon active coordinates, which is absent in traditional second-order

TN ψ (r; R) χ (R) =

∑ Tαψ (r; R) χ(R) α

= −∑ α

∂ 2ψ (r; R) ∂ψ (r; R) ∂χ (R) ℏ ⎧ ⎨χ (R) +2 2μα ⎩ ∂R α ∂R α ∂R α2 2

+ ψ (r; R) 13963

∂ 2χ (R) ⎫ ⎬ ∂R α2 ⎭

(A1)

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

where α is an index for nuclear coordinates. IC transitions can be induced by either of the first two terms on the right-hand side of eq A1, whereas the third term gives rise to a set of vibrational eigenstates when combined with a potential energy surface that is generated by solving the electronic part of the Schrödinger equation. In the present work, the second term in eq A1 is treated as the perturbation and the first term is neglected under the assumption that it possesses a relatively small magnitude.66 For a system with a single promoting mode, the matrix element between electronic states b and c (with vibrational levels β and ξ) is then written as Vbβ , cξ = −

ℏ2 ∂ ∂ ⟨b| |c⟩⟨β| |ξ ⟩ μp ∂R p ∂R p

where the time arguments indicate action of the propagator, ⟨β(t)| = ⟨β|exp(iHct/ℏ). Expressions for the dynamic overlap integrals in eq B4 can be found elsewhere.60,61 In Section IIC, we make use of the following unitless quantity: Φδβ (t3) =

■

APPENDIX C

For the terms indicated in Figure 3, the general forms of the line shape functions f1j(t1,t2,t3) and f 2j(t1,t2,t3) are given by59 * (t 2 ) f1j (t1 , t 2 , t3) = gbb , j(t1 + t 2 + t3) + gcc*, j(t3) + gbb ,j

(A2)

− gbc , j(t1 + t 2 + t3) − gbc*, j(t3) + gbc , j(t1 + t 2) * (t3) + g (t1) − g (t1 + t 2) − g * (t 2 + t3) + gbb bb , j bb , j bb , j ,j − gbc*, j(t3) − gbc*, j(t 2) + gbc*, j(t 2 + t3)

* (t1 + t 2) f2j (t1 , t 2 , t3) = gbb , j(t 2 + t3) + gcc*, j(t3) + gbb ,j * (t3) − gbc , j(t 2 + t3) − gbc*, j(t3) + gbc , j(t 2) + gbb ,j

Evaluating Matrix Elements in a Basis of Harmonic Oscillators

* (t1) − g (t 2) − g * (t1 + t 2 + t3) − g * (t3) + gbb ,j bb , j bb , j bc , j

Equations used to evaluate matrix elements between vibronic eigenstates are outlined below. First, the momentum operator is written in terms of creation and annihilation operators for the promoting mode as ℏμp ωp ℏ ∂ =i (Bp† − Bp) i ∂R p 2

− gbc*, j(t1 + t 2) + gbc*, j(t1 + t 2 + t3)

μp ωp(β + 1) 2ℏ

(B1)

|β + 1⟩ +

⎤ |β − 1⟩⎥ ⎥ 2ℏ ⎦

μp ωpβ

(B2)

The vibrational part of the non-adiabatic coupling that arises in the second-order description is 2

∂ ⟨ξ| |β ⟩ ∂R p μp ωp = |− β + 1 ⟨ξ|β + 1⟩ + β ⟨ξ|β − 1⟩|2 (B3) 2ℏ The matrix elements in the fourth-order rate formula are readily evaluated using similar rules. For example, the matrix element in eq 17, which includes a propagator in t3, is given by

* (t3) + g * (t 2) f1j (t1 , t 2 , t3) = gbb , j(t1) + (1 − η)[2gbb ,j bb , j * (t 2 + t3) + g (t1 + t 2 + t3)] − gbb , j(t1 + t 2) −gbb ,j bb , j (C3)

and * (t1) + (1 − η)[2g * (t3) − g (t 2) f2j (t1 , t 2 , t3) = gbb ,j bb , j bb , j * (t1 + t 2) +g (t 2 + t3) − g * (t1 + t 2 + t3)] +gbb ,j bb , j bb , j

⟨δ|∂/∂R p exp(iHct3/ℏ) ∂/∂R p|β⟩ μp ωp = {− (δ + 1)(β + 1) ⟨δ + 1(t3)|β + 1⟩ 2ℏ

(C4)

The form of gbb,j(t) depends on the structure of the bath. Both underdamped intramolecular modes and overdamped collective modes can be included. In α-TP, the symmetric C C stretch is known to dominate the vibronic progression associated with the HOMO-to-LUMO transition. Therefore, this mode is treated explicitly as part of the system in the model

(δ + 1)β ⟨δ + 1(t3)|β − 1⟩

+ δ(β + 1) ⟨δ − 1(t3)|β + 1⟩ − − 1⟩}

(C2)

The number of parameters can be reduced by simple physical arguments to obtain insights into the relaxation processes in αTP. First, because state c is optically dark, the line shape function gcc,j(t) cannot be estimated by inspection of the linear absorbance spectrum. However, in α-TP, it is thought that both the radiative and non-radiative transitions predominantly involve promotion of an electron from the highest occupied molecular orbital (HOMO) to the LUMO orbital. Therefore, we set the functions gbb,j * (t3) and gcc,j * (t3) equal by assuming that these two transitions possess similar spectral densities. Correlation effects are conveniently incorporated by introducing a parameter η that interpolates between the fully anticorrelated (η = −1) and fully correlated (η = 1) regimes.67 By setting gbc,j(t) = ηgbb,j(t), the line shape functions are rewritten as

It follows that the gradient operator transforms the vibronic eigenstate |β⟩ according to

+

(C1)

and

APPENDIX B

⎡ ∂ |β ⟩ = ⎢ − ⎢ ∂R p ⎣

(B5)

Line Shape Functions

The part of the matrix element involving the vibrational quantum states is readily evaluated in a basis of harmonic oscillators, whereas parameterization of the electronic part, Velbc = ⟨b|∂Rp|c⟩, presents a greater challenge. In the present empirical approach, Velbc enters the rate formula as a constant and can therefore be adjusted to reproduce measured IC time scales.

■

2ℏ ⟨δ|∂/∂R p exp(iHct3/ℏ) ∂/∂R p|β⟩ μp ωp

δβ ⟨δ − 1(t3)|β (B4) 13964

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

Excitation Processes in Polyatomic Molecules. Proc. Natl. Acad. Sci. U.S.A. 1969, 63, 31−35. (11) Tully, J. C.; Preston, R. K. Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H+ with D2. J. Chem. Phys. 1971, 55, 562−572. (12) Miller, W. H.; George, T. F. Semiclassical Theory of Electronic Transitions in Low Energy Atomic and Molecular Collisions Involving Several Nuclear Degrees of Freedom. J. Chem. Phys. 1972, 56, 5637− 5652. (13) Freed, K. F. The Theory of Radiationless Processes in Polyatomic Molecules. Top. Curr. Chem. 1972, 31, 105−139. (14) Coalson, R. D. A Spin−Boson Model for Spectroscopy Involving Nonadiabatically Coupled Potential Energy Surfaces. J. Chem. Phys. 1987, 86, 995−1009. (15) Neria, E.; Nitzan, A. Semiclassical Evaluation of Nonadiabatic Rates in Condensed Phases. J. Chem. Phys. 1993, 99, 1109. (16) Stock, G. A Semiclassical Self-Consistent-Field Approach to Dissipative Dynamics. II. Internal Conversion Processes. J. Chem. Phys. 1995, 103, 2888−2902. (17) Donoso, A.; Martens, C. C. Semiclassical Multistate Liouville Dynamics in the Adiabatic Representation. J. Chem. Phys. 2000, 112, 3980−3989. (18) Farrow, D. A.; Qian, W.; Smith, E. R.; Ferro, A. A.; Jonas, D. M. Polarized Pump-Probe Measurements of Electronic Motion via a Conical Intersection. J. Chem. Phys. 2008, 128, 144510. (19) Kosma, K.; Trushin, S. A.; FuB, W.; Schmid, W. E. Cyclohexadiene Ring Opening Observed with 13 fs Resolution: Coherent Oscillations Confirm the Reaction Path. Phys. Chem. Chem. Phys. 2009, 11, 172−181. (20) Khurmi, C.; Berg, M. A. Dispersed Kinetics without Rate Heterogeneity in an Ionic Liquid Measured with Multiple PopulationPeriod Transient Spectroscopy. J. Phys. Chem. Lett. 2010, 1, 161−164. (21) Kotur, M.; Weinacht, T. C.; Zhou, C.; Matsika, S. Strong-Field Molecular Ionization from Multiple Orbitals. Phys. Rev. X 2011, 1, 021010. (22) Ward, C. L.; Elles, C. G. Controlling the Excited-State Reaction Dynamics of a Photochromic Molecular Switch with Sequential TwoPhoton Excitation. J. Phys. Chem. Lett. 2012, 3, 2995−3000. (23) Cresp-Otero, R.; Barbatti, M.; Yu, H.; Evans, N. L.; Ullrich, S. Ultrafast Dynamics of UV-Excited Imidazole. ChemPhysChem 2011, 12, 3365−3375. (24) Fuss, W. Where Does the Energy Go in fs-Relaxation? Chem. Phys. 2013, 425, 96−103. (25) Zhang, Y.; Oliver, T. A. A.; Ashfold, M. N. R.; Bradforth, S. E. Contrasting the Excited State Reaction Pathways of Phenol and ParaMethylthiophenol in the Gas and Liquid Phases. Faraday Discuss. 2012, 157, 141−163. (26) Krebs, N.; Pugliesi, I.; Riedle, E. Pulse Compression of Ultrashort UV Pulses by Self-Phase Modulation in Bulk Material. Appl. Sci. 2013, 3, 153−167. (27) Lee, J.; Challa, J. R.; McCamant, D. W. Pump Power Dependence in Resonance Femtosecond Stimulated Raman Spectroscopy. J. Raman Spectrosc. 2013, 44, 1263−1272. (28) Yarkony, D. R. Conical Intersections: Diabolical and Often Misunderstood. Acc. Chem. Res. 1998, 31, 511−518. (29) Piryatinski, A.; Stepanov, M.; Tretiak, S.; Chernyak, V. Semiclassical Scattering on Conical Intersections. Phys. Rev. Lett. 2005, 95, 223001. (30) Worth, G. A.; Cederbaum, L. S. Beyond Born-Oppenheimer: Molecular Dynamics Through a Conical Intersection. Annu. Rev. Phys. Chem. 2004, 55, 127−158. (31) Jasper, A. W.; Nangia, S.; Zhu, C.; Truhlar, D. G. Non-Born− Oppenheimer Molecular Dynamics. Acc. Chem. Res. 2006, 39, 101− 108. (32) Levine, B. G.; Martínez, T. J. Isomerization Through Conical Intersections. Annu. Rev. Phys. Chem. 2007, 58, 613−634. (33) Matsika, S.; Krause, P. Nonadiabatic Events and Conical Intersections. Annu. Rev. Phys. Chem. 2011, 62, 621−643.

calculations. As in ref 53, all other photoinduced nuclear motions are incorporated using a collective overdamped Brownian oscillator coordinate, which is related to the collective bath coordinates in eq 4 by gab , j(t ) =

1 ℏ2

∫0

t

dτ2

∫0

τ2

dτ1 ⟨Q aj(τ1)Q bj(0)⟩

(C5)

The line shape function is written as ⎛ 2λ k T λ ⎞ gbb , j(t ) = ⎜⎜ 2bb 2B − i bb ⎟⎟[exp(−Λbb , jt ) + Λbb , jt − 1] ℏΛbb , j ⎠ ⎝ ℏ Λbb , j (C6)

where λbb is the reorganization energy for the HOMO-toLUMO excitation and Λbb,j is the rate at which a nonequilibrium state relaxes to equilibrium.

■

ASSOCIATED CONTENT

S Supporting Information *

A derivation of eq 26, auxiliary equations for the underdamped wavepacket in eq 27, an analytical expression for the fourthorder rate constant, discussion of the normalized excitation probability (eq 23), calculations of the internal conversion rate at second-order, and a movie that illustrates the dynamics shown in Figure 9. This material is available free of charge via the Internet at http://pubs.acs.org.

■ ■ ■

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.

ACKNOWLEDGMENTS This work is supported by the National Science Foundation under CHE-0952439. REFERENCES

(1) Frank, H. A.; Cogdell, R. J. Carotenoids in Photosynthesis. Photochem. Photobiol. 1996, 63, 257−264. (2) Gustavsson, T.; Bányász, A.; Lazzarotto, E.; Markovitsi, D.; Scalmani, G.; Frisch, M. J.; Barone, V.; Improta, R. Singlet ExcitedState Behavior of Uracil and Thymine in Aqueous Solution: A Combined Experimental and Computational Study of 11 Uracil Derivatives. J. Am. Chem. Soc. 2006, 128, 607−619. (3) Middleton, C. T.; de La Harpe, K.; Su, C.; Law, Y. K.; CrespoHernández, C. E.; Kohler, B. DNA Excited-State Dynamics: From Single Bases to the Double Helix. Annu. Rev. Phys. Chem. 2009, 60, 217−239. (4) Polívka, T.; Sundström, V. Ultrafast Dynamics of Carotenoid Excited States−From Solution to Natural and Artificial Systems. Chem. Rev. 2004, 104, 2021−2072. (5) Meredith, P.; Riesz, J. Radiative Relaxation Quantum Yields for Synthetic Eumelanin. Photochem. Photobiol. 2004, 79, 211−216. (6) Piletic, I. R.; Matthews, T. E.; Warren, W. S. Probing nearinfrared photorelaxation pathways in eumelanins and pheomelanins. J. Phys. Chem. A 2010, 114, 11483−11491. (7) Arruda, B. C.; Smith, B.; Spears, K. G.; Sension, R. J. Ultrafast ring-opening reactions: A comparison of α-terpinene, α-phellandrene, and 7-dehydrocholesterol with 1,3-cyclohexadiene. Faraday Discuss. 2013, 163, 159−171. (8) Bixon, M.; Jortner, J. Intramolecular Radiationless Transitions. J. Chem. Phys. 1968, 48, 715−726. (9) Chock, D. P.; Jortner, J.; Rice, S. A. Theory of Radiationless Transitions in an Isolated Molecule. J. Chem. Phys. 1968, 49, 610−621. (10) Rhodes, W.; Henry, B. R.; Kasha, M. A Stationary State Approach to Radiationless Transitions: Radiation Bandwidth Effect on 13965

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966

The Journal of Physical Chemistry A

Article

(34) Nelson, T.; Fernandez-Alberti, S.; Chernyak, V.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics Modeling of Photoinduced Dynamics in Conjugated Molecules. J. Phys. Chem. B 2011, 115, 5402−5414. (35) Tapavica, E.; Meyer, A. M.; Furche, F. Unravelling the Details of Vitamin D Photosynthesis by Non-Adiabatic Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2011, 13, 20986−20998. (36) Domcke, W.; Yarkony, D. R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. Annu. Rev. Phys. Chem. 2012, 63, 325−352. (37) Alguire, E.; Subotnik, J. E. Optimal Diabatic States Based on Solvation Parameters. J. Chem. Phys. 2012, 137, 194108. (38) Wittig, C. Geometric Phase and Gauge Connection in Polyatomic Molecules. Phys. Chem. Chem. Phys. 2012, 14, 6409−6432. (39) Shu, Y.; Levine, B. G. Communication: Non-Radiative Recombination via Conical Intersection at a Semiconductor Defect. J. Chem. Phys. 2013, 139, 081102. (40) Myers, A. B. Resonance Raman Intensity Analysis of ExcitedState Dynamics. Acc. Chem. Res. 1997, 30, 519−527. (41) Kelley, A. M. Resonance Raman Intensity Analysis of Vibrational and Solvent Reorganization in Photoinduced Charge Transfer. J. Phys. Chem. A 1999, 103, 6891−6903. (42) Mukamel, S. Principles of Nonlinear Optical Spectroscopy; Oxford University Press: New York, 1995. (43) Nitzan, A. Chemical Dynamics in Condensed Phases; Oxford University Press: Oxford, 2006. (44) Barbara, P. F.; Meyer, T. J.; Ratner, M. A. Contemporary Issues in Electron Transfer Research. J. Phys. Chem. 1996, 100, 13148− 13168. (45) Abramavicius, D.; Palmieri, B.; Voronine, D. V.; Sanda, F.; Mukamel, S. Coherent Multidimensional Optical Spectroscopy of Excitons in Molecular Aggregates; Quasiparticle versus Supermolecule Perspectives. Chem. Rev. 2009, 109, 2350−2408. (46) Sperling, J.; Nemeth, A.; Hauer, J.; Abramavicius, D.; Mukamel, S.; Kauffmann, H. F.; Milota, F. Excitons and Disorder in Molecular Nanotubes: A 2D Electronic Spectroscopy Study and First Comparison to a Microscopic Model. J. Phys. Chem. A 2010, 114, 8179−8189. (47) Womick, J. M.; Moran, A. M. Vibronic Enhancement of Exciton Sizes and Energy Transport in Photosynthetic Complexes. J. Phys. Chem. B 2011, 115, 1347−1356. (48) Chenu, A.; Christensson, N.; Kauffmann, H. F.; Mančal, T. Enhancement of Vibronic and Ground-State Vibrational Coherences in 2D Spectra of Photosynthetic Complexes. Sci. Rep. 2013, 3, 2029. (49) Valkunas, L.; Abramavicius, D.; Mančal, T. Molecular Excitation Dynamics and Relaxation: Quantum Theory and Spectroscopy; WileyVCH: Weinheim, Germany, 2013. (50) Garavelli, M.; Page, C. S.; Celani, P.; Olivucci, M.; Schmid, W. E.; Trushin, S. A.; Fuss, W. Reaction Path of a sub-200 fs Photochemical Electrocyclic Reaction. J. Phys. Chem. A 2001, 105, 4458−4469. (51) West, B. A.; Womick, J. M.; Moran, A. M. Influence of Temperature on Thymine-to-Solvent Vibrational Energy Transfer. J. Chem. Phys. 2011, 135, 114505. (52) West, B. A.; Moran, A. M. Two-Dimensional Electronic Spectroscopy in the Ultraviolet Wavelength Range. J. Phys. Chem. Lett. 2012, 3, 2575−2581. (53) West, B. A.; Molesky, B. P.; Montoni, N. P.; Moran, A. M. Nonlinear Optical Signatures of Ultraviolet Light-Induced Ring Opening in α-Terpinene. New. J. Phys. 2013, 15, 025007. (54) West, B. A.; Giokas, P. G.; Molesky, B. P.; Ross, A. D.; Moran, A. M. Toward Two-Dimensional Photon Echo Spectroscopy with 200 nm Laser Pulses. Opt. Express 2013, 21, 2118−2125. (55) West, B. A.; Molesky, B. P.; Giokas, P. G. Uncovering Molecular Relaxation Processes with Nonlinear Spectroscopies in the Deep UV. Chem. Phys. 2013, 423, 92−104. (56) West, B. A.; Womick, J. M.; Moran, A. M. Interplay between Vibrational Energy Transfer and Excited State Deactivation in DNA Components. J. Phys. Chem. A 2013, 117, 5865−5874.

(57) Tang, K.-C.; Rury, A.; Orozco, M. B.; Egendorf, J.; Spears, K. G.; Sension, R. J. Ultrafast Electrocyclic Ring Opening of 7-Dehydrocholesterol in Solution: The Influence of Solvent on Excited State Dynamics. J. Chem. Phys. 2011, 134, 104503. (58) Arruda, B. C.; Peng, J.; Smith, B.; Spears, K. G.; Sension, R. J. Photochemical Ring-Opening and Ground State Relaxation in αTerpinene with Comparison to Provitamin D3. J. Phys. Chem. B 2013, 117, 4696−4704. (59) Abramavicius, D.; Mukamel, S. Many-Body Approaches for Simulating Coherent Nonlinear Spectroscopies of Electronic and Vibrational Excitons. Chem. Rev. 2004, 104, 2073−2098. (60) Myers, A. B.; Mathies, R. A.; Tannor, D. J.; Heller, E. J. Excited State Geometry Changes from Preresonance Raman Intensities: Isoprene and Dexatriene. J. Chem. Phys. 1982, 77, 3857−3866. (61) Yan, Y.-J.; Mukamel, S. Eigenstate-Free, Green Function, Calculation of Molecular Absorption and Fluorescence Line Shapes. J. Chem. Phys. 1986, 85, 5908−5923. (62) Zhang, Y.; Berg, M. A. Ultrafast Dichroism Spectroscopy of Anthracene in Solution. II. Solvation Dynamics from a OneDimensional Experiment. J. Chem. Phys. 2001, 115, 4223−4230. (63) Kwac, K.; Cho, M. Two-Color Pump−Probe Spectroscopies of Two- and Three-Level Systems: 2-Dimensional Line Shapes and Solvation Dynamics. J. Phys. Chem. A 2003, 107, 5903−5912. (64) Lee, M.; Holtom, G. R.; Hochstrasser, R. M. Observation of the Kramers Turnover Region in the Isomerism of Trans-Stilbene in Fluid Ethane. Chem. Phys. Lett. 1985, 118, 359−363. (65) Lee, H.; Cheng, Y.-C.; Fleming, G. R. Coherence Dynamics in Photosynthesis: Protein Protection of Excitonic Coherence. Science 2007, 316, 1462−1465. (66) Kelley, A. M. Condensed-Phase Molecular Spectroscopy and Photophysics; John Wiley & Sons: Hoboken, NJ, 2012. (67) Venkatramani, R.; Mukamel, S. Correlated Line Broadening in Multidimensional Vibrational Spectroscopy. J. Chem. Phys. 2002, 117, 11089−11101.

13966

dx.doi.org/10.1021/jp4079162 | J. Phys. Chem. A 2013, 117, 13954−13966