Dynamics of Coupled ElectronâBoson Systems with the Multiple

Subscriber access provided by UNIV OF MISSOURI ST LOUIS

Article

Dynamics of Coupled Electron-Boson Systems With the Multiple Davydov D Ansatz and the Generalized Coherent State 1

Lipeng Chen, Raffaele Borrelli, and Yang Zhao J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b07069 • Publication Date (Web): 24 Oct 2017 Downloaded from http://pubs.acs.org on October 30, 2017

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 A 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 38

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

Dynamics of coupled electron-boson systems with the multiple Davydov D1 Ansatz and the generalized coherent state Lipeng Chen,† Raffaele Borrelli,‡ and Yang Zhao∗,† †Division of Materials Science, Nanyang Technological University, Singapore 639798, Singapore ‡Department of Agricultural, Forestry and Food Science, Universit` a di Torino, I-10095 Grugliasco, TO, Italy E-mail: [email protected] Phone: +65 6513 7990. Fax: +65 6790 9081 Abstract Dynamics of a coupled electron-boson system is investigated by employing a multitude of the Davydov D1 trial states, also known as the multi-D1 Ansatz, and a second trial state based on a superposition of the time-dependent generalized coherent state (GCS Ansatz ). The two Ans¨ atze are applied to study population dynamics in the spin-boson model and the Holstein molecular crystal model, and a detailed comparison with numerically exact results obtained by the (multi-layer) multiconfiguration timedependent Hartree method and the hierarchy equations of motion approach is drawn. It is found that the two methodologies proposed here have significantly improved over that with the single D1 Ansatz, yielding quantitatively accurate results even in the critical cases of large energy biases and large transfer integrals. The two methodologies provide new effective tools for accurate, efficient simulation of many-body quantum

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

dynamics thanks to a relatively small number of parameters which characterize the electron-nuclear wave functions. The wave function-based approaches are capable to track explicitly detailed bosonic dynamics, which is absent by construct in approaches based on the reduced density matrix. The efficiency and flexibility of our methods is also one of the advantages as compared with numerically exact approaches such as QUAPI and HEOM, especially at low temperatures and in the strong coupling regime.

Introduction Understanding the intricate interplay between electronic and vibrational degrees of freedom (DOFs) is of particular importance for deciphering eletronic dephasing and energy relaxation processes in organic solar cells 1–4 . However, fully quantum mechanical calculations of those dynamic processes remain computationally prohibitive due to a huge number of DOFs. A common technique alleviating the computational bottleneck is to provide a reduced description of the system by tracing out the bath DOFs. The bath impact on the system is then encoded into the correlation function of the collective energy gap coordinate, which can be exactly solved by the cumulant expansion for harmonic baths 5 . In the language of the theory of open quantum systems, this is often referred to as Quantum Master Equations (QMEs) approach 6,7 . The hierarchical equations of motion (HEOM) 8–10 and the quasi-adiabatic path integral (QUAPI) 11 are among the most successful QMEs approaches that propagate the reduced density matrix in a numerically exact manner. However, both methods are computationally expensive in the strong electron-phonon coupling regime and at low temperatures due to the non-Markovian and the non-perturbative nature. Furthermore, the application of HEOM is currently restricted to certain forms of phonon spectral densities, albeit several extensions have been proposed to circumvent this difficulty 12–15 . A large number of efficient, approximate QMEs methods based on a polaron transformation have been formulated to tackle both incoherent and coherent dynamics 16–20 . In the aforementioned QMEs methods, all information on bosonic dynamics is implicitly 2


Page 2 of 38

Page 3 of 38

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


integrated into system observables such as electronic populations. However, recent ultrafast spectroscopic experiments have revealed the potential role of discrete high frequency intramolecular vibrational modes on long-lived quantum beating in 2D electronic spectra of light harvesting complexes 21–23 . It is also suggested that high frequency intramolecular vibrational modes may facilitate efficient singlet fission processes in organic molecular crystals 24–26 . In this regard, wave function propagation methods treating electronic and vibrational DOFs on an equal footing are more promising alternatives for describing finer details of coupled electronic-vibrational dynamics. The multiconfiguration time-dependent Hartree (MCTDH) method, and its multilayer variant (ML-MCTDH), 27–30 are typical examples of wave function propagation methods, which provide powerful numerical tools for accurate simulations of quantum dynamics with many DOFs at very low temperatures. Recently, the time-dependent Davydov Ans¨ atze 31,32 and their multistate extensions 33–35 (i.e., superpositions of the Davydov Ans¨ atze) have been shown to provide highly accurate, yet efficient descriptions for dynamical properties of coupled electronic-vibrational systems with many DOFs. By employing the Dirac-Frenkel time-dependent variational principle, these Ans¨ atze have been widely used to study the dynamical behavior of the Holstein polaron and the spin-boson model (SBM) 31–38 , as well as energy transfer processes in various light harvesting complexes 39–43 . It is noted that the bosonic wave functions of the Davydov Ansatz are described by coherent states, i.e., Gaussian wave packets, or by a linear combination of them, which are also the basis of many other quantum dynamics descriptions of many-body systems 44–51 . A natural extension of the Gaussian wave packets representation of the bosonic dynamics is the use of the generalized Gauss-Hermite functions 52,53 , also known as generalized coherent states in the quantum optics literature 54–56 . Borrelli, et al. recently proposed a new Ansatz based on the superposition of time-dependent Gauss-Hermite wave packets to handle the dynamics of several electronic states coupled to a bath of harmonic oscillators 57,58 . The favorable scaling property of the method renders it a promising numerical tool for the simulation of dynamical processes in condensed phases.

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 38

We first compare in this work the theoretical aspects of the multiple Davydov D1 Ansatz and the generalized coherent state, and then show how accurate simulations of quantum dynamics of coupled electronic and vibrational systems can be achieved by these trial states. The rest of the paper is organized as follows. The multiple Davydov D1 Ansatz and the generalized coherent state are introduced in Sec. II. Population dynamics is calculated for two representative models, i.e., the SBM and the Holstein molecular crystal model, using the Davydov D1 Ansatz and the generalized coherent state, and results are presented in Sec. III. Conclusions are drawn in Sec. IV.

METHODOLOGY The multiple Davydov D1 Ansatz ˆ the time evolution of the wave function Ψ(t) satisfies the For a general Hamiltonian H, Schrödinger equation i

∂ ˆ |Ψ(t)i = H|Ψ(t)i. ∂t

(1)

In this work, a multitude of the Davydov D1 trial states, also known as the multi-D1 Ansatz, will be employed as an approximate solution to the above Schrödinger equation. As an improved Ansatz to the single D1 Ansatz, the multi-D1 Ansatz can be constructed as follows

|DM 1 (t)i

=

N M X X i

Ain (t)|ni exp

n

X q

[finq (t)ˆb†q − H.c.] |0iph

(2)

Here we consider a model Hamiltonian with several electronic states coupled to a bath of harmonic oscillators. |ni specifies electronic states (n=1,2,· · · ,N), ˆb†q (ˆbq ) is the boson creation (annihilation) operator with the momentum q, and |0iph is the vacuum state of the boson bath. Ain (t) and finq (t) are the amplitudes of the electronic state |ni and phonon displacements of the boson mode q, respectively, and i denotes the coherent superposition state. When M = 1, the multi-D1 Ansatz reduces to the standard Davydov D1 Ansatz. It is 4


Page 5 of 38

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


readily seen that in the multi-D1 Ansatz the wave function is parametrised using M N (1 + d) complex numbers, with d denoting the number of the bath DOFs. The equations of motion for the variational parameters Ain (t) and finq (t) can be obtained by the Dirac-Frenkel time-dependent variational principle 59 . For this purpose, we first formulate the Lagrangian associated with the multi-D1 Ansatz " # → − ← − i~ ∂ ∂ L = hD1M (t)| |D1M (t)i − hD1M (t)| |D1M (t)i 2 ∂t ∂t ˆ M (t)i −hD1M (t)|H|D 1

(3)

The Dirac-Frenkel time-dependent variational principle then yields the equations of motion for the variational parameters Ain (t) and finq (t) ∂L d ∂L ( =0 ∗) − dt ∂ A˙in ∂A∗in ∂L d ∂L ( =0 ∗) − ∗ ˙ dt ∂ finq ∂finq

(4)

Detailed derivations of the equations of the motion for the variational parameters can be found in Refs. 34,35 . Expectation values of any physical observables of interest can be then calculated in terms of the multi-D1 Ansatz as ˆ ˆ M hO(t)i = hDM 1 (t)|O|D1 (t)i

(5)

ˆ is a time independent operator. where O

The generalized coherent state Ansatz Generalized coherent states (GCSs) were first proposed by Senitzky 60 and then extended by Nieto 54 to the generalized squeezed coherent state. The GCS of a harmonic oscillator is

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 38

defined as 54,55,61,62 ˆ |k, αi = D(α)|ki

(6)

ˆ where D(α) is the displacement operator given by ˆ D(α) = exp(αˆ a† − a∗ a ˆ)

(7)

Here a ˆ and a ˆ† are boson annihilation and creation operators, and |ki is the eigenstate of a harmonic oscillator. The k = 0 case reduces to the celebrated coherent state 63 Without loss of generality, we can write the wave function for a physical system of N interacting electronic states coupled to d vibrational DOFs as a direct product of electronic and vibrational states 57,58

|Ψ(t)i =

X

ClK (t)|K, αl (t)i|li

(8)

lK

where |li represents the lth electronic state (l=1,2,· · · ,N), K is a multi-index K = (k1 , k2 , · · · , kd ) ∈ N0d , and |K, αl (t)i are a complete orthonormal set for the boson DOFs which can be expressed as the product of d time-dependent GCSs

|K, αl (t)i = |k1 , · · · , kd , αl1 (t), · · · , αld (t)i =

d Y

|ki , αli (t)i

(9)

i=1

It should be noted that the size of the bosonic Hilbert space increases exponentially with d and thus needs to be truncated for practical calculations. The truncation scheme employed here is based on partitioning the entire Hilbert space into a set of subspaces where only c vibrational modes can be excited simultaneously (c=0,1,· · · ,d). Since there exists dc

distinct combinations with c excited modes, the wave function of Eq. 8 can be reformulated

6


Page 7 of 38

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


as

Ψ(t) =

" d (dc) N X XX

wi1 ···wic

X

Clki1 ···kic (t)|ki1 , αli1 (t)i

c=0 i1 ···ic ki1 ···kic =1

l=1

#

· · · |kic , αlic (t)i |li =

N X l=1

"

Cl (t)|αl (t)i +

′

Clki (t)|ki , αli (t)i|αl (t)i

i=1 ki =1

() X wj wi X X d 2

+

wi d X X

′′

Clki kj (t)|ki , αli (t)i|kj , αlj (t)i|αl (t)i

i,j=1 ki =1 kj =1

+

#

· · · |li

(10)

′

′′

where αl is partitioned as αl = (αli , αl ) in the case of 1 excited mode, and αl = (αli , αlj , αl ) in the case of 2 excited modes. For simplicity, all indices of expansion coefficients with zero value are omitted, and wi is the maximum quantum number for the ith excited mode. This kind of partitioning reduces the size of the bosonic Hilbert space dramatically, and has been applied to simulate molecular electronic lineshapes 64 and electron superexchange in molecular chains 65 . It should be noted that the partition scheme used in Eq. 10 leads to wave functions similar to the configuration interaction (CI) expansion in electronic structure calculations. In the CI method, the wave function is expanded as a linear combination of configuration state functions (CSFs) constructed from spin orbitals (SO)

Ψ=

X

SO SO SO pI ΦSO I = p 0 Φ0 + p 1 Φ1 + p 2 Φ2 + · · ·

(11)

I=0

SO where the first term (p0 ΦSO 0 ) is normally the Hartree-Fock determinant, and pI ΦI (I =

1, 2, · · · ) denote CSFs where I spin orbitals that are swapped with virtual orbitals from the Hartree-Fock determinant. For this truncation scheme, we can test numerical convergence by varying c until there

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 38

are no significant changes for the expectation value of physical observables of interest. Approximations with c up to 2 are most promising for studying quantum dynamics in systems with many DOFs. In the following, we refer to the Ansatz (Eq. 10) with excitation levels c = 1, 2 as ΨGCS1 (t) and ΨGCS2 (t), respectively. In the GCS1 Ansatz, the wave function is parametrised using N (1 + d + dw) complex numbers, whereas in the GCS2 Ansatz the number of parameters increases to N (1 + d + dw +

d(d−1) 2 w ), 2

w being the size of the basis

set. It should be also mentioned that the Ansatz (Eq. 10) with c = 0 corresponds to the single Davydov D1 Ansatz (Eq. 2 with M = 1). A further flexibility of the basis set is to have different values of the excitation level c for different electronic states (see infra). The equations of motion for the parameters ClK (t) and αli (t) can be derived by applying the standard Dirac-Frenkel variational principle (see Eq. 4) with the Lagrangian " # → − ← − ∂ i~ ∂ hΨ(t)| |Ψ(t)i − hΨ(t)| |Ψ(t)i L = 2 ∂t ∂t ˆ −hΨ(t)|H|Ψ(t)i

(12)

The readers are referred to Refs. 57,58 for a detailed derivation of the equations of the motion for the variational parameters.

Results and Discussions In this section, we present results of population dynamics based on the multi-D1 Ansatz and on the GCS Ansatz, and compare them with those from the numerically exact MCTDH (or ML-MCTDH) and HEOM approaches. In Subsection IIIA, we focus on spin dynamics of the SBM. In Subsection IIIB, we treat dynamics properties of a one-dimensional Holstein polaron.

8


Page 9 of 38

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


Spin boson model As an archetype model for studying environment-induced dissipation, the SBM 66,67 finds wide spread applications in condensed phase physics and physical chemistry processes, such as quantum entanglement and computation 68–70 , quantum phase transition 66,71,72 and excitation energy transfer in biological systems 73 . The SBM describes a two level system linearly coupled to an ensemble of harmonic oscillators, its Hamiltonian can be written as

H = ǫσz + ∆σx +

X ωj j

2

(p2j + qj2 ) + σz

X

g j qj

(13)

j

where σx and σz are the Pauli matrices, ǫ and ∆ are the spin bias and the tunnelling constant, respectively, and pj , qj represent the dimensionless momentum and coordinate of the jth boson mode, gj is the coupling strength between the spin system and the jth boson mode. The coupling between the spin system and the harmonic bath is captured by the spectral density function J(ω) =

πX 2 g δ(ω − ωj ) 2 j j

(14)

In this paper, we use an Ohmic spectral density

J(ω) =

π αω exp(−ω/ωc ) 2

(15)

where ωc is the cutoff frequency, and α is the so-called Kondo parameter characterizing the system-bath coupling strength. For practical calculations, a discretization procedure of the spectral density is required in order to obtain gj . Following Refs. 74,75 ,the spectral density is divided into Nb modes in the frequency domain [0, ωmax ], where ωmax = 8ωc is the frequency upper bound. The frequency of the jth mode is then given by ω0 ωj = −ωc ln 1 − j ωc

9

!


(16)


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 38

and its linear coupling strength is

gj =

p 2αωc ωj ω0

(17)

with ω0 = ωc (1 − e−ωmax /ωc )/Nb . In the discussion to follow, we will present results for cases with an effective number of bath modes Nb = 15 and 60. We restrict our discussion to a relatively short time (up to 500 fs) as numerical errors are expected to accumulate at long times. The long-time limit (equilibrium) of the system dynamics, which is beyond the scope of this paper, can be calculated using imaginary time path integrals 76–78 or the non-Markovian dynamical maps approach 79 . In the non-Markovian dynamical maps approach, information on dynamical correlations can be retrieved from the initial trajectories generated experimentally or numerically to formulate a collection of non-Markovian dynamical maps, which are then transformed to a set of non-Markovian transfer tensors. Under the assumption of time-translational invariance, these transfer tensors can be used to propagate system dynamics to arbitrarily long time scales 79 . In the electron transfer theory 7,80 , the electron transfer reactions can be classified as adiabatic (∆/ωc > 1) or nonadiabatic (∆/ωc < 1), where ∆ and ωc are electronic coupling and characteristic frequency of the bath, respectively. The characteristic timescale of the bath is slow (fast) compared to that of the electronic tunneling in the adiabatic (nonadiabatic) reactions. We will check the extent of applicability of the GCS and multi-D1 Ansätze in the nonadiabatic, intermediate and adiabatic parameter regimes using benchmarking by the numerically exact MCTDH and ML-MCTDH methdologies 27–30 . All MCTDH (or MLMCTDH) results are obtained using the MCTDH Heidelberg package version 8.5 81 . We first consider the nonadiabatic regime. Fig. 1 displays the time evolution of the down state population as calculated by the GCS Ansatz (left panel) and the multi-D1 Ansatz (right panel) for three values of the bias ǫ = 0 (Figs. 1(a,d)), ǫ = ∆ (Figs. 1(b,e)), and ǫ = 2∆ (Figs. 1(c,f)). The other model parameters are the tunnelling constant ∆ = 100cm−1 , the

10


Page 11 of 38

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


cutoff frequency ωc = 10∆, the Kondo constant α = 0.3, and the effective number of bath modes Nb = 15. The spin is assumed to initially reside in the up state, while the factorized initial condition corresponding to the phonon vacuum state is applied to the boson DOFs. As shown in Fig. 1(a), the GCS Ansatz with excitation level c = 0 (the single Davydov D1 Ansatz ) already captures the oscillatory features of population dynamics in the absence of bias (ǫ = 0), and the exact results are well reproduced with the GCS Ansatz with c = 2. It is found that P (t) obtained by the 2-D1 Ansatz are already in agreement with those by the MCTDH method (see Fig. 1(d)), and further increasing multiplicity M does not change the dynamics in the time regime considered. For ǫ = ∆, the GCS Ansatz with c = 0 completely fails to depict the population dynamics, while the GCS Ansatz with c = 2 describes the correct dynamical behaviors. Similar to Fig. 1(d), convergent results have been achieved with 4-D1 Ansatz for the case of ǫ = ∆ (see Fig. 1(e)). Upon increasing the bias to ǫ = 2∆, a challenging case for most numerical methods 82,83 , one can readily observe that the GCS Ansatz with c = 2 predicts only qualitatively the overall beating behavior of the population dynamics. The amplitudes of these beatings are found to be about 5% lower than those from MCTDH method. On the other hand, the 7-D1 Ansatz yields the amplitudes as well as the frequencies of these beatings in perfect agreement with those predicted by MCTDH method (see Fig. 1(f)). For α = 0.5 the correct time evolution of population dynamics under zero bias and ǫ = ∆ is nicely predicted by the GCS2, whereas the performance of the GCS Ansatz deteriorates for the case of a large energy bias ǫ = 2∆, as illustrated in Fig. 2. It is interesting to note that the accurate time evolution of population dynamics for all three values of ǫ, i.e. ǫ = 0, ∆, and 2∆, are obtained by the multi-D1 Ansatz with multiplicity M = 6, 6, and 7, respectively. The time evolution of population dynamics for the strong coupling case of α = 1.0 is depicted in Fig. 3. As compared with results in Figs. 1 and 2, the overall performance of the GCS Ansatz worsens for α = 1.0 due to the fact that the multiphonon excitation with c > 2 should have been included in the Ansatz in order to correctly characterize the population

11



1

1 (a)

(d)

0.6

0.6

P(t)

0.8

P(t)

0.8

0.4

0.4

0

2D1

c=0 c=1 c=2 MCTDH

0.2

0

100

200

300

400

4D1 6D1

0.2

0

500

MCTDH 0

100

200

t(fs) 1

300

400

500

t(fs) 1

(b)

(e)

0.6

0.6

P(t)

0.8

P(t)

0.8

0.4

0.4

0.2

0

2D1

c=0 c=1 c=2 MCTDH 0

100

200

300

400

4D1 6D1

0.2

MCTDH 0

500

0

100

200

t(fs)

300

400

500

t(fs)

1

1 (c)

(f)

0.8

0.6

0.6

P(t)

0.8

P(t)

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 12 of 38

0.4

0.4

0.2

0

2D1

c=0 c=1 c=2 MCTDH 0

100

200

300

400

6D1 7D

0.2

1

0

500

MCTDH 0

100

t(fs)

200

300

400

500

t(fs)

Figure 1: The time evolution of the down state population calculated by the GCS Ansatz (left panel) and the multi-D1 Ansatz (right panel). The parameters of the model are ∆ = 100cm−1 , ωc = 10∆, α = 0.3; (a,d) ǫ = 0, (b,e) ǫ = ∆, (c,f) ǫ = 2∆. The number of bath modes is Nb = 15.

12


Page 13 of 38

1

1

(a)

(d)

0.6

0.6

P(t)

0.8

P(t)

0.8

0.4

0.4

0

2D1

c=0 c=1 c=2 MCTDH

0.2

0

100

200

300

400

4D1 6D1

0.2

0

500

0

100

200

t(fs) 1

300

MCTDH 400 500

t(fs) 1

(b)

(e)

0.6

0.6

P(t)

0.8

P(t)

0.8

0.4

0.4 2D1 c=0 c=1 c=2 MCTDH

0.2

0

0

100

200

300

400

4D1

0.2

6D1 MCTDH

0

500

0

100

200

t(fs)

300

400

500

t(fs)

1

1 (c)

(f)

0.8

0.6

0.6

P(t)

0.8

P(t)

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


0.4

0.4

0.2

0

2D1

c=0 c=1 c=2 MCTDH 0

100

200

300

400

6D 1 7D

0.2

1

0

500

MCTDH 0

100

t(fs)

200

300

400

500

t(fs)


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

dynamics for the large bias case of ǫ = 2∆. It is found that a slightly larger multiplicity of M=10 for the multi-D1 Ansatz is needed to capture the accurate time evolution of population dynamics (see the right panel of Fig. 3). The multi-D1 Ansatz is thus valid for a wide range of coupling strengths and energy biases. We then consider the intermediate regime with ∆/ωc = 1 as shown in Fig. 4. For zero bias, the GCS1 Ansatz already yields results in good agreement with those of the MCTDH method. Increasing the bias ǫ to ∆, the GCS2 Ansatz starts to deviate from the MCTDH results after 100fs. In particular, its long-time value of about 0.7 is much lower than the exact value of 0.9. For even a larger bias ǫ = 2∆, the GCS2 Ansatz fails to describe the population dynamics, predicting a steady state value of about 0.19 after 200fs, whereas the MCTDH results continue to increase with time. On the other hand, the multi-D1 Ansatz reproduces the accurate population dynamics quite well, albeit there is a slightly deviation from the exact results at long times for the case of ǫ = 2∆. We finally consider the adiabatic regime with ∆/ωc = 4 as displayed in Fig. 5. The population dynamics exhibits weakly damped coherent oscillations because the reorganization energy, Er = 2αωc = 50cm−1 , is smaller than the electronic coupling ∆ = 100cm−1 . For zero bias, the GCS1 Ansatz reproduces the exact population dynamics. The GCS results show undamped phase-shifted population dynamics as compared to the exact results for the bias cases (see Fig. 5(b,c)). For all three biases, the multi-D1 results are in perfect agreement with the numerically exact results of the MCTDH method. Figures 6 and 7 show the population dynamics for the model with a larger number of modes Nb = 60. Two values of Kondo parameter (α = 0.3 and 0.5) and an energy bias ǫ = ∆ and an cutoff frequency ωc = 10∆ are considered respectively. In Fig. 6(a), the population dynamics obtained by the GCS Ansatz with c1 = c2 = 0 exhibit oscillatory patterns with an average value close to 0.5, an artefact induced by an insufficient number of excited states as the system should be in a localized state in the presence of a finite bias. With increasing the excitation level to c1 = c2 = 1, the steady

14


Page 14 of 38

Page 15 of 38

0.5

0.5

(a)

0.3

0.3

(d)

P(t)

0.4

P(t)

0.4

0.2

0.2

c=0 c=1 c=2 MCTDH

0.1

0

0

100

200

300

400

6D1

0.1

10D

0

500

MCTDH

0

100

200

300

400

500

t(fs) 0.6

(b)

0.5

0.5

0.4

0.4

P(t)

P(t)

0.6

0.3 0.2

(e)

0.3 0.2 8D

c=0 c=1 c=2 MCTDH

0.1 0

0

100

200

300

400

0.8

1

10D1

0.1

12D

0

500

1

MCTDH

0

100

200

t(fs)

300

400

500

t(fs) 0.8

(c)

0.7

0.6

0.6

0.5

0.5

P(t)

0.7

0.4 0.3

(f)

0.4 0.3

0.2

0

100

200

300

400

6D

0.2

c=0 c=1 c=2 MCTDH

0.1 0

1

12D1

t(fs)

P(t)

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


0

500

1

8D1

0.1

10D1 MCTDH

0

100

t(fs)

200

300

400

500

t(fs)


15



1

1

(a)

0.6

0.6

(d)

P(t)

0.8

P(t)

0.8

0.4

0.4

0

2D1

c=0 c=1 c=2 MCTDH

0.2

0

100

200

300

400

0.2

4D

1

1

6D1

0

500

MCTDH

0

100

200

t(fs)

300

400

500

t(fs) 1

(b)

0.6

0.6

(e)

P(t)

0.8

P(t)

0.8

0.4

0.4 4D1 c=0 c=1 c=2 MCTDH

0.2

0

0

100

200

300

400

0.2

6D1 8D1

0

500

MCTDH

0

100

200

t(fs) 0.8

0.8

(c)

400

500

(f)

P(t)

0.6

0.4

0.2

0

300

t(fs)

0.6

P(t)

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 16 of 38

100

200

300

400

6D1

0.2

c=0 c=1 c=2 MCTDH 0

0.4

8D

1

10D1

0

500

MCTDH

0

100

t(fs)

200

300

400

500

t(fs)

Figure 4: The time evolution of the down state population calculated by the GCS Ansatz (left panel) and the multi-D1 Ansatz (right panel). The parameters of the model are ∆ = 100cm−1 , ωc = ∆, α = 1.0; (a,d) ǫ = 0, (b,e) ǫ = ∆, (c,f) ǫ = 2∆. The number of bath modes is Nb = 15.

16


Page 17 of 38

1

1

(a)

0.6

0.6

(d)

P(t)

0.8

P(t)

0.8

0.4

0.4

0.2

0

0.2

c=0 c=1 c=2 MCTDH

0

100

200

300

400

2D1 0

500

MCTDH 0

100

200

t(fs) 0.6

(b)

0.5

0.5

0.4

0.4

P(t)

P(t)

0.6

0.3 0.2

400

500

(e)

0.3 0.2 2D1

c=0 c=1 c=2 MCTDH

0.1 0

300

t(fs)

0

100

200

300

400

4D1 6D1

0.1 0

500

MCTDH 0

100

200

t(fs) 0.25

300

400

500

t(fs) 0.25

(c)

0.2

0.15

0.15

(f)

P(t)

0.2

P(t)

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


0.1

0.1

0.05

0

0

100

200

300

400

2D1

0.05

c=0 c=1 c=2 MCTDH

4D

1

6D1

0

500

MCTDH

0

100

t(fs)

200

300

400

500

t(fs)

Figure 5: The time evolution of the down state population calculated by the GCS Ansatz (left panel) and the multi-D1 Ansatz (right panel). The parameters of the model are ∆ = 100cm−1 , ωc = ∆/4, α = 1.0; (a,d) ǫ = 0, (b,e) ǫ = ∆, (c,f) ǫ = 2∆. The number of bath modes is Nb = 15.

17



1

1

(a)

0.8

0.6

0.6

(b)

P(t)

0.8

P(t)

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 18 of 38

0.4

0.4 c1=c 2=0

2D1

c1=c 2=1 c1=1,c 2=2

0.2

4D1 5D1

0.2

ML-MCTDH 0

0

100

200

300

400

0

500

ML-MCTDH 0

100

t(fs)

200

300

400

500

t(fs)

Figure 6: The time evolution of the down state population calculated by the GCS Ansatz (left panel) and the multi-D1 Ansatz (right panel). The parameters of the model are ∆ = 100cm−1 , ωc = 10∆, α = 0.3, and ǫ = ∆. The number of bath modes is Nb = 60. state value of the down state population approaches a maximum of about 0.82. When states with two excited vibrations are included in the GCS Ansatz, it is found that the population transfer completely to the down state after 500 fs. The performance of the GCS Ansatz slightly worsens for α = 0.5 (see Fig. 7(a)), and the GCS Ansatz with c1 = 1, c2 = 2 yields a maximum transition probability of 0.85. Comparisons between results obtained from the GCS Ansatz and from the numerically exact ML-MCTDH method clearly shows that the GCS Ansatz predicts a qualitatively correct behavior of the population dynamics. On the other hand, it is shown that the 2-D1 Ansatz significantly improves the description of population dynamics over the single Davydov D1 counterpart, and exact population dynamics can be reproduced by the multi-D1 Ansatz with a multiplicity of 5 (see Figs. 6(b) and 7(b)). The results presented in the left panel of Figs. 1-7 suggest that the GCS Ansatz performs better in case of small energy bias ǫ and system-bath coupling strength α. Increasing ǫ and α, vibrational states with higher excitation levels should be included into the Ansatz in order to characterize the correct population dynamics. On the other hand, the right panels of Figs. 1-7, clearly show that the multi-D1 Ansatz provides excellent population dynamics over a wide range of parameter regimes (ǫ, α, ωc and Nb ). Even for the challenging cases of large energy biases and strong system-bath coupling, one can readily observe that a medium multiplicity can yield accurate population dynamics at long times. Furthermore, our multi18


Page 19 of 38

1

1

(a)

0.8

0.6

0.6

(b)

P(t)

0.8

P(t)

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


0.4

0.4 2D1

c1=c2=0 c1=c 2=1 c1=1,c 2=2

0.2

4D1 5D1

0.2

ML-MCTDH

ML-MCTDH 0

0

100

200

300

400

0

500

0

100

t(fs)

200

300

400

500

t(fs)

Figure 7: The time evolution of the down state population calculated by the GCS Ansatz (left panel) and the multi-D1 Ansatz (right panel). The parameters of the model are ∆ = 100cm−1 , ωc = 10∆, α = 0.5, and ǫ = ∆. The number of bath modes is Nb = 60. D1 Ansatz provides a unified treatment of population dynamics for the adiabatic (slow bath), nonadiabatic (fast bath), as well as the intermediate regime.

Holstein model As one of the most important theoretical models to elucidate the complex interplay between the electron and its surrounding phonons, the so-called electron-phonon cornucopia, the Holstein molecular crystal model has been widely applied to study charge and energy transport in organic solar cells and natural light harvesting systems 4,7 . The Holstein model describes a number of singly excited molecule coupled to intra- and inter-molecular vibrational degrees of freedom. The Hamiltonian of the one-dimensional Holstein molecular crystal model with a periodic boundary condition can be written as ˆ =H ˆ ex + H ˆ ph + H ˆ ex−ph H

19


(18)


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 20 of 38

with ˆ ex = −J H

X

a ˆ†n (ˆ an+1 + a ˆn−1 )

n

ˆ ph = H

X

ωqˆb†qˆbq

q

ˆ ex−ph = −g H

N X

ωq a ˆ†n a ˆn (ˆbq eiqn + ˆb†q e−iqn )

(19)

q,n=1

ˆ ex is the Frenkel-exciton Hamiltonian, and a where H ˆ†n (ˆ an ) is the creation (annihilation) operator of an exciton at the nth site, J is the nearest-neighbor inter-molecular coupling constant. ˆ ph describes the vibrational modes where ˆb† (ˆbq ) is the creation (annihilation) operator of H q ˆ ex−ph is the intra-molecule excitona phonon mode with momentum q and frequency ωq . H phonon coupling in a site-diagonal form, and g is the corresponding coupling strength. N is the number of sites in the Holstein model. For the sake of simplicity, a linear phonon dispersion band ωq = ω0 [1 + W (2|q| − 1)] with q =

2π l(l N

= − N2 + 1, · · · , −1, 0, 1, · · · , N2 ) is

assumed in the calculation. ω0 is the central frequency of the phonon band which is taken as the energy unit, i.e., ω0 = 1, and phonon bandwidth is 2W ω0 . In this subsection, we mainly focus on the time evolution of the exciton population Pex (t, n), which is defined as follows Pex (t, n) = hΦ(t)|ˆ a†n a ˆn |Φ(t)i

(20)

where Φ(t) refers to the multi-D1 Ansatz (Eq. 2) or the GCS Ansatz (Eq. 10). In the simulation, we assume that the system initially has one exciton at site n = 0. The time evolution of the exciton probability Pex (t, n) obtained by the multi-D1 Ansatz and the GCS Ansatz for two values of transfer integral J = 0.5, 1.0 are displayed in Figs. 8 and 9. As a comparison, we also include results calculated by the numerically exact HEOM method for each case 14 . As shown in Figs. 8(a,d), quite obvious difference is found in the exciton propagation calculated by the single D1 Ansatz and the HEOM method when 20


Page 21 of 38

6

1 c) GCS1

0.9

5

0.9

5

0.8 0.7

4

0.4

3

0.5 0.4

2

0.3 0.2

1

0.6

)

0.5

2

0.7

0

t(2 /

) 3

0.8

4

0.6

0

t(2 /

6

1 a) D1

0.3 0.2

1

0.1

0 -5

-4

-3

-2

-1

0

1

2

3

4

5

0.1

0 -5

0

-4

-3

-2

site index n 6

-1

0

1

2

3

4

5

6

1

b) 6D1

1 d) HEOM 0.9

0.9

5

0.8 0.7

4

0.6

) t(2 /

0.4

3

0.5 0.4

2

0.3 0.2

1

0.7

0

) 0

0.5

2

0.8

4

0.6

3

0.3 0.2

1

0.1

0.1

0 -5

-4

0

site index n

5

t(2

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


-3

-2

-1

0

1

2

3

4

5

0 -5

0

-4

-3

-2

-1

0

1

2

3

4

5

0

site index n

site index n

Figure 8: Time evolution of the exciton probability Pex (t, n) calculated by (a) the single D1 Ansatz, (b) the DM=6 Ansatz, (c) the GCS1 Ansatz with a quantum number of 4 for each 1 mode, and (d) the HEOM method, for the case of J = 0.5, W = 0.5, g = 0.1, and N = 10 t/(2π/ω0 ) > 3. While the exciton probability is nearly trapped at the opposite site of the initial excitation of the ring for the former case, it continues to propagate in two opposite directions for the latter case. The multi-D1 Ansatz and the GCS Ansatz, on the other hand, improve significantly the description of exciton propagation dynamics. It is found that the =6 DM Ansatz and the GCS1 Ansatz provide almost exact results as compared to converged 1

HEOM calculations [see Figs. 8(b-d)], demonstrating their validity as efficient and robust methodologies for studying polaron dynamics(see Appendix for a clearer view on population dynamics after t = 0.4π/ω0 ). Fig. 9 illustrates the time evolution of the exciton probability at a larger transfer integral J = 1.0. As the transfer integral J is increased from 0.5 to 1.0, one can readily observes

21



6

1

c) GCS2

0.9

0.9

5

0.8 0.7

4

0.4

3

0.5 0.4

2

0.3 0.2

1

0.6

) t(2 /

0.5

2

0.7

0

) 3

0.8

4

0.6

0

t(2 /

6

1

a) D1

5

0.3 0.2

1

0.1

0 -5

0.1

0

-4

-3

-2

-1

0

1

2

3

4

0 -5

5

0

-4

-3

-2

site index n 6

-1

0

1

2

3

4

5

site index n 6

1

b) 8D1

1 d) HEOM 0.9

0.9

5

5

0.8 0.7

4

0.4

3

0.5 0.4

2

0.3 0.2

1

0.6

)

0.5

2

0.7

0

t(2 /

) 3

0.8

4

0.6

0

t(2 /

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 22 of 38

0.3 0.2

1

0.1

0.1

0 -5

-4

-3

-2

-1

0

1

2

3

4

5

0 -5

0

0

-4

-3

-2

-1

0

1

2

3

4

5

site index n

site index n

Figure 9: Time evolution of the exciton probability Pex (t, n) calculated by (a) the single D1 Ansatz, (b) the DM=8 Ansatz, (c) the GCS2 Ansatz with a quantum number of 42 for each 1 mode, and (d) the HEOM method, for the case of J = 1.0, W = 0.5, g = 0.1, and N = 10 that the exciton wave packets propagate from the initial site n = 0 to the entire chain with a larger speed. In order to accurately characterize the acceleration of the exciton propagation, it is found that the multi-D1 Ansatz with a much larger multiplicity M and the GCS Ansatz with more mode excitations are required to achieve good agreements with exact results [see Figs. 9(b-d)]. This can be understood by recalling that a large transfer integral induces a delocalized exciton, which in turn triggers phonon distortion along the exciton path, leading to a plane wave-like feature for the phonon wave function. In Fig. 10, the exciton probabilities Pex (t, n) at site n = 0 are plotted for the case of J = 0.5 and J = 1.0, respectively. For J = 0.5, accurate exciton dynamics at both short =6 and long times is reproduced by the DM Ansatz and the GCS1 Ansatz. When increasing 1 =8 J to 1.0, it is found that the exciton probability calculated by the DM Ansatz starts to 1

22


Page 23 of 38

deviates slightly from exact results after t/(2π/ω0 ) > 3, and GCS2 results are in quantitative agreement with numerically exact HEOM. 1

a)

HEOM D1 6D1

0.8

P ex (t,0)

GCS1 0.6

0.4

0.2

0

0

1

2

3

4

5

6

t(2π/ω0) 1 b)

HEOM D1 8D1

0.8

GCS2

P ex (t,0)

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


0.6

0.4

0.2

0

0

1

2

3

4

5

6

t(2π/ω0)

Figure 10: Time evolution of the exciton probability Pex (t, n) at n = 0 with model parameters corresponding to (a) Fig. 8 and (b) Fig. 9

Conclusions In this study, we have applied the multi-D1 Ansatz and the GCS Ansatz to the population dynamics of the Ohmic SBM and the Holstein molecular crystal model. Both variational Ans¨ atze have been validated by comparing the computed population dynamics with those from the numerically exact MCTDH (or ML-MCTDH) and the HEOM methods. In general, the performance of multi-D1 and GCS approaches improves significantly over the single

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

Davydov D1 Ansatz, especially for the SBM with a large energy bias and for the Holstein model with a large transfer integral. Specifically, the GCS Ansatz performs better for small energy bias and weak spin-bath coupling. It has been shown that at least two simultaneously excited vibrational modes are required to obtain accurate dynamics for increased energy bias and spin-bath coupling. Furthermore, in some limiting cases moving from GCS1 to GCS2 Ansatz does not provide a significant improvement (see panel c of Figures 4 and 5). This is very likely due to the nature of the GCS Ansatz, which, as discussed, is similar to a CI expansion: single and double excitations might not be enough to reproduce the dynamical features of the system even if the primitive basis set is time dependent. Increasing the size of the CI expansion would probably lead to better results but is computationally too demanding. On the other hand, the multi-D1 Ansatz provides almost exact spin dynamics over a wide range of energy bias and spin-bath coupling strength. For the Holstein model, exciton probabilities calculated by the GCS2 Ansatz are found to agree nicely with the exact HEOM results for both intermediate and large transfer integrals. As the transfer integral increases, a multiplicity of at least 8 is required to obtain numerically convergent results for the multi-D1 Ansatz. It is worth mentioning that the multi-D1 Ansatz depends on M N (1 + d) parameters, where M is the multiplicity, N the number of electronic states, and d the number of nuclear DOFs. In the GCS1 Ansatz the number of parameters is given by N (1 + d + dw), where w is the size of the basis set, while in the GCS2 Ansatz the number of variables is significantly larger due to the binomial factor. In applications where w ≪ d, GCS1 has better scaling properties than multi-D1 . However, we have shown that in systems with large energy gaps and/or large transfer integrals increasing the basis set in the GCS1 Ansatz does not lead to any improvement, and one has to resort to either GCS2 or the multi-D1 Ansatz. In this case, the multi-D1 Ansatz is extremely advantageous and provides numerically exact results if a sufficiently large multiplicity is used. The advantage of wave function-based methods is to be able to track explicitly the bosonic

24


Page 24 of 38

Page 25 of 38

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


dynamics, which is usually hidden in the reduced density matrix descriptions. Recent developments of quantum optical techniques 84 and ultrafast nonlinear spectroscopy 4,5,85–87 have shed light on the importance of the quantum mixing between electronic and vibrational (photonic) DOFs. The concept of the polaron states (quantum mechanically mixed states between the exciton and intramolecular vibrational modes) has been proposed to explain the long-lived quantum beating in 2D electronic spectra of light harvesting complexes 21–23 . Polaronic states are also important to characterize singlet fission processes in organic crystals 26,88–90 . K¨ uhn and co-workers have recently investigated the coupled exciton-vibrational dynamics of the Fenna-Matthews-Olson complex using the ML-MCTDH approach, which demonstrated the exceptional role of the vibrational and vibronic excitation in photosynthetic energy transfer 91,92 . In addition, an accurate description of a quantum mechanically mixed excitonic and photonic state (a polariton) is essential to study cavity quantum electrodynamics (QED) in the strong light-matter coupling regime 93 . Those processes can be readily investigated by applying the variational approaches based on the multi-D1 Ansatz and the GCS Ansatz. Furthermore, with the above methodologies it is possible to study charge transport mechanisms in realistic organic materials in the presence of an external field by employing a Holstein model with a phase-factor modified transfer integral 94,95 . Finally, the extension to finite temperature cases can be implemented by adopting the Monte Carlo importance sampling method 36 or the thermo field dynamics approach 96 . Work in these directions is in progress.

Acknowledgement The authors thank Maxim Gelin for useful discussion. Support from the Singapore National Research Foundation through the Competitive Research Programme (CRP) under Project No. NRF-CRP5-2009-04, and the Singapore Ministry of Education Academic Research Fund Tier 1 (Grant No. RG106/15) are gratefully acknowledged.

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

Figure 11: Time evolution of the exciton probability Pex (t, n) after t = 0.4π/ω0 calculated by (a) the single D1 Ansatz, (b) the DM=6 Ansatz, (c) the GCS1 Ansatz with a quantum 1 number of 4 for each mode, and (d) the HEOM method, for the case of J = 0.5, W = 0.5, g = 0.1, and N = 10

Appendix: A clearer view on Holstein polaron dynamics In this Appendix, a clearer view is presented on the exciton probability Pex (t, n) after t = 0.4π/ω0 for the case of J = 0.5, W = 0.5, g = 0.1, and N = 10 as shown in Fig. 11. Both 6D1 and GCS1 results are in good agreement with the exact results of HEOM method. As compared to the computationally expensive HEOM method, the variational Ansätze provide a numerically accurate and computationally efficient description of polaron dynamics .

26


Page 26 of 38

Page 27 of 38

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


References (1) G¨ unes, S.; Neugebauer, H.; Sariciftci, N. S. Conjugated polymer-based organic solar cells. Chem. Rev 2007, 107, 1324-1338. (2) Br´ edas, J. L.; Norton, J. E.; Cornil, J.; Coropceanu, V. Molecular understanding of organic solar cells: The challenges. Acc. Chem. Res 2009, 42, 1691-1699. (3) Clarke, T. M.; Durrant, J. R. Charge photogeneration in organic solar cells. Chem. Rev 2010, 110, 6736-6767. (4) Chen, L. P.; Shenai, P.; Zheng, F. L.; Somoza, A.; Zhao, Y. Optimal energy transfer in light-harvesting systems. Molecules 2015, 20, 15224-15272. (5) Mukamel, S. Principles of Nonlinear Optical Spectroscopy; Oxford University Press: New York, 1995. (6) Breuer, H. P.; Petruccione, F. The Theory of Open Quantum Systems; Oxford University Press: New York, 2002. (7) May, V.; K¨ uhn, O. Charge and Energy Transfer Dynamics in Molecular Systems; WileyVCH: Berlin, 2004. (8) Tanimura, Y.; Kubo, R. Time evolution of a quantum system in contact with a nearly gaussian-markoffian noise bath. J. Phys. Soc. Jpn 1989, 58, 101-114. (9) Tanimura, Y. Nonperturbative expansion method for a quantum system coupled to a harmonic oscillator bath. Phys. Rev. A 1990, 41, 6676. (10) Tanimura, Y. Stochastic Liouville, Langevin, Fokker-Planck, and master equation approaches to quantum dissipative systems. J. Phys. Soc. Jpn 2006, 75, 082001.

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

(11) Makri, N.; Makarov, D. E. Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory. J. Chem. Phys 1995, 102, 4600. Makri, N. J. Math. Phys 1995, 36, 2430. (12) Tang, Z.; Yang, X. Q.; Gong, Z.; Wang, H.; Wu, J. L. Extended hierarchy equations of motion for the spin-boson model. J. Chem. Phys 2015, 143, 224112. (13) Moix, J. M.; Cao, J. S. A hybrid stochastic hierarchy equations of motion approach to treat the low temperature dynamics of non-Markovian open quantum systems. J. Chem. Phys 2013, 139, 134106. (14) Chen, L. P.; Zhao, Y.; Tanimura, Y. Dynamics of a one-dimensional Holstein polaron with the hierarchical equations of motion approach. J. Phys. Chem. Lett 2015, 6, 31103115. (15) Ding, J. J.; Xu, J.; Hu, J.; Xu, R. X.; Yan, Y. J. Optimized hierarchical equations of motion theory for drude dissipation and efficient implementation to nonlinear spectroscopies. J. Chem. Phys 2011, 135, 164107. (16) L¨ u, Z. G.; Zheng, H. Quantum dynamics of the dissipative two-state system coupled with a sub-Ohmic bath. Phys. Rev. B 2007, 75, 054302. (17) McCutcheon, D. P. S.; Nazir, A. Consistent treatment of coherent and incoherent energy transfer dynamics using a variational master equation. J. Chem. Phys 2011, 135, 114501. (18) McCutcheon, D. P. S.; Dattani, N. S.; Gauger, E. M.; Lovett, B. W.; Nazir, A. A general approach to quantum dynamics using a variational master equation: application to phonon-damped Rabi rotations in quantum dots. Phys. Rev. B 2011, 84, 081305(R). (19) Fujihashi, Y.; Kimura, A. Improved variational master equation theory for the excitation energy transfer. J. Phys. Soc. Jpn 2014, 83, 014801.

28


Page 28 of 38

Page 29 of 38

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


(20) Sun, K. W.; Fujihashi, Y.; Ishizaki, A.; Zhao, Y. A variational master equation approach to quantum dynamics with off-diagonal coupling in a sub-Ohmic environment. J. Chem. Phys 2016, 144, 204106. (21) Romero, E.; Augulis, R.; Novoderezhkin, V. I.; Ferretti, M.; Thieme, J.; Zigmantas, D.; van Grondelle, R. Quantum coherence in photosynthesis for efficient solar-energy conversion. Nat. Phys 2014, 10, 676-682. (22) Fuller, F. D.; Pan, J.; Gelzinis, A.; Butkus, V.; Senlik, S. S.; Wilcox, D. E.; Ogilvie, J. P. Vibronic coherence in oxygenic photosynthesis. Nat. Chem 2014, 6, 706-711. (23) Oliver, T. A.; Lewis, N. H.; Fleming, G. R. Correlating the motion of electrons and nuclei with two-dimensional electronic-vibrational spectroscopy. Proc. Natl. Acad. Sci. U. S. A 2014, 111, 10061-10066. (24) Musser, A. J.; Liebel, M.; Schnedermann, C.; Wende, T.; Kehoe, T. B.; Rao, A.; Kukura, P. Evidence for conical intersection dynamics mediating ultrafast singlet exciton fission. Nat. Phys 2015, 11, 352-357. (25) Ito, S.; Nagami, T.; Nakano, M. Density analysis of intra- and intermolecular vibronic couplings toward bath engineering for singlet fission. J. Phys. Chem. Lett 2015, 6, 49724977. (26) Fujihashi, Y.; Chen, L. P.; Ishizaki, A.; Wang, J. L.; Zhao, Y. Effect of high frequency modes on singlet fission dynamics. J. Chem. Phys 2017, 146, 044101. (27) Beck, M. H.; Jackle, A. Worth, G. A. Meyer, H. D. The multiconfiguration timedependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets. Phys. Rep 2000, 324, 1-105. (28) Wang, H.; Thoss, M. Multilayer formulation of the multiconfiguration time-dependent Hartree theory. J. Chem. Phys 2003, 119, 1289-1299. 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

(29) Wang, H.; Thoss, M. Quantum-mechanical evaluation of the Bolzmann operator in correlation functions for large molecular systems: A multilayer multiconfiguration timedependent Hartree approach. J. Chem. Phys 2006, 124, 034114. (30) Manthe, U. A multilayer multiconfigurational time-dependent Hartree approach for quantum dynamics on general potential energy surfaces. J. Chem. Phys 2008, 128, 164116. (31) Luo, B.; Ye, J.; Guan, C.; Zhao, Y. Validity of time-dependent trial states for the Holstein polaron. Phys. Chem. Chem. Phys 2010, 12, 15073. (32) Sun, J.; Luo, B.; Zhao, Y. Dynamics of a one-dimensional Holstein polaron with the Davydov Ansatz. Phys. Rev. B 2010, 82, 014305. (33) Zhou, N. J.; Huang, Z. K.; Zhu, J. F.; Chernyak, V.; Zhao, Y. Polaron dynamics with a multitude of Davydov D2 trial states. J. Chem. Phys 2015, 143, 014113. (34) Zhou, N. J.; Chen, L. P.; Huang, Z. K.; Sun, K. W.; Tanimura, Y.; Zhao, Y. Fast, accurate simulation of polaron dynamics and multidimensional spectroscopy by multiple Davydov trial states. J. Phys. Chem. A 2016, 120, 1562-1576. (35) Wang, L.; Chen, L. P.; Zhou, N. J.; Zhao, Y. Variational dynamics of the sub-Ohmic spin-boson model on the basis of multiple Davydov D1 states. J. Chem. Phys 2016, 144, 024101. (36) Wang, L.; Fujihashi, Y.; Chen, L.; Zhao, Y. Finite-temperature time-dependent variation with multiple Davydov states. J. Chem. Phys 2017, 146, 124127. (37) Huang, Z.; Wang, L.; Wu, C.; Chen, L.; Grossmann, F.; Zhao, Y. Polaron dynamics with off-diagonal coupling: beyond the Ehrenfest approximation. Phys. Chem. Chem. Phys 2017, 19, 1655-1668.

30


Page 30 of 38

Page 31 of 38

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


(38) Huang, Z.; Chen, L.; Zhou, N.; Zhao, Y. Transient dynamics of a one-dimensional Holstein polaron under the influence of an external electric field. Ann. Phys 2017, 529, 1600367. (39) Ye, J.; Sun, K.; Zhao, Y.; Yu, Y.; Lee, C. K.; Cao, J. Excitonic energy transfer in light-harvesting complexes in purple bacteria. J. Chem. Phys 2012, 136, 245104. (40) Sun, K.; Ye, J.; Zhao, Y. Path induced coherent energy transfer in light-harvesting complexes in purple bacteria. J. Chem. Phys 2014, 141, 124103. (41) Somoza, A.; Chen, L. P.; Sun, K. W.; Zhao, Y. Probing ultrafast excitation energy transfer of the chlorosome with exciton-phonon variational dynamics. Phys. Chem. Chem. Phys 2016, 18, 20298. (42) Chen, L.; Gelin, M. F.; Domcke, W.; Zhao, Y. Theory of femtosecond coherent doublepump single-molecule spectroscopy: application to light harvesting complexes. J. Chem. Phys 2015, 142, 164106. (43) Huynh, T. D.; Sun, K.; Gelin, M. F.; Zhao, Y. Polaron dynamics in two-dimensional photon-echo spectroscopy of molecular rings. J. Chem. Phys 2013, 139, 104103. (44) Martinazzo, R.; Nest, M.; Saalfrank, P.; Tantardini, G. F. A local coherent-state approximation to system-bath quantum dynamics. J. Chem. Phys 2006, 125, 194102. (45) Mendiv-Tapia, D.; Lasorne, B.; Worth, G. A.; Robb, M. A.; Bearpark, M. J. Towards converging non-adiabatic direct dynamics calculations using frozen-width variational Gaussian product basis functions. J. Chem. Phys 2012, 137, 22A548. (46) Lasorne, B.; Worth, G. A.; Robb, M. A. Excited-state dynamics. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 460. (47) Martinez, T. J.; Levine, R. D. First principles molecular dynamics on multiple electronic states: A case study of NaI. J. Chem. Phys 1996, 105, 6334. 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

(48) Burghardt, I.; Giri, K.; Worth, G. A. Multimode quantum dynamics using Gaussian wavepackets: The Gaussian-based multiconfiguration time-dependent Hartree (GMCTDH) method applied to the absorption spectrum of pyrazine. J. Chem. Phys 2008, 129, 174104. (49) Richings, G. W.; Polyak, I.; Spinlove, K. E.; Worth, G. A.; Burghardt, I.; Lasorne, B. Quantum dynamics simulations using Gaussian wavepackets: the vMCG method. Int. Rev. Phys. Chem 2015, 34, 269-308. (50) Shalashilin, D. V.; Child, M. S. The phase space CCS approach to quantum and semiclassical molecular dynamics for high-dimensional systems. Chem. Phys 2004, 304, 103120. (51) Shalashilin, D. V. Quantum mechanics with the basis set guided by Ehrenfest trajectories: Theory and application to spin-boson model. J. Chem. Phys 2009, 130, 244101. (52) Hagedorn, G. A. Semiclassical quantum mechanics. I. The ~ → 0 limit for coherent states. Commun. Math. Phys 1980, 71, 77-93. (53) Hagedorn, G. A. Semiclassical quantum mechanics. III. The large order asymptotics and more general states. Ann. Phys 1981, 135, 58-70. (54) Nieto, M. M. Displaced and squeezed number states. Phys. Lett. A 1997, 229, 135-143. (55) de Oliveira, F. A. M.; Kim, M. S.; Knight, P. L.; Buˇ z ek, V. Properties of displaced number states. Phys. Rev. A 1990, 41, 2645. (56) Combescure, M. The squeezed state approach of the semiclassical limit of the timedependent Schrödinger equation. J. Math. Phys 1992, 33, 3870. (57) Borrelli, R.; Peluso, A. Quantum dynamics of electronic transitions with Gauss-Hermit wave packets. J. Chem. Phys 2016, 144, 114102.

32


Page 32 of 38

Page 33 of 38

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


(58) Borrelli, R.; Maxim, F. G. The generalized coherent state ansatz: application to quantum electron-vibrational dynamics. Chem. Phys 2016, 481, 91-98. (59) Dirac, P. A. M. Note on exchange phenomena in the Thomas atom. Math. Proc. Cambridge Philos. Soc 1930, 26, 376. Frenkel, J. Wave Mechanics; Oxford University Press: New York, 1934. (60) Senitzky, I. R. Harmonic oscillator wave functions. Phys. Rev 1954, 95, 1115-1116. (61) Roy, S. M.; Singh, V. Generalized coherent states and the uncertainty principle. Phys. Rev. D 1982, 25, 3413-3416. (62) Satyanarayana, M. V. Generalized coherent states and generalized squeezed coherent states. Phys. Rev. D 1985, 32, 400. (63) Schrödinger, E. Der stetige bergang von der Mikro- zur Makromechanik. Naturwiss 1926, 44, 664-666. (64) Jankowiak, H. C.; Stuber, J. L.; Berger, R. Vibronic transitions in large molecular systems: rigorous prescreening conditions for Franck-Condon factors. J. Chem. Phys 2007, 127, 234101. (65) Borrelli, R.; Capobianco, A.; Landi, A.; Peluso, A. Vibronic couplings and coherent electron transfer in bridged systems. Phys. Chem. Chem. Phys 2015, 17, 30937-30945. (66) Leggett, A. J.; Chakravary, S.; Dorsey, A.; Fisher, M.; Garg, A.; Zwerger, W. Dynamics of the dissipative two-state system. Rev. Mod. Phys 1987, 59, 1. (67) Weiss, U. Quantum Dissipative Systems; World Scientific: Singapore, 2007. (68) Costi, T. A.; Mckenzie, R. H. Entanglement between a qubit and the environment in the spin-boson model. Phys. Rev. A 2003, 68, 034301. Khveshchenko, D. V. Quantum impurity models of noisy qubits. Phys. Rev. B 2004, 69, 153311. 33



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

(69) Hur, K. L.; Beaupr` a, P. C.; Hofstetter, W. Entanglement and criticality in quantum impurity systems. Phys. Rev. Lett 2007, 99,126801. (70) Yao, Y.; Duan, L. W.; L¨ u, Z. G.; Wu, C. Q.; Zhao, Y. Dynamics of the sub-Ohmic spin-boson model: A comparison of three numerical approaches. Phys. Rev. E 2013, 88, 023303. Duan, L. W.; Wang, H.; Chen, Q. H.; Zhao, Y. Entanglement dynamics of two qubits coupled individually to Ohmic baths. J. Chem. Phys 2013, 139, 044115. Deng, T. R.; Yan, Y. Y.; Chen, L. P.; Zhao, Y. Dynamics of the two-spin spin-boson model with a common bath. J. Chem. Phys 2016, 144, 144102. (71) Bulla, R.; Tong, N. H.; Vojta, M. Numerical renormalization group for bosonic systems and application to the sub-Ohmic spin-boson model. Phys. Rev. Lett 2003, 91, 170601. Vojta, M.; Tong, N. H.; Bulla, R. Quantum phase transitions in the sub-Ohmic spin-boson model: failure of the quantum-classical mapping. 2005, 94, 070604. (72) Zhao, Y.; Yao, Y.; Chernyak, V.; Zhao, Y. Spin-boson model with diagonal and offdiagonal coupling to two independent baths: ground-state phase transition in the deep sub-Ohmic regime. J. Chem. Phys 2014, 140, 161105. Zhou, N. J.; Chen, L. P.; Zhao, Y. Mozyrsky, Chernyak, V.; Zhao, Y. Ground-state properties of sub-Ohmic spin-boson model with simultaneous diagonal and off-diagonal coupling. Phys. Rev. B 2014, 90, 155135. Zhou, N. J.; Chen, L. P.; Xu, D. Z.; CHernyak, V.; Zhao, Y. Symmetry and the critical phase of the two-bath spin-boson model: Ground-state properties. Phys. Rev. B 2015, 91, 195129. (73) Hu, X.; Ritz, T.; Damjanovi´ c, A.; Schulten, K. Pigment organization and transfer of electronic excitation in the photosynthetic unit of purple bacteria. J. Phys. Chem. B 1997, 101, 3854-3871. (74) Makri, N. The linear response approximation and its lowest order corrections: An influence functional approach. J. Phys. Chem. B 1999, 103, 2823-2829.

