Theoretical Analysis of the Relaxation Dynamics in Perylene Bisimide

Article pubs.acs.org/JPCA

Theoretical Analysis of the Relaxation Dynamics in Perylene Bisimide Dimers Excited by Femtosecond Laser Pulses Alexander Schubert,† Mirjam Falge,† Martin Kess,† Volker Settels,† Stefan Lochbrunner,§ Walter T. Strunz,∥ Frank Würthner,‡ Bernd Engels,*,† and Volker Engel*,† †

Institut für Physikalische und Theoretische Chemie, Universität Würzburg, Hubland Campus Nord, Emil-Fischer-Strasse 42, 97074 Würzburg, Germany ‡ Institut für Organische Chemie, Universität Würzburg, Am Hubland, 97074 Würzburg, Germany § Institut für Physik, Universität Rostock, Universitätsplatz 3, 18055 Rostock, Germany ∥ Institut für Theoretische Physik, TU Dresden, 01062 Dresden, Germany ABSTRACT: We present a model for the relaxation dynamics in perylene bisimide dimers, which is based on ab initio electronic structure and quantum dynamics calculations including effects of dissipation. The excited-state dynamics proceeds via a mixing of electronic states of local Frenkel and charge-transfer characters, which becomes effective upon a small distortion of the dimer geometry. In this way, it is possible to explain the fast depopulation of the photoexcited state, which we characterize by femtosecond transient absorption measurements. The combined theoretical and experimental analysis hints at a trapping mechanism, which involves nonadiabatic and dissipative dynamics in an excited-state vibronic manifold and provides an atomistic picture that might prove valuable for future design of photovoltaic materials.

I. INTRODUCTION The fast and efficient transport of energy absorbed from light is of central importance in photosynthesis as well as photovoltaic processes. In this context, molecular π-aggregates play a fundamental role as, for example, coupled chromophores in light-harvesting complexes or dye molecules in organic semiconductor devices.1−3 What restricts the energy transfer are de-excitation mechanisms leading to a trapping of energy on local sites, which, in turn, is the reason for the rather limited energy transfer in most organic materials.4 While this is wellknown as a bottleneck for organic photovoltaics, the underlying processes are essentially not understood, and it is of great interest to obtain an atomistic understanding of such quenching processes. The molecular dimer, consisting of two monomers (M), is an important subunit of aggregates. It has been shown that the optical properties of perylene bisimide aggregates can excellently be explained using dimer models.5−7 Also, after excitation, the energy tends to localize very fast in these subunits.8 Therefore, in what follows, we restrict our discussion to elementary processes in dimers, taking place after photo absorption. We recently investigated experimentally and theoretically7 aggregates of 3,4,9,10-perylene tetracarboxylic acid bisimide (PBI) to understand the optical properties of such important building blocks for molecular electronics.9,10 The monomer structure and the dimer geometry, as derived from our © 2014 American Chemical Society

quantum chemistry calculations (see section III.A), are illustrated in Figure 1. The picture that evolved from our studies can be summarized as follows. According to the dimer transition-dipole geometry with monomer moments at a

Figure 1. Upper part: Structure of the PBI-monomer. For the experimental work, R = 3,4,5-dodecylphenyl substituents were employed. Lower part: Ground-state configuration of the dimer. In the equilibrium geometry, two monomers are arranged at a torsional angle of φ = 30°. Received: December 12, 2013 Revised: January 30, 2014 Published: January 30, 2014 1403

dx.doi.org/10.1021/jp412166a | J. Phys. Chem. A 2014, 118, 1403−1412

The Journal of Physical Chemistry A

Article

relative angle of 30°, the upper of the two local Frenkel states accessible in the absorption band centered around 480 nm6 possesses the higher transition rate from the electronic ground state. Excitation of this band is followed by relaxation to the lowest excited state. Thereby, the two monomers orient and assume a parallel transition-dipole geometry. Radiative decay to the ground state, which is forbidden for the exact parallel configuration, is only possible if nuclear degrees-of-freedom are incorporated, and, as a consequence, a long radiative lifetime (∼32 ns) and a low fluorescence yield as compared to the monomer case are found.6 Because of the flat excited-state torsional potential, the red-shifted emission band is broadened significantly. Up to this point, a consistent picture of the absorption and emission properties could be established. In a recent letter, we reported results of femtosecond transient absorption measurements on PBI.11 It was found that the photoexcited state is depopulated on a 200 fs time scale. A quantum model based on ab initio electronic structure and nuclear wave packet calculations was introduced, which was able to interpret the data. For the present system, it is essential to use quantum chemical methods that provide potential energy curves as accurately as possible, because the underlying mechanisms are initiated by couplings between different electronic states. This forbids the use of, for example, combined classical/semiempirical quantum methods.12 In this Article, we give all details of the model and investigate several aspects being connected to the laser-excitation and dissipation processes. In particular, we replace the formerly used too simple dissipation model by a more accurate approach based on a stochastic Schrödinger equation. This Article is organized as follows. The experiment is briefly described in section II. The electronic structure calculations and the construction of the diabatic Hamiltonian for the nuclear motion are described in section III, which also contains details about the wave packet propagation. The numerical results are given in section IV, which also presents a comparison to experiment and a conclusion.

Figure 2. The time-dependent absorption change (at 520 nm) (○) is decomposed into three components (lower curves) representing an oscillatory, a decay, and a step contribution, as indicated. Also shown is the signal from the pure solvent (upper fine line). Reprinted with permission from ref 11. Copyright 2013 American Chemical Society.

damping time of 190 fs. The combination of the oscillatory and the exponential component results in an almost complete decay of the stimulated emission within the first 200 fs.

