Ultrafast Nonadiabatic Dynamics of Singlet Fission: Quantum

Subscriber access provided by NORTH DAKOTA STATE UNIV

Article

Ultrafast Nonadiabatic Dynamics of Singlet Fission: Quantum Dynamics with the Multilayer Multiconfigurational Time-dependent Hartree (ML-MCTDH) Method Jie Zheng, Yu Xie, Shengshi Jiang, and Zhenggang Lan J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.5b09921 • Publication Date (Web): 04 Jan 2016 Downloaded from http://pubs.acs.org on January 8, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Ultrafast Nonadiabatic Dynamics of Singlet Fission: Quantum Dynamics with the Multilayer Multiconfigurational Time-dependent Hartree (ML-MCTDH) Method Jie Zheng,ab Yu Xie,*a Shengshi Jiangac and Zhenggang Lan*ab a

Key Laboratory of Biobased Materials, Qingdao Institute of Bioenergy and Bioprocess Technology, Chinese Academy of Sciences, 266101 Qingdao, Shandong, People’s Republic of China

b

University of Chinese Academy of Sciences, 100049 Beijing, People’s Republic of China

c

College of Physics, Qingdao University, 266071 Qingdao, Shandong, People’s Republic of China

* E-mail: [email protected], [email protected]

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT Singlet fission (SF) is supposed to potentially improve the efficiency of solar energy conversion in organic photovoltaic systems. The multilayer multiconfigurational time-dependent hartree (ML-MCTDH) method was employed to describe the singlet fission of the pentacene system with a three-state model. The ML-MCTDH result agrees well with the previous simulations using the Redfield theory, the hierarchical equation of motion (HEOM) and the symmetrical quasi-classical (SQC) theory. We carefully investigated the role of vibrational modes with different frequencies in singlet fission dynamics. Interestingly, we observed the important contribution of a few modes with frequency resonance to electronic transition. Such a finding can be understood by revisiting the super-exchange mechanism within the framework of Fermi’s golden rule. As a numerically exact method, ML-MCTDH not only provides an accurate description of the microscopy insight of the SF dynamics but also provides benchmark results to examine the performance of other approximated dynamical methods.

KEYWORDS ML-MCTDH, singlet fission, super-exchange mechanism

2


Page 2 of 41

Page 3 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION Singlet fission (SF) refers to a photoinduced excited-state process that converts one excited singlet exciton into two triplet excitons, namely a correlated pair of triplet states.1-5 In recent years, many works have presumed that organic photovoltaic materials with efficient SF may suppress the Shockley−Queisser limit,6 resulting in a large improvement in the efficiency of organic solar cells. Thus, the understanding of the physical insight of the SF processes at the molecular/atomic level is critical to the proposal of useful ideas for the rational design of photovoltaic materials.2, 7 Normally, a simplified kinetic model2 is used to illustrate the SF process in molecular aggregates composed of several chromophores. One widely accepted mechanism is one in which the excited singlet state (S1) first converts into a correlated triplet state pair (TT) on adjacent chromophores. Previous work noted that such a TT state is in fact a multi-excitonic (ME) state3, 8-11, and thus, the energy of the TT state is not exactly equivalent to the sum of two triplet states. Such a TT state (or ME state) subsequently dissociates into two separated triplet states (T+T). Needless to say, the efficient S1TT internal conversion is a prerequisite for effective SF processes; therefore, the TT formation step has been the subject of much research. Please refer to several interesting reviews for a discussion on the SF mechanism.1-5, 10, 12-13

Many experimental works employed the time-resolved spectra to study the SF

dynamics of various systems,8, 1,3-diphenylisobenzofuran.3,

14-41

for example, tetracene, pentacene and

8, 14-23, 25-27

Extensive theoretical efforts were also

undertaken to understand the SF mechanism from different perspectives—for instance, the construction of the effective model Hamiltonian from electronic-structure calculations methods

9-10, 42-53

, the simulation of the SF processes with different dynamical

10, 41-43, 47-48, 53-58

and the examination of the morphology effects 53, 59-62. In the

process of these works, some basic knowledge on SF dynamics has been accumulated. Two mechanisms were proposed to explain the short-time SF dynamics 3, 10—namely, the direct mechanism driven by the S1TT electronic couplings and the mediated mechanism by an intermediate charge transfer (CT) state. Although extensive debates exist on the importance of these two mechanisms for specific systems5, 10, 46-48, 53-54, 57, 63-65

, it is also widely accepted that the dominant mechanism may be highly dependent

on the systems under study.48, 53, 66 For instance, Tamura et al.53 once pointed out that 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 41

the ultrafast SF dynamics in stacked 6,13-bis(triisopropyl-silylethynyl)-pentacene (TIPS-pentacene) complex occurs mainly through the CT-mediated mechanism, although both CT-mediated and direct pathways are possible. However, the direct mechanism via a conical intersection becomes dominant for the C2h stacked rubrene complex, while the symmetry-breaking vibrations must be involved. Yost et al.66 demonstrated that the strong mixture of the CT and local excited state is responsible for the SF dynamics of pentacene, TIPS-pentacene and 6,13-di(2’-thienyl)pentacene (DTP), while the involvements of the CT state is not possible for the 6,13-di-benzothiophenepentacene (DBTP), 6,13-di-biphenyl-4-yl-pentacene (DBP) and 6,13-diphenylpentacene (DPP) systems. Berkelbach et al.48 once figured out that the mixture of the local excitation and CT transition is highly dependent on the intermolecular distance, electronic couplings, and dielectric environments. Among all reported organic systems, the pentacene system displays rather high SF yield and thus has received a high level of research interest3,

10

. Recent

time-resolved spectroscopy experiments clarified that the formation of the TT state in crystalline pentacene, pentacene or pentacene/C60 thin films falls into the ultrafast time scale of 80–200 fs10, 20, 67-68. In addition, a few experimental works suggested that the SF dynamics of pentacene is governed by the coherent mechanism, i.e. coherent superposition of S1 and multiexcitionic state.3, 15 Recent experimental work of the ultrafast two-dimensional electronic photon echo (2DPE) spectroscopy pointed out that vibronic couplings, instead of electronic couplings, play an essential role in the internal conversion of S1 to TT. 69 Zimmerman et al. employed the restricted active space double spin-flip (RAS-2SF) method to examine the excited state of pentacene dimer5, 63. Because the CT state is high-lying with respect to S1 and TT in their calculations, they ruled out the possible involvement of the CT state in the SF dynamics. They also emphasized that the direct SF dynamics may be efficiently achieved via S1-TT avoided crossing with the change of the intermolecular distance. However, several research groups noted that the ultrafast SF dynamics via the CT-mediated mechanism is possible even if the CT state lies higher than the S1 and TT states. Particularly, the involvement of the CT-mediated mechanism in the ultrafast SF dynamics of pentacene and other systems was studied within different dynamical frameworks, such as the quantum 4