34


Page 34 of 38

Page 35 of 38

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


(75) Wang, H.; Thoss, M.; Miller, W. H. Systematic convergence in the dynamical hybrid approach for complex systems: A numerically exact methodology. J. Chem. Phys 2001, 115, 2979. (76) Moix, J. M.; Zhao, Y.; Cao, J. S. Equilibrium-reduced density matrix formulation: Influence of noise, disorder, and temperature on localization in excitonic systems. Phys. Rev. B 2012, 85, 115412. (77) Tanimura, Y. Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities. J. Chem. Phys 2014, 141, 044114. (78) Tanimura, Y. Real-time and imaginary-time quantum hierarchal Fokker-Planck equations. J. Chem. Phys 2015, 142, 144110. (79) Cerrillo, J.; Cao, J. S. Non-Markovian dynamical maps: Numerical processing of open quantum trajectories. Phys. Rev. Lett 2014, 112, 110401. (80) Thoss, M.; Wang, H.; Miller, W. H. Self-consistent hybrid approach for complex systems: application to the spin-boson model with Debye spectral density. J. Chem. Phys 2001, 115, 2991. (81) Vendrell, O.; Meyer, H. D. Multilayer multiconfiguration time-dependent Hartree method: Implementation and application to a Henon-Heiles Hamiltonian and to pyrazine. J. Chem. Phys 2011, 134, 044135. (82) Borrelli, R.; Donato, M. D.; Peluso, A. Quantum dynamics of electron transfer from bacteriochlorophyll to pheophytin in bacterial reaction centers. J. Chem. Theory. Comput 2007, 3, 673-680. (83) Borrelli, R.; Donato, M. D.; Peluso, A. Role of intramolecular vibrations in long-range