III. THEORY AND MODEL III.A. Electronic Structure Calculations. The model we develop to explain the experiments is based on PBI dimer computations. In a local picture, the lowest excitation of a dimer gives rise to two localized Frenkel (F) states |φF1 ⟩, |φF2 ⟩ with configurations M−M*, M*−M and two localized chargetransfer states (CT) |φ3CT⟩, |φ4CT⟩ (M−−M+, M+−M−), respectively. Model Hamiltonians incorporating coupled Frenkel- and CT-states have been used to describe aggregate properties; for recent work, see, for example, refs 15−20. The adiabatic electronic states obtained from quantum chemistry calculations, which look upon the dimer as a supermolecule, are mixtures of these configurations where the mixing coefficients strongly depend on geometry. This leads to four adiabatic excited states |φan⟩ (n = 1−4) with geometry-dependent energies. The supermolecule ansatz describes the interaction between the monomers on a full quantum chemical level, which allows predicting the relative energies of the adiabatic states as a function of chosen coordinates. These considerations show that dimers properly represent all electronic interactions including Frenkel and CT states in particular. Semiempirical approaches that would allow one to include various nuclear degrees of freedom turned out to be too inaccurate for the present example. Even time-dependent density functional theory (TD-DFT), which provides quite accurate results for many organic compounds, is not suited for the calculations because Frenkel and CT states might be involved in the described processes. Such approaches drastically underestimate the energetic position of CT states even for long-range corrected functionals.21,22 Hence, to obtain reliable information about ground and excited states, we use the approximate coupled cluster second-order method using spincomponent-scaling (SCS-CC2) for the computation of the excited states, while the spin-component-scaling Møller− Plesset perturbation theory to second order (SCS-MP2) with

II. EXPERIMENT The experiment has been described in detail before.11,13 Briefly, ultrafast pump−probe absorption measurements are performed on an aggregating PBI-derivative14 in methylcyclohexane. The sample is excited at 480 nm (which is at the respective absorption maximum) with 30 fs long pump pulses generated by a noncollinear optical parametric amplifier. The weak probe pulses (40 fs) are centered at 520 nm. The measured transientabsorption traces are shown in Figure 2 (upper panel). The figure also contains the signal from the pure solvent. We decompose the signal into three contributions: (a) one component is a step function that starts with the excitation and afterward remains constant; (b) a (negative) exponentially decaying component from which a time constant of 215 fs can be extracted; and (c) an additional function that exhibits oscillations with a period of 381 fs. We interpret these findings as follows. The step function reflects the bleach of the electronic ground state, which is switched on with the excitation and does not change thereafter. The exponential decay, which has a negative amplitude, resembles stimulated emission. It sets in with the excitation and decays with a time constant of 215 fs. We note that this decay is not seen in the case of nonaggregated PBIs.13 The oscillatory contribution has a frequency of 86 cm−1 and a 1404



Article

the resolution of identity approximation23−25 is employed for the ground state. SCS-CC2 typically provides excitation energies with an accuracy around 0.1 eV.26 Similar error bars can also be expected for the energy differences between the predominantly neutral and CT states because error compensation will be small due to the different electronic structure of the states. For the SCS-CC2 computations, we employ the SVP (split valence plus polarization) basis set.24,27,28 The quantum chemical ab initio calculations are carried out using the Turbomole 6.1 program package. The analysis of the character of the excited states is performed by a method that was proposed by Liu.29 This approach is based on a projection of the delocalized dimer MOs to the corresponding Löwdin orthogonalized monomer ones. The projection coefficients of the configurations can be used to analyze the states by their amount of CT versus Frenkel character. Because the localization of an exciton is unclear, the use of trimers or even larger aggregates as basic unit would be desirable, but the SCS-CC approach is so expensive that computations of larger aggregates (trimer, tetramers, etc.) are impossible. Furthermore, the use of aggregates larger than dimers seems to be unnecessary because recent model calculations proved the dimer ansatz to be able to describe the main effects in the aggregate spectra. They show that the extension to trimers and tetramers only leads to subtle changes.30 This further justifies the use of dimers instead of larger aggregates. The experiments indicate that the depopulation of the laser excited upper Frenkel state occurs within a few hundred femtoseconds.11 The question is which nuclear motion can efficiently drive such a fast nonadiabatic transition from the laser excited state to lower electronic states. A comprehensive answer could be obtained if one computes the gradients of the laser excited state starting at the Franck−Condon region and relaxing the system to its next minimum. A fast depopulation could then, for example, happen if along a given path crossings or avoided crossings with the lower states take place. Even strong nonadiabatic couplings to a lower state would already be sufficient for a fast depopulation. Because such computations are far too expensive (see above), we considered various possible nuclear motions and checked if they could lead to the depopulation. For all considerations, we started from the Franck−Condon region, that is, the ground-state equilibrium geometry of the PBI dimer. In the equilibrium geometry, the monomers are arranged with a center of mass separation of 3.3 Å and a torsional angle of φ = 30°7 (see Figure 1). The internal geometries of the monomers are determined with the dispersion corrected density functional theory (DFT-D2) method31,32 using the BLYP33−35 functional in combination with the TZVP basis set at the non-hydrogen atoms and the TZV36 basis set at the hydrogen atoms. For the subsequent computations, these internal coordinates were kept fixed. Because of the large moment of inertia of the monomer− monomer motion, the respective dynamics along the torsional coordinate is slow so that the fast decay seen in our experiment cannot be explained by such intermolecular geometry changes. This argument also holds for intermolecular motions, which shift both monomers with respect to each other. In the equilibrium geometry of the PBI dimer, both monomers form a kind of H-aggregate. If due to such a shift motion the dimer starts to change into J-orientations, both Frenkel states would cross, giving rise to the fast depopulation. However, such motions also have to be excluded because they are too slow.

