# Exciton Transfer in Organic Photovoltaic Cells: A Role of Local and Nonlocal
Electron-Phonon Interactions in a Donor Domain

Mauro Cainelli<EMAIL_ADDRESS>Department of Chemistry, Graduate
School of Science, Kyoto University, Kyoto 606-8502, Japan  Yoshitaka Tanimura
<EMAIL_ADDRESS>Department of Chemistry, Graduate School of
Science, Kyoto University, Kyoto 606-8502, Japan

###### Abstract

We theoretically investigate an exciton transfer process in a donor domain of
organic photovoltaic cells focusing on the roles of local and nonlocal
electron-phonon interactions. Our model consists of a three-level system
described by the Holstein-Peierls Hamiltonian coupled to multiple heat baths
for local and nonlocal molecular modes characterized by Brownian spectral
distribution functions. We chose tetracene as a reference donor molecule,
where the spectral distribution functions of the local and nonlocal modes
exist. We then employ the reduced hierarchy equations of motion (HEOM)
approach to simulating the dynamics of the system under the influence of the
environment as a function of the electron-phonon coupling strength and
temperature. We rigorously calculate the reduced density matrix elements to
explain the timescale of dynamics under the influence of the dissipative local
and nonlocal modes. The results indicate that the strong nonlocal electron-
phonon interaction under high temperature conditions favors the exciton
transfer process and enhances the efficiency of organic photovoltaic
materials, while the lifetime of the exciton becomes shorter due to a low
frequency local mode.

††preprint: AIP/123-QED

##  I Introduction

The development of better photovoltaic devices is a matter of great interest
in the investigation of renewable energy sources for a wide variety of systems
including biology- and material-oriented venues. Organic solar cells are
promising materials because of its low fabrication cost, lightweight and the
flexibility of design in comparison to silicon-based materials exist in the
market.[1, 2] Typical composition consists of polymers of small molecules as
donors and fullerene derivatives as acceptors, although non-fullerene
derivatives have also exhibited promising results.[3] However, the power
conversion efficiency of organic solar cells remains with maximum values of
11-17% that are still much lower than the inorganic counterparts.[4, 5]
Moreover, the mechanism on the organic devices is more complicated, mainly
because of the strong electron-phonon interaction due to the low dielectric
constant, where the quantum coherence between the electron and phonon plays an
essential role. Thus, for the investigation to improve the efficiency of such
systems, we have to employ fully quantum mechanical descriptions of broad
validity.[6]

During the process, sunlight is absorbed predominantly by the molecular donor
with the creation of strongly bounded electron-hole pairs (excitons). The
excitons then diffuse to the heterojunction, in which they dissociate into
electrons and holes, and are then collected at the corresponding
electrodes.[7] While most of the electronic dynamics occur in the active
layer, recent studies indicate that timescales and length scales of these
processes depend on the local geometry and electrostatic field. [8, 9]
Therefore, for the design of photovoltaic devices, understanding and
controlling the morphology is significant.[10, 11, 12] Although charge
transport and charge separation processes in organic solar cells have been
studied intensively, investigations concerning the role of incoherent
(hopping) [13] and coherent (delocalization) [14, 15, 16] transfer mechanisms
in different structure and physical conditions have not been well
explored.[17, 18, 19]

As recent theoretical analysis indicated, the effect of the nonlocal
(intermolecular) electron-phonon coupling on charge transport processes, in
addition to that of the local (intramolecular) one, is essential to account
for the energy transfer mechanism on such devices.[20] In heterojunction
blends, the entire process is ultrafast and takes place within 100fs [21, 22]
due to the presence of short lived delocalized excited species,[23, 24, 25]
which associate with the formation of hot charge transfer (CT) states.[26] The
hot CT states are assumed to be the main contribution to the generated
photocurrent, because they allow us to create a charge separation before the
thermalization process occur, while incoherent charge separation mechanisms
are characterized by slower separation of charges via multiple hopping steps
between localized states. In organic solar cells, both mechanisms are assumed
to be present.

It should be noted that the predominant dynamical processes in solar cells are
a system dependent, in which the molecules used as donor and acceptor,
geometry disposition, and electron-phonon coupling strengths determine the
power conversion efficiency. In order to account for the problem of this kind,
processes such as light absorption, energy- and charge-transports in organic
semiconductors must be taken into account by simulating dynamics in both
electronic and molecular degrees of freedom simultaneously.[27] A widely used
approach in order to properly treat vibrational modes in the framework of
quantum mechanics is based on the Holstein-Peierls model,[28] where the
Holstein (local) part of the Hamiltonian describes the variation of the site
energies[29] and the Peierls (nonlocal) part of it represents the modulation
of the transfer integrals.[30] Although the importance of nonlocal electron-
phonon coupling on the charge transport properties has been realized
recently,[31] the Hamiltonian of this model is not diagonalizable and the
analysis of the system is not easy. Thus, the effects of the local and
nonlocal couplings are usually treated collectively under the mean field
polaron (MFP) approximation,[32] ignoring the interaction between the local
and nonlocal modes. Moreover, because the electronic processes occur in
condensed phases, where the surrounding molecules provide the thermal
fluctuation and dissipation, we have to adapt open quantum dynamics treatment
for the investigation of time-irreversible reaction dynamics.[33, 34]

In order to study electron transfer problems in connection to open quantum
dynamics theory, a commonly used model considers the electronic states coupled
to an intermediate harmonic oscillator, which is further coupled to a heat
bath.[35] Then, the harmonic mode can be included in the bath by carrying out
a canonical transformation, which leads to a multilevel system coupled to the
heat bath with the Brownian spectral distribution (BSD) function.[36] Such a
system can be treated using the HEOM formalism, in a numerically rigorous
manner,[37] even in the low temperature case.[38, 39, 40] The HEOM formalism
can handle not only the strong system-bath coupling but also quantum coherence
between the system and bath that is important for the investigation of the
organic materials, as demonstrated in photosynthesis antenna systems[41, 42,
43, 44] and DNA.[45] Furthermore, this method does not require any
approximation, most notably the rotating wave approximation (RWA) or the MFP
approximation.

The HEOM approach has been applied to the case of solid state materials
described by a deformation potential [46] and the Holstein Hamiltonian.[47] In
the Holstein case, because we can study only a small system with the finite
number of phonon modes associated with the system site, special treatment is
necessary to maintain the stability of the equations.[48] The situation is the
same in the present case, and thus we further introduce the heat-bath into the
regular Holstein-Peierls model.

In this paper, we present a study of the nonlocal electron-phonon coupling
effect and the relation with a local coupling effect on the charge
transference process in organic photovoltaic devices, particularly in the
donor layer. We employ the HEOM approach for a BSD function to simulate the
non-Markovian dynamics of the Holstein-Peierls model that is further coupled
to the heat-bath.

The organization of this paper is as follows: In Sec. II, we introduce the
Holstein-Peierls + bath (HPB) model. Then, the HEOM for the HPB model is
explained. In Sec. III we present the details of the calculation and the
results of the simulation for various conditions. Sec. IV is devoted to
concluding remarks.