35



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

electron transfer between pheophytin and ubiquinone in bacterial photosynthetic reaction centers. Biophys. J 2005, 89, 830-841. (84) Haroche, S.; Raimond, J. M. Exploring the Quantum: Atoms, Cavities, and Photons; Oxford University Press, New York, 2006. (85) Chenu, A.; Scholes, G. D. Coherence in energy transfer and photosynthesis. Annu. Rev. Phys. Chem 2015, 66, 69-96. (86) Romero, E.; Novoderezhkin, V. I.; van Grondelle, R. Quantum design of photosynthesis for bio-inspired solar-energy conversion. Nature 2017, 543, 355-365. (87) Scholes, G. D.; Fleming, G. R.; Chen, L. X.; Aspuru-Guzik, A.; Buchleitner, A.; Coker, D. F.; Engel, G. S.; van Grondelle, R.; Ishizaki, A.; Jonas, D. M. et al. Using coherence to enhance function in chemical and biophysical systems. Nature 2017, 543, 647-656. (88) Bakulin, A. A.; Morgan, S. E.; Kehoe, T. B.; Wilson, M. W.; Chin, A. W.; Zigmantas, D.; Egorova, D.; Rao, A. Real-time observation of multiexcitonic states in ultrafast singlet fission using coherent 2D electronic spectroscopy. Nat. Chem 2016, 8, 16-23. (89) Tempelaar, R.; Reichman, D. R. Vibronic exciton theory of singlet fission. I. Linear absorption and the anatomy of the correlated triplet pair state. R. J. Chem. Phys 2017, 146, 174703. (90) Tempelaar, R.; Reichman, D. R. Vibronic exciton theory of singlet fission. II. Twodimensional spectroscopic detection of the correlated triplet pair state. J. Chem. Phys 2017, 146, 174704. (91) Schulze, J.; K¨ u, O. Explicit correlated exciton-vibrational dynamics of the FMO complex. J. Phys. Chem. B 2015, 119, 6211-6216.