Photoinduced vibrations of the single monomers, on the other hand, might lead to population transfer between the two Frenkel states if a distortion along the respective nuclear coordinates enhances the coupling strength and/or the energy difference between them. One intuitive depopulation mechanism can be envisioned by the following argument: In the case of a perfect localization of the excitation onto one monomer, the (Davydov) splitting disappears and the Frenkel states become degenerated. If, afterward, localization is lost, the splitting reappears and the exciton can relax into the lower Frenkel state. This motion can be realized by relaxing one monomer to a structure M* and keeping the other fixed at its ground-state structure M. To do so, the PBI monomer structure is optimized for the optically bright 11B2u state using the SCS-MP2/TZV(P) level of theory. The dimer is then arranged at the intermonomer configuration as described above. A reaction coordinate is introduced, which linearly relates the initial ground state M−M to the final M−M* configuration; see below. For this motion, the ground state is calculated by SCS-MP2/SVP23−25 and excitation energies by SCS-CC2/SVP.23,37 We find that all states are stabilized during the motion.38 The effect is moderate for the ground state, whereas it is larger in the excited states. As an important fact, one finds that along the potential energy curve the Davydov splitting between the Frenkel states does not decrease but increases slightly. This indicates that the excitation does not localize during the proposed motion. Interestingly, the lowest CT state shows the largest stabilization. Consequently, a motion stabilizing this CT state could be more promising as the one discussed above to depopulate the bright Frenkel state. At the vertical excitation geometry, the upper predominantly Frenkel-state and the lower predominantly CT-state differ only by about 0.04 eV. Taking into account that our geometry is optimized for neutral monomers, it might be that a relaxation to a corresponding CT-geometry induces a crossing of both states. To clarify this situation, we characterize the influence of a deformational motion on the relative energies. Because a full optimization is not feasible, we compute the ground-state geometry of a monomeric cation and anion and assemble them to the corresponding dimer structure. Because of problems with spin contamination in HF calculations, the monomer structures of the anion and cation are optimized using B3LYP-D/ TZV(P).33,34,39,40 We define a reaction coordinate q, which linearly connects the initial ground state (M−M) configuration with coordinates Rα(q = 0) to the final (M+−M−) configuration with coordinates Rα(q = 1). With the numbering of the coordinates of the N atoms as {Rα, α = 1 − 3N}, they are parametrized as R α(q) = R α(q = 0) + q[R α(q = 1) − R α(q = 0)]

(1)

The molecular deformation, which takes place in going from q = 0 to q = 1, is, in fact, only minor. In the cation, the imide groups are slightly stretched, while they are compressed in the anion. The maximal deviation in the atomic positions is found to be on the order of 0.01 Å.13 Despite the small structural changes, the effect on the excited-state energies is quite large. The four adiabatic potential curves Van(q) are displayed in Figure 3, upper panel. The potential of the electronic ground state has a shallow minimum at q = 0 (where it assumes the value of zero) and is not included in the figure. In our experiment, femtosecondexcitation from the Franck−Condon region (q = 0) populates the second excited state. 1405



Article

potentials coincide with the diabatic ones. We use the parametrization: 3

V nd(q) =

∑ ai ,nqi

(2)

i=0

where the coefficients are collected in Table 1. From the adiabatic potentials and CT-characters (Figure 3), it is Table 1. Parameters for the Diabatic Potentials and Coupling Elements n an,0 [eV] an,1 [eV] an,2 [eV] an,3 [meV] Jn [eV] βn qc,n

Figure 3. Upper panel: Adiabatic potential energy curves Van(q) as a function of the reaction coordinate q. The lower panel shows the quantum chemically determined CT-characters for the different electronic states (straight lines). Also shown are the CT-characters derived from the diabatization procedure (dashed lines).

The potential curves show (at least) two avoided crossings. The first one occurs at a value of q = 0 (between Va3 and Va4), and the other one is expected around q = 0.37 (between Va2 and Va3). The fact that at regions of avoided crossings the electronic character of the adiabatic states changes fast as a function of the nuclear coordinates is reflected in rapid changes of molecular properties. This is manifested, for example, in the oscillator strengths and also the CT-characters, which are shown in Figure 3, lower panel. It is seen that the laser excited state |φa2⟩ at the (M−M)* configuration (q = 0) is of predominantly Frenkel-character. The next higher state (|φa3⟩), however, that is close in energy, has almost exclusively CT-character. These characters change rapidly in passing the avoided crossing at a value of q = 0.37 with the result that, for larger values of the reaction coordinate, the nature of the two states interchanges. For larger values of q, the states |φa1⟩ and |φa2⟩ also start to interchange their character, which hints at a strong and delocalized coupling between these states. Finally, the upper state |φa4⟩ is separated from the other ones and only plays a minor role in the relaxation dynamics, see section IV, whereas it is of great significance in the following diabatization scheme. The question arises if further nuclear motions can contribute to the depopulation, that is, if a one-dimensional model is sufficient. This is difficult to answer because the necessary quantum chemical calculations are prohibitively expensive. Nevertheless, the coupling along the motion described above seems to be sufficiently strong to explain the fast depopulation (see below). The question remains if the distortion identified above has to be coupled to important intermolecular motions, for example, the torsional motion around φ, which was found to be responsible for the strong red shift of the emission spectra of PBI aggregates. In our opinion, this is not necessary. The coupling between both is expected to be very weak because the associated dynamics proceed on very different time scales. Hence, we remain with a one-dimensional model. III.B. Diabatization Scheme. For an efficient treatment of the quantum dynamics, we set up a diabatic model for the manifolds of excited electronic states. The ground state is far separated energetically and for that reason needs not to be incorporated. The diagonal elements of the diabatic potential matrix (Vdn(q),n = 1−4) are numbered according to their energetic order at a value of q = −0.5, for which the adiabatic

1

2

3

4

2.64 −0.193 0.177 −3.20 0.1 0.5 2.15

2.95 −0.158 0.182 −2.10 0.007 20.0 0.37

3.02 −0.0875 0.173 −2.45 0.028 0.1 0.02

3.02 −0.333 0.160 4.66

reasonable to assume a coupling between the pairs of adiabatic states (n, n′) = (2,3), (3,4) and, for larger distances, (1,2), respectively. We translate this into a diabatic picture and take only couplings V14(q), V24(q), and V34(q) into account. Note that the numbering of the adiabatic and diabatic curves (which are identical at q = −0.5) differs in energy at larger positive

Figure 4. Upper panel: Diabatic potentials Vdn(q) as a function of the reaction coordinate q. The eigenenergies in the four states are indicated as horizontal lines in the respective potentials. The lower panel contains the diabatic coupling elements Vdn4.

distances of the reaction coordinate (see Figure 4). The coupling elements are of the form: 2

V nd4(q) = V 4dn(q) = Jn e−βn(q − qc,n)

(3)

where the maxima occur at the crossing points qc,n of the respective diabatic potentials. The coupling strengths Jn are taken as one-half the splitting between the adiabatic potentials at these positions. What is left to be specified are the values of the parameters βn, which determine the extent of the coupling regions. To fix their values, we first impose the condition that the diagonalization of the diabatic potential matrix V̂ d(q) leads to the diagonal adiabatic potential matrix V̂ a(q), that is: 1406