Page 5 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


master equation method based on the reduced density matrix by several groups—Reichman1, 47-48, 54, Grozema60, 70, Ratner

43, 71

, Eaves58, etc.—and SQC by

Tao55-57, etc. In their views, the ultrafast TT formation in the experimental observation implies the involvement of the CT-mediated mechanism for pentacene because the S1-TT direct electronic coupling is rather weak.2, 47, 54 When the pentacene crystal is considered, several theoretical works clarified that the extensive involvement of the CT excitation in the lowest adiabatic excited states.46, 48, 72-74 The initial single-excited states prepared by the photoexcitation should already contains the significant contribution of the CT characters, thus the low-lying TT state should be quickly formed by the one-electron transfer mechanism. In this sense, the involvement of the CT excitation is still responsible for the ultrafast SF dynamics of pentacene crystal, while here the CT character and other electronic configurations should be directly mixed instead of virtually mixed.46, 48 However, many open questions still remain to clarify the details of the SF dynamics. For example, the stacking configurations may have potential effects on the SF dynamics.5, 50, 53, 60, 63, 70, 75-76 Alternatively, the energies of the involved state seem to be relevant to the electronic structure methods in the calculations. For instance, recent work by Thoss and co-workers75 noted that the CT state may become the lowest excited state and the TT state may be located in the high energy domain at some dimer or trimer configuration at the CASPT2 level. Thus, the accurate description of the involved electronic states is still quite challenging. Thus far a number of theoretical methodologies have been applied to the ultrafast singlet fission dynamics in pentacene systems, such as the Redfield theory10, 47-48, 54, SQC55-57, the trajectory surface-hopping (TSH)77 and the HEOM method47. However, the Redfield theory is valid only in the case of weak system–bath couplings.78 Although the HEOM method within the framework of the reduced density matrix approach is assumed to be an accurate way to address the dynamics governed by the intermediated system–bath coupling Hamiltonian47, 79-80, it is very computationally demanding. Both SQC and TSH methods are quite approximated. Thus, more theoretical work is needed to examine the microscopic picture of the SF dynamics, which are essentially electronically nonadiabatic dynamics.

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 41

In recent years, it is also possible to study the nonadiabatic dynamics of complex systems using a numerically exact full quantum dynamics method—for instance, the multilayer multiconfigurational time-dependent Hartree (ML-MCTDH) theory first formulated by Wang and Thoss.81-85 Meyer and Worth also independently provided a formally identical formulation termed as “cascading MCTDH” but did not provide any implementation owing to coding difficulties86. Manthe found a recursive way to set up the ML-MCTDH expansion to an infinite number of layers, which also allows us to easily derive the corresponding working equation of motion (EOM) for each layer.87-88 Later, Meyer and co-workers employed this recursive idea to successfully implement ML-MCTDH into the Heidelberg MCTDH package89, which also allows large-scale ML-MCTDH calculation of the quantum evolution of extremely complex systems.90-92 Both MCTDH and ML-MCTDH have also been applied to check the photoinduced dynamics of various photovoltaic systems53,

80, 82-85, 93-111

, such as

dye-sensitized solar cells82-85 and organic solar cells80, 93-97, 99-106, 108-111. The research topics range from the photoinduced charge transfer96-97,

101-106, 108, 111

to the

excited-state energy transfer80, 95, 109-110 and even the singlet fissions53. For example, Tamura, Burghardt and co-workers recently employed the standard MCTDH to study the SF dynamics of acene derivatives.53 They found that the stacked status in acene crystal may have a strong impact on the interstate couplings and the SF mechanism. Here, we applied the ML-MCTDH method to understand the ultrafast nonadiabatic dynamics in a model system for singlet fission process in pentacene system. The purpose of this work is mainly twofold. First, we wished to compare the ML-MCTDH study to recent theoretical works by using the Redfield, HEOM and SQC methods. Thus, we employed the three-state Hamiltonian proposed by Gao-Zhu-Reichman and co-workers for comparison.10, 47 Second, instead of obtaining only the overall SF dynamics (such as population dynamics with respect to time), we wished to study more dynamical details, such as the influence of vibrational modes (phonon modes) with different frequencies on SF dynamics and the electronic coherence. The ML-MCTDH calculations explicitly figured out the roles of modes with different frequencies in the SF dynamics, and then the corresponding key features were qualitatively understood by Fermi’s golden rule.9,

78, 112

The

investigation of the effects of different modes also helped us speed up the 6


Page 7 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


convergence of ML-MCTDH calculation by the optimal and reasonable construction of ML-MCTDH expansion (tree). This work should provide a complementary view of the mechanism of the SF dynamics.

2. COMPUTATIONAL DETAILS 2.1 Hamiltonian The current work employs a three-state model Hamiltonian to describe the SF dynamics of a pentacene pair, which was first developed and employed by Gao-Zhu-Reichman and co-workers in the dynamics study using the Redfield and HEOM theories10, 47. This model was also employed by Tao in his SQC dynamics simulation.57 The three-state model consists of a single excited state S1 of the molecular pair, a charge transfer state CT, and a coupled triplet pair state TT. A system–bath Hamiltonian was used to describe the SF dynamics. The Hamiltonian includes three parts: the electronic Hamiltonian

H el (system),

environmental phonon part H ph (bath), and the electronic-phonon couplings H el − ph (system–bath couplings); i.e., H = H el + H ph + H el − ph H el = ∑ φk Ek φk + ∑ φk Vkl φl l ≠k

k Nb

1 H ph = ∑ ω j  Pj2 + Q 2j  j 2 H el − ph = ∑ φk φk k

(1)

Nb

∑ ( −c

kj

Qj )

j

Here, φk represents the electronic state and three electronic states are included in this work—namely, S1, TT and CT. Ek is the energy of each electronic state, and Vkl is their interstate coupling. Qj, Pj represent the position and momentum of the jth bath mode. ωj, and ckj are the frequency and the electronic-phonon coupling constant of the corresponding mode, respectively. Nb is the number of bath modes coupled to each electronic state. Similar to previous work, we neglected the system–bath couplings for the off-diagonal elements of the electronic Hamiltonian.

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 41