##  II The HEOM for a Holstein-Peierls + bath system

The Holstein-Peierls model employs a tight-binding picture of an electronic
system that is linearly coupled to an optical phonon mode at each site. In the
present study, we further introduce a harmonic oscillator heat-bath.

The total Hamiltonian is expressed in terms of the electron Hamiltonian
($\hat{H}_{el}$), phonon Hamiltonian ($\hat{H}_{ph}$) and electron-phonon
interaction with the heat-bath ($\hat{H}_{ph-B}$) as

$\displaystyle\hat{H}_{tot}=\hat{H}_{el}+\hat{H}_{el-ph}+\hat{H}_{ph-B},$ (1)

where

$\displaystyle\hat{H}_{el}=\sum_{i}\varepsilon_{i}^{0}\hat{a}_{i}^{\dagger}\hat{a}_{i}+\sum_{i\neq
j}t_{ij}^{0}\hat{a}_{i}^{\dagger}\hat{a}_{j},$ (2) $\displaystyle\hat{H}_{el-
ph}$ $\displaystyle=$
$\displaystyle\sum_{i}\sum_{{\alpha}}{\hat{V}}_{\alpha}^{i}\left(\hat{b}_{\alpha}^{\dagger}+\hat{b}_{\alpha}\right)+\sum_{ij}\sum_{{\alpha}}{\hat{V}}_{\alpha}^{ij}\left(\hat{b}_{\alpha}^{\dagger}+\hat{b}_{\alpha}\right)$
(3)

with

$\displaystyle{\hat{V}}_{\alpha}^{i}=\frac{1}{\sqrt{{\cal
N}_{C}}}\hbar\Omega_{\alpha}g_{i}^{\alpha}\hat{a}_{i}^{\dagger}\hat{a}_{i},$
$\displaystyle{\hat{V}}_{\alpha}^{ij}=\frac{1}{\sqrt{{\cal
N}_{C}}}\hbar\Omega_{\alpha}g_{ij}^{\alpha}\hat{a}_{i}^{\dagger}\hat{a}_{j},$
(4)

and

$\displaystyle\hat{H}_{ph-B}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\hbar\Omega_{\alpha}\left(\hat{b}_{\alpha}^{\dagger}\hat{b}_{\alpha}+\frac{1}{2}\right)+\sum_{\alpha}\sum_{k=1}^{N_{\alpha}}\hbar\omega_{k}\left(\hat{b}_{k}^{\dagger}\hat{b}_{k}+\frac{1}{2}\right)$
(5)
$\displaystyle+\sum_{\alpha}\left(\hat{b}_{\alpha}^{\dagger}+\hat{b}_{\alpha}\right)\sum_{k=1}^{N_{\alpha}}c_{k}^{\alpha}\left(\hat{b}_{k}^{\dagger}+\hat{b}_{k}\right).$

Here, $\hat{a}_{i}^{\dagger}$ ($\hat{a}_{i}$) is the creation (annihilation)
operator of the electronic quasi-particles (electrons, holes, or Frenkel
excitons) and $\varepsilon_{i}^{0}$ and $t_{ij}^{0}$ are the on-site
electronic energy and the ampliture of the transfer integral for the $i-j$th
electronic states that correspond to the diagonal and off-diagonal elements of
the on-site Hamiltonian, respectively. The creation (annihilation) operator of
the phonon (or vibron) mode ${\alpha}$ with the frequency $\Omega_{\alpha}$ is
expressed as $\hat{b}_{\alpha}^{\dagger}$ ($\hat{b}_{\alpha}$) . The local and
nonlocal system-bath interactions are expressed as ${\hat{V}}_{\alpha}^{i}$
and ${\hat{V}}_{\alpha}^{ij}$ with the dimensionless local and nonlocal
electron-phonon couplings strength $g_{i}^{\alpha}$ and $g_{ij}^{\alpha}$,
respectively. The constant ${\cal N}_{C}$ is the number of unitary cells
considered. The operator $\hat{b}_{k}^{\dagger}$ ($\hat{b}_{k}$) is the
creation (annihilation) operator of the heat-bath that can be characterized by
the oscillator-bath coupling strength and the inverse temperature $\beta\equiv
1/k_{\mathrm{B}}T$, where $k_{\mathrm{B}}$ is the Boltzmann constant. The bath
system is typically modeled by the spectral distribution function (SDF),
defined by

$\displaystyle
J_{\alpha}(\omega)\equiv\sum_{\alpha}\sum_{k=1}^{N_{\alpha}}\left(c_{k}^{\alpha}\right)^{2}\delta(\omega-\omega_{k}).$
(6)

For the heat bath to be an unlimited heat source possessing an infinite heat
capacity, the number of heat bath oscillators $N_{\alpha}$ is effectively made
infinitely large by replacing $J_{\alpha}(\omega)$ with a continuous
distribution. This Hamiltonian can be regarded as either empirical or first-
principles model, whether the system parameters are obtained from experimental
data or quantum chemistry calculations. The electron-phonon interactions
produce a time-dependent variation of the transport parameters and thus
introduce a dynamic disorder: Local and nonlocal electron-phonon couplings
correspond to diagonal and off-diagonal dynamic disorder mechanisms.

The challenges in using the Holstein-Peierls model are the following.[49]
First, the Hamiltonian is not exactly solvable. Second, the commonly used
approximations to solve this Hamiltonian are not fully satisfied for the
investigation of organic semiconductors. Third, it is difficult to find
efficient and reliable schemes to evaluate the system parameters
($\varepsilon_{i}^{0}$, $t_{ij}^{0}$, $\omega_{\alpha}$, $g_{i}^{\alpha}$,and
$g_{ij}^{\alpha}$).

In general, the second point in the above arises as a consequence of the first
and, in order to study the model, several approximations have been introduced:
First, the electron-electron interaction is treated in a framework of mean
field theory. Second, the amplitude $t_{ij}^{0}$ must be small in comparison
with the excitation energy of the intramolecular mode. Third, the electron-
phonon coupling is considered to be a linear form. Four, the nonlocal coupling
is treated under the band renormalization theory, assuming the collective
effect of low-frequency modes (MFP approximation). Five, phonon modes are
treated under rigid-body approximation that ignores the internal molecular
vibrations with environmental degrees of freedom that consist of the
intermolecular translational and librational molecular modes, in addition to
the acoustic modes. Moreover, in order to understand the structural effects on
the charge dynamics and the role of mixing morphology on the efficient pathway
to free charge photogeneration, it is necessary to account for the individual
effects of the couplings at different phonon frequencies and the interaction
between intra- and intermolecular vibrational modes.