36


Page 36 of 38

Page 37 of 38

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


(92) Schulze, J.; K¨ u, O. Multi-layer multi-configuration time-dependent Hartree (MLMCTDH) approach to the correlated exciton-vibrational dynamics in the FMO complex. J. Chem. Phys 2016, 144, 185101. (93) Díaz-Camacho, G.; Bermudez, A.; Carcía-Ripoll, J. J. Dynamical polaron Ansatz. A theoretical tool for the ultrastrong-coupling regime of circuit QED. Phys. Rev. A 2016, 93, 043843. (94) Yu, J. F.; Wu, C. Q.; Sun, X.; Nasu, K. Breather in the motion of a polaron in an electric field. Phys. Rev. B 2004, 70, 064303. (95) Vidmar, L.; Bonca, J.; Mierzejewski, M.; Prelovsek, P.; Trugman, S. A. Nonequilibrium dynamics of the Holstein polaron driven by an external electric field. Phys. Rev. B 2011, 83, 134301. (96) Borrelli, R.; Gelin, M. F. Quantum electron-vibrational dynamics at finite temperature: Thermo field dynamics approach. J. Chem. Phys 2016, 145, 224101.

37



Graphical TOC Entry 1

1

(a)

0.8

0.6

0.6

(b)

P(t)

0.8

P(t)

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 38 of 38

0.4

0.4 c1=c 2=0

2D1

c1=c 2=1 c1=1,c 2=2

0.2

4D1

0.2

5D1

ML-MCTDH 0

0

100

200

300

400

0

500

t(fs)

ML-MCTDH 0

100

200

300

t(fs)

38


400

500

Dynamics of Coupled ElectronâBoson Systems with the Multiple

Recommend Documents

Dynamics of Coupled ElectronâBoson Systems with the Multiple