Article

The numbers for βn are taken as collected in Table 1. The adiabatic potential energy curves obtained by diagonalization of the diabatic potential matrix are in perfect agreement with the ab initio results and are not shown here. The CTcharacters obtained from the diabatization procedure are compared to the quantum chemically determined ones (Â a(q)) in Figure 3, lower panel. Although the agreement is not perfect, it is sufficiently accurate to be used in our dynamics calculation. A better agreement could be obtained by a finetuning of the parameters. This, however, seems to be too ambitious regarding the simplicity of our one-dimensional model. At this point, we also have to consider the accuracy of the ab initio data. As stated in section III.A, the accuracy of the potential curves is in the order of 0.1 eV. To illustrate the effect of this uncertainty, we have modified the set of curves Van such that the potential with quantum number n = 2 is shifted down in energy. In Figure 5, we compare the ab initio potentials,

a d −1 V̂ (q) =D̂ (q)V̂ (q)D̂ (q)

⎛ V d(q) 0 0 ⎜ 1 ⎜ 0 V 2d(q) 0 −1 = D̂ (q)⎜ ⎜ 0 V 3d(q) ⎜ 0 ⎜ d d ⎝ V 41(q) V42(q) V 43

d V14 (q)⎞ ⎟ d V 24(q)⎟⎟ D̂ (q) ⎟ d V 34 (q)⎟ ⎟ V 4d(q) ⎠

⎛ V1a(q) 0 0 0 ⎞ ⎜ ⎟ ⎜ 0 V 2a(q) 0 0 ⎟ ⎟ =⎜ ⎜ 0 0 V 3a(q) 0 ⎟ ⎜⎜ ⎟⎟ 0 0 V 4a(q)⎠ ⎝ 0 (4)

where D̂ (q) is a unitary matrix. This procedure, however, is not unique, and we proceed and use a “property-based” approach,41 which takes the character of the electronic states into account. Here, we use the CT-character as an input, which is readily available from our electronic structure calculation (see Figure 3). Within a diabatic representation, the character of the electronic basic functions {|φdn⟩} should only weakly (and ideally not) depend on the reaction coordinate, whereas the adiabatic basis {|φan⟩} (built by solutions of the electronic Schrödinger equation) depends on the coordinate q. We define localized basis functions of pure Frenkel- (|φF1 ⟩,|φF2 ⟩) and CTCT character (|φCT 3 ⟩, |φ4 ⟩). The projector onto the CT-states then is Â = |φ3CT⟩⟨φ3CT| + |φ4CT⟩⟨φ4CT|

Figure 5. Adiabatic potentials as obtained from the ab initio calculations are shown as solid lines. Here, the potential with quantum number n = 2 is shifted downward in energy by an amount of 0.05 eV, which is within the accuracy of the quantum chemical calculation. Also shown are the potentials Van(dia.) derived from the diagonalization of the diabatic potential matrix. To achieve agreement of the two sets of curves, the coupling element Vd24(q) has to be modified; see text.

(5)

having matrix elements in the adiabatic (j = a) and diabatic (j = d) basis: j Â nm = ⟨φnj|Â |φmj⟩

where Va2(q) is shifted by an amount of −0.05 eV (solid lines), to those obtained from the diabatization (dashed lines). It is seen that only minor deviations occur. To achieve this agreement, the diabatization needs a stronger coupling element J24 = 0.0325 eV and, because the coupling becomes much more delocalized, a smaller width parameter β2 = 1. If we take a shift of −0.075 eV (and β2 = 1), which is still within the error bars of the ab initio approach, the coupling becomes even stronger and takes a value of J24 = 0.0445 eV. III.C. Kinetic Energy Operator. The kinetic energy operator along the linearized reaction-path can be derived in a straightforward manner.42,43 Therefore, we first introduce the mass-scaled (atomic masses mα) displacement vectors:

(6)

and the CT-character in the electronic state (n) given in one or the other representation is readily calculated as: j Â nn = |⟨φ3CT|φnj⟩|2 + |⟨φ4CT|φnj⟩|2

(7)

Upon a change of basis, the operator Â transforms like the potential operators: a −1 d Â (q) = D̂ (q)Â D̂ (q)

(8)

For a diabatic basis, which does not depend on the nuclear coordinate, the operator Â d is as well independent of q. Thus, we can vary its constant diagonal matrix elements together with the coupling parameters βn, such that the transformation eq 8 yields the matrix Â anm with diagonal elements, which correspond to the calculated CT-characters. In the calculation, we used the adiabatic CT characters at the geometry corresponding to q = −0.5, where interactions are weak, to set up the constant diagonal diabatic CT-matrix. Neglecting off-diagonal elements, one obtains the following CT matrix: ⎛ 0.16 0 0 0 ⎞ ⎟ ⎜ 0 0.06 0 0 ⎟ d Â = ⎜ ⎜ 0 0 0.88 0 ⎟ ⎟ ⎜ ⎝ 0 0 0 0.86 ⎠

Δ0α =

mα [R α(1) − R α(0)]

(10)

where the index runs from α = 1 to α = 3N, and Rα(0), Rα(1) are the initial and final molecular configurations as defined in section III.A. Note that we treat every degree-of-freedom α separately, while N is the total number of atoms. The squared norm of these displacement vectors is M (0) =

∑ mα(R α(1) − R α(0))2 α

(11)

Along the reaction path, the atomic positions can be parametrized as:

(9) 1407


The Journal of Physical Chemistry A R α(q0) = R α(0) +

Article

Δα0 q mα 0

p2̂ γ 1 1 + ⟨pn̂ ⟩n pn̂ + i⟨qn̂ ⟩n pn̂ − pn2̂ − ⟨pn̂ ⟩n2 2 2 2 2 γ γ {p ̂ − ⟨pn̂ ⟩n }χ (t ) + +i (21) 2 n 4

(12)

xα(q0) = xα(0) + Δα0 q0

(13)

(14)

xα(q0 , ..., q3N − 1) = xα(0) +

+

∑ i=1

(15)

1 M

(i)

(16)

Using the chain rule, we obtain: ∂ = ∂xα ′

∑ i