In order to overcome some of the limitation mentioned above, here we employ
the HEOM formalism. The MFP approximation can be avoided because the
Hamiltonian does not require to be diagonalized. Thus, the phonon modes are
treated dynamically employing the system-bath model and large values of
$t_{ij}^{0}$ can be treated by choosing an appropriate hierarchy depth.
Although in this work we limited our analysis to the linear electron-phonon
couplings case, the extension to study nonlinear coupling is also possible.

By using the canonical transformation, the Hamiltonian in Eq.(1) can be
expressed in terms of the creation and annihilation operators of the heat-bath
+ primary oscillation ($\hat{b}_{k}^{\prime\dagger}$ and
$\hat{b}_{k}^{\prime}$, respectively) in the following general form

$\displaystyle\hat{H}_{tot}$ $\displaystyle=$
$\displaystyle\hat{H}_{el}+\sum_{\alpha}\sum_{k=1}^{N_{\alpha}}\sum_{i}d_{k}^{ii\alpha}\hat{a}_{i}^{\dagger}\hat{a}_{i}(\hat{b}_{k}^{\prime\dagger}+\hat{b}_{k}^{\prime})$
(7)
$\displaystyle+\sum_{\alpha}\sum_{k=1}^{N_{\alpha}}\sum_{ij}f_{k}^{ij\alpha}\hat{a}_{i}^{\dagger}\hat{a}_{j}(\hat{b}_{k}^{\prime\dagger}+\hat{b}_{k}^{\prime})+\sum_{\alpha}\sum_{k=1}^{N_{\alpha}}\hbar\omega_{k}\left(\hat{b}_{k}^{\prime\dagger}\hat{b}_{k}^{\prime}+\frac{1}{2}\right),$

where $d_{k}^{ii\alpha}$ and $f_{k}^{ij\alpha}$ are the dimensionless
constants that represent the local and nonlocal interactions of the system
with the oscillator-bath. Assuming that the SDF of the later to be
$J_{\alpha}(\omega)=\gamma_{\alpha}\omega$, we have the BSD defined as

$\displaystyle
J_{\alpha}^{\prime}(\omega)=\frac{\hbar\lambda_{\alpha}}{2\pi}\frac{\gamma_{\alpha}\Omega_{\alpha}^{2}\omega}{(\Omega_{\alpha}^{2}-\omega^{2})^{2}+\gamma^{2}\omega^{2}},$
(8)

where $\lambda_{\alpha}$ relates to the reorganization energy and nonlocal
relaxation energy which accounts for the displacement of the excited state in
a relation to the ground state in coordinate space as a consequence of the
local and nonlocal electron-phonon interactions, respectively. The parameter
$\gamma_{\alpha}$ is the coupling strength between the oscillator and the
bath, which is related to the peak width of $J_{\alpha}^{\prime}(\omega)$.

The reduced hierarchy equations of motion for the Brownian distribution can be
expressed as (see Appendix A)

$\displaystyle\frac{\partial}{\partial t}$
$\displaystyle\hat{\rho}_{\\{n^{\alpha},m^{\alpha};j_{1}^{\alpha}...j_{K^{\alpha}}^{\alpha}\\}}(t)$
(9) $\displaystyle=$
$\displaystyle-\left[\frac{i}{\hbar}\hat{H}_{el}^{\times}+\sum_{\alpha}\left\\{\frac{(n^{\alpha}+m^{\alpha})}{2}\gamma_{\alpha}-i(n^{\alpha}-m^{\alpha})\zeta_{\alpha}+\sum_{k=1}^{K^{\alpha}}j_{k}^{\alpha}\nu_{k}^{\alpha}-\hat{\Xi}_{\alpha}\right\\}\right]\hat{\rho}_{\\{n^{\alpha},m^{\alpha};j_{1}^{\alpha}...j_{K^{\alpha}}^{\alpha}\\}}(t)$
$\displaystyle+\sum_{\alpha}\hat{V}^{\times}_{\alpha}\bigg{[}\hat{\rho}_{\\{n^{\alpha}+1,m^{\alpha};j_{1}^{\alpha}...j_{K^{\alpha}}^{\alpha}\\}}(t)+\hat{\rho}_{\\{n^{\alpha},m^{\alpha}+1;j_{1}^{\alpha}...j_{K^{\alpha}}^{\alpha}\\}}(t)\bigg{]}$
$\displaystyle+\sum_{\alpha}n^{\alpha}\hat{\Theta}_{-}^{\alpha}\hat{\rho}_{\\{n^{\alpha}-1,m^{\alpha};j_{1}^{\alpha}...j_{K^{\alpha}}^{\alpha}\\}}(t)+\sum_{\alpha}m^{\alpha}\hat{\Theta}_{+}^{\alpha}\hat{\rho}_{\\{n^{\alpha},m^{\alpha}-1;j_{1}^{\alpha}...j_{K^{\alpha}}^{\alpha}\\}}(t)$
$\displaystyle+\sum_{\alpha}\sum_{k=1}^{K^{\alpha}}\hat{V}_{\alpha}^{\times}\hat{\rho}_{\\{n^{\alpha},m^{\alpha};j_{1}^{\alpha}...j_{k+1}^{\alpha},...j_{K^{\alpha}}^{\alpha}\\}}(t)+\sum_{\alpha}\sum_{k=1}^{K^{\alpha}}j_{k}^{\alpha}\nu_{k}^{\alpha}\hat{\Psi}_{k}^{\alpha}\hat{\rho}_{\\{n^{\alpha},m^{\alpha};j_{1}^{\alpha}...j_{k-1}^{\alpha},...j_{K^{\alpha}}^{\alpha}\\}}(t),$

where the first term in brackets is the Liouvillian
($\hat{\mathcal{L}}^{(n,m)}$) that involves the temperature correction terms.
For a high number of hierarchy elements, the term
$(n^{\alpha}+m^{\alpha})\gamma_{\alpha}/2+\sum_{k=1}^{K^{\alpha}}j_{k}^{\alpha}\nu_{k}^{\alpha}$
becomes much larger than the characteristic time of the system. In this case,
the hierarchy can be truncated by the terminator.[50]

##  III Results and Discussion

Figure 1:  Schematic view of the model that is described by the Hamiltonian
given in Eqs.(1)-(5). The ground state and the two excited states are coupled
to a local mode that is further coupled to the local bath ($H_{Bl}$),
respectively. The excited states are coupled to the nonlocal mode that is
further coupled to the nonlocal bath ($H_{Bnl}$). Figure 2: Two configurations
of the excited state pathways. Each pair represents the dimer molecule
illustrated in Fig. 1.