The electronic state energies and the interstate couplings in the above Hamiltonian were taken from previous works10, 57; see Table 1. These values are good representatives of transition energy and electronic couplings obtained from experimental113-116 and theoretical54, 63, 73, 117 results.

Table 1. Elements of the electronic Hamiltonian matrix10, 57—i.e., energies Ek of the three states and their couplings Vkl (unit: eV). The initial energy of TT state (ETT) is set to 0 eV. S1

TT

CT

S1

0.2

0

-0.05

TT

0

0

-0.05

CT

-0.05

-0.05

0.3

2.2 Spectral density Following previous studies, the frequencies and coupling constants of bath modes (photon modes) are characterized by the continuous Debye-type spectral density, J (ω ) =

2λωωc , ω 2 + ωc2

(2)

where λ is the reorganization energy that represents the coupling strength between the electronic degrees of freedom and the bath, and ωc is the characteristic frequency of the bath. In this work, ωc = 0.18 eV and λ = 0.1 eV, according to previous works2, 10, 54, 57, 118

. Each electronic state couples to its corresponding bath. Once the continuous spectral density is generated, it can be discretized to an

arbitrary number of bath modes according to

J k (ω ) =

π 2

N

∑ c δ (ω − ω ) 2 ki

i

(3)

i =1

The coupling coefficient cki of each mode is obtained by Eq. (4), given a sampling 8


Page 9 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interval ∆ω98 2  cki =  J k (ωi ) ∆ω  π 

1

2

(4)

The sampling interval defines the Poincare recurrence period τ p = 2π / ∆ω .96 For the case of t < τ p , the observed dynamics are effectively irreversible. In the current study, the long-time propagation of the quantum dynamics was employed to examine whether the recurrence appears. The result shows that such effects can be safely neglected within our propagation time duration. This electron-phonon coupling constant can also be viewed as the vibronic coupling strengths if we transformed the system-plus-bath Hamiltonian into the equivalent vibronic coupling Hamiltonian. Thus, in the below discussion, we always use the picture of the vibronic coupling Hamiltonian, which treats the phonon modes as vibrational modes and electron-phonon coupling as vibronic couplings.

2.3 ML-MCTDH ML-MCTDH81, 87-88, 90-92, 119, which can address thousands degrees of freedom, is a powerful extension of the standard MCTDH86 approach. The MCTDH method uses the time-dependent basis to represent the wave function n1

nf

f

j1 =1

j f =1

κ =1

Ψ ( Q1 ,L, Q f , t ) = ∑L ∑ Aj1L j f ( t )∏ ϕ (jκκ ) ( Qκ , t )

(5)

where Q1 ,L , Q f are coordinates describing degrees of freedom (DOFs), Aj1L j f denote the time-dependent expansion coefficients, and {

ϕ (jκκ ) } denote the

time-dependent basis functions, known as single particle functions (SPFs). In the ML-MCTDH scheme, the wave function is expanded by SPFs recursively to form a multilayer tree structure, which is in nature a tensor decomposition method. In other words, an SPF in the upper layer is an expansion of SPFs in the lower layer; i.e.,

ϕ

l −1;κ1Lκ l −1 m

(Q

l −1;κ1Lκl −2

κl −1

)

n1

nκl

, t = ∑L ∑ A j1 =1

j pκ =1 l

l ;κ1Lκl −1 m; j1L j pκ

l

pκl

( t )∏ϕ lj;κ Lκ 1

κl =1

κl

9


l

(Q

l ;κ1Lκl −1 κl

,t

)

(6)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 41

where l denotes the layer depth, and κ1, . . ., κl−1 denote the indices of the logical degrees of freedom starting from each node on the top layer down to a particular primary coordinate. According to this idea, a tree (or hierarchy) structure was constructed to represent the wave function. Applying the variational principle and the recursive algorithm provided by Manthe87-88, the equations of motion (EOMs) and all quantities entering the EOMs can be obtained. The final ML-MCTDH expansion can be represented as a multilayer tree; see Figure 1. In this paper, all ML-MCTDH calculations were performed using the Heidelberg MCTDH package89.

Figure 1. Tree structure for the ML-MCTDH wavefunction. Here, circles denote SPF bases on each layer, squares on the bottom layer denote primitive bases and qelec denotes electronic degree of freedom.

2.3.1 Initial condition and physical observables Within the Condon approximation, the initial wave packet ψ (t = 0) is obtained by the vertical excitation of the ground vibrational level of the electronic ground state to the S1 state. Thus, at time zero, the nuclear wavefunction is assumed to be the ground vibrational level of the ground electronic state at T=0 K. The temperature effect was also examined in a rather approximated way; see Appendix I for details. The diabatic population and electronic coherence were measured by the diagonal elements Pi ,i (t ) and the off-diagonal elements Pi, j (t ) of the reduced electronic density matrix, which are given as

{

∧

}

Pi , j (t ) = Trb φ j φi ρ (t ) = ∫ ρi , j ( Q, t ) dQ 10


(7)

Page 11 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Here, the density operator is defined as ρˆ (t ) = ψ (t ) ψ (t ) , and the trace (Tr) is over all bath degrees of freedom. To give an alternative theoretical view in the description of the SF dynamics, the rate equation based on Fermi’s golden rule9, 78, 112 was also addressed. All details can be found in the discussions below and the Supporting Information (SI).

3. RESULTS and DISCUSSION

3.1 ML-MCTDH study of SF dynamics In principle, after the discretization of bath degrees of freedom, all bath modes with distinguished frequencies should be considered in the completed treatment of the SF dynamics. However, instead of including all modes directly in the ML-MCTDH calculations using a brute force method, we decided to first check the contribution of the vibrational modes within different frequency domains and then slowly access the overall SF dynamics by adding increasing vibrational degrees of freedom. In this treatment, we first considered several “important” active modes to define a reduced dimensional Hamiltonian, examined the SF dynamics governed by these important modes, then added “bystander” modes step by step until the reasonable convergent result is obtained. This way not only clearly identifies which modes play major roles in the SF dynamics and which modes are bystander inactive modes but also largely reduces our efforts to build the reasonable ML-MCTDH tree in the dynamical calculations of the high-dimensional model.