∂qi ∂ = ∂xα ′ ∂qi

Δiα ′ ∂ M (0) ∂qi

∑ i

(22)

1 M (0)ωn

(23)

p̂ (24)

and the expectation values ⟨q̂n⟩n, ⟨p̂n⟩n are evaluated with the diabatic wave functions ψdn(qn,t) in the various states. We have introduced the function χ(t) = dW(t)/dt, where W(t) is the complex Wiener increment, which fulfills the usual requirements.46 Furthermore, γ is a dissipation constant, which is a parameter in our calculation. We then arrive at the (SplitOperator47) propagation scheme:

∑ Δiα[xα − xα(0)] α

}

M (0)ωn (q ̂ − qe, n)

pn̂ =

By projection, the latter equation can be inverted as (i ≠ 0): qi =

1 ⟨q ̂ ⟩n2 2 n

qn̂ =

3N − 1

Δiα qi

{

Here, we use the definitions:

The (mass-scaled) atomic coordinates are expressed as: Δ0α q0

γ 1 ⟨qn̂ ⟩n qn̂ − i⟨pn̂ ⟩n qn̂ − qn̂ 2 2 2 γ + {q ̂ − ⟨qn̂ ⟩n }χ (t ) 2 n

Bn̂ (q)̂ = Vn(q) + −

α

}

and

where the last equation employs mass scaled coordinates xα(q0) = (mα)1/2Rα(q0), and, for now, the formerly defined reaction coordinate q is named q0. We define (3N − 1) additional orthogonal displacement vectors with components Δiα, so that:

∑ ΔiαΔαj = M (i)δij