We chose tetracene as a reference donor molecule under study where the SDFs of
the electron-phonon coupling have been known.[51] In order to simulate the
exciton transfer process between the molecules, we consider a multi-state
system in which each state represents a dimer in the donor layer. While we
might simulate a molecular chain in the framework of the present model by
increasing the number of excited states, here we consider a simple three-state
system in order to investigate the role of local and nonlocal modes. Thus, the
system is described by a single ground state and two excited states expressed
as $\big{|}1\big{\rangle}$, $\big{|}2\big{\rangle}$, and
$\big{|}3\big{\rangle}$, respectively. As illustrated in Fig. 1, the excited
states $\big{|}2\big{\rangle}$ and $\big{|}3\big{\rangle}$ interact with the
nonlocal and local modes, while the ground state $\big{|}1\big{\rangle}$ only
interacts with the local mode. The local and nonlocal modes are further
coupled to their own heat-bath. We then consider the transfer between the two
excited states in two directions, represented by two different pair of
geometry dispositions of the dimer (Fig. 2), whose difference is described by
the values of the transfer integral and the nonlocal coupling strength.

The exciton Hamiltonian is then expressed in the matrix form as

$\displaystyle\hat{H}_{el}=\hbar\begin{bmatrix}\omega_{1}&\Delta_{12}&0\\\
\Delta_{21}&\omega_{2}&\Delta_{23}\\\ 0&\Delta_{32}&\omega_{3}\end{bmatrix},$
(10)

where $\hbar\omega_{i}$ is the Frenkel exciton energy of the $i$th state
($\varepsilon_{i}$) and $\hbar\Delta_{ij}$ is the coupling strength between
the $i$th and $j$th states (transfer integral $t_{ij}$). For all of our
computations, we fixed the eigenenergy of the system as $\omega_{1}=0$cm-1 and
$\omega_{2}=\omega_{3}=500$cm-1. We then use the excitation energy as the
characteristic frequency of the system $\omega=500$cm-1 with the time scale
$1/{\omega}=66.71$fs. Although the band gap between the ground and singlet
state in oligoacenes is of the order of 2-2.5eV ($16000-20000$cm-1), here we
chose a much smaller value for $\omega$ to qualitatively investigate the
effect of phonon interactions upon ultrafast dynamics with suppressing
numerical costs. The off-diagonal elements representing the interaction
between the ground and excited state are $\Delta_{12}=\Delta_{21}=0.32\omega$.
It should be noted that, while similar single crystals, which include
rubrene[52, 53] and pentacene[54] have been modeled and investigated, the
parameter values of the off-diagonal element of tetracene has not been well
explored. Here we employ the value of pentacene as a representative case,
while the interactions between the excited states are chosen to be
$\Delta_{23}=\Delta_{32}=1.23\omega$ for the pair 1 and
$\Delta_{23}=\Delta_{32}=0.12\omega$ for the pair 2, respectively. We consider
the room temperature case ($T=300$K) with the value of
$\beta\hbar\omega=0.382$ and the low temperature case (T=10K) with the value
of $\beta\hbar\omega=11.45$.

The system-bath interactions for the local and nonlocal modes are defined in
the matrix form as

$\displaystyle\hat{V}^{i}=\begin{bmatrix}-0.5&0&0\\\ 0&0.5&0\\\
0&0&0.5\end{bmatrix},$ (11)

and

$\displaystyle\hat{V}^{ij}=\begin{bmatrix}0&0&0\\\ 0&0&1\\\
0&1&0\end{bmatrix}.$ (12)

The diagonal and off-diagonal elements of the system-bath coupling relate to
the local and nonlocal dynamic disorder mechanisms. Although the SDFs of local
and nonlocal baths consist of two or three oscillator modes, here we consider
one primary mode to simplify the model and analysis. Then, the coupling
strength between the local/nonlocal oscillator and the bath, which acts as the
inverse noise correlation time of the BSD, is fixed to $\gamma=0.05\omega$.
The initial condition of the system at $t=0$ is set by $\rho_{22}(0)=1$ and
$\rho_{11}(0)=\rho_{33}(0)=0$. We simulate the time evolution of the reduced
density matrix by numerically integrating the HEOM with respect to the time
with a time step of ${0.001}/{\omega}$, using the fourth-order Runge-Kutta
method. The depth and the truncation number of hierarchy are chosen to $N=65$
with $K=0$ for the high temperature case and $N=12$ with $K=3$ for the low
temperature case, respectively.

Figure 3:  Time-evolution of the reduced density matrix elements
$\rho_{11}(t)$ (blue), $\rho_{22}(t)$ (orange), and $\rho_{33}(t)$ (green) at
$T=300K$ for different parameter values of local and nonlocal oscillator bath.
Here, (a) and (b) represent the pair 1 case and (c) represents the pair 2
case. The very early stage of the time evolution in each figure is presented
in Figs. 5(a)-(c) in Appendix B.

The time-evolution of the reduced density matrix for different sets of bath
modes at room temperature is shown in Fig. 3. For the local/nonlocal baths, we
consider (a)
$(\Omega_{l}^{1},~{}\lambda_{l}^{1})=(160,~{}0.032\omega)$/$(\Omega_{nl}^{1},~{}\lambda_{nl}^{1})=(50,~{}0.19\omega)$
and (b)
$(\Omega_{l}^{1},~{}\lambda_{l}^{1})=(475,~{}0.005\omega)$/$(\Omega_{nl}^{1},~{}\lambda_{nl}^{1})=(50,~{}0.19\omega)$
for the pair 1, while (c)
$(\Omega_{l}^{2},~{}\lambda_{l}^{2})=(475,~{}0.005\omega)$/$(\Omega_{nl}^{2},~{}\lambda_{nl}^{2})=(50,~{}0.137\omega)$
for the pair 2, with the unit of cm-1.

As illustrated in Fig. 3(a), the local mode of lower frequency slightly
enhances the exciton relaxation at a longer time period in comparison to that
of higher frequency presented in Fig. 3(b), because the coupling strength of
the local mode $\lambda_{l}$ is larger in the former case. The population
relaxation in Fig. 3(b) is also slower than that in Fig. 3(c), because, in the
case of pair 1, the exciton transfer between $\big{|}2\big{\rangle}$ and
$\big{|}3\big{\rangle}$ is favored due to the large nonlocal electron-phonon
coupling $\lambda_{nl}$ and transfer integral (excited state interaction)
$\Delta_{23}$. The very early stage of the time-evolution in each figure is
presented in Figs. 5 (a)-(c) in Appendix B, respectively, and it is shown that
this transfer is characterized by an ultrafast coherent oscillation in the
time period of $t$=1 [1/$\omega$] for pair 1 (Figs. 5(a) and 5(b)) and $t$=1.5
[1/$\omega$] for pair 2 (Fig. 5(c)). These facts indicate that the local and
nonlocal modes play an opposite role: The nonlocal mode promotes ultrafast
exciton transfer as a function of the coupling strength of the nonlocal mode,
$\lambda_{nl}$, that suppresses the exciton relaxation to the ground sate,
whereas the local mode enhances the relaxation through the localization of the
exciton as a function of the coupling strength of the local mode,
$\lambda_{l}$.

It should be mentioned that the small oscillatory motion of the excited state
population observed at approximately $t$=20 [1/$\omega$] in Figs. 3(a) and
3(b) are numerical errors that arise from an insufficient depth of hierarchy
elements $N$ due to the computational limitation of the HEOM calculation. This
can be suppressed by increasing $N$, although the calculations become
computationally more expensive.