3.1.1 SF dynamics including modes within different frequency regions We divided all discretized vibrational modes into several groups, each of which contains the modes belonging to a particular frequency domain. It is well known that ultrafast nuclear-electronic coupled dynamics may be strongly modified by vibrational modes with strong vibronic couplings (system–bath couplings) or modes resonant with the electronic transitions. These basic principles provide us with some clues to determine the frequency range of each domain. 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The initial test calculations indicate that the vibrational modes with frequencies larger than 0.4 eV do not have any influence on the SF dynamics, as shown in Figure S1 in SI. This observation is consistent with the fact that the vibrational frequencies of these bath modes are too far from any electronic transition; thus, these modes do not play any role in the SF dynamics, and we did not consider them in the following calculations. After some trial calculations, the frequency domain (0–0.4 eV) was divided into eight regions (R1–R8)—namely, 0–0.075 eV (R1), 0.075–0.108 eV (R2), 0.108–0.125 eV (R3), 0.125–0.165 eV (R4), 0.165–0.2 eV (R5), 0.2–0.3 eV (R6), 0.3–0.35 eV (R7), and 0.35–0.4 eV (R8). In each frequency domain, several modes were constructed, and their vibronic coupling strengths were computed from the spectral density. To provide a reasonable representation of the system–bath coupling, different numbers of bath modes were employed in the discretization of each frequency domain. This approach allows us to ensure that the number of modes is sufficient for the correct results of the ML-MCTDH calculations. The corresponding convergent results of dynamical calculations are shown in Figure 2.

12


Page 12 of 41

Page 13 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Population dynamics when the vibrational modes in different frequency regions are considered. (a) 0–0.075 eV (R1), (b) 0.075–0.108 eV (R2), (c) 0.108– 0.125 eV (R3), (d) 0.125–0.165 eV (R4), (e) 0.165–0.2 eV (R5), (f) 0.2–0.3 eV (R6), (g) 0.3–0.35 eV (R7), (h) 0.35–0.4 eV (R8). Different numbers of bath modes were employed in each frequency domain for the convergence test.

After the inclusion of the vibrational modes within frequency domain R3 (0.108– 0.125 eV), R6 (0.2–0.3 eV), or R8 (0.35–0.4 eV), the resulting SF dynamics are 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

almost identical to those of pure electronic evolution (Appendix II). In one word, these bath modes are bystanders that are not involved in the SF dynamics. The further ML-MCTDH calculations with high-dimensional Hamiltonians support this conclusion (see discussions below). Thus, we may neglect them in the SF dynamics calculations of high-dimensional cases. Alternatively, we may employ a limited number of modes to represent their existence, and the minimum mode number is determined by extensive test calculations. This helps us largely reduce the number of modes (or dimensionality) in the overall ML-MCTDH calculations. In Figure 2, the most significant finding is that the vibrational modes within frequency domain R5 (0.165–0.2 eV) induce the ultrafast SF dynamics even without the attendance of other modes. Although the CT state lies higher than the S1 and TT states, the ultrafast S1 decay and the TT formation are observed. Through all SF dynamics, the CT state is weakly populated. The current reduced-dimensional model calculations supports that the ultrafast SF dynamics are possible via the CT-mediated super-exchange mechanism. This conclusion is consistent with previous studies of dynamics10, 47, 54-57. Most importantly, even if only a few modes in R5 were included, the current dynamics results—for instance, the timescale of the TT formation—are very close to that in previous work. In this sense, our calculations clearly show that the modes within the R5 domain play an essential role in the ultrafast SF dynamics governed by the super-exchange mechanism. The key impact of these modes on the SF dynamics can easily be explained by the fact that the frequencies of these modes (0.165–0.2 eV) are nearly resonance with the S1-TT energy gap. The similar resonance effects have been mentioned by previous studies.54, 60 More discussions on the resonance issue and the super-exchange mechanism are found in the discussion below and in SI. Certainly, visible differences exist between the dynamics, including the modes in R5 (Figure 2(e)) and the dynamics with more modes from a board frequency domain (see discussions below). The time-dependent S1 and CT populations display obvious oscillation in Figure 2 (e), whereas such oscillation patterns almost disappear in the previous theoretical work10, 47, 54-57 employing a large number of bath modes with a rather board frequency distribution. The roles of the vibrational modes belonging to two other frequency domains (R4 and R7) were also investigated, as shown in Figure 2(d) and (g). Although both of 14


Page 14 of 41

Page 15 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


them do not induce the obvious population transfer to the TT states, they tend to weaken the oscillation pattern of the time-dependent S1 and CT population, which becomes very minor at 400 fs. Possibly because the low-frequency modes (R1 and R2) display the visible vibronic coupling strength, they show some effects on the SF dynamics, but they seem to play a minor role here. The similar fact was noticed by Reichman and coworkers.54 Their minor effects may result from the fact that their slow motion should be well separated from the fast SF dynamics. Renaud et al. once argued that the intermolecular modes with low frequency might have a potential influence on the SF dynamics60. This view holds true when the SF dynamics occur within a similar time scale to the vibrational motions, the resonance condition is satisfied or the different stacking statuses (static disorders) created by intermolecular vibrations are considered. Thanks to the weak impact of the low-frequency modes on the ultrafast SF dynamics, the convergence difficulties of the ML-MCTDH calculations are largely alleviated when the multidimensional models were considered. In principle, the large number of basis functions should be employed in the ML-MCTDH calculation to expand the single-particle function of each low-frequency mode in the case of visible vibronic coupling. In this situation, the convergence becomes very difficult, particularly for the long-time propagation. However, the SF dynamics are very fast, and only the short-time propagation is required. Thus, it is still possible to achieve convergence in the ML-MCTDH calculations when the multidimensional models were considered.

3.1.2 The SF dynamics of the high-dimensional Hamiltonian

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. SF dynamics governed by the multidimensional Hamiltonian. (a) Modes in the R4, R5 and R7 regions are included in the ML-MCTDH calculations. Two different Ml-MCTDH trees give the same results. (b) The overall dynamics, including all frequency regions (R1–R8, 183 modes), compared to the dynamics of basic_tree (90 modes). Two different ML-MCTDH trees—namely, tree1 and tree2—were built to check the convergence. The coupling parameters ckl of the modes in all frequency regions are given in Table S1.