{

Â n (p ̂) =

(17)

a

∂2 = ∂xα2′

∑ i,j

−1 ̂ ̂ ̂ ̂ t /2 ψ̲ (p , t + dt ) = e−iA(p̂)dt /2F ̂ e−iB(q)d D̂ (q)̂ e−idtV (q)̂ D̂ (q)̂

Δiα ′Δαj ′

∂2 M (i)M (j) ∂qi ∂qj

̂

e−iB(p̂)dt /2F ̂

(18)

1 2

∑ α

∂2 1 =− 2 2 ∂xα

∑ i

∂2 2 M (i) ∂qi 1

(19)

Thus, T̂ consists of the sum of the operators, which act along the direction of the different displacement vectors. The reaction-path kinetic operator then can be separated:

T̂(q) = −

∂2 2 2M (0) ∂q

e

ψ (p , t )

(25)

where F̂ denotes the Fourier-transform from coordinate to momentum space and F̂ −1 its inverse transform. The propagation conserves the norm to first order, but due to a finite time step dt in the numerical implementation the norm fluctuates slightly. Therefore, the total wave function is normalized after every time step. For recent application of the numerical procedure, see refs 48−50. By calculating Nr realizations of the wave function, the density matrix can be determined as:

Because of the orthogonality of the displacement vectors, the kinetic energy operator takes the form: T̂ = −

−1 −iÂ (p ̂)dt /2

1

ρ̲ (q , q′, t ) = (20)

where we have renamed the reaction coordinate as q = q0. Apparently, the squared norm M(0) given by eq 11 represents an effective mass, which, in our case, takes the numerical value of M(0) = 973 with the unit (Hartree)−1. III.D. Quantum State Diffusion. We use a version of a stochastic Schrödinger equation44 to incorporate system−bath interactions. In more detail, we employ a quantum state diffusion equation in the Itô form45 for harmonic oscillators coupled to a zero temperature bath. It leads to a stochastic Schrödinger equation, which is equivalent to the Lindblad master equation.45 In our case, the wave function has four diabatic components. Each diagonal element Ĥ dn of the diabatic Hamiltonian is, to a very good approximation, harmonic because the cubic components in the expansion eq 2 are small. The oscillator frequencies ωn are derived from the second derivative of the respective potential curve Vdn(q) at the equilibrium distance qe,n. Because the frequencies and equilibrium distances are different in the four states, the propagation involves diagonal matrix operators Â , B̂ with diagonal matrix elements:

1 Nr

Nr

∑ ψ̲ jd(q , t )( ψ̲ jd(q′, t ))* j=1

(26)

The populations in the different states are then obtained as the trace over respective components of the density matrix, that is: Pn(t ) =

1 Nr

Nr

∑ ∫ dq |ψnd,j(q , t )|2 j=1

(27)

We neglect the effect of dissipation on the diabatic coupling (diabatic damping approximation (DDA)). It has been shown that this approximation eventually breaks down51,52 and leads to trapping of population in excited electronic states so that a complete relaxation to the ground state does not occur. Following the analysis of Kleinekathöfer et al., we regard the level structure of our diabatic states. The respective eigenenergies are marked in Figure 4, upper panel. Because the diabatic state |φd3⟩ plays only a minor role in the dynamics, it is sufficient to discuss the validity of the DDA for the (2,4) and (4,1) transitions. From the figure, it can be seen that the lowest eigenenergy of the laser excited state |φd2⟩ lies above the crossing point of the potentials Vd2 and Vd4 and that the levels of the two respective states are not in resonance. Under these 1408



Article

conditions, the DDA was shown to give only minor deviations from the exact results.51 The same applies for the (4,1) transition where the potentials represent a situation of a Marcus inverted region53 and the levels of the two states are nonresonant.52 Thus, to develop an overall picture of the nonadiabatic dissipation process, the DDA, as adopted in this work, offers a reasonable first approach.

IV. NUMERICAL RESULTS We regard the population dynamics employing the quantumstate diffusion method to describe the dissipation. In Figure 6,

Figure 7. Population dynamics obtained for a stronger coupling element Vd24(q).

populated transiently, and its population reaches values up to 20%. Note that asymptotically about 3% of the population is trapped in the excited state |φd4⟩. This hints at errors due to the DDA, which, however, are not severe. We note that the same decay dynamics is obtained if the diabatization uses the shifted potential curve Va2, see Figure 5, which means that, within the accuracy of the ab initio data, the modification of the coupling is reasonable. An analysis of the wave packet motion in the initial state (|φd2⟩) exhibits a vibrational period of Tq = 41 fs.11 This is too fast to be detected with the experimental cross-correlation width of 50 fs. On the other hand, the absorption trace in Figure 2 shows an oscillatory pattern with a period of Tb = 381 fs. This suggests that the signal reflects a quantum beat structure, which stems from the coherent superposition of several vibrational modes. The nature of these modes is unclear. To identify them, it is necessary to explore the multidimensional potential energy surface, which, for reasons explained above, is out of range. From the spectral width of the excitation pulse (which is about ΓE = 0.12 eV), vibrational wave packets can only be prepared if states with a vibrational energy spacing less than the pulse-width are coherently excited. Thus, for example, a coherent superposition of states for the C−C stretching vibration at 0.175 eV, which is prominent in the monomer absorption spectrum, and also influences the dimer spectrum,5,6 cannot be excited. However, there are other monomer modes that fall into the spectral range of our laser pulse and have non-negligible transition strengths.54 Because there are many vibrational levels seen in aggregate spectra,55 there might be also other deformational motions of the dimer with appropriate vibrational frequencies. Finally, the results of our numerical calculations are compared to experiment. Figure 8 shows the decaying signal, which corresponds to stimulated emission and thus represents the depopulation of the excited state. This curve is compared to the population Pdecay(t) = 1 − Pd2(t), where Pd2(t) is the diabatic population in the laser-excited state |φd2⟩; see Figure 7. Note that in the calculation, the excitation and detection process is not included. In particular, the calculated population results from a δ-pulse excitation, which explains the deviations at early times. One has to keep in mind that the decay constant γ is adjusted to give best agreement with the measurement. We use a value of γ = 0.005. For this choice, the internal energy decays exponentially and reaches its final value within 1.5 ps. Noticing the excellent agreement with experiment, the limitations of our one-dimensional model have to be kept in mind. Nonadiabatic,

Figure 6. Decay of the laser-excited state |φd2⟩ determined for different values of the dissipation constant γ, as indicated. The decay is exponential and, in the present range of values, proceeds faster with increasing value of γ.

we show the population P2d(t) of the laser-excited state calculated for various values of the damping constant γ. Here, we use the ground-state vibrational wave function placed in state |φd2⟩ as the initial state. We note that for a resonant excitation with a 30 fs pulse (as in the experiment), an identical population dynamics is obtained. An exponential decay of the population is seen in each case. Furthermore, small oscillations exist that stem from a population back- and forth-transfer between the different states reflecting the vibrational wave packet motion. For the smallest value γ = 0.0001, only about 3% of the population decays in the first 2 ps. The decay is accelerated for larger values of γ. For example, increasing γ by a factor of 10 (γ = 0.001) results in a loss of 25% of the population. The fastest decay is seen for γ = 0.005. We found that the time-scale for the decay is much more sensitive to the diabatic coupling element than to the choice of γ. A faster decay is obtained for a larger coupling between the states |φd2⟩ and |φd4⟩. The ab initio data provide us with a coupling strength of J2 = 0.007 eV. However, due to the expected errors in the quantum chemical calculation, it can also be chosen stronger; see the discussion in section III.B. Furthermore, the choice of a larger coupling element is, in general, reasonable because here a one-dimensional dynamics is described. Each motion along the (orthogonal) coordinates qn is likely to couple the two electronic states, leading to an enhanced population transfer (although for some modes there might be a vanishing coupling). We parametrize the coupling element Vd24(q) with a strength of J2 = 0.04 eV and employ a value of γ = 0.005. This leads to the curves displayed in Figure 7. It is seen that the decay of the initially populated state now proceeds much faster. As before, the intermediate state |φd4⟩ is 1409



Article

Furthermore, the subsequent relaxation proceeds in the lower state on the potential energy curve of the intermonomer torsional motion. It takes place on a much longer time-scale caused by a rather flat potential energy curve and a large moment of inertia corresponding to the rotation of the two monomers with respect to each other. Following this torsional motion to the minimal energy leads to a final state geometry given by a face-to-face orientation, which acts as an efficient trap for the exciton.7 In addition to the clarification of the trapping mechanism, the here presented model allows predictions on how to avoid trapping effects. Actually, for the given morphology (stacked Haggregate in which the monomers are rotationally displaced with respect to each other), we come to the conclusion that a suppression of trapping is not possible because the accompanied relaxational motion involves only tiny distortion and proceeds on an ultrafast time-scale. Suppression can only be achieved in geometrical arrangements, which exhibit larger energy gaps between Frenkel- and CT-states.

Figure 8. Comparison of theory and experiment. Shown is the decaying component extracted from the measured signal (see Figure 2). The theoretical curve corresponds to the depletion of population of the initially excited state Pdecay(t) = 1 − Pd2(t).

multidimensional wave packet dynamics is certainly the way to go but is not feasible for the present system. Our proposal here is that the motion along the well-defined reaction coordinate is the dominant relaxation path driving the relaxation process. Along this reaction path, the system Hamiltonian is derived from ab initio calculations. In section III.A, we have given arguments for the choice of the reaction coodinate and that along this path the excited states are mixed, whereas other paths do not show this mixing. Of course, we cannot exclude that still other kinds of deformational motions are important, which could not be identified yet. Nevertheless, due to the tremendous computational costs of the high-level electronic structure calculations, we restrict ourselves to the onedimensional approach. To conclude, we have introduced a model that aims at the characterization of the decay-dynamics of PBI-dimers after femtosecond excitation. In a first step, a reaction coordinate is defined, which relates Frenkel and charge-transfer configurations of the dimer. The slight deformation of geometry yields a mixing and strong coupling between different electronic states of the dimer, which are of partly Frenkel- and CTcharacter. Afterward, we set up a diabatic Hamiltonian for the nuclear motion on coupled electronic states. The timedependent Schrödinger equation for this motion is solved including effects of dissipation. Taking the accuracy of the ab initio data into account, it is then possible to describe the measured decay. In the future, one could hope that computational technology and quantum chemical methodology will make it possible to provide more-dimensional potential surfaces. It would then be possible to extend the here proposed model and study the influence of further degrees of freedom on the exciton trapping process. With the findings of this work, we are now in the position to discuss the multistep trapping mechanism in rotationally displaced π-stacked dimers. After femtosecond excitation of the energetically higher Frenkel state, a transition of the system to the lower lying Frenkel state is induced. In this process, the system assumes a transient CT character but is passed on directly to the lower lying Frenkel state. The transfer takes place on the subpicosecond time-scale. Thus, considering larger aggregates, this time-scale is comparable to the one of excitation energy transfer (exciton transfer) from the dimer unit to adjacent subunits. However, the accompanying loss of energy during the transfer already diminishes the probability for an exciton transfer.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Phone: +49931-31-86376. Fax: +49-931-31-85331. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support by the German Science Foundation (FOR 1809, GRK 1221, SFB 652) and the VW-Stiftung is gratefully acknowledged. We are thankful to Stefan Schindlbeck for perfoming the transient absorption measurements and Ulrich Kleinekathöfer for helpful discussions.

■

REFERENCES

(1) van Amerongen, H.; Valkunas, L.; van Grondelle, R. Photosynthetic Excitons; World Scientific: Singapore, 2000. (2) Brédas, J.-L.; Beljonne, D.; Coropceanu, V.; Cornil, J. Chargetransfer and Energy-transfer Processes in Pi-conjugated Oligomers and Polymers: A Molecular Picture. Chem. Rev. 2004, 104, 4971−5003. (3) Wasielewski, M. R. Self-Assembly Strategies for Integrating Light Harvesting and Charge Separation in Artificial Photosynthetic Systems. Acc. Chem. Res. 2009, 42, 1910−1921. (4) Lunt, R. R.; Benzinger, J. B.; Forrest, S. R. Relationship between Crystalline Order and Exciton Diffusion Length in Molecular Organic Semiconductors. Adv. Mater. 2010, 22, 1233−1236. (5) Seibt, J.; Marquetand, P.; Engel, V.; Chen, Z.; Dehm, V.; Würthner, F. On the Geometry Dependence of Molecular Dimer Spectra With an Application to Aggregates of Perylene Bisimide. Chem. Phys. 2006, 328, 354−362. (6) Chen, Z.; Stepanenko, V.; Dehm, V.; Prins, P.; Siebbeles, L.; Seibt, J.; Marquetand, P.; Engel, V.; Würthner, F. Photoluminescence and Conductivity of Self-assembled Pi-Pi Stacks of Perylene Bisimide Dyes. Chem.Eur. J. 2007, 13, 436−449. (7) Fink, R. F.; Seibt, J.; Engel, V.; Renz, M.; Kaupp, M.; Lochbrunner, S.; Zhao, H.-M.; Pfister, J.; Würthner, F.; Engels, B. Exciton Trapping in Pi-Conjugated Materials: A Quantum-ChemistryBased Protocol Applied to Perylene Bisimide Dye Aggregates. J. Am. Chem. Soc. 2008, 130, 12858−12859. (8) Howard, I. A.; Laquai, F.; Keivanidis, P. E.; Friend, R. H.; Greenham, N. C. Perylene Tetracarboxydiimide As an Electron Acceptor in Organic Solar Cells: A Study of Charge Generation and Recombination. J. Phys. Chem. C 2009, 113, 21225−21232. 1410



Article

(9) Zhan, X.; Facchetti, A.; Barlow, S.; Marks, T. J.; Ratner, M. A.; Wasielewski, M. R.; Marder, S. R. Rylene and Related Diimides for Organic Electronics. Adv. Mater. 2011, 23, 268−284. (10) Würthner, F.; Stolte, M. Naphthalene and Perylene Diimides for Organic Transistors. Chem. Commun. 2011, 47, 5109−5115. (11) Schubert, A.; Settels, V.; Liu, W.; Würthner, F.; Meier, C.; Fink, R.; Schindlbeck, S.; Lochbrunner, S.; Engels, B.; Engel, V. Ultrafast Exciton Self-Trapping upon Geometry Deformation in Perylene-Based Molecular Aggregates. J. Phys. Chem. Lett. 2013, 4, 792−796. (12) 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. (13) See Supporting Information from ref 11. (14) Würthner, F.; Chen, Z.; Dehm, V.; Stepanenko, V. Onedimensional Luminescent Nanoaggregates of Perylene Bisimides. Chem. Commun. 2006, 105, 1188−1190. (15) Tamura, H.; Ramon, J. G. S.; Bittner, E. R.; Burghardt, I. Phonon-Driven Exciton Dissociation at Donor−Acceptor Polymer Heterojunctions: Direct versus Bridge-Mediated Vibronic Coupling Pathways. J. Phys. Chem. B 2008, 112, 495−506. (16) Tamura, H.; Ramon, J. G. S.; Bittner, E. R.; Burghardt, I. Phonon-Driven Ultrafast Exciton Dissociation at Donor-Acceptor Polymer Heterojunctions. Phys. Rev. Lett. 2008, 100, 107402. (17) Lalov, I. J.; Warns, C.; Reineker, P. Model of Mixed Frenkel and Charge-Transfer Excitons in Donor−Acceptor Molecular Crystals: Investigation of Vibronic Spectra. New J. Phys. 2008, 10, 085006. (18) Heinemeyer, U.; Scholz, R.; Gisslén, L.; Alonso, M. I.; Ossó, J. O.; Garriga, M.; Hinderhofer, A.; Kytka, M.; Kowarik, S.; Gerlach, A.; Schreiber, F. Exciton-Phonon Coupling in Diindenoperylene Thin Films. Phys. Rev. B 2008, 78, 085210. (19) Gisslén, L.; Scholz, R. Crystallochromy of Perylene Pigments: Interference Between Frenkel Excitons and Charge-Transfer States. Phys. Rev. B 2009, 80, 115309. (20) Gangilenka, V. R.; Titova, L. V.; Smith, L. M.; Wagner, H. P.; DeSilva, L. A. A.; Gisslén, L.; Scholz, R. Selective Excitation of Exciton Transitions in PTCDA Crystals and Films. Phys. Rev. B 2010, 81, 155208. (21) Liu, W.; Settels, V.; Harbach, P. H. P.; Dreuw, A.; Fink, R. F.; Engels, B. Assessment of TD-DFT- and TD-HF-based Approaches for the Prediction of Exciton Coupling Parameters, Potential Energy Curves, and Electronic Characters of Electronically Excited Aggregates. J. Comput. Chem. 2011, 32, 1971−1981. (22) Zhao, H. M.; Pfister, J.; Settels, V.; Renz, M.; Kaupp, M.; Dehm, V. C.; Würthner, F.; Fink, R. F.; Engels, B. Understanding Groundand Excited-State Properties of Perylene Tetracarboxylic Acid Bisimide Crystals by Means of Quantum Chemical Computations. J. Am. Chem. Soc. 2009, 131, 15660−15668. (23) Grimme, S. Improved second-order Møller−Plesset Perturbation Theory by Separate Scaling of Parallel- and Antiparallel-Spin Pair Correlation Energies. J. Chem. Phys. 2003, 118, 9095−9102. (24) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: Optimized Auxiliary Basis Sets and Demonstration of Efficiency. Chem. Phys. Lett. 1998, 294, 143−152. (25) Weigend, F.; Häser, M. RI-MP2: First Derivatives and Global Consistency. Theor. Chem. Acc. 1997, 97, 331−340. (26) Hättig, C.; Köhn, A. Transition Moments and Excited-State First-Order Properties in the Coupled-Cluster Model CC2 Using the Resolution-of-the-Identity Approximation. J. Chem. Phys. 2002, 117, 6939−6951. (27) Eichkorn, K.; Treutler, O.; Ö hm, H.; Häser, M.; Ahlrichs, R. Auxiliary Basis Sets to Approximate Coulomb Potentials. Chem. Phys. Lett. 1995, 242, 652−660. (28) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary Basis Sets for Main Row Atoms and Transition Metals and Their Use to Approximate Coulomb Potentials. Theor. Chem. Acc. 1997, 97, 119−124. (29) Liu, W. Dissertation; University of Würzburg: Würzburg, 2010.

(30) Seibt, J.; Winkler, T.; Renziehausen, K.; Dehm, V.; Würthner, F.; Engel, V. Vibronic Transitions and Quantum Dynamics in Molecular Oligomers: A Theoretical Analysis with an Application to Aggregates of Perylene Bisimides. J. Phys. Chem. A 2009, 113, 13475− 13482. (31) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (32) Grimme, S. Semiempirical GGA-type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (33) Dirac, P. A. M. Quantum Mechanics of Many-Electron Systems. Proc. R. Soc. A 1929, 123, 714−733. (34) Slater, J. C. A Simplification of the Hartree-Fock Method. Phys. Rev. 1951, 81, 385. (35) Lee, C. T.; Yang, W. T.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (36) Schäfer, A.; Huber, C.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets of Triple Zeta Valence Quality for Atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5835. (37) Hellweg, A.; Grün, S. A.; Hättig, C. Benchmarking the Performance of Spin-Component Scaled CC2 in Ground and Electronically Excited States. Phys. Chem. Chem. Phys. 2008, 10, 4119−4127. (38) Settels, V. Dissertation; University of Würzburg: Würzburg, 2012. (39) Vosko, S. H.; Wilk, L.; Nusaire, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J. Phys. 1980, 58, 1200−1211. (40) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5642. (41) Werner, H. J.; Meyer, W. MCSCF Study of the Avoided Curve Crossing of the Two Lowest 1+ States of LiF. J. Chem. Phys. 1981, 74, 5802−5807. (42) Podolsky, B. Quantum-Mechanically Correct Form of Hamiltonian Function for Conservative Systems. Phys. Rev. 1928, 32, 812−816. (43) Lauvergnat, D.; Nauts, A. Exact Numerical Computation of a Kinetic Energy Operator in Curvilinear Coordinates. J. Chem. Phys. 2002, 116, 8560−8570. (44) van Kampen, N. G. Stochastic Processes in Physics and Chemistry; Elsevier: Amsterdam, 1992. (45) Gisin, N.; Percival, I. C. The Quantum-State Diffusion Model Applied to Open Systems. J. Phys. A 1992, 25, 5677−5691. (46) Gardiner, C. W. Handbook of Stochastic Methods; Springer: Berlin, 2002. (47) Feit, M. D.; Fleck, J. A.; Steiger, A. Solution of the Schrödinger Equation by a Spectral Method. J. Comput. Phys. 1982, 47, 412−433. (48) Schlesinger, M.; Mudrich, M.; Stienkemeier, F.; Strunz, W. T. Dissipative Vibrational Wave Packet Dynamics of Alkali Dimers Attached to Helium Nanodroplets. Chem. Phys. Lett. 2010, 490, 245− 248. (49) Grüner, B.; Schlesinger, M.; Heister, P.; Strunz, W. T.; Stienkemeier, F. Vibrational Relaxation and Dephasing of Rb2 Attached to Helium Nanodroplets. Phys. Chem. Chem. Phys. 2011, 13, 6816−6826. (50) Schlesinger, M.; Strunz, W. T. Detailed Model Study of Dissipative Quantum Dynamics of K2 Attached to Helium Nanodroplets. New J. Phys. 2012, 14, 013029. (51) Kleinekathöfer, U.; Kondov, I.; Schreiber, M. Perturbative Treatment of Intercenter Coupling in the Framework of Redfield Theory. Chem. Phys. 2001, 268, 121−130. (52) Egorova, D.; Kühl, A.; Domcke, W. Modeling of Ultrafast Electron-Transfer Dynamics: Multi-Level Redfield Theory and Validity of Approximations. Chem. Phys. 2001, 268, 105−120. (53) May, V.; Kühn, O. Charge and Energy Transfer Dynamics in Molecular Systems, 3rd ed.; Wiley-VCH: Weinheim, 2011. 1411



Article

(54) Guthmuller, J.; Zutterman, F.; Champagne, B. Multimode Simulation of Dimer Absorption Spectra from First Principles Calculations: Application to the 3,4,9,10-Perylenetetracarboxylic Diimide Dimer. J. Chem. Phys. 2009, 131, 154302. (55) Roden, J.; Eisfeld, A.; Dvorak, M.; Bünemann, O.; Stienkemeier, F. Vibronic Line Shapes of PTCDA Oligomers in Helium Nanodroplets. J. Chem. Phys. 2011, 134, 054907.

1412


Theoretical Analysis of the Relaxation Dynamics in Perylene Bisimide

Recommend Documents