Figure 4:  Time evolution of the reduced density matrix elements
$\rho_{11}(t)$ (blue), $\rho_{22}(t)$ (orange), and $\rho_{33}(t)$ (green) at
$T=10K$ for (a) pair 1 and (b) pair 2 with the same set of the bath parameters
presented in Figs. 3(b) and 3(c), respectively. The very early stage of the
time-evolution in each figure is presented in Figs. 6(a) and 6(b) in Appendix
B.

Next we discuss the low temperature cases ($T$=10K) of Figs. 3(b) and 3(c) in
Figs. 4(a) and 4(b), respectively. The initial stage of Figs. 4(a) and 4(b)
are separately plotted in Figs. 6(a) and 6(b) in Appendix B, respectively. We
observe ultrafast coherent oscillations between the excited states in the time
period of $t$=1.2 [1/$\omega$] for the pair 1 and $t$=7.5 [1/$\omega$] for the
pair 2, respectively. The frequency of the ultrafast oscillation in the low
temperature case is slower than that in the room temperature case, and thus
the population exchange between $\big{|}2\big{\rangle}$ and
$\big{|}3\big{\rangle}$ toward the equilibrium state takes longer time. In
fact, we observe that, in the case of the pair 2, the exciton relaxation
starts at time $t$=2.5 [1/$\omega$], which is shorter than the time period of
the nonlocal oscillation. This indicates that, under the low temperature
condition, the effect of the bath through the local and nonlocal couplings is
suppressed, while the system dynamics is determined from the interplay of the
ground and the excited states coupling $\Delta_{12}$ and the excitation
coupling $\Delta_{23}$. For the pair 1 with larger $\Delta_{23}$, the exciton
transfer between $\big{|}2\big{\rangle}$ and $\big{|}3\big{\rangle}$ is
favored and the relaxation to the ground state $\big{|}1\big{\rangle}$ becomes
slower than that for the pair 2.

In this regards, although the ground state population of the thermal
equilibrium state becomes larger for longer time period at lower temperature,
the equilibration process that arises from the local and nonlocal baths is
faster in the high temperature case.In any case, when the coupling strength
between the excited states and the nonlocal coupling strength are larger, the
relaxation to the ground state become slower. Thus, large nonlocal electron-
phonon coupling seems to favors the mechanism of the exciton transfer among
the excited states even in the weak excitation coupling strength condition,
while the lifetime of the exciton is suppressed by low frequency local modes.
Because the configuration of pair 1 satisfies both the strong excitation
coupling and strong nonlocal coupling conditions, it is more favorable than
that of pair 2.

Although it was found from investigations of polaron transport process that
the nonlocal electron-phonon coupling promotes not only delocalization, but
also localization due to a bandwidth narrowing effect, in particular for the
strong coupling case at high temperatures,[55, 56, 57] the model system
employed in this work is too small and further study in a larger system needs
to be done for the verification.

##  IV Conclusion

In this work, we introduced Holstein-Peierls + bath Hamiltonian in order to
investigate the effects of local and nonlocal electron-phonon interactions for
the exciton dynamics of organic molecules in photovoltaic cells. We employed
the HEOM formalism to simulate the dynamics of excitons described by this
model under various physical conditions for local and nonlocal environmental
modes without employing approximations, most notably the RWA and MFP
approximations. It should be noted that the situations that we considered here
are classified as non-perturbative and non-Markovian system-bath interaction
case in open quantum dynamics theory.

The calculated results demonstrate the effects of dynamic disorder on the
order of the reported time scale of the exciton processes that are
characterized by the presence of ultrafast coherent oscillations, in
particular at low temperatures. We found that a stronger nonlocal electron-
phonon coupling favors the exciton transfer between excited states and may be
a key to enhancing the efficiency of exciton transfer processes even in
materials with low excitation couplings. Moreover, the effect seems to be
favored at high temperature, while the number of low frequency local bath
modes should be suppressed.

In the present work, we limited our analysis to a small system with a specific
parameter set focusing on a role of local and nonlocal system-phonon
interactions. Because the exciton transfer is a long range effect, an
extension of the present investigation to a larger system should be realized
to provide better insight in this topic, specially regarding further
discussions between the delocalization/localization mechanisms. We plan to
investigate this direction in a future study.

###### Acknowledgements.

The financial support from The Kyoto University Foundation is acknowledged.
M.C. is supported by the Japanese Government (MEXT) Postgraduate Scholarships.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Reduced Hierarchy Equations of Motion for Holstein-Peierls
Model

We consider the density matrix of the total system with the factorized initial
condition. The reduced density matrix element is then expressed in the path
integral form as

$\displaystyle\rho(z,z^{\prime},t)$ $\displaystyle=$
$\displaystyle\iint_{z(t)=z}^{z^{\prime}(t)=z^{\prime}}\mathcal{D}[z,z^{\prime}]e^{iS_{0}[z,z^{\prime}]/\hbar}$
(13)
$\displaystyle\times\rho(z_{0},z^{\prime}_{0},t_{0})e^{-iS_{0}[z,z^{\prime}]/\hbar}\mathcal{F}[z,z^{\prime},t],$

where $\iint_{z(t)=z}^{z^{\prime}(t)=z^{\prime}}\mathcal{D}[z,z^{\prime}]$
represents the functional integral of a set of Grassmann variables, which
describe the system states, $\rho(z_{0},z^{\prime}_{0},t_{0})$ is the
factorized system part density, $S_{0}$ is the action for the system
Hamiltonian (Lagrangian), and $\mathcal{F}[z,z^{\prime},t]$ is the called
Feynman-Vernon influence functional that includes all the information
regarding the bath and the system bath-interaction.[58, 59, 60, 61]

Considering the time dependent kernels $iL_{1}(t)=i\int_{0}^{\infty}d\omega
J(\omega)\sin(\omega t)$ and $L_{2}(t)=\int_{0}^{\infty}d\omega
J(\omega)\coth(\beta\hbar\omega/2)\cos(\omega t)$ that correspond to the
fluctuation and dissipation effects, respectively, introduced by the bath,[58,
62, 63] the influence functional can be written, evaluating the contour
integrals for the Brownian spectral distribution as[39, 40, 64, 65]

$\displaystyle\mathcal{F}[z,z^{\prime},t]$ $\displaystyle=$ $\displaystyle
exp\Bigg{\\{}-\int_{t_{0}}^{t}ds\int_{t_{0}}^{s}duV^{\times}(s)$ (14)
$\displaystyle\times\bigg{[}\Theta_{-}(u)e^{-\big{(}\frac{\gamma}{2}-i\zeta\big{)}(s-u)}+\Theta_{+}(u)e^{-\big{(}\frac{\gamma}{2}+i\zeta\big{)}(s-u)}\bigg{]}\Bigg{\\}}$
$\displaystyle\times\exp\Bigg{[}\int_{t_{0}}^{t}ds\int_{t_{0}}^{s}duV^{\times}(s)\sum_{k=1}^{\infty}\Psi_{k}(u)\nu_{k}e^{-\nu_{k}(s-u)}\Bigg{]}$
$\displaystyle\times\exp\bigg{[}\int_{t_{0}}^{t}ds\Xi(s)\bigg{]},$