We included the modes belonging to three frequency domains (R4, R5 and R7) to perform the ML-MCTDH calculations. The convergence of the ML-MCTDH calculations was examined by employing different numbers of modes in each frequency domain, as shown in Figure 3(a). As discussed in Figure 2, the modes belonging to frequency domain R5 induce the pronounced population transfer to the TT state, whereas the populations of S1 and CT states display the large oscillation. By adding the modes within the R4 and R7 frequency range, the oscillation is largely suppressed. Although the weak oscillation patterns do not fully vanish, both the increase in the TT population and the decrease in the S1 population become almost monotonic, except that the strong S1 population recurrence appears in the early-stage dynamics (< 50 fs). Thus far, the feature of the SF dynamic has been almost fully consistent with the result of previous studies10, 47, 54-57, but here we included only a small number of the vibrational modes from three different frequency domains (R4, R5 and R7). The above calculations with the inclusion of the modes belonging to three different frequency domains (R4, R5 and R7) defined a reasonable ML-MCTDH tree. This “basic” tree provides us with a foundation for the following ML-MCTDH 16


Page 16 of 41

Page 17 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculations with higher dimensionality. In practice, this “basic” tree serves as a “key branch” in the bigger ML-MCTDH tree used for a higher dimensional calculation. When more modes are included, we simply added a few new branches into the whole ML-MCTDH tree and do not change the topology of the “key branch”. This largely reduces our efforts to build the ML-MCTDH tree for efficient calculation. When the modes from regions R2, R3, R6 or R8 are added, the overall dynamics remain almost unchanged, indicating that these modes do not strongly couple with the SF dynamics. The inclusion of the modes with extremely low frequencies (R1) slightly slows down the population transfer, whereas the basic features of the SF dynamics are kept; see Figure S2 in SI. Finally, the modes from R1–R8 were all introduced to construct the high-dimensional multimode Hamiltonian to access the continuum bath limit. Two different ML-MCTDH trees were designed, and their results are identical, providing a double check for the convergence; see Figure 3(b). Compared to the results of the basic tree in Figure 3(a), the time-dependent populations of all electronic states become very smooth, except that the S1 and CT populations display oscillation in the very early stage of the SF dynamics. The calculated time scale of the TT formation is less than 300 fs, which agrees very well with available experimental

3, 10, 38

and

theoretical results10, 47, 54-57.

Figure 4. The SF dynamics governed by the Hamiltonian including the modes of R4, R5 and R7. (a) The vibronic couplings (k) of all modes are set to the 1/2 and 3/2 of their original values. (b) The influence of the energy of CT with respect to TT on the SF dynamics.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

It is well known that the vibronic coupling strength should have a visible influence on the electronic-nuclear coupled dynamics. To examine this effect, we modified the vibronic coupling constant in the ML-MCTDH calculation that employed the basic tree, as shown in Figure 4(a). When the vibronic coupling strengths of included modes decrease to one-half of their original values, the SF dynamics become much slower, and the TT population only reaches ~30% at 400 fs. At the same time, the initial oscillation pattern of the S1 population lasts longer. When the vibronic coupling strengths increase to 3/2 of their original values, a dramatic speedup of the SF dynamics is observed, and the TT population quickly rises to ~0.85 at 400 fs. The current ML-MCTDH studies show that the modes within the frequency domain of 0.165–0.2 eV play a major role in the SF dynamics. Next, we attempted to provide a more detailed investigation into the contribution of these modes. Starting from the case in which the vibronic coupling strengths increase to 3/2 of their original values for the basic tree (R4+R5+R7 in Figure 4(a)), we designed two additional situations for comparison. (1) Only the vibronic coupling strengths for the modes within the R5 frequency domain (0.165–0.2 eV) increase to 3/2 of their original values, whereas the coupling parameter of all modes in R4 and R7 remains unchanged. (2) For only the modes within the R4 and R7 domains, we increased their vibronic coupling strengths to 3/2. The first situation results in faster SF dynamics, as shown in Figure S3(a) in SI, similar to the results of increasing the vibronic coupling for all involved modes (R4+R5+R7). However, the second case does not result in significantly faster SF dynamics. This finding again demonstrates the importance of the modes belonging to the R5 frequency domain (0.165–0.2 eV). Among all modes in the R5 region, we successfully identified that only four modes with frequencies 0.1895 eV, 0.193 eV, 0.1965 eV, and 0.2 eV (in the basic tree) play major roles in the SF dynamics. If we considered only the above four modes, which couple with the TT state and set up three-state and four-mode Hamiltonian models, the SF dynamics based on this model are easily obtained; see Figure S3(b) in SI. Although this model results in slightly slower SF dynamics, the overall time scale of the TT formation does not differ dramatically with the result of high-dimensional models. This difference becomes even less when we increase the vibronic coupling 18


Page 18 of 41

Page 19 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


strengths of the above four modes by a factor of 3/2. We also assessed the influence of the energy the CT state on the SF dynamic in Figure 4(b). When the CT state energy rises to 0.5 eV (S1-CT energy gap: 0.3 eV), the decay of the S1 state and the increase in the TT state become significantly slower. In this sense, the S1-CT energy gap is critical to determine the timescale of the SF dynamics. The underlying reason can be understood by the super-exchange mechanism; see the discussion below for more details.

3.1.3 Electronic coherence

Figure 5. The time-dependent S1-CT electronic coherence (real part (a) and imaginary part (b)). Four models were considered here: the pure electronic Hamiltonian, the Hamiltonian including the modes belonging to the R5 domain, the Hamiltonian including the modes belonging to R4+R5+R7, and the Hamiltonian including all involved modes (R1–R8).

When increasing vibrational modes are included in the dynamical calculations, the oscillation pattern of the electronic population is largely suppressed, and the S1 population decay becomes almost monotonic. In principle, such an oscillation pattern should be relevant to the electronic coherent motions, which are characterized by the off-diagonal elements of the reduced electronic density matrix. As shown in Figure 5 and Figure S4 in SI, the oscillation pattern of the electronic population (Appendix II) in the pure electronic dynamics is the result of the electronic coherence because the

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

oscillation periods of both quantities are identical. Significant S1-CT coherent motion was observed, whereas the coherences between other states are much weaker; see Figure S4 in SI. The electronic coherence lasts forever in the pure electronic dynamics, but it vanishes quickly when more modes are included. In the high-dimensional cases, both the real and imaginary parts of ρ S1 ,CT ( t ) display obvious oscillatory patterns in the very early state of the SF dynamics, whose periods are also consistent with the appearance of the population recurrence. However, after only one or two periods, the electronic coherence disappears quickly. This indicates that the early stage of the SF dynamics is governed by coherent electronic motions, whereas the decoherence occurs rapidly when a sufficient number of bath modes are introduced. In the end, the imaginary term of ρ S1 ,CT ( t ) decays completely to zero, whereas the real part of

ρ S ,CT ( t ) becomes a small constant value. A similar phenomena is also obtained for 1

the S1-TT and CT-TT coherence terms; see Figure S4 in SI, although their values are much smaller. In fact, the presence of the significant electronic coherence in the early stage of the ultrafast dynamics and the rapid decoherent effects was noticed by previous studies on the photoinduced reactions of organic photovoltaic molecules111, 120

.

3.2 SF mechanism revisited 3.2.1 Super-exchange mechanism Both previous calculations and the current work demonstrated that the ultrafast SF dynamics are possible via the super-exchange mechanism10, 43, 47, 54-57, even if the CT state lies higher than the initial S1 state. Previous work also provided the qualitative description of the SF dynamics governed by the super-exchange mechanism in the framework of the Marcus-like theory, including the CT state as a “bridge”54-55. In fact, the super-exchange mechanism may also be responsible for the ultrafast dynamics of the bridged-mediated electron transfer, which was widely studied by various theoretical methods78,

112, 121-124

, ranging from approximated

perturbative treatments78, 112, 121-122 to rigorous quantum dynamics123-124. 20


Page 20 of 41

Page 21 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Schematic diagram of super-exchange mechanism. Three electronic states (S1, CT and TT) are included, and each state supports many vibrational levels. The minimum energies of three electronic states are E10, E20 and E30. In the ML-MCTDH dynamics, the modes belonging to the particular frequency domain (0.165–0.2 eV) are mainly responsible for the ultrafast TT formation. It is possible to employ a more approximated theory to understand this observation. In past decades, intensive efforts have been made to develop an approximated theory to describe the electron transfer rate by considering the quantum vibrational motions explicitly. For instance, Bixon and Jortner estimated the transition rate between different vibronic doorway states in the electron transfer via a bridge within the framework of Fermi’s golden rule.9, 78, 112 Troisi and coworker once used this approach to analyze the SF process of linear chain molecules9. We borrowed this idea for the discussions below. As shown in Figure 6, three electronic states (S1, CT and TT) are included, and each state supports many vibrational levels. The minimum energies of the three electronic states are labeled E10, E20 and E30, respectively. The interstate couplings are VS1,CT, VS1,TT and VCT,TT, which refers to V12, V13, and V23 in Eq. (8)–(12). Here, the initial vibronic doorway state |α (with energy ) on the S1 state is coupled with the bridged vibronic doorway state |β (with energy ) on the CT state. The latter doorway state |β further interacts with the final target doorway state |γ (with energy ) on the TT state. When the |α |β coupling is considered, the first-order perturbation theory gives the initial vibronic doorway state (labeled | ):

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

α% = α + ∑ β β

where

χ β χα

V12 χ β χα Eβ − Eα

Page 22 of 41

(8)

denotes the overlap integral between two different vibrational levels.

Here, we use V12 to represent the electronic S1-CT couplings to be consistent with Eq. (1). The | |γ coupling is easily obtained as

α% Hˆ γ = V12V23 ∑ β

χα χ β χ β χγ Eβ − Eα

(9)

Here, V23 represents the electronic CT-TT couplings. The super-exchange rate of the SF process starting from the single initial doorway state |α is now given by

ksup er =

2π h

∑γ

2

α% Hˆ γ δ ( Eα − Eγ )

(10)

This equation shows that the rate of electronic transition between the initial and final doorway states is governed by several factors—i.e., the electronic S1-CT and CT-TT couplings, the overlap between different vibrational levels, the |α |β energy gap and the possible resonance condition between the initial and final doorway states ( ). The current ML-MCTDH studies showed that the ultrafast SF dynamics are mainly governed by the vibrational modes within the 0.165–0.20 eV frequency domain, particularly by a few modes (with frequencies 0.1895 eV, 0.193 eV, 0.1965 eV, and 0.2 eV) displaying vibronic couplings with the TT state. These observations can be easily explained by a close examination of Eq. (10). We noticed that the frequencies of these modes are very close to the electronic energy difference between the S1 and TT states. Thus, after the inclusion of the one-quanta vibrational excitation of those modes explicitly, it is highly possible that the possible resonance condition between the initial and final doorway states is nearly satisfied (, , ≈ 0). However, great attention should be paid when drawing such a simple institutive conclusion before the overlaps among different vibrational states are examined carefully. 22


Page 23 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In the current model Hamiltonian, the overall reorganization energy is very small, implying rather weak vibronic couplings for most modes except the modes with very low frequencies. Thus, along most normal coordinates, the shift from the ground-state minimum to the excited-state minimum should be very small. This indicates that the initial doorway state is very close to the lowest vibrational level of the S1 state; i.e.,

α ≈ S1, 0ν

(except the very low-frequency modes). In addition, the overlap

between the vibrational states with the same quanta (for instance O000 = 0 | 0 0 | 0 ) should be rather large (even close to one), whereas the overlap terms may become much smaller when the vibrational states with different quanta are involved (for instance, O001 = 0 | 0 0 |1

should be much smaller—namely, O000 >> O001 ). A

similar analysis also predicts the minor contribution of other overlap terms, such as

O010 = 0 |1 1| 0

or O011 = 0 |1 1|1 .

When we consider the modes away from the resonance conditions, the

S1, 0ν → TT , 0ν

transition certainly becomes important. When a few modes with

frequencies close to the resonance condition (, , ≈ 0) are included, the

S1, 0ν → TT ,1ν

transition also becomes possible because the large resonance term

may compensate for the small O001 term. In other words, both S1, 0ν → TT , 0ν and S1, 0ν → TT ,1ν

transitions may be involved, whereas their ratio should be

determined by the values of the overlap of different vibronic doorway states, which are controlled by the vibronic coupling strength of each mode. The approximated estimations of the contribution of the overlap factor and resonance condition are found in Figure S5 in SI. Thus far, we have clear shown that the efficient SF dynamics may be achieved by the involvement of particular vibrational modes if such modes not only satisfy the resonance condition , , ≈ 0 but also provide enough overlap of different doorway states. This explains why only a few modes (0.1895 eV, 0.193 eV, 0.1965 eV, and 0.2 eV) with vibronic couplings to the TT state play an essential role in the ultrafast SF dynamics obtained by the ML-MCTDH calculations. In the current applied spectral density, the vibronic coupling strength increases with the lowering of 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the mode frequency in the current model Hamiltonian, resulting in larger values of O001 . Thus, the modes with the largest contribution to the SF dynamics always have a slightly smaller frequency than the exact resonance condition—i.e., ν V13 → ksup er > kdirect E2 − E1 V12V23 ≈ V13 → ksup er ≈ kdirect E2 − E1