where

$\displaystyle\Theta_{\pm}=\frac{\lambda\omega_{0}^{2}}{2\zeta}\Bigg{\\{}\pm\coth\bigg{[}\frac{\beta\hbar}{2}\Big{(}i\frac{\gamma}{2}\mp\zeta\Big{)}\bigg{]}V^{\times}(u)\mp
V^{\circ}(u)\Bigg{\\}},$ (15)
$\displaystyle\Psi_{k}(u)=\frac{4\lambda\gamma\omega_{0}^{2}}{\beta\hbar}\frac{\nu_{k}}{(\omega_{0}^{2}+\nu_{k}^{2})^{2}-\gamma^{2}\nu_{k}^{2}}V^{\times}(u),$
(16) $\displaystyle\Xi(s)=V^{\times}(s)\sum_{k=K+1}^{\infty}\Psi_{k}(s).$ (17)

Here, $\zeta=\sqrt{\omega_{0}^{2}-\gamma^{2}/4}$ for $\gamma<2\omega_{0}$ and
$\nu_{k}=(2\pi/\beta h)k$ are the Matsubara frequencies, $\beta={1}/{k_{B}T}$
is the inverse temperature, and $k_{B}$ is the Boltzmann constant.

The commutator and anticommutator of the interaction operator $V$ are
expressed as $V^{\times}(t)=V(z(t))-V(z^{\prime}(t))$ and
$V^{\circ}(t)=V(z(t))+V(z^{\prime}(t))$, and the electron-phonon couplings are
introduced in the spectral distribution. The local and nonlocal transition
rates are determined from the diagonal and off-diagonal elements of $V$,
respectively.

Thanks to the exponential function form of the kernels, any higher-order time
derivative of the correlation functions is also expressed in terms of
exponential functions. These terms are called the auxiliary density matrices
(ADMs) in the HEOM formalism and are necessary to describe the non-Markovian
effects of the environment. The HEOM for the Brownian distribution can be
derived by calculating the time derivative on the left and right side of the
factorized system density and the influence functional.[64, 65] Summing over
all bath modes results in Eq. (9).

For a high number of hierarchy elements, the hierarchy can be truncated by the
terminator

$\displaystyle\frac{\partial}{\partial
t}\hat{\rho}_{j_{1}…j_{k}}^{(n,m)}(t)=-\bigg{[}\frac{i}{\hbar}\hat{H}_{el}^{\times}-i(n-m)\zeta-\hat{\Xi}\bigg{]}\hat{\rho}_{j_{1}…j_{k}}^{(n,m)}(t).$
(18)

The total number of hierarchy elements can be evaluated as
$L=(N+M+K+1)/(K+1)!(N+M)!$ with the total number of terminator elements
$L_{term}=(N+M+K)/K!(N+M)!$, where $N$ is the hierarchy depth for $n,m_{+}$
and $m_{-}$, and $M$ is the increase of the hierarchy for each spectral mode
($M=2$ for $m_{+}$ and $m_{-}$ in the case of the Brownian oscillator). In the
HEOM formalism, each member of the hierarchy is coupled to the lower and
higher terms, and the value for the initial density element is the exact
solution of the total Hamiltonian defined by $\rho_{0…0}^{(0,0)}(t)$, which
includes all the system-bath interactions.

##  Appendix B The ultrafast oscillations of the density matrix elements in
the very early stage