(12)

V12V23 < V13 → ksup er < kdirect E2 − E1 Going back to the current model for the pentacene pair, V12V23 ( E2 − E1 ) −1 is approximately 0.025 eV, whereas the S1-TT direct coupling V13 given by previous studies is approximately 0.005 eV.63 Thus, within the current model, the super-exchange mechanism seems to lead to a much faster SF rate and become dominant here, which also explains the ultrafast dynamics (80–200 fs) of the pentacene cluster in experimental studies10, 20, 67. Zimmerman et al. once noted that S1 and TT may cross each other by alternating the distance of two pentacenes.5, 63-64 Normally, such a crossing point may not be accessed rapidly by the intermolecular vibration on the excited states because such motion should be very slow, at least slower than the ultrafast nonadiabatic SF dynamics in experimental observations. However, the current conclusion is based on the fact that the SF dynamics occur within the extremely ultrafast time domain. If the time scale of the SF dynamics becomes comparable to the period of the intermolecular vibrations, the distance alternation may drive the system assessing the S1-TT crossing region, and the direct mechanism may become possible. In other realistic systems, two mechanisms may compete with each other. Several factors may potentially influence the alternation of the reaction pathways. For example, if the S1-CT energy gap increases dramatically, the CT state may not be

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

involved. Additionally, as discussed in the last paragraph, the SF dynamics may be controlled by the stacked statues5,

50, 53, 60, 63, 70, 75-76

, which strongly modify the

electronic couplings and resonance condition. If the S1-TT degeneracy becomes easily accessible, the ultrafast SF may then also be possible via the direct mechanism.

3.3 Discussion Although the current calculations employed a rather approximated model, some valuable comments can still be derived for the understanding of the SF dynamics of pentacene pairs and other systems, at least at the qualitative level. We clearly showed that the super-exchange mechanism may result in the ultrafast SF dynamics. In addition, a few modes with frequency resonance and the S1-TT energy gap play essential roles in inducing the large population transfer. Such resonance effects were also noticed by previous works.54, 60 This effect can be well understood by Fermi’s golden rule9,

78, 112

. Although most other modes seem not to result in efficient

population transfer, they rush out the electronic coherence significantly. For the modes with very low or very high frequencies, they are basically bystanders.54 Possibly because only low-frequency modes can be excited at room temperature, the SF dynamics are almost independent of the temperature. However, when the SF dynamics occur within a much longer time-scale or the resonance condition is satisfied, the contribution of the low-frequency modes becomes important60 and the temperature effect may become pronounced53. As a summary, the SF dynamics may be modified by several factors, such as vibrational relaxation, influence of the static and dynamical disorder, and the energy levels of involved electronic states and their couplings. We expect that the discussion of these issues may arouse more interesting debates. The spectral density in the current calculations includes the large contribution from the high-frequency domain. At first glance, this is slightly counterintuitive because the highest frequency for molecular vibrations is generally below 0.5 eV. This counterintuitive issue was clearly addressed in the previous work by Reichman and coworkers48 . The nonlinear effects should appear in the realistic systems, and only the bilinear system–bath coupling is considered within the current model Hamiltonian. 26


Page 26 of 41

Page 27 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


If the spectral density is simply obtained by the Fourier transformation of the autocorrelation function of the bath operator, the high-order term results in the contribution in the high-frequency domains. For practical reasons, we can consider these modes or neglect them completely because they do not play any role in the current SF dynamics (see Figure S1 in SI). One idea to expand the current ML-MCTDH study is to employ more advanced or realistic models constructed from the electronic structure calculations. We noted that in recent years, extensive efforts were made to compute the parameters in the diabatic Hamiltonian for SF9-10, 42-48, 51-53 . Along this research line, we may obtain more precise energy levels of the involved electronic states. After the evaluation of electronic energies, the vibronic couplings may be determined by the calculations of the excited-state energy surface along different normal mode coordinates. A similar idea was widely used in the study of excited-state excitation energy transfer and electron transfer80, 83-84, 96-97, 109-111. Following this procedure, we may construct the diabatic Hamiltonian for more realistic systems. The employment of ML-MCTDH to study these models is also possible, and more electronic states and more vibrational degrees of freedom should be considered. The setup of the ML-MCTDH tree to achieve efficient convergence then becomes even more challenging. We noted that the dynamical calculations based on current model give similar results no matter which dynamics theory is employed. These dynamical theories include ML-MCTDH, Redfield, HEOM and SQC dynamics47-48, 54-57, 126. This strongly indicates that the system–bath coupling in the current model is very weak. When strong system–bath interaction is involved, the perturbative Redfield approach may not work properly78. In this situation, the more precise approaches, HEOM47, 54 or ML-MCTDH126, become essential for the correct treatment. In addition, we also noticed that the super-exchange mechanism also exists in the ultrafast ET dynamics, which were also widely studied by various dynamics methods—for instance, real-time path integral123-124 , which is also supposed to be a rather precise method. It is also very challenging to employ such an approach to study SF dynamics. Thus, the combination of different theoretical approaches and the comparison of their outputs are still critical because such works help us not only understand more physical insight of the SF process but also examine the performance of different dynamical methods. 27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The current work employed the diabatic model to study the SF processes. Krylov and co-workers once argued that the diabatic model may not provide the precise description of all involved electronic states and their couplings65. Thus, they recommended the use of an adiabatic rather than diabatic representation. However, this idea brings more difficulties for dynamics studies. Many dynamical methods are not applicable in adiabatic representation127, and only a few trajectory-based methods become suitable128. However, the employment of these trajectory-based methods to describe the SF dynamics also faces some problems. For instance, one of widely used trajectory-based methods, the surface-hopping approach, suffers from many deficiencies129-133, such as over-coherence. Other advanced methods, such as the quantum-classical Liouville equation128, 137

value representation

134-136

, semiclassical dynamics with initial

and the ab Initio multiply spawning method138, may require

huge numerical efforts in the treatment of the SF dynamics. The recently developed SQC method seems to provide a promising way to solve the nonadiabatic dynamics of large systems

55-57, 139

, whereas more test calculations in more general problems still

must be performed. Thus, the accurate treatment of the SF process, particularly in the adiabatic representation, is still beyond the current computational ability. However, it is still possible to construct a reasonable diabatic Hamiltonian if the suitable quasidiabatic process is taken. Thus, an understanding on the SF dynamics is still possible based on the approximated diabatic model.

4. CONCLUSIONS This paper represents our first effort to employ the ML-MCTDH method to investigate the SF dynamics in polyacenes. We used the electron-phonon coupling Hamiltonian constructed in several previous works, which describes the SF process of pentacene systems. The Hamiltonian includes three electronic states (S1, CT and TT) and a phonon bath. The continuous Debye-type spectral density used by previous work was taken to represent the bath degrees of freedom, which is re-discretized to construct the finite number of bath vibrational modes. The ML-MCTDH method is then employed to study the quantum evolution of the SF dynamics. Instead of placing all modes in the dynamics calculations using the brute force 28


Page 28 of 41

Page 29 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


method, we first investigated the influence of bath modes with different frequencies on the SF dynamics in detail. The results show that the ultrafast SF dynamics are mainly relevant to the vibrational modes with the proper frequency resonance with the S1-TT energy gap, whereas other modes provide minor effects, such as inducing visible electronic decoherence. Previous studies have already pointed out that the mode with frequency resonance with energy gap basically drives the efficient SF dynamics.54, 60 Here we provide a more detailed study on such topic. Our dynamical study also supports the theory that the ultrafast SF dynamics of pentacene are possible via the CT-mediated mechanism, consistent with previous theoretical studies.47, 54, 57 Through a quantitative analysis in the basis of Fermi’s golden rule9, 78, 112, we clarified that the resonance condition achieved by including vibrational levels may also be critical for efficient SF dynamics. To satisfy this condition, both vibrational wavefunction overlap and energy gap should be considered. When the frequency of a mode is resonant with the S1-TT energy gap, the larger overlap between vibronic doorway levels significantly improves the SF rate. If we treat electron-phonon coupling in the way of continuum spectral density, an efficient SF is achieved in the situation in which the large amplitude of the spectral density coincides with the S1-TT energy gap. Similar observations were made in previous studies. This analysis provides a useful idea to conduct the computational design of the novel photovoltaic compounds with high SF efficiency. The results of the current ML-MCTDH study are very similar to those of previous works obtained by the Redfield, HEOM and SQC methods10, 47-48, 54-57. This indicates that the vibronic couplings in the current model Hamiltonian are very weak. In a general molecular system conducting the efficient SF, it is unclear whether the vibronic couplings are weak. This issue can be clarified only by the construction of the realistic vibronic coupling Hamiltonian from the ab initio calculations. We notice that extensive efforts have been recently made to build the Hamiltonian for the SF processes, and some of them have already tried to access the accurate vibronic couplings of excited states. It should be highly interesting to employ different dynamical methods to re-examine more accurate models for a semi-quantitative or quantitative study of the SF dynamics in other realistic systems. This represents a great theoretical challenge for the future. 29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 41

APPENDIX I. Temperature effect The temperature effect on the current ET process was assessed in a rather approximated way. Various initial conditions were created by the Monte-Carlo sampling of the vibration levels for each normal mode according to the Boltzmann distribution at a certain temperature:

pi =

exp ( − Ei / kT )

∑ exp ( − E

j

/ kT )

(13)

j

where Pi is the probability, Ei and Ej are the energy levels for the particular vibrational mode, k is the Boltzmann constant, and T is the temperature. The quantum dynamics were then performed for each initial condition. The ensemble average of the quantum evolution over all initial conditions finally gives the quantum SF dynamics at a certain temperature. As shown in Figure 7, the SF dynamics at 300 K were compared to the previous results at 0 K. Although the overall SF dynamics may be slightly sped up by raising the temperature, the minor difference indicates that the temperature effect on the current ET process is rather weak. Similar conclusions were also drawn in previous experimental and theoretical studies.97, 140

Figure 7. Effect of temperature: the SF dynamics at 300 K were compared to the results at 0 K. 30


Page 31 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


II. Pure electronic dynamics The pure electronic dynamics governed only by the system Hamiltonian without inclusion of the bath degrees of freedom is investigated. As expected, the population dynamics in such a case displays the typical Rabi-type oscillation, as shown in Figure 8. Owing to the presence of the S1-CT couplings, the obvious population interchange occurs between the S1 and CT states. The oscillation period of the S1 population is approximately 28.3 fs, and the amplitude in the population oscillation is approximately 0.4. Although visible couplings exist between the CT and TT states, the population of the TT state remains very low (less than 0.027) during the dynamical evolution.

Figure 8.Pure electronic dynamics III. The optimal construction of the ML-MCTDH tree The method for constructing a reasonable ML-MCTDH expansion is critical for the computational efficiency and numerical convergence in the dynamical study of complex systems. Previous studies provided some useful principles for this purpose78, 111

, but it is still challenging to construct a detailed and practical procedure that

provides a step-by-step tutorial to guide the construction of the rational ML-MCTDH tree for general problems. Thus, we wish to provide some personal suggestions on this topic based on our empirical observations. Although these views are mainly based on our employment of ML-MCTDH calculations in the Heidelberg package,89 we believe that a few of principles may be valid for broad applications that are more general. 31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

When many degrees of freedom are involved in the ML-MCTDH calculations, it is not wise to simply put all of them together and build a whole ML-MCTDH tree with physical intuition. Such a construction may easily result in an arbitrary tree with very low computational efficiency and convergent difficulty. It is wiser to divide all degrees of freedom into different groups and perform several initial test calculations. For vibronic coupling problems, this method also provided us with a rough idea of the roles of vibrational modes in the nonadiabatic dynamics. In other words, the importance of each degree of freedom was determined by the test calculations, and the corresponding setups were taken as the reasonable initial guess for the large ML-MCTDH calculations. It should be more effective to try the above protocol iteratively and increase the expansion size of the ML-MCTDH tree step-by-step by including increasingly more modes. Sometimes, it is even necessary to build different ML-MCTDH trees to check the final convergence and computational efficiency. At each layer, similar basis numbers should be used for different branches to achieve the “balance” of the whole ML-MCTDH tree. According to this rule, it should be better to put several “inactive” (unimportant) modes together to form a branch. In contrast, less “active” (important) modes, even a single one, should be grouped to define a branch. The reason is obvious: more bases are required for each important mode. The mode combination in the bottom layer may potentially impact the computational cost. Similar to the above rule, more inactive modes should be combined even in the bottom layer, whereas less active modes should be combined or each active mode should be treated separately. For systems with rather few degrees of freedom (

Ultrafast Nonadiabatic Dynamics of Singlet Fission: Quantum

Recommend Documents