Here, we replot the very early stage of the time-evolution ($t\leq 10$
[1/$\omega$] of the reduced density matrix elements presented in Figs. 3 and
4, respectively. In each figure, we observe ultrafast coherent oscillations
between $\big{|}2\big{\rangle}$ and $\big{|}3\big{\rangle}$. As it is
illustrated in Figs. 5 (b) and 5 (c) and Figs. 6 (a) and 6 (b), the ultrafast
oscillation is enhanced for larger $\lambda_{nl}$, which indicates that
nonlocal electron-phonon interaction that relates to the dynamics disorder
enhances the transition between the excited states.

Figure 5: Very early stage of the time-evolution of the reduced density matrix
elements $\rho_{11}(t)$ (blue), $\rho_{22}(t)$ (orange), and $\rho_{33}(t)$
(green) at $T=300$K presented in Figs. 3 (a)-(c) as Figs. (a)-(c),
respectively, are illustrated. Here, (a) and (b) correspond to the case of
pair 1 and (c) corresponds to the case of pair 2. Figure 6:  Very early stage
of the time-evolution of the reduced density matrix elements $\rho_{11}(t)$
(blue), $\rho_{22}(t)$ (orange), and $\rho_{33}(t)$ (green) at $T=10$K
presented in Figs. 4 (a) and (b) as Figs. (a) and (b), respectively, are
illustrated.

## References

  * Horowitz [1998] G. Horowitz, “Organic field-effect transistors,” Adv. Mater. 10, 365–377 (1998).
  * Dimitrakopoulos and Malenfant [2002] C. D. Dimitrakopoulos and P. R. L. Malenfant, “Organic thin film transistors for large area electronics,” Adv. Mater. 14, 99–117 (2002).
  * Zhang  _et al._ [2018] J. Zhang, H. S. Tan, X. Guo, A. Facchetti, and H. Yan., “Material insights and challenges for non-fullerene organic solar cells based on small molecular acceptors,” Nature Energy 3, 720–731 (2018).
  * Scharber  _et al._ [2006] M. C. Scharber, D. Mühlbacher, M. Koppe, P. Denk, C. Waldauf, A. J. Heeger, and C. J. Brabec, “Design rules for donors in bulk-heterojunction solar cells—towards 10% energy-conversion efficiency,” Adv. Mater. 18, 789–794 (2006).
  * Meng  _et al._ [2018] L. Meng, Y. Zhang, X. Wan, C. Li, X. Zhang, Y. Wanga, X. Ke, Z. Xiao, L. Ding, R. Xia, H.-L. Yip, Y. Cao, and Y. Chen, “Organic and solution-processed tandem solar cells with 17.3% efficiency,” Science 361, 1094–1098 (2018).
  * Brédas  _et al._ [2009] J.-L. Brédas, J. E. Norton, J. Cornil, and V. Coropceanu, “Molecular understanding of organic solar cells: The challenges,” Acc. Chem. Res. 42, 1691–1699 (2009).
  * Heeger [2014] A. J. Heeger, “25th anniversary article: Bulk heterojunction solar cells: Understanding the mechanism of operation,” Adv. Mater. 26, 10–28 (2014).
  * Devizi  _et al._ [2010] A. Devizi, K. Meerholz, D. Hertel, and V. Gulbinas, “Hierarchical charge carrier motion in conjugated polymers,” Chem. Phys. Lett. 498, 302–306 (2010).
  * Oviedo-Casado, Urbina, and Prior [2017] S. Oviedo-Casado, A. Urbina, and J. Prior, “Magnetic field enhancement of organic photovoltaic cells performance,” Sci. Rep. 7, 4297 (2017).
  * Botiz and Stigelin [2014] I. Botiz and N. Stigelin, “Influence of molecular conformations and microstructure on the optoelectronic properties of conjugated polymers,” Materials 7, 2273–2300 (2014).
  * Vandewal, Himmelberger, and Salleo [2013] K. Vandewal, S. Himmelberger, and A. Salleo, “Structural factors that affect the performance of organic bulk heterojunction solar cells,” Macromolecules 46, 6379–6387 (2013).
  * Beaujuge and Frechet [2011] P. M. Beaujuge and J. M. J. Frechet, “Molecular design and ordering effects in $\pi$-functional materials for transistor and solar cell applications,” J. Am. Chem. Soc. 133, 20009–20029 (2011).
  * Athanasopoulos  _et al._ [2017] S. Athanasopoulos, S. Tscheuschner, H. Bässler, and A. Köhler, “Efficient charge separation of cold charge-transfer states in organic solar cells through incoherent hopping,” J. Phys. Chem. Lett. 8, 2093–2098 (2017).
  * Picon, Bussac, and Zuppiroli [2007] J.-D. Picon, M. N. Bussac, and L. Zuppiroli, “Quantum coherence and carriers mobility in organic semiconductors,” Phys. Rev. B 75, 235106 (2007).
  * Birech  _et al._ [2014] Z. Birech, M. Schwoerer, T. Schmeiler, J. Pflaum, and H. Schwoerer, “Ultrafast dynamics of excitons in tetracene single crystals,” J. Chem. Phys. 140, 114501 (2014).
  * Marciniak  _et al._ [2007] H. Marciniak, M. Fiebig, M. Huth, S. Schiefer, B. Nickel, F. Selmaier, and S. Lochbrunner, “Ultrafast exciton relaxation in microcrystalline pentacene films,” Phys. Rev. Lett. 99, 176402 (2007).
  * Elton and Wang [2016] J. G. S. Elton and W. L. Wang, “Ultrafast charge-transfer in organic photovoltaic interfaces: geometrical and functionalization effects,” Nanoscale 8, 15902–15910 (2016).
  * Zheng  _et al._ [2017] Z. Zheng, N. R. Tummala, Y.-T. Fu, V. Coropceanu, and J.-L. Brédas, “Charge-transfer states in organic solar cells: Understanding the impact of polarization, delocalization, and disorder,” ACS Appl. Mater. Interfaces 9, 18095–18102 (2017).
  * Chen, Coropceanu, and Brédas [2018] X.-K. Chen, V. Coropceanu, and J.-L. Brédas, “Assessing the nature of the charge-transfer electronic states in organic solar cells,” Nat. Commun. 9, 5295 (2018).
  * Li, Coropceanua, and Brédas [2013] Y. Li, V. Coropceanua, and J.-L. Brédas, “Nonlocal electron-phonon coupling in organic semiconductor crystals: The role of acoustic lattice vibrations,” J. Chem. Phys. 138, 204713 (2013).
  * Bakulin  _et al._ [2012] A. A. Bakulin, A. Rao, V. G. Pavelyev, P. H. M. van Loosdrecht, M. S. Pshenichnikov, D. Niedzialek, J. Cornil, D. Beljonne, and R. H. Friend, “The role of driving energy and delocalized states for charge separation in organic semiconductors,” Science 335, 1340–1344 (2012).
  * Gélinas  _et al._ [2014] S. Gélinas, A. Rao, A. Kumar, S. L. Smith, A. W. Chin, J. Clark, T. S. van der Poll, G. C. Bazan, and R. H. Friend, “Ultrafast long-range charge separation in organic semiconductor photovoltaic diodes,” Science 343, 512–516 (2014).
  * Collini and Scholes [2009a] E. Collini and G. D. Scholes, “Coherent intrachain energy migration in a conjugated polymer at room temperature,” Science 323, 369–373 (2009a).
  * Collini and Scholes [2009b] E. Collini and G. D. Scholes, “Electronic and vibrational coherences in resonance energy transfer along meh-ppv chains at room temperature,” J. Phys. Chem. A 113, 4223–4241 (2009b).
  * Sio and Lienau [2017] A. D. Sio and C. Lienau, “Vibronic coupling in organic semiconductors for photovoltaics,” Phys. Chem. Chem. Phys. 19, 18813–18830 (2017).
  * Grancini  _et al._ [2013] G. Grancini, M. Maiuri, D. Frazzil, A. Petrozzal, H. Egelhaaf, D. Brida, G. Cerullo, and G. Lanzani, “Hot exciton dissociation in polymer solar cells,” Nat. Mater. 12, 29–33 (2013).
  * Tully [2012] J. C. Tully, “Perspective: Nonadiabatic dynamics theory,” J. Chem. Phys. 137, 22A301 (2012).
  * Hannewald  _et al._ [2004] K. Hannewald, V. M. Stojanović, J. M. T. Schellekens, P. A. Bobbert, G. Kresse, and J. Hafner, “Theory of polaron bandwidth narrowing in organic molecular crystals,” Phys. Rev. B 69, 075211 (2004).
  * Holstein [1959] T. Holstein, “Studies of polaron motion: Part i. the molecular-crystal model,” Ann. Phys. 8, 325–342 (1959).
  * Peierls [1930] R. E. Peierls, “Zur theorie der elektrischen und thermischen leitfähigkeit von metallen.” Ann. Phys. 396, 121–148 (1930).
  * Duan  _et al._ [2019] H.-G. Duan, P. Nalbach, R. J. D. Miller, and M. Thorwart, “Ultrafast energy transfer in excitonically coupled molecules induced by a nonlocal peierls phonon,” J. Phys. Chem. Lett. 10, 1206–1211 (2019).
  * Florencio Jr. and Goodman [1986] J. Florencio Jr. and B. Goodman, “Mean-field approximation for the molecular polaron with q2 coupling,” Phys. Rev. B 34, 3645–3648 (1986).
  * Zusman [1980] L. D. Zusman, “Outer-sphere electron transfer in polar solvents,” Chem. Phys. 49, 295–304 (1980).
  * Tachiya [1989] M. Tachiya, “Relation between the electron-transfer rate and the free energy change of reaction,” J. Phys. Chem. 93, 7050–7052 (1989).
  * Wolynes [1987] P. Wolynes, “Dissipation, tunneling, and adiabaticity criteria for curve crossing problems in the condensed phase,” J. Chem. Phys. 86, 1957–1966 (1987).
  * Garg, Onuchic, and Ambegaokar [1985] A. Garg, J. N. Onuchic, and V. Ambegaokar, “Effect of friction on electron transfer in biomolecules,” J. Chem. Phys. 83, 4491–4503 (1985).
  * Tanimura [2020] Y. Tanimura, “Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (heom),” J. Chem. Phys. 153, 020901 (2020).
  * Ishizaki and Tanimura [2005] A. Ishizaki and Y. Tanimura, “Quantum dynamics of system strongly coupled to low-temperature colored noise bath: Reduced hierarchy equations approach,” J. Phys. Soc. Jpn. 74, 3131–3134 (2005).
  * Tanimura and Mukamel [1994] Y. Tanimura and S. Mukamel, “Optical stark spectroscopy of a brownian oscillator in intense fields,” J. Phys. Soc. Jpn. 63, 66–77 (1994).
  * Tanaka and Tanimura [2009] M. Tanaka and Y. Tanimura, “Quantum dissipative dynamics of electron transfer reaction system: Nonperturbative hierarchy equations approach,” J. Phys. Soc. Jpn. 78, 073802 (2009).
  * Ishizaki and Fleming [2009] A. Ishizaki and G. R. Fleming, “Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature,” Proc. Natl. Acad. Sci. U.S.A. 106, 17255 (2009).
  * Kreisbeck and Kramer [2012] C. Kreisbeck and T. Kramer, “Long-lived electronic coherence in dissipative exciton dynamics of light-harvesting complexes,” J. Phys. Chem. Lett. 3, 2828 (2012).
  * Fujihashi, Fleming, and Ishizaki [2015] Y. Fujihashi, G. R. Fleming, and A. Ishizaki, “Impact of environmentally induced fluctuations on quantum mechanically mixed electronic and vibrational pigment states in photosynthetic energy transfer and 2d electronic spectra,” J. Chem. Phys. 142, 212403 (2015).
  * Sakamoto and Tanimura [2017] S. Sakamoto and Y. Tanimura, “Long-lived electronic coherence in dissipative exciton dynamics of light-harvesting complexes,” J. Phys. Chem. Lett. 8, 5390 (2017).
  * Dijkstra and Tanimura [2010] A. G. Dijkstra and Y. Tanimura, “Correlated fluctuations in the exciton dynamics and spectroscopy of dna,” New J. Phys. 12, 055005 (2010).
  * Tanimura [1990] Y. Tanimura, “Nonperturbative expansion method for a quantum system coupled to a harmonic-oscillator bath,” Phys. Rev. A 41, 6676–6687 (1990).
  * Chen, Zhao, and Tanimura [2015] L. Chen, Y. Zhao, and Y. Tanimura, “Dynamics of a one-dimensional holstein polaron with the hierarchical equations of motion approach,” J. Phys. Chem. Lett. 6, 3110–3115 (2015).
  * Dunn, Tempelaar, and Reichman [2019] I. S. Dunn, R. Tempelaar, and D. R. Reichman, “Removing instabilities in the hierarchical equations of motion: Exact and approximate projection approaches,” J. Chem. Phys. 150, 184109 (2019).
  * Zhugayevych and Tretiak [2015] A. Zhugayevych and S. Tretiak, “Theoretical description of structural and electronic properties of organic photovoltaic materials,” Rev. Phys. Chem. 66, 305–330 (2015).
  * Tanimura [2006] Y. Tanimura, “Stochastic liouville, langevin, fokker–planck, and master equation approaches to quantum dissipative systems,” J. Phys. Soc. Jpn. 75, 082001 (2006).
  * Xie, Santana-Bonilla, and Troisi [2018] X. Xie, A. Santana-Bonilla, and A. Troisi, “Nonlocal electron–phonon coupling in prototypical molecular semiconductors from first principles,” J. Chem. Theory Comput. 14, 3752–3762 (2018).
  * Hu  _et al._ [2020] W. Hu, K. Sun, Q. Xu, L. Chen, and Y. Zhao, “Ultrafast dynamics in rubrene and its spectroscopic manifestation,” J. Chem. Phys. 153, 174105 (2020).
  * Sun  _et al._ [2020] K. Sun, Q. Xu, L. Chen, M. F. Gelin, and Y. Zhao, “Temperature effects on singlet fission dynamics mediated by a conical intersection,” J. Chem. Phys. 153, 194106 (2020).
  * Fuemmeler  _et al._ [2016] E. G. Fuemmeler, S. N. Sanders, A. B. Pun, E. Kumarasamy, T. Zeng, K. Miyata, M. L. Steigerwald, X.-Y. Zhu, M. Y. Sfeir, L. M. Campos, and N. Ananth, “A direct mechanism of ultrafast intramolecular singlet fission in pentacene dimers,” ACS Cent. Sci. 2, 316–324 (2016).
  * Chen  _et al._ [2011] D. Chen, J. Ye, H. Zhang, and Y. Zhao, “On the munn-silbey approach to polaron transport with off-diagonal coupling and temperature-dependent canonical transformations,” J. Phys. Chem. B 115, 5312–5321 (2011).
  * Zhao, Brown, and Lindenberg [1994] Y. Zhao, D. W. Brown, and K. Lindenberg, “On the munn–silbey approach to nonlocal exciton–phonon coupling,” J. Chem. Phys. 100, 2335–2345 (1994).
  * Zhou  _et al._ [2016] N. Zhou, L. Chen, Z. Huang, K. Sun, Y. Tanimura, and Y. Zhao, “Fast, accurate simulation of polaron dynamics and multidimensional spectroscopy by multiple davydov trial states,” J. Phys. Chem. A 120, 1562–1576 (2016).
  * Feynman and Vernon [1963] R. P. Feynman and F. L. Vernon, “The theory of a general quantum system interacting with a linear dissipative system,” Ann. Phys. 271, 547–607 (1963).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, “Dynamics of the dissipative two-state system,” Rev. Mod. Phys. 59, 1–85 (1987).
  * Tanimura [2014] Y. Tanimura, “Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities,” J. Chem. Phys. 141, 044114 (2014).
  * Tanimura [2015] Y. Tanimura, “Real-time and imaginary-time quantum hierarchal fokker-planck equations,” J. Chem. Phys. 142, 144110 (2015).
  * Tanimura and Kubo [1989] Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly gaussian-markoffian noise bath,” J. Phys. Soc. Jpn. 58, 101–114 (1989).
  * Tanimura and Wolynes [1991] Y. Tanimura and P. G. Wolynes, “Quantum and classical fokker-planck equations for a gaussian-markovian noise bath,” Phys. Rev. A 43, 4131–4142 (1991).
  * Tanaka and Tanimura [2010] M. Tanaka and Y. Tanimura, “Multistate electron transfer dynamics in the condensed phase: Exact calculations from the reduced hierarchy equations of motion approach,” J. Chem. Phys. 132, 214502 (2010).
  * Tanimura [2012] Y. Tanimura, “Reduced hierarchy equations of motion approach with drude plus brownian spectral distribution: Probing electron transfer processes by means of two-dimensional correlation spectroscopy,” J. Chem. Phys. 137, 22A550 (2012).