# Driving Force and Nonequilibrium Vibronic Dynamics
in Charge Separation of Strongly Bound Electron-Hole Pairs

Alejandro D. Somoza<EMAIL_ADDRESS>Insitut für Theoretische Physik
and IQST, Albert-Einstein-Allee 11, Universität Ulm, D-89069 Ulm, Germany
German Aerospace Center (DLR), Institute of Engineering Thermodynamics,
Wilhelm-Runge Straße 10, 89081 Ulm  Nicola Lorenzoni  Insitut für Theoretische
Physik and IQST, Albert-Einstein-Allee 11, Universität Ulm, D-89069 Ulm,
Germany  James Lim  Insitut für Theoretische Physik and IQST, Albert-Einstein-
Allee 11, Universität Ulm, D-89069 Ulm, Germany  Susana F. Huelga  Insitut für
Theoretische Physik and IQST, Albert-Einstein-Allee 11, Universität Ulm,
D-89069 Ulm, Germany  Martin B. Plenio<EMAIL_ADDRESS>Insitut für
Theoretische Physik and IQST, Albert-Einstein-Allee 11, Universität Ulm,
D-89069 Ulm, Germany

###### Abstract

Electron-hole pairs in organic photovoltaics dissociate efficiently despite
their Coulomb-binding energy exceeding thermal energy at room temperature. The
electronic states involved in charge separation couple to structured
vibrational environments containing multiple underdamped modes. The non-
perturbative simulations of such large, spatially extended electronic-
vibrational (vibronic) systems remains an outstanding challenge. Current
methods bypass this difficulty by considering effective one-dimensional
Coulomb potentials or unstructured environments. Here we extend and apply a
recently developed method for the non-perturbative simulation of open quantum
systems to the dynamics of charge separation in one, two and three-
dimensional donor-acceptor networks. This allows us to identify the precise
conditions in which underdamped vibrational motion induces efficient long-
range charge separation. Our analysis provides a comprehensive picture of
ultrafast charge separation by showing how different mechanisms driven either
by electronic or vibronic couplings are well differentiated for a wide range
of driving forces and how entropic effects become apparent in large vibronic
systems. These results allow us to quantify the relative importance of
electronic and vibronic contributions in organic photovoltaics and provide a
toolbox for the design of efficient charge separation pathways in artificial
nanostructures.

††preprint: PREPRINT

##  I Introduction

When a solar cell made out of an inorganic semiconductor like silicon is
exposed to light, electrons can be readily extracted from the valence band to
the conduction band and then captured at the electrodes. If, however, light is
absorbed by carbon-based materials, photons produce strongly bound electron-
hole pairs called excitons, which are collective optical excitations that may
be delocalized across several molecular units [1]. Excitons are charge
neutral, namely the electron and the hole occupy, respectively, the lowest
unoccupied molecular orbital (LUMO) and the highest occupied molecular orbital
(HOMO) of the same molecular unit, and require dissociation in order to
produce a current [2]. In contrast, the charge-transfer (CT) states describe
partially separated electron-hole pairs where an electron and a hole occupy,
respectively, the LUMO and HOMO levels that belong to different molecular
sites. The transfer from exciton to CT states is thus suitable to describe the
dynamics of electron transfer [3, 4]. In photosynthetic organisms excitons are
split in pigment-protein complexes called reaction centers [5, 6, 7]. In
organic photovoltaics (OPV), blends of materials with different electron
affinities are used to provide an energetic landscape that is favourable to
charge separation at the interface [8]. These devices exhibit ultrafast, long-
range charge separation with high quantum efficiencies [9, 10, 11]. This means
that a large proportion of absorbed photons produces excitons or strongly
bound CT states that are successfully dissociated. Some of these electron-hole
pairs however thermalize towards the lowest-energy CT state localized at the
interface, which is for this reason considered an energetic trap that leads to
non-radiative electron-hole recombination [12, 13, 14, 15], as schematically
shown in Fig.1(a). This localization process is predominantly mediated by
high-frequency vibrational modes that can bridge the energy gap between high-
lying exciton/CT states and the lowest-energy interfacial CT state. The energy
loss associated to this process is typically larger than $0.6\,$eV per photon
[16, 17, 18], leading to a low power conversion efficiency in OPV with respect
to their inorganic counterparts that results in a small open circuit voltage
[19].

Figure 1: Coulomb potential and vibrational environments. (a) Schematic
representation of a one-dimensional chain consisting of a donor and $(N-1)$
acceptors. The Coulomb-binding energy of electron and hole is modelled by
$\Omega_{k}=-V/k$ with $V=0.3\,$eV for $k\geq 1$. The energy-gap between
exciton and interfacial CT states is defined as driving force
$\Delta=\Omega_{0}-\Omega_{1}$. (b) Vibrational environments consist of low-
frequency phonon baths with room temperature energy scales ($k_{B}T\approx
200\,{\rm cm}^{-1}$) and high-frequency vibrational modes. In this work, the
low-frequency phonons are modelled by an Ohmic spectral density ${\cal
J}_{l}(\omega)$ with an exponential cutoff, while the high-frequency modes are
described by Lorentzian spectral densities ${\cal J}_{h}(\omega)$ centered at
vibrational frequency $\omega_{h}=1200\,{\rm cm}^{-1}$ (or $1500\,{\rm
cm}^{-1}$). The vibrational damping rate of the high-frequency modes is taken
to be $\gamma=(50\,{\rm fs})^{-1}$ or $(500\,{\rm fs})^{-1}$, as shown in red
and blue, respectively, to investigate the role of non-equilbrium vibrational
motion in long-range charge separation.

Although energetically costly, dissociation of strongly bound electron-hole
pairs [10, 20, 11] takes place despite the much lower thermal energy at room
temperature. The energy of bound CT states is largely dependent on the offset
between the LUMO of the acceptor and the HOMO of the donor [21, 22, 23].
Fixing the acceptor and employing different donor materials (or vice versa) is
a popular strategy to investigate the energetics at the interface and achieve
a high voltage, small energy losses and sufficient photocurrent density [24,
25, 26, 27, 28, 19]. Surprisingly, some of these blends show ultrafast and
efficient exciton dissociation despite having small or no apparent driving
force [29, 30, 31, 32, 15, 18, 33, 34, 35]. The driving force is a crucial
parameter in charge separation and refers to the energy difference between
exciton and interfacial CT state (see $\Delta$ in Fig.1(a)). Hybridization
between exciton and CT states has been thought to be behind the successful
ultrafast charge separation of these promising materials, which are often
based on small molecules (oligomers) with acceptor-donor-acceptor structures
that have reached power conversion efficiencies of up to $17\,\%$ [36, 31, 15,
37]. This represents an astonishing 50$\,\%$ increase in the state-of-the-art
performance of organic photovoltaics in less than a decade. From finite
molecular clusters to periodic molecular solids, ultrafast long-range charge
separation has appeared across a wide variety of photovoltaic platforms, but
the underlying mechanism has not been understood fully, leading some to
advocate for a deeper analysis of charge separation processes [38, 39].

Some experimental studies rule out thermal activation as an important
mechanism for charge separation in a large number of photovoltaic devices [40,
41, 42]. In contrast, the vibronic coupling to underdamped vibrational modes
is presumed to enable coherent charge separation [43, 44, 45, 46, 47, 48, 49],
which requires non-perturbative simulation tools for a reliable description of
the vibronic interaction between exciton/CT states and molecular vibrations.
However, in many theoretical studies on the charge separation in extended
systems, a broad and unstructured environmental spectral density has been
considered [50, 51, 52] to reduce simulation costs, neglecting the ubiquitous
presence of underdamped vibrational modes in organic molecules and their role
in charge separation. In addition, the non-Markovian vibronic effects proposed
to suppress the localization of electron-hole pairs at the interfaces, e.g.
suggested in Ref.[52], are found to be well described by a Markovian quantum
master equation, as shown in Appendix A, due to weak vibronic coupling
strength and no underdamped modes considered in simulations. This indicates
that a vibronic mechanism inferred solely on the basis of non-perturbative
numerical results without the subsequent formulation of an accurate physical
mechanism may lead to ambiguities in the interpretation of the underlying
mechanism. Some first-principles numerical methods have been employed to
simulate vibronic charge separation [53, 54, 55, 56, 57, 58, 59, 60], where
underdamped vibrational modes are considered. However, the interpretation of
simulated results is a non-trivial issue here. For instance, in Ref.[53], an
effective one-dimensional Coulomb potential is considered where electron-hole
binding energy is assumed to be reduced by instantaneous electron
delocalization in three-dimensional acceptor aggregates and as a result the
electronic coupling being responsible for a hole transfer becomes larger in
magnitude than the detunings in energy levels of the effective potential. In
Appendix B, we show how, in this case, deactivating completely the vibrational
environment has little impact in charge separation dynamics. This leads us to
conclude that the ultrafast long-range charge separation observed in Ref.[53]
is not necessarily enhanced by vibronic couplings, but merely induced by the
weak Coulomb-binding energy. Other theoretical studies have focused on
intermolecular modes as the relevant vibrations behind charge separation [61],
while sometimes, intramolecular modes are attributed a hampering role [62].
This is, as we will demonstrate, in sharp contrast to our findings, as
intramolecular modes can induce both effects.

Given the heterogeneity of donor-acceptor materials and model parameters
employed across the literature, we aim to discern the underlying mechanisms of
charge separation dynamics as a function of the driving force and the
structure of vibrational environments based on non-perturbative simulations
and detailed reduced model analysis. We determine under what conditions
underdamped vibrational motion induces efficient long-range charge separation
in the presence of strong Coulomb-binding energy $V\sim 0.3\,$eV. To this end
we consider one, two and three-dimensional donor-acceptor networks, instead of
effective one-dimensional Coulomb potentials, by using our non-perturbative
simulation method called dissipation assisted matrix product factorization
(DAMPF) [63, 64, 65, 66], to investigate how coherent vibronic couplings
promote long-range charge separation in high-dimensional multi-site systems.
We show that there are two available mechanisms for ultrafast long-range
charge separation in donor-acceptor interfaces. For low driving forces
$\Delta\sim 0.15\,$eV, the transitions between near-resonant exciton and
delocalised CT states occur on a sub-ps time scale even if vibronic couplings
are not considered. For high driving forces $\Delta\sim 0.3\,$eV, the vibronic
coupling of underdamped high-frequency vibrational modes with frequencies
$\omega_{h}\sim 0.15\,$eV induces the transitions between exciton and CT
states delocalised over multiple acceptors. Here a vibrationally cold exciton
can interact resonantly with vibrationally hot lower-energy CT states and,
subsequently, also with vibrationally cold high-energy CT states. The charge
separation process becomes significantly inefficient in this case when
vibronic couplings are ignored in simulations, hinting the genuine vibronic
effects induced by underdamped vibrational modes. For both low and high
driving forces, we demonstrate that the time scale of the charge localization
towards the donor-acceptor interfaces is determined by the lifetime of the
high-frequency vibrational modes, as strongly damped modes promote the
transitions to the lowest-energy interfacial CT state. These results
demonstrate that experimentally measured long-lived vibrational and vibronic
coherences in OPV [43, 44, 45, 46, 47, 48, 49] may have a functional relevance
in charge separation processes.

Figure 2: Charge separation under different vibrational environments and
electronic eigenstate structures. (a) Time-averaged mean electron-hole
distance $\langle x\rangle_{t\leq 400\,{\rm fs}}$ is displayed as a function
of driving force $\Delta$ when vibrational environments are absent (${\cal
J}(\omega)=0)$ or only low-frequency phonon baths are present (${\cal
J}(\omega)={\cal J}_{l}(\omega)$), shown in yellow and gray, respectively.
Here we consider a linear chain consisting of a donor and nine acceptors
($N=10$ sites). (b-c) Time-averaged mean electron-hole distance when
electronic states are coupled to (b) high-frequency vibrational modes only
(${\cal J}(\omega)={\cal J}_{h}(\omega)$) or to (c) the total vibrational
environments (${\cal J}(\omega)={\cal J}_{l}(\omega)+{\cal J}_{h}(\omega)$).
Here the vibrational damping rate of the high-frequency modes with frequency
$\omega_{h}=1200\,{\rm cm}^{-1}$ is taken to be $\gamma=(50\,{\rm fs})^{-1}$
or $(500\,{\rm fs})^{-1}$, shown in red and blue, respectively. (d) Electronic
eigenstates in the absence of exciton-CT couplings where probability
distributions for finding an electron at the $k$-th acceptor are vertically
shifted depending on electronic energy-levels $E_{\alpha}$. (e-f) Electronic
eigenstates in the presence of exciton-CT couplings where driving force is
taken to be (e) $\Delta_{e}=0.15\,$eV or (f) $\Delta_{v}=0.3\,$eV. In (e),
hybrid exciton-CT states contributing to initial charge separation dynamics
are colored in blue/green. In (f), the exciton and delocalised CT states,
governing initial charge separation via a vibronic mixing, are highlighted in
blue and green, respectively. In (d-f), the probabilities for finding an
electron at donor/acceptor interface are shown in red.

##  II Model

To investigate the influence of underdamped vibrational motion on the charge
separation of strongly bound electron-hole pairs, we consider a one-
dimensional chain consisting of $N$ sites, composed of an electron donor in
contact with a chain of $(N-1)$ electron acceptors, as schematically shown in
Fig.1(a). Two and three-dimensional donor-acceptor networks will be considered
later. The electronic Hamiltonian is modeled by

$H_{e}=\sum_{k=0}^{N-1}\Omega_{k}\left|k\right>\left<k\right|+\sum_{k=0}^{N-2}J_{k,k+1}(\left|k\right>\left<k+1\right|+h.c.),$
(1)

where $h.c.$ denotes the Hermitian conjugate. Here $\left|0\right>$ denotes an
exciton state localized at the donor, while $\left|k\right>$ with $k\geq 1$ is
a CT state with an electron localized at the $k$-th acceptor. For simplicity,
we assume that the hole is fixed at the donor within the time scale of our
simulations due to its lower mobility with respect to the electron [67, 52,
50]. The energy-levels of the CT states take into account the Coulomb
attraction between electron and hole, modelled by $\Omega_{k}=-V/k$ with
$V=0.3\,{\rm eV}$ for $k\geq 1$. We take $J_{k,k+1}=500\,{\rm cm}^{-1}\approx
0.06\,{\rm eV}$ for the electronic coupling being responsible for an electron
transfer, a common value found in acceptor aggregates such as fullerene
derivatives [53, 54, 68]. The exciton energy $\Omega_{0}$ depends on the
molecular properties of the donor [21, 69, 22, 10], which will be considered a
free variable parametrized by the driving force
$\Delta=\Omega_{0}-\Omega_{1}$, as shown in Fig.1(a).

For simplicity, we assume that each electronic state $\left|k\right>$ is
coupled to an independent vibrational environment that is initially in a
thermal state at room temperature. The vibrational Hamiltonian is written as

$H_{v}=\sum_{k=0}^{N-1}\sum_{q}\omega_{q}b_{k,q}^{\dagger}b_{k,q},$ (2)

with $b_{k,q}$ ($b_{k,q}^{\dagger}$) describing the annihilation (creation)
operator of a vibrational mode with frequency $\omega_{q}$ that is locally
coupled to the electronic state $\left|k\right>$. The vibronic interaction is
modeled by

$H_{e-v}=\sum_{k=0}^{N-1}\left|k\right>\left<k\right|\sum_{q}\omega_{q}\sqrt{s_{q}}(b_{k,q}+b_{k,q}^{\dagger}),$
(3)

where the vibronic coupling strength is quantified by the Huang-Rhys (HR)
factors $s_{q}$. The vibrational environments are fully characterized by a
phonon spectral density ${\cal
J}(\omega)=\sum_{q}\omega_{q}^{2}s_{q}\delta(\omega-\omega_{q})$ with
$\delta(\omega)$ denoting the Dirac delta function. According to first-
principles calculations of functionalized fullerene electron acceptors, the
vibrational environment consists of multiple low-frequency modes, with
vibrational frequencies smaller than the thermal energy at room temperature
($k_{B}T\approx 200\,{\rm cm}^{-1}\approx 0.025\,{\rm eV}$), and a few
discrete modes with high vibrational frequencies of the order of $\sim
1000\,{\rm cm}^{-1}$ and HR factors $\lesssim 0.1$ [70, 71, 72, 73]. Motivated
by these observations, we consider a phonon spectral density ${\cal
J}(\omega)={\cal J}_{l}(\omega)+{\cal J}_{h}(\omega)$ where ${\cal
J}_{l}(\omega)=\frac{\lambda_{l}}{\omega_{l}}\omega e^{-\omega/\omega_{l}}$,
with $\omega_{l}=80\,{\rm cm}^{-1}$ and $\lambda_{l}=50\,{\rm cm}^{-1}$,
describing a low-frequency phonon spectrum (see gray curve in Fig.1(b)). The
high-frequency vibrational modes are modeled by a Lorentzian function ${\cal
J}_{h}(\omega)=\frac{4\omega_{h}s_{h}\gamma(\omega_{h}^{2}+\gamma^{2})}{\pi}\omega((\omega+\omega_{h})^{2}+\gamma^{2})^{-1}((\omega-\omega_{h})^{2}+\gamma^{2})^{-1}$
with vibrational frequency $\omega_{h}=1200\,{\rm cm}^{-1}\approx
V/2=0.15\,$eV and HR factor $s_{h}=0.1$. Here the reorganization energy of the
high-frequency mode, defined by $\int_{0}^{\infty}d\omega{\cal
J}_{h}(\omega)\omega^{-1}=w_{h}s_{h}$, is independent of its vibrational
damping rate $\gamma$ (see red and blue curves in Fig.1(b)).

In order to tackle the problem of simulating large vibronic systems, we have
extended DAMPF [64], where a continuous vibrational environment is described
by a finite number of oscillators undergoing Markovian dissipation
(pseudomodes) and a tensor network formalism is used. With DAMPF the reduced
electronic system dynamics can be simulated in a numerically accurate manner
for highly structured phonon spectral densities by fitting the corresponding
bath correlation functions via an optimal set of parameters of either coupled
or uncoupled pseudomodes [63, 64, 65, 66]. The extended DAMPF-method opens the
door to non-perturbative simulations of many body systems consisting of
several tens of sites coupled to structured environments in one, two- and
three spatial dimensions, as will be demonstrated in this work. More details
about the method and the explicit equation of motion in terms of pseudomodes
can be found in Appendix C.

##  III Results

###  III.1 Driving Force and Vibrational Environments

Here we investigate the charge separation dynamics on a sub-ps time scale
simulated by DAMPF. For simplicity, we consider a linear chain consisting of a
donor and nine acceptors ($N=10$). Longer one-dimensional chains and higher-
dimensional donor/acceptor networks will be considered later. We assume that
an exciton state $\left|0\right>$ localised at the donor site is created at
the initial time $t=0$ and then an electron transfer through the acceptors
induces the transitions from the exciton to the CT states $\left|k\right>$
with $k\geq 1$. The mean distance between electron and hole is considered a
figure of merit for charge separation, defined by $\langle
x(t)\rangle=\sum_{k=0}^{N-1}kP_{k}(t)$ with $P_{k}(t)$ representing the
populations of the exciton and CT states $\left|k\right>$ at time $t$, with
the assumption that the distance between nearby sites is uniform. To
investigate how the initial charge separation dynamics depends on the exciton
energy $\Omega_{0}$ and the structure of vibrational environments, we analyse
the time-averaged electron-hole distance, defined by $\langle x\rangle_{t\leq
T}=\frac{1}{T}\int_{0}^{T}dt\langle x(t)\rangle$ with $T=400\,{\rm fs}$, as a
function of the driving force $\Delta=\Omega_{0}+V$ for various environmental
structures. The role of high-frequency vibrational modes and their non-
equilibrium dynamics in charge separation processes is identified by
considering (i) no environments (${\cal J}(\omega)=0$), (ii) low-frequency
phonon baths (${\cal J}(\omega)={\cal J}_{l}(\omega)$, see gray curve in
Fig.1(b)), (iii) high-frequency vibrational modes with controlled damping
rates $\gamma\in\\{(50\,{\rm fs})^{-1},(500\,{\rm fs})^{-1}\\}$ (${\cal
J}(\omega)={\cal J}_{h}(\omega)$, see red and blue curves in Fig.1(b)), and
(iv) the total vibrational environments including both low-frequency phonon
baths and high-frequency vibrational modes (${\cal J}(\omega)={\cal
J}_{l}(\omega)+{\cal J}_{h}(\omega)$).

Figure 3: Long-range charge separation dynamics. (a-b) Total probability for
separating an electron-hole pair more than four molecular units,
$\sum_{k=5}^{9}P_{k}(t)$, is shown as a function of time $t$ and driving force
$\Delta$. Here a linear chain consisting of $N=10$ sites is considered where
electronic states are coupled to total vibrational environments (${\cal
J}(\omega)={\cal J}_{l}(\omega)+{\cal J}_{h}(\omega)$). The damping rate of
the high-frequency vibrational modes with $\omega_{h}=1200\,{\rm cm}^{-1}$ is
taken to be (a) $\gamma=(500\,{\rm fs})^{-1}$ or (b) $\gamma=(50\,{\rm
fs})^{-1}$.

In Fig.2(a), the time-averaged electron-hole distance is shown as a function
of the driving force $\Delta$ when vibrational environments are not considered
($\mathcal{J}(\omega)=0$). In this case, the charge separation dynamics is
purely electronic and the mean electron-hole distance shows multiple peaks for
$\Delta\lesssim 0.3\,{\rm eV}$. When electronic states are only coupled to
low-frequency phonon baths (${\cal J}(\omega)={\cal J}_{l}(\omega)$), these
peaks are smeared out, resulting in a smooth, broad single peak centered at
$\Delta_{e}\approx 0.15\;\text{eV}$. In Fig.2(b) where the electronic states
are coupled to high-frequency vibrational modes (${\cal J}(\omega)={\cal
J}_{h}(\omega)$), the time-averaged electron-hole distance is displayed for
different vibrational damping rates $\gamma=(50\,{\rm fs})^{-1}$ and
$\gamma=(500\,{\rm fs})^{-1}$, shown in red and blue, respectively. With
$\omega_{h}$ denoting the vibrational frequency of the high-frequency modes,
the electron-hole distance is maximized at $\Delta_{e}\approx 0.15\,{\rm eV}$,
$\Delta_{e}+\omega_{h}\approx 0.3\,{\rm eV}$, $\Delta_{e}+2\omega_{h}\approx
0.45\,{\rm eV}$, making the charge separation process efficient for a broader
range of the driving force $\Delta$ when compared to the cases that the high-
frequency modes are ignored (see Fig.2(a)). It is notable that the electron-
hole distance is larger for the lower damping rate $\gamma=(500\,{\rm
fs})^{-1}$ of the high-frequency vibrational modes than for the higher damping
rate $\gamma=(50\,{\rm fs})^{-1}$. These results imply that non-equilibrium
vibrational dynamics can promote long-range charge separation. This
observation still holds even if the low-frequency phonon baths are considered
in addition to the high-frequency vibrational modes (${\cal J}(\omega)={\cal
J}_{l}(\omega)+{\cal J}_{h}(\omega)$), as shown in Fig.2(c), where the
electron-hole distance is maximized at $\Delta_{e}\approx 0.15\;\text{eV}$ and
$\Delta_{v}=\Delta_{e}+\omega_{h}\approx 0.3\;\text{eV}$. We note that the
electron-hole distance at low driving forces $\Delta\sim\Delta_{e}$ is
insensitive to the presence of vibrational environments, while at high driving
forces $\Delta\sim\Delta_{v}$, the charge separation process becomes
significantly inefficient when the high-frequency vibrational modes are
ignored. These results suggest that vibrational environments may play an
essential role in the long-range charge separation at high driving forces,
while the exciton dissociation at low driving forces may be governed by
electronic interactions.

So far the time-averaged mean electron-hole distance has been considered to
identify under what conditions the charge separation on a sub-ps time scale
becomes efficient. However, it does not show how much populations of the CT
states with well-separated electron-hole pairs are generated and how quickly
the long-range electron-hole separation takes place. In Fig.3, we show the
population dynamics of the CT states where electron and hole are separated
more than four molecular units, defined by $\sum_{k=5}^{9}P_{k}(t)$, for the
case that electronic states are coupled to the total vibrational environments
(${\cal J}(\omega)={\cal J}_{l}(\omega)+{\cal J}_{h}(\omega)$). When the high-
frequency vibrational modes are weakly damped with $\gamma=(500\,{\rm
fs})^{-1}$, the electron is transferred to the second half of the acceptor
chain within 100 fs and then the long-range electron-hole separation is
sustained on a sub-ps time scale for a wide range of the driving forces
$\Delta$, as shown in Fig.3(a). When the high-frequency modes are strongly
damped with $\gamma=(50\,{\rm fs})^{-1}$, for low driving forces around
$\Delta_{e}\approx 0.15\,$eV the long-range charge separation occurs within
100 fs, but the electron is quickly transferred back to the donor-acceptor
interface, as shown in Fig.3(b). For high driving forces around
$\Delta_{v}\approx 0.3\,$eV, the long-range charge separation and subsequent
localization towards the interface take place on a slower time scale when
compared to the case of the low driving forces. These results demonstrate that
underdamped vibrational motion can promote long-range charge separation when
the excess energy $\Delta-V$, defined by the energy difference between exciton
state and fully separated free charge carriers, is negative or close to zero
[74, 34, 33, 75, 35].

Figure 4: Vibration-assisted exciton dissociation and charge localization
towards donor/acceptor interfaces. (a) With a hole fixed at donor, the
probability distribution for finding an electron at the donor ($k=0$,
corresponding to exciton) or at the $k$-th acceptor ($k\geq 1$) is displayed
as a function of time $t$, with mean electron-hole distance $\langle
x(t)\rangle$ shown in red. With high driving force $\Delta_{v}=0.3\,$eV, here
we consider a linear chain consisting of $N=10$ sites and total vibrational
environments (${\cal J}(\omega)={\cal J}_{l}(\omega)+{\cal J}_{h}(\omega)$)
including strongly damped high-frequency modes with $\omega_{h}=1200\,{\rm
cm}^{-1}$ and $\gamma=(50\,{\rm fs})^{-1}$. (b) Charge separation dynamics
when the high-frequency vibrational modes are weakly damped with
$\gamma=(500\,{\rm fs})^{-1}$. (c) Population dynamics of interfacial CT state
$\left|1\right>$. (d) Probability for separating an electron-hole pair more
than four molecular units, $\sum_{k=5}^{9}P_{k}(t)$, shown as a function of
time $t$. In both (c) and (d), the strongly (weakly) damped case with
$\gamma=(50\,{\rm fs})^{-1}$ ($\gamma=(500\,{\rm fs})^{-1}$ is shown in red
(blue).

###  III.2 Electronic mixing at low driving forces

The long-range charge separation observed in DAMPF simulations can be
rationalized by analysing the energy-levels and delocalization lengths of the
exciton and CT states. In Fig.2(d), we consider the eigenstates
$|E_{\alpha}^{\rm(CT)}\rangle$ of the electronic Hamiltonian where the exciton
state $\left|0\right>$ and its coupling $J_{0,1}$ to the CT states are
ignored, namely $H_{\rm
CT}=\sum_{k=1}^{N-1}\Omega_{k}\left|k\right>\left<k\right|+\sum_{k=1}^{N-2}J_{k,k+1}(\left|k\right>\left<k+1\right|+h.c.)$.
With a hole fixed at the donor site, the probability distributions $|\langle
k|E_{\alpha}^{\rm(CT)}\rangle|^{2}$ for finding an electron at the $k$-th
acceptor site is displayed, which are vertically shifted by
$E_{\alpha}^{\rm(CT)}+V$ with $E_{\alpha}^{\rm(CT)}$ representing the
eigenvalues of $H_{\rm CT}$. The lowest-energy CT eigenstate is mainly
localised at the interface due to the strong Coulomb-binding energy considered
in simulations ($\Omega_{2}-\Omega_{1}=V/2=0.15\,{\rm eV}>J_{1,2}\approx
0.06\,{\rm eV}$). The other higher-energy CT eigenstates are significantly
delocalised in the acceptor domain with smaller populations $|\langle
1|E_{\alpha}^{\rm(CT)}\rangle|^{2}$ at the interface for higher energies
$E_{\alpha}^{\rm(CT)}$, as highlighted in red.

We now consider the full electronic Hamiltonian $H_{e}$ including the exciton
state at the low driving force $\Delta_{e}\approx 0.15\,{\rm eV}$ where
efficient long-range charge separation occurs even in the absence of
vibrational environments (see Fig.2(a)). The exciton state $\left|0\right>$ is
coupled to the eigenstates $|E_{\alpha}^{\rm(CT)}\rangle$ of $H_{\rm CT}$ via
the electronic coupling $H_{i}=J_{0,1}(\left|0\right>\left<1\right|+h.c.)$ at
the interface, leading to the exciton-CT couplings in the form
$\left<0\right|H_{i}|E_{\alpha}^{\rm(CT)}\rangle=J_{0,1}\langle
1|E_{\alpha}^{\rm(CT)}\rangle$. This implies that the transition between
exciton and CT state $|E_{\alpha}^{\rm(CT)}\rangle$ is enhanced when the
exciton energy $\Omega_{0}$ is near-resonant with the CT energy
$E_{\alpha}^{\rm(CT)}$ and the CT state has sufficiently high population
$|\langle 1|E_{\alpha}^{\rm(CT)}|^{2}$ at the interface (see red bars in
Fig.2(d)). For $\Delta_{e}=\Omega_{0}+V=0.15\,{\rm eV}$, the exciton state can
be strongly mixed with a near-resonant CT state delocalised over multiple
acceptor sites (see Fig.2(d)), leading to two hybrid exciton-CT eigenstates of
the total electronic Hamiltonian $H_{e}$, described by the superpositions of
$\left|0\right>$ and multiple $\left|k\right>$ with $k\geq 1$ (see Fig.2(e)).
This indicates that the multiple peaks in the time-averaged electron-hole
distance $\langle x\rangle_{t\leq 400\,{\rm fs}}$ shown in Fig.2(a) originate
from the resonances between exciton and CT states
$|E_{\alpha}^{\rm(CT)}\rangle$. Here the high-lying CT states with energies
$E_{\alpha}^{\rm(CT)}+V\gtrsim 0.3\,{\rm eV}$ do not show long-range electron-
hole separation, as the interfacial electronic couplings $J_{0,1}\langle
1|E_{\alpha}^{\rm(CT)}\rangle$ are not strong enough to induce notable
transitions between exciton and CT states within the time scale $T=400\,{\rm
fs}$ considered in Fig.2(a). These high-energy CT states can be populated via
a near-resonant exciton state, but the corresponding purely electronic charge
separation occurs on a slower ps time scale, as shown in Appendix D, and
therefore this process can be significantly affected by low-frequency phonon
baths. This is contrary to the charge separation at the low driving force
$\Delta_{e}\approx 0.15\,{\rm eV}$, which takes place within 100 fs and
therefore the early electronic dynamics is weakly affected by vibrational
environments. We note that when this analysis is applied to the charge
separation model in Ref.[53], it can be shown that an exciton state is
strongly mixed with near-resonant CT states delocalised in an effective one-
dimensional Coulomb potential and as a result the ultrafast long-range charge
separation reported in Ref.[53] can be well described by a purely electronic
model where vibrational environments are ignored (see Appendix B).

###  III.3 Vibronic mixing at high driving forces

Contrary to the case of $\Delta_{e}=0.15\,{\rm eV}$, the eigenstates of the
full electronic Hamiltonian $H_{e}$ with $\Delta_{v}=0.3\,{\rm eV}$ show a
weak mixing between exciton and CT states, as displayed in Fig.2(f), where the
eigenstate $\left|E_{\rm XT}\right>$ with the most strong excitonic character
$|\langle 0|E_{\rm XT}\rangle|\approx 1$, marked in blue, has small amplitudes
around the interface, $|\langle k|E_{\rm XT}\rangle|\ll 1$ for $k\geq 1$. Here
the energy-gaps between the exciton state $\left|E_{\rm XT}\right>$, shown in
blue, and lower-energy eigenstates $\left|E_{\rm CT}\right>$ with strong CT
characters, shown in green, are near-resonant with the vibrational frequency
of the high-frequency modes, $E_{\rm XT}-E_{\rm CT}\approx\omega_{h}$.
Therefore, the vibrationally cold exciton state $\left|E_{\rm
XT},0_{v}\right>$ can resonantly interact with vibrationally hot CT states
$\left|E_{\rm CT},1_{v}\right>$ where one of the high-frequency modes is
singly excited. Here the CT states are delocalised in the acceptor domain, but
have non-negligible amplitudes around the interface, leading to a moderate
vibronic coupling to the exciton state, $\left<E_{\rm
XT}\right|H_{e-v}\left|E_{\rm CT}\right>=\sum_{k=0}^{N-1}\langle E_{\rm
XT}|k\rangle\langle k|E_{\rm
CT}\rangle\omega_{h}\sqrt{s_{h}}(b_{k,h}+b_{k,h}^{\dagger})$ with $b_{k,h}$
($b_{k,h}^{\dagger}$) denoting the annihilation (creation) operator of the
high-frequency vibrational mode locally coupled to electronic state
$\left|k\right>$. The other high-lying CT states $\left|E_{\rm
CT}^{\prime}\right>$ near-resonant with the exciton state, $E_{\rm
CT}^{\prime}\approx E_{\rm XT}$, may have relatively small amplitudes around
the interface, so the direct vibronic coupling to the exciton state could be
small. However, the transitions from the exciton $\left|E_{\rm
XT},0_{v}\right>$ to the vibrationally hot low-lying CT states $\left|E_{\rm
CT},1_{v}\right>$ can allow subsequent transitions to vibrationally cold high-
lying CT states $\left|E_{\rm CT}^{\prime},0_{v}\right>$, as the delocalised
CT states $\left|E_{\rm CT}\right>$ and $\left|E_{\rm CT}^{\prime}\right>$ are
spatially overlapped. Such consecutive transitions are mediated by vibrational
excitations and can delay the process of charge localization at donor-acceptor
interfaces if the damping rate of the high-frequency vibrational modes is
sufficiently lower than the transition rates amongst exciton and CT states.
This picture is in line with the vibronic eigenstate analysis where the high-
frequency modes are included as a part of system Hamiltonian in addition to
the electronic states, as summarised in Appendix E.

Figure 5: Charge separation in large vibronic systems in one, two and three
spatial dimensions. (a) With a hole fixed at donor, the probability
distribution for finding an electron at the donor ($k=0$) or at the $k$-th
acceptor ($k\geq 1$) is shown as a function of time $t$ for a linear chain
consisting of $N=20$ sites. Here electronic states are coupled to total
vibrational environments (${\cal J}(\omega)={\cal J}_{l}(\omega)+{\cal
J}_{h}(\omega)$) including weakly damped high-frequqency modes with
$\omega_{h}=1500\,{\rm cm}^{-1}$ and $\gamma=(500\,{\rm fs})^{-1}$. (b) Mean
electron-hole distances $\langle x(t)\rangle$ of linear chains composed of
$N\in\\{10,20,30\\}$ sites, shown in $\\{$red,blue,black$\\}$, respectively.
Here we compare the case of the full environments, shown in solid lines, with
that of no vibrational environments (${\cal J}(\omega)=0$), shown in dashed
lines. (c) Schematic representation of one-, two- and three-dimensional donor-
acceptor networks where a donor is coupled to acceptor aggregates. (d) Time-
averaged mean electron-hole distances $\langle x\rangle_{t\leq 400\,{\rm fs}}$
of one- and two-dimensional networks consisting of $L=4$ acceptor layers,
displayed as a function of driving force $\Delta$ and computed by DAMPF [64].
Here the full and no environment cases are shown in blue and yellow,
respectively. Note that the overlapped regions are shown in darker yellow,
indicating that the no environment case is well covered by the full
environment case. (e) Time-averaged mean electron-hole distances of one-, two-
and three-dimensional networks consisting of $L=4$ or $L=9$ acceptor layers,
simulated by a reduced vibronic model (see Appendix E for more details). For
the three-dimensional network with $L=4$, $\langle x\rangle_{t\leq 400\,{\rm
fs}}$ is maximized at high driving forces $\Delta\approx 0.35\,$eV. (f) For
$\Delta=0.35\,$eV, time evolution of mean electron-hole distance $\langle
x(t)\rangle$ is shown for one-, two- and three-dimensional networks with
$L=4$, computed by DAMPF.

###  III.4 Functional Relevance of Long-lived Vibrational Motion

So far we have discussed the underlying mechanisms behind long-range charge
separation on a sub-ps time scale. We now investigate how subsequent charge
localization towards the donor-acceptor interface depends on the lifetimes of
high-frequency vibrational modes to demonstrate that non-equilibrium
vibrational dynamics can maintain long-range electron-hole separation.

In Fig.4(a) and (b), where the high-frequency modes are strongly and weakly
damped, respectively, with $\gamma=(50\,{\rm fs})^{-1}$ and $\gamma=(500\,{\rm
fs})^{-1}$, the population dynamics $P_{k}(t)$ of the exciton $\left|0\right>$
and CT states $\left|k\right>$ with $k\geq 1$ is shown as a function of time
$t$ up to 1.5 ps in addition to the mean electron-hole distance $\langle
x(t)\rangle$. Here we consider the high driving force $\Delta_{v}=0.3\,{\rm
eV}$ where the vibronic transition from exciton $\left|E_{\rm
XT},0_{v}\right>$ to delocalised CT states $\left|E_{\rm CT},1_{v}\right>$
takes place. When the high-frequency modes are strongly damped, the
vibrationally hot CT states $\left|E_{\rm CT},1_{v}\right>$ quickly dissipate
to $\left|E_{\rm CT},0_{v}\right>$, leading to subsequent vibronic transitions
to vibrationally hot interfacial CT states $\left|E_{\rm ICT},1_{v}\right>$
(see Fig.2(f)). After that, the vibrational damping of the high-frequency
modes generates the population of the lowest-energy interfacial CT state
$\left|E_{\rm ICT},0_{v}\right>$ and makes the electron-hole pair trapped at
the interface, as shown in Fig.4(a). When the high-frequency vibrational modes
are weakly damped, the mean electron-hole distance is maximized at $\sim
700\,{\rm fs}$, as shown in Fig.4(b), and then the population $P_{1}(t)$ of
the CT state $\left|1\right>$ localised around the interface starts to be
increased. This localized interfacial state $\left|1\right>$ has been
considered an energetic trap that leads to non-radiative losses [19]. In
Fig.4(c), the population dynamics of $P_{1}(t)$ is shown in red and blue,
respectively, for $\gamma=(50\,{\rm fs})^{-1}$ and $\gamma=(500\,{\rm
fs})^{-1}$. In the strongly damped case, $P_{1}(t)$ rapidly increases in time
and then saturates at $\sim 0.9$ on a picosecond time scale. This is contrary
to the weakly damped case where $P_{1}(t)$ is quickly saturated at $\sim 0.1$
within 100 fs and then does not increase until $\sim 500\,{\rm fs}$,
demonstrating that the charge localization towards the interface can be
delayed by the underdamped nature of the high-frequency vibrational modes. The
delayed charge localization makes long-range electron-hole separation
maintained on a picosecond time scale, as shown in Fig.4(d) where
$\sum_{k=5}^{9}P_{k}(t)$ is plotted. These results suggest that long-lived
vibrational and vibronic coherences observed in nonlinear optical spectra of
organic solar cells [44, 46, 48] may have a functional relevance in long-range
charge separation.

###  III.5 Large Vibronic Systems

So far we have considered a one-dimensional chain consisting of $N=10$ sites.
Here we investigate the charge separation dynamics in larger multi-site
systems, including longer linear chains, and donor-acceptor networks in two
and three spatial dimensions.

For the linear chains consisting of a donor and $(N-1)$ acceptors, we consider
the total vibrational environments including low-frequency phonon baths and
high-frequency vibrational modes with $\gamma=(500\,{\rm fs})^{-1}$ (${\cal
J}(\omega)={\cal J}_{l}(\omega)+{\cal J}_{h}(\omega)$). The driving force is
taken to be $\Delta_{v}=0.3\,{\rm eV}$, for which long-range charge separation
occurs mediated by vibronic couplings in the case of $N=10$ sites. In
Fig.5(a), a longer linear chain is considered with $N=20$ and the population
dynamics $P_{k}(t)$ of the exciton and CT states $\left|k\right>$ is shown. It
is notable that an electron-hole pair is separated more than ten molecular
units within $\sim 200$ fs. Interestingly, with a hole fixed at the donor
site, the probability distributions $P_{k}(t)$ for finding an electron at the
$k$-th acceptor are strongly delocalised over the entire acceptor chain, which
are maximized at $k\approx 6$ and locally minimized at $k\approx 3$. This
implies that an exciton state is vibronically mixed with strongly delocalised
CT states, as the detunings $\Omega_{k+1}-\Omega_{k}=V(k(k+1))^{-1}$ in the
energy-levels of the Coulomb potential become smaller in magnitude than the
electronic coupling $J_{k,k+1}=500\,{\rm cm}^{-1}$ being responsible for an
electron transfer when $V=0.3\,{\rm eV}$ and $k>1$. This is in line with the
dynamics of the mean electron-hole distance $\langle x(t)\rangle$ of the
linear chains consisting of $N\in\\{10,15,20\\}$ sites, shown in solid lines
in Fig.5(b), where the exciton dissociation becomes more efficient for longer
acceptor chains. The mean electron-hole distance is decreased when the energy-
levels $\Omega_{k}$ of the Coulomb potential are randomly generated based on
independent Gaussian distributions, as the delocalization lengths of the CT
states are reduced on average (see Appendix F). Importantly, for the high
driving force $\Delta_{v}=0.3\,{\rm eV}$, the charge separation process
becomes significantly less efficient when vibrational environments are not
considered (${\cal J}(\omega)=0$), as shown in dashed lines in Fig.5(b). These
results suggest that non-equilibrium vibrational dynamics in ordered
donor/acceptor aggregates can promote long-range charge separation.

From the perspective of the microcanonical ensemble, the number of charge-
separated states becomes much larger than that of interfacial CT states as the
dimension of donor-acceptor aggregates is increased [76]. The statistical
advantage results in an entropic drive that further promotes charge separation
[77] and is relevant in two- and three-dimensional donor-acceptor networks in
the thermodynamic limit. To corroborate these ideas, we consider a variety of
donor-acceptor networks with different sizes and dimensions. In Fig.5(c), the
schematic representations of one-, two- and three-dimensional donor-acceptor
networks considered in our simulations are displayed where the size of each
network is quantified by the number $L$ of acceptor layers. In the one-
dimensional chains, the number of acceptors in each layer is unity, while in
the two-dimensional triangular (three-dimensional pyramidal) structures, the
number of acceptors in each layer increases linearly (quadratically) as a
function of the minimum distance to the donor site. We assume that the
distances between nearby sites are uniform and the corresponding nearest-
neighbour electron-transfer couplings are taken to be $500\,{\rm cm}^{-1}$.
The electronic Hamiltonian is described by the exciton and CT states
$\left|k\right>$ where a hole is fixed at the donor while an electron is
localized at the $k$-th acceptor. The corresponding CT energy is modelled by
$\Omega_{k}=-V/|{\bf r}_{0}-{\bf r}_{k}|$ with $V=0.3\,{\rm eV}$ where ${\bf
r}_{0}$ and ${\bf r}_{k}$ denote, respectively, the positions of the donor and
$k$-th acceptor with the distance between nearby sites taken to be unity and
dimensionless. To increase the size of the donor-acceptor networks that can be
considered in simulations, we only consider the high-frequency vibrational
modes (${\cal J}(\omega)={\cal J}_{h}(\omega)$) with $\omega_{h}=1500\,{\rm
cm}^{-1}$, $s_{h}=0.1$ and $\gamma=(500\,{\rm fs})^{-1}$.

In Fig.5(d), the time-averaged electron-hole distance $\langle x\rangle_{t\leq
400\,{\rm fs}}$ simulated by DAMPF is shown as a function of the driving force
$\Delta$ for one- and two-dimensional networks with $L=4$. Here we consider
the minimum distance between donor and each acceptor layer in the computation
of the mean electron-hole distance, instead of the distances between donor and
individual acceptors. We compare the case that the high-frequency vibrational
modes are coupled to electronic states (${\cal J}(\omega)={\cal
J}_{h}(\omega)$), shown in blue, with that of no vibrational environments
(${\cal J}(\omega)=0$), shown in yellow (slightly darker in the overlapped
regions). Note that vibronic couplings make charge separation efficient for a
broader range of the driving force $\Delta$ in both one- and two-dimensional
networks, and that long-range charge separation is further enhanced in the
higher-dimensional network. To simulate larger vibronic systems, in Fig.5(e),
we consider a reduced vibronic model constructed within vibrational subspaces
describing up to four vibrational excitations distributed amongst the high-
frequency vibrational modes in the polaron basis (see Appendix E for more
details). For $L=4$, the simulated results obtained by the reduced models of
one- and two-dimensional networks are qualitatively similar to the numerically
exact DAMPF results shown in Fig.5(d). The reduced model results demonstrate
that long-range charge separation can be enhanced by considering a three-
dimensional donor-acceptor network with $L=4$, or by increasing the number of
layers to $L=9$ in the one- and two-dimensional cases. In Fig.5(f), the
dynamics of the mean electron-hole distance $\langle x(t)\rangle$ of the one-,
two- and three-dimensional systems with $L=4$, computed by DAMPF, is shown for
a high driving force $\Delta=0.35\,{\rm eV}$ where the time-averaged electron-
hole distance of the three-dimensional system shown in Fig.5(e) is maximized.
These results demonstrate that long-range charge separation can be enhanced by
considering higher-dimensional multi-site systems with vibronic couplings.

##  IV Conclusions

We have extended the non-perturbative simulation method DAMPF to provide
access to charge separation dynamics of a strongly bound electron-hole pair in
one-, two- and three-dimensional donor-acceptor networks where a donor is
coupled to acceptor aggregates. By controlling the driving force and the
structure of vibrational environments, we identified two distinct mechanisms
for long-range charge separation. The first mechanism, activated at low
driving forces, is characterized by hybrid exciton-CT states where long-range
exciton dissociation takes place on a sub-100 fs time scale, which is not
assisted by underdamped high-frequency vibrational modes. In the second
mechanism, dominating charge separation at high driving forces, the exciton-CT
hybridization occurs and it is mediated by vibronic interaction with
underdamped high-frequency vibrational modes, leading to efficient charge
separation for a broad range of driving forces. For both mechanisms, we have
demonstrated that long-range charge separation is significantly suppressed
when the high-frequency vibrational modes are strongly damped or
delocalization lengths of the CT states are reduced by static disorder in the
energy-levels of Coulomb potentials. These results suggest that non-
equilibrium vibrational motion can promote long-range charge separation in
ordered donor-acceptor aggregates.

The formulation and analysis of a reduced model whose validity became
accessible to numerical corroboration thanks to the extension of the
numerically exact simulation tool DAMPF allows us to identify unambiguously
the mechanisms that underlie charge separation dynamics. The methods employed
here can be applied to more realistic models where multiple donors are coupled
to acceptor aggregates, without introducing effective one-dimensional Coulomb
potentials, and vibrational environments are highly structured, which deserves
a separate investigation. We expect our findings to help open up the
engineering of vibrational environments for efficient long-range charge
separation in organic solar cells and the identification of charge separation
processes in other systems such as photosynthetic reaction centers and other
biological processes driven by electron transfer.

## Acknowledgments

This work was supported by the ERC Synergy grant HyperQ (grant no. 856432),
the BMBF project PhoQuant (grant no. 13N16110) under funding program quantum
technologies - from basic research to market, and an IQST PhD fellowship. The
authors acknowledge support by the state of Baden-Württemberg through bwHPC
and the German Research Foundation (DFG) through grant no. INST 40/575-1 FUGG
(JUSTUS 2 cluster)

##  Appendix A Markovian quantum master equation

In Ref.[52], a linear chain consisting of a donor and ten acceptors ($N=11$
sites) was considered where each electronic state $\left|k\right>$ is coupled
to a broad phonon spectrum, modelled by an Ohmic spectral density with a
Lorentz-Drude cutoff function, ${\cal
J}(\omega)=\frac{2}{\pi}\lambda\gamma\omega(\omega^{2}+\gamma^{2})^{-1}$,
where $\lambda=\int_{0}^{\infty}d\omega{\cal J}(\omega)\omega^{-1}$ is the
reoganization energy quantifying overall environmental coupling strength and
$\gamma$ determines the width of the phonon spectral density. Here underdamped
high-frequency vibrational modes are not considered, and as a result, quantum
correlations between electronic states and vibrational environments decay
within a 100 fs time scale for the values of $\gamma$ considered in Ref.[52].
The energy-levels of the CT states $\left|k\right>$ with $k\geq 1$ were
modelled by $\Omega_{k}=-V/k$ with $V=0.3$ eV, as is the case of our work. To
take into account non-Markovian environmental effects, hierarchical equations
of motion (HEOM) [78, 79] were used in Ref.[52]. However, electronic couplings
($J_{0,1}=0.15\,$eV and $J_{k,k+1}=0.10\,$eV for $k\geq 1$) were assumed to be
several times stronger than the overall environmental coupling strength
($\lambda=0.02$ eV), already suggesting that the influence of the broad phonon
spectral density considered in Ref.[52] on reduced electronic dynamics could
be well described by a Markovian quantum master equation where non-Markovian
environmental effects are ignored [80]. Fig.6 shows that the mean electron-
hole distances $\langle x(t)\rangle$ computed by HEOM, shown in dashed lines,
are well matched to the perturbative electronic dynamics computed by the
Redfield equation, shown in solid lines, for several values of $\gamma$
considered in Ref.[52]. These results demonstrate that the electronic dynamics
reported in Ref.[52] can be well reproduced by a Markovian noise model and
therefore non-Markovian effects may not play an important role in charge
separation processes for the model parameters considered in Ref.[52].

Figure 6: Charge separation dynamics in the presence of weakly coupled
unstructured vibrational environments. The mean electron-hole distances
$\langle x(t)\rangle$, computed by a non-perturbative method (HEOM), are shown
in dashed lines, while the predictions of a perturbative method (Redfield
equation) are shown in solid lines. Here we consider several values of bath
relaxation times $\gamma^{-1}\in\\{10,20,50,100,200\\}\,$fs (see blue, red,
yellow, purple, green lines, respectively), considered in Ref.[52].

##  Appendix B Weak Coulomb-binding energy

Figure 7: Electronic charge separation of weakly bound electron-hole pair. (a)
Time-averaged mean electron-hole distance $\langle x\rangle_{t\leq 400\,{\rm
fs}}$ is shown as a function of driving force $\Delta$, where a linear chain
consisting of $N=20$ sites is considered with electronic parameters from
Ref.[53]. Here ${\cal J}(\omega)=0$ (no vibrational modes), and $\Delta\approx
0.09\,{\rm eV}$ is the value considered in Ref.[53]. (b) The population
dynamics of $P_{0}(t)$, $P_{1}(t)$ and $\sum_{k=2}^{19}P_{k}(t)$, in the
absence of vibrational modes, are shown in blue, red and grey, respectively.
The transient of $\sum_{k=2}^{19}P_{k}(t)$ in the presence of vibrational
environments estimated by first-principles methods is shown in black with data
taken from Ref.[53].

Here we show that long-range charge separation can take place even in the
absence of vibrational environments when electronic couplings are sufficiently
stronger than detunings in the energy-levels of Coulomb potentials. In this
case, vibronic enhancement of charge separation can be negligible even if
vibrational environments are modelled based on first-principles methods.

In Ref.[53], where vibronic effects in charge separation process were
investigated by using first-principles methods, a hole transfer through a
linear chain of donors was simulated with the assumption that an electron is
delocalized in acceptor aggregates. It was suggested that the electron
delocalisation can reduce the Coulomb attraction between electron and hole.
This leads to an effective one-dimensional model where the detunings
$|\Omega_{k}-\Omega_{k+1}|$ in the energy-levels of the CT states
$\left|k\right>$ are smaller in magnitude than electronic couplings,
$J_{0,1}=0.2\,{\rm eV}$ and $J_{k,k+1}=0.12\,{\rm eV}$ for $1\leq k\leq 18$.
For these electronic parameters, the time-averaged mean electron-hole distance
$\langle x\rangle_{t\leq 400\,{\rm fs}}$ in the absence of vibrational
environments (${\cal J}(\omega)=0$) is shown as a function of driving force
$\Delta$ in Fig.7(a), which is maximised at $\Delta\approx 0.09\,{\rm eV}$,
being close to the driving force considered in numerical simulations of
Ref.[53]. For $\Delta\approx 0.09\,{\rm eV}$, the population dynamics of
$P_{0}(t)$, $P_{1}(t)$ and $\sum_{k=2}^{19}P_{k}(t)$ in the absence of
vibrational environments are shown in blue, red and grey lines, respectively,
in Fig.7(b), which are qualitatively similar to the results of Ref.[53] where
electronic-vibrational couplings were considered in simulations (see a black
line in Fig.7(b) displaying the transient of $\sum_{k=2}^{19}P_{k}(t)$ shown
in Ref.[53]). These results demonstrate that ultrafast long-range charge
separation can take place even if it is not assisted by vibronic couplings,
contrary to the vibronic mechanism proposed in Ref.[53], when electronic
couplings are sufficiently strong.

##  Appendix C DAMPF

In this work, non-perturbative simulations of charge separation dynamics have
been carried out by using dissipation-assisted matrix product factorization
(DAMPF) method [64] to take into account non-Markovian effects induced by a
continuous spectrum of vibrational environments. The DAMPF method relies on a
tensor network formalism and pseudomode theory [63, 66], which enables one to
consider a finite number of oscillators under Markovian noise (pseudomodes) to
mimic in a numerically accurate way the action of the original continuous
vibrational environments on reduced electronic system dynamics. In this work,
vibrational environments consist of low-frequency phonon baths and high-
frequency vibrational modes. The low-frequency phonon bath, modelled by an
Ohmic spectral density with an exponential cutoff function, has been described
by five pseudomodes coupled to each other. On the contrary, every high-
frequency vibrational mode has been described by an independent pseudomode,
which can well describe a narrow Lorentzian spectral density. In this case,
the vibrational Hamiltonians of the coupled and uncoupled pseudomodes,
respectively, are described by

$\displaystyle H_{v}^{(l)}=$
$\displaystyle\sum_{k=0}^{N-1}\sum_{q=1}^{5}(\omega_{q}a^{\dagger}_{k,q}a_{k,q}+g_{q}a_{k,q}a^{\dagger}_{k,q+1}+h.c.),$
(4) $\displaystyle H_{v}^{(h)}=$
$\displaystyle\sum_{k=0}^{N-1}\omega_{h}b^{\dagger}_{k,h}b_{k,h},$ (5)

where $a_{k,q}$ and $b_{k,h}$ represent pseudomodes coupled to an electronic
state $\left|k\right>$, and $g_{q}$ is the coupling between pseudomodes
$a_{k,q}$ and $a_{k,q+1}$ with $q+1$ reduced to 1 when $q=5$. The vibronic
couplings to the pseudomodes are described by

$\displaystyle H_{e-v}$
$\displaystyle=\sum_{k=0}^{N-1}\left|k\right>\left<k\right|\sum_{q=1}^{5}(c_{q}a_{k,q}+c_{q}^{*}a_{k,q}^{\dagger})$
(6)
$\displaystyle\quad+\sum_{k=0}^{N-1}\left|k\right>\left<k\right|\omega_{h}\sqrt{s_{h}}(b_{k,h}+b_{k,h}^{\dagger}).$
(7)

The time evolution of a vibronic density matrix $\rho$ is governed by a
Lindblad equation in the form

$\displaystyle\dot{\rho}$
$\displaystyle=-i[H_{e}+H_{v}^{(l)}+H_{v}^{(h)}+H_{e-v},\rho]$ (8)
$\displaystyle\quad+\sum_{k=0}^{N-1}\sum_{q=1}^{5}\Big{[}\gamma_{q}(n(\omega_{q})+1)(a_{k,q}\rho
a_{k,q}^{\dagger}-\frac{1}{2}\\{a_{k,q}^{\dagger}a_{k,q},\rho\\})$
$\displaystyle\quad\qquad\qquad\qquad\,+\gamma_{q}n(\omega_{q})(a_{k,q}^{\dagger}\rho
a_{k,q}-\frac{1}{2}\\{a_{k,q}a_{k,q}^{\dagger},\rho\\})\Big{]}$
$\displaystyle\quad+\sum_{k=0}^{N-1}\Big{[}\gamma_{h}(n(\omega_{h})+1)(b_{k,h}\rho
b_{k,h}^{\dagger}-\frac{1}{2}\\{b_{k,h}^{\dagger}b_{k,h},\rho\\})$
$\displaystyle\quad\qquad\qquad\,\,\,+\gamma_{h}n(\omega_{h})(b_{k,h}^{\dagger}\rho
b_{k,h}-\frac{1}{2}\\{b_{k,h}b_{k,h}^{\dagger},\rho\\})\Big{]},$

with $n(\omega)=(\exp(\hbar\omega/k_{B}T)-1)^{-1}$. The parameters
$\omega_{q}$, $g_{q}$, $c_{q}$ and $\gamma_{q}$ of the coupled pseudomodes can
be found in Table IV of Ref.[66] with a center frequency $\omega_{l}=80$ cm-1.
The vibrational frequency $\omega_{h}$ of the high-frequency modes is
approximately an order of magnitude larger than the thermal energy $k_{B}T$ at
room temperature, for which $n(\omega_{h})\ll 1$ and $\gamma_{h}$ is reduced
to the vibrational damping rate $\gamma\in\\{(50\,{\rm fs})^{-1},(500\,{\rm
fs})^{-1}\\}$ considered in the main text.

In DAMPF, where the number of uncoupled (coupled) pseudomodes is denoted by
$Q_{u}=N$ ($Q_{c}=5N$), we consider a vibronic density matrix $\rho$
represented in the form

$\displaystyle\rho=$
$\displaystyle\sum_{k,k^{\prime}=0}^{N-1}\left|k\right>\left<k^{\prime}\right|\otimes\sum_{i_{1}}\sum_{i_{2}}\cdots\sum_{i_{Q_{c}}}\sum_{i_{Q_{u}}}$
(9) $\displaystyle
A^{(k,k^{\prime})}_{(1;i_{1})}A^{(k,k^{\prime})}_{(2;i_{2})}\dots
A^{(k,k^{\prime})}_{(Q_{c};i_{Q_{c}})}A^{(k,k^{\prime})}_{(Q_{u};i_{Q_{u}})}x_{i_{1}}\otimes\dots\otimes
x_{i_{Q_{u}}},$

describing a collection of $N^{2}$ matrix product states (MPSs), via a
vectorization of the density matrix, conditional to the populations
$\left|k\right>\left<k\right|$ or coherences
$\left|k\right>\left<k^{\prime}\right|$ with $k\neq k^{\prime}$ of electronic
states $\left|k\right>$. This representation retains thus all quantum
correlations that have a clear electronic origin, restricting the
factorization at the level of different vibrational environments. We found
that the simulation costs of DAMPF required to achieve convergence in reduced
electronic dynamics are significantly higher for genuine vibronic dynamics
occurring at high driving forces $\Delta_{v}\approx 0.3$ eV than for
electronic dynamics taking place at low driving forces $\Delta_{e}\approx
0.15$ eV. This implies that strong correlations between electronic states and
vibrational modes, quantified by bond dimensions, occur when a vibronic mixing
induces long-range charge separation.

##  Appendix D Purely coherent electronic dynamics

In Fig.8, the probability for separating an electron-hole pair more than four
molecular units, defined by $\sum_{k=5}^{9}P_{k}(t)$, is shown as a function
of time $t$ and driving force $\Delta$ for a linear chain composed of a donor
and nine acceptors ($N=10$) where vibrational environments are not considered
(${\cal J}(\omega)=0$). Long-range charge separation takes place at discrete
values of the driving force $\Delta$ when exciton $\left|0\right>$ is near-
resonant with a delocalised CT state $|E_{\alpha}^{({\rm CT})}\rangle$. Here
oscillatory features of $\sum_{k=5}^{9}P_{k}(t)$ on a picosecond time scale
are due to the absence of noise induced by vibrational environments (${\cal
J}(\omega)=0$). Note that the period of oscillations is longer at higher
driving force $\Delta$, demonstrating that the coupling
$\left<0\right|H_{i}|E_{\alpha}^{\rm(CT)}\rangle=J_{0,1}\langle
1|E_{\alpha}^{\rm(CT)}\rangle$ between exciton and CT states is weaker for
higher-energy CT state. This implies that long-range charge separation at
higher driving force $\Delta$ requires a smaller energy-gap between exciton
and CT states since the degree of exciton-CT mixing is determined by the ratio
of the coupling and detuning between $\left|0\right>$ and $|E_{\alpha}^{({\rm
CT})}\rangle$, leading to narrower vertical line shapes at higher driving
forces in Fig.8.

Figure 8: Charge separation in the absence of vibrational environments. For a
linear chain consisting of $N=10$ sites, the probability for separating an
electron-hole pair more than four molecular units, $\sum_{k=5}^{9}P_{k}(t)$,
is shown as a function of time $t$ and driving force $\Delta$.

##  Appendix E Vibronic eigenstate analysis

To analyse vibronic mixing of exciton and delocalized CT states at high
driving force $\Delta_{v}=0.3\,$eV, we consider a reduced vibronic model where
electronic states are coupled to high-frequency vibrational modes. For
simplicity, the damping of the high-frequency modes and the vibronic coupling
to low-frequency phonon baths are not considered. The reduced model
Hamiltonian consists of three parts,
$H_{r}=H_{e}+H_{v}^{\prime}+H_{e-v}^{\prime}$, defined by

$\displaystyle
H_{e}=\sum_{k}\Omega_{k}\left|k\right>\left<k\right|+\sum_{k,k^{\prime}}J_{k,k^{\prime}}\left|k\right>\left<k^{\prime}\right|,$
(10) $\displaystyle
H_{v}^{\prime}=\sum_{k}\omega_{h}b_{k,h}^{\dagger}b_{k,h},$ (11)
$\displaystyle
H_{e-v}^{\prime}=\sum_{k}\left|k\right>\left<k\right|\omega_{h}\sqrt{s_{h}}(b_{k,h}+b_{k,h}^{\dagger}).$
(12)

To minimize the number of vibrational states required to achieve the numerical
convergence in reduced electronic dynamics, we consider a displaced
vibrational basis defined by unitary displacement operator conditional to
electronic states (polaron transformation),
$U=\sum_{k}\left|k\right>\left<k\right|D_{k}$ with
$D_{k}=\exp(\sqrt{s_{h}}(b_{k,h}^{\dagger}-b_{k,h}))$. In the polaron basis,
the reduced vibronic Hamiltonian is expressed as

$\displaystyle UH_{r}U^{\dagger}$
$\displaystyle=\sum_{k}(\Omega_{k}-\omega_{h}s_{h})\left|k\right>\left<k\right|+\sum_{k}\omega_{h}b_{k,h}^{\dagger}b_{k,h}$
(13)
$\displaystyle\quad+\sum_{k,k^{\prime}}J_{k,k^{\prime}}\left|k\right>\left<k^{\prime}\right|D_{k}D_{k^{\prime}}^{\dagger}.$

In simulations, we consider a vibrational subspace spanned by up to $N_{v}$
vibrational excitations distributed amongst multiple high-frequency modes
$b_{k,h}$. When $N_{v}=0$, only the global vibrational ground state
$\left|0_{v}\right>$ is considered where all the high-frequency modes are in
their vibrational ground states. When $N_{v}=1$, singly excited vibrational
states $|1_{v}^{(k)}\rangle$ are considered, in addition to
$\left|0_{v}\right>$, where only the $k$-th high-frequency mode $b_{k,h}$ is
singly excited, while all the other modes are in their vibrational ground
states. When $N_{v}=2$, doubly excited vibrational states are considered, in
addition to $\left|0_{v}\right>$ and $|1_{v}^{(k)}\rangle$, where only a
single vibrational mode is doubly excited ($|2_{v}^{(k)}\rangle$) or two
different modes are singly excited at the same time
($|1_{v}^{(k)},1_{v}^{(k^{\prime})}\rangle$).

Figure 9: Vibronic eigenstate analysis. (a) Transient of mean electron-hole
distance $\langle x(t)\rangle$ of a linear chain consisting of a donor and
nine acceptors where $\Delta_{v}=0.3\,{\rm eV}$. Here we consider a reduced
vibronic model where up to $N_{v}$ vibrational excitations are considered in
total. (b) Electron-hole distance $\langle x_{j}\rangle$ of vibronic
eigenstate $\left|\psi_{j}\right>$ is displayed as a function of the overlap
$|\langle\psi_{j}|\psi_{\rm initial}\rangle|^{2}$ with vibrationally cold
initial exciton state $\left|\psi_{\rm initial}\right>$. (c) Electronic
eigenstate structure. (d) For vibronic eigenstates that have large overlaps
with the initial state, denoted by $\left|\psi_{j=1,2,3}\right>$ in (b), the
populations $P_{j}(E_{\alpha})$ in electronic eigenbasis are shown, which are
dominated by exciton state $\left|E_{\rm XT}\right>=\left|E_{8}\right>$, low-
lying CT states $\left|E_{\rm
CT}\right>\in\\{\left|E_{3}\right>,\left|E_{4}\right>\\}$ and high-lying CT
state $\left|E_{\rm CT}^{\prime}\right>=\left|E_{9}\right>$.

In Fig.9, we consider a linear chain consisting of $N=10$ sites and
$\Delta_{v}=0.3\,{\rm eV}$ where the vibronic mixing of exciton and
delocalized CT states is induced by high-frequency vibrational modes with
$\omega_{h}=1200\,{\rm cm}^{-1}$ and $s_{h}=0.1$. In Fig.9(a), the mean
electron-hole distance $\langle x(t)\rangle$ computed by the reduced
Hamiltonian is displayed where the total number $N_{v}$ of vibrational
excitations is increased from 0 to 4. When vibrational excitations are not
considered in simulations ($N_{v}=0$), the mean electron-hole distance remains
below 1 up to $2\,{\rm ps}$, as shown in red. When vibrational excitations are
considered in simulations ($N_{v}\geq 1$), charge separation process is
significantly enhanced by vibronic couplings. The electronic dynamics shows
convergence when $N_{v}\geq 3$ where a blue line for $N_{v}=3$ is well
overlapped with a black line for $N_{v}=4$. The qualitative features of the
numerically converged electronic dynamics can be well reproduced by
approximate results for $N_{v}=1$ and $N_{v}=2$, as shown in orange and green,
respectively.

Table 1: Vibronic eigenstate analysis. For vibronic eigenstates $\left|\psi_{j=1,2,3}\right>$ shown in Fig.9(b), populations $P_{j}(E_{\alpha},N_{v})$ in electronic eigenbasis are summarized where $N_{v}$ denotes the total number of vibrational excitations distributed amongst high-frequency modes. $P_{j}(E_{\alpha},N_{v})<0.01$ are not shown here. $|\psi_{1}\rangle$ |  | $|\psi_{2}\rangle$ |
---|---|---|---
$|E_{\alpha},N_{v}\rangle$ | $P_{1}(E_{\alpha},N_{v})$ | $|E_{\alpha},N_{v}\rangle$ | $P_{2}(E_{\alpha},N_{v})$
$|E_{8},0_{v}\rangle$ | 0.55 | $|E_{4},1_{v}\rangle$ | 0.41
$|E_{4},1_{v}\rangle$ | 0.26 | $|E_{8},0_{v}\rangle$ | 0.20
$|E_{3},1_{v}\rangle$ | 0.06 | $|E_{3},1_{v}\rangle$ | 0.15
$|E_{9},0_{v}\rangle$ | 0.05 | $|E_{9},0_{v}\rangle$ | 0.12
$|E_{1},2_{v}\rangle$ | 0.02 | $|E_{1},2_{v}\rangle$ | 0.04
$|E_{2},1_{v}\rangle$ | 0.02 | $|E_{10},0_{v}\rangle$ | 0.03
$|E_{1},1_{v}\rangle$ | 0.01 | $|E_{2},1_{v}\rangle$ | 0.02
|  | $|E_{5},1\rangle$ | 0.01
$|\psi_{3}\rangle$ |  |  |
$|E_{\alpha},N_{v}\rangle$ | $P_{3}(E_{\alpha},N_{v})$ |  |
$|E_{4},1_{v}\rangle$ | 0.81 |  |
$|E_{8},0_{v}\rangle$ | 0.09 |  |
$|E_{1},2_{v}\rangle$ | 0.05 |  |

To identify the origin of the vibronic enhancement of charge separation, we
consider vibronic eigenstates of the polaron-transformed Hamiltonian
$UH_{r}U^{\dagger}$, represented by

$\left|\psi_{j}\right>=\sum_{\alpha=1}^{N}\sum_{\vec{n}_{v}}\psi_{j}(E_{\alpha},\vec{n}_{v})\left|E_{\alpha},\vec{n}_{v}\right>,$
(14)

where $\left|E_{\alpha}\right>$ denote the eigenstates of electronic
Hamiltonian $H_{e}$, while $\vec{n}_{v}$ describe the distributions of
vibrational excitations of the high-frequency modes. The populations of the
vibronic eigenstates $\left|\psi_{j}\right>$ in the electronic basis
$\left|k\right>$, where a hole is fixed at donor and an electron is localised
at the donor ($k=0$) or at the $k$-th acceptor ($k\geq 1$), are defined by
$P_{j}(k)={\rm Tr}_{v}[\langle k|\psi_{j}\rangle\langle\psi_{j}|k\rangle]$
where ${\rm Tr}_{v}$ denotes the partial trace over vibrational degrees of
freedom. Similarly we consider the populations of the vibronic eigenstates in
the electronic eigenbasis $\left|E_{\alpha}\right>$, defined by
$P_{j}(E_{\alpha})={\rm Tr}_{v}[\langle
E_{\alpha}|\psi_{j}\rangle\langle\psi_{j}|E_{\alpha}\rangle]$. In Fig.9(b),
the overlap $|\langle\psi_{j}|\psi_{\rm initial}\rangle|^{2}$ between vibronic
eigenstates $\left|\psi_{j}\right>$ and vibrationally cold initial exciton
state $\left|\psi_{\rm initial}\right>=U|0,0_{v}\rangle$ in the polaron basis
is shown as a function of electron-hole distance of the vibronic eigenstates,
defined by $\langle x_{j}\rangle=\sum_{k=0}^{N-1}kP_{j}(k)$. It is notable
that the overlap with the initial state is dominated by three vibronic
eigenstates $\left|\psi_{j=1,2,3}\right>$ with comparable electron-hole
distances of $\langle x_{j}\rangle\approx 4$, which is close to the maximum
electron-hole distance $\langle x(t)\rangle$ shown in Fig.9(a). The
populations $P_{j}(E_{\alpha})$ of $\left|\psi_{j=1,2,3}\right>$ in the
electronic eigenbasis are shown in Fig.9(d) as a function of the energy-levels
$E_{\alpha}$ of electronic eigenstates. It is notable that the vibronic
eigenstates are dominated by the populations of the exciton state
$\left|E_{\rm XT}\right>=\left|E_{8}\right>$ and lower-energy CT states
$\left|E_{\rm CT}\right>\in\\{\left|E_{4}\right>,\left|E_{5}\right>\\}$ with a
relatively small contribution of a high-lying CT state $\left|E_{\rm
CT}^{\prime}\right>=\left|E_{9}\right>$ near-resonant with the exciton state.
These results demonstrate that the exciton state is vibronically mixed with
delocalised CT states.

The presence of the vibronic mixing can be demonstrated more clearly by
analysing both electronic and vibrational states of the vibronic eigenstates
$\left|\psi_{j}\right>$. For a linear chain consisting of $N=10$ sites, the
electronic states are coupled to ten high-frequency vibrational modes in
total. To simplify analysis, for each vibronic eigenstate
$\left|\psi_{j}\right>$, we investigate the overlap with electronic
eigenstates $\left|E_{\alpha}\right>$ in the presence of $N_{v}$ vibrational
excitations, defined by $P_{j}(E_{\alpha},N_{v})=\sum_{{\rm
sum}(\vec{n}_{v})=N_{v}}|\psi_{j}(E_{\alpha},\vec{n}_{v})|^{2}$ where the
summation runs over all possible vibrational states with $N_{v}$ excitations
distributed amongst the ten high-frequency modes. In Table 1, all
$P_{j}(E_{\alpha},N_{v})$ being larger than 0.01 are shown, demonstrating that
the vibronic eigenstates $\left|\psi_{j=1,2,3}\right>$ governing initial
charge separation dynamics are well described by the superpositions of
vibrationally cold exciton state $\left|E_{\rm XT},0_{v}\right>$,
vibrationally hot $\left|E_{\rm CT},1_{v}\right>$ and cold CT states
$\left|E_{\rm CT}^{\prime},0_{v}\right>$ delocalised in the acceptor domain.

Figure 10: Charge separation dynamics in the presence of static disorder.
(a-c) Ensemble-averaged mean electron-hole distance $\langle x(t)\rangle$ of
one-dimensional chains consisting of (a) $N=10$, (b) $N=20$, (c) $N=30$ sites.
The degree of static disorder in the energy-levels $\Omega_{k}$ of Coulomb
potentials is quantified by the standard deviation $\sigma$ of independent
Gaussian distributions, which is taken to be $\sigma=50$, $100$, $200$, $500$,
$1000\,{\rm cm}^{-1}$, respectively, shown in red, orange, green, blue, black.
For each ensemble-averaged transient, 1000 randomly generated sets of energy-
levels $\Omega_{k}$ were considered.

##  Appendix F Static disorder

To investigate how static disorder in the energy-levels $\Omega_{k}$ of
Coulomb potentials affects charge separation dynamics, we consider one-
dimensional chains consisting of a donor and $(N-1)$ acceptors. The energy-
levels $\Omega_{k}$ of electronic states $\left|k\right>$ are modelled by
independent Gaussian distributions with mean values taken to be
$\langle\Delta\rangle=\langle\Omega_{0}-\Omega_{1}\rangle=0.3\,$eV and
$\langle\Omega_{k}\rangle=-V/k$ for $k\geq 1$. For simplicity, we employ the
reduced vibronic model in Appendix E.

In Fig.10(a), where a linear chain consisting of $N=10$ sites is considered,
the ensemble-averaged mean electron-hole distance $\langle x(t)\rangle$ is
shown for several values of the standard deviation $\sigma$ of the Gaussian
distributions. For each value of $\sigma$, 1000 randomly generated sets of the
energy-levels $\Omega_{k}$ were considered. It is notable that the electron-
hole distance is suppressed as the standard deviation $\sigma$ of the Gaussian
distributions is increased from 50, via 100, 200, 500, to $1000\,{\rm
cm}^{-1}$. This implies that charge separation is suppressed as delocalization
lengths of the CT states are decreased by static disorder. In Fig.10(b) and
(c), longer linear chains with $N=15$ and $N=20$ are considered, respectively,
where vibronic charge separation is also found to be suppressed by static
disorder. Note that as the degree of static disorder, quantified by $\sigma$,
is increased, the charge separation dynamics becomes almost independent of the
size $(N-1)$ of acceptor aggregates (see blue lines in Fig.10(a)-(c) where
$\sigma=J_{k,k+1}=500\,{\rm cm}^{-1}$).

## References

  * Scholes and Rumbles [2006] G. D. Scholes and G. Rumbles, Excitons in Nanoscale Systems, Nat. Mater. 5, 920 (2006).
  * Hedley  _et al._ [2017] G. J. Hedley, A. Ruseckas, and I. D. W. Samuel, Light Harvesting for Organic Photovoltaics, Chem. Rev. 117, 796 (2017).
  * Vandewal [2016] K. Vandewal, Interfacial Charge Transfer States in Condensed Phase Systems, Annu. Rev. Phys. Chem. 67, 113 (2016).
  * Coropceanu  _et al._ [2019] V. Coropceanu, X.-K. Chen, T. Wang, Z. Zheng, and J.-L. Brédas, Charge-transfer Electronic States in Organic Solar Cells, Nat. Rev. Mater. 4, 689 (2019).
  * van Amerongen  _et al._ [2000] H. van Amerongen, R. van Grondelle, and L. Valkunas,  _Photosynthetic Excitons_ (World Scientific, 2000).
  * Huelga and Plenio [2013] S. F. Huelga and M. B. Plenio, Vibrations, Quanta and Biology, Contemp. Phys. 54, 181 (2013).
  * Romero  _et al._ [2017] E. Romero, V. I. Novoderezhkin, and R. van Grondelle, Quantum Design of Photosynthesis for Bio-Inspired Solar-Energy Conversion, Nature 543, 355 (2017).
  * 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 (2009).
  * Park  _et al._ [2009] S. H. Park, A. Roy, S. Beaupré, S. Cho, N. Coates, J. S. Moon, D. Moses, M. Leclerc, K. Lee, and A. J. Heeger, Bulk Heterojunction Solar Cells with Internal Quantum Efficiency Approaching 100%, Nat. Photonics 3, 297 (2009).
  * Clarke and Durrant [2010] T. M. Clarke and J. R. Durrant, Charge Photogeneration in Organic Solar Cells, Chem. Rev. 110, 6736 (2010).
  * Gao and Inganäs [2014] F. Gao and O. Inganäs, Charge Generation in Polymer–Fullerene Bulk-Heterojunction Solar Cells, Phys. Chem. Chem. Phys. 16, 20291 (2014).
  * Faist  _et al._ [2011] M. A. Faist, T. Kirchartz, W. Gong, R. S. Ashraf, I. McCulloch, J. C. de Mello, N. J. Ekins-Daukes, D. D. C. Bradley, and J. Nelson, Competition Between the Charge Transfer State and the Singlet States of Donor or Acceptor Limiting the Efficiency in Polymer:Fullerene Solar Cells, J. Am. Chem. Soc. 134, 685 (2011).
  * Nan  _et al._ [2016] G. Nan, X. Zhang, and G. Lu, The Lowest-Energy Charge-Transfer State and its Role in Charge Separation in Organic Photovoltaics, Phys. Chem. Chem. Phys. 18, 17546 (2016).
  * Azzouzi  _et al._ [2018] M. Azzouzi, J. Yan, T. Kirchartz, K. Liu, J. Wang, H. Wu, and J. Nelson, Nonradiative Energy Losses in Bulk-Heterojunction Organic Photovoltaics, Phys. Rev. X 8, 31055 (2018).
  * Eisner  _et al._ [2019] F. D. Eisner, M. Azzouzi, Z. Fei, X. Hou, T. D. Anthopoulos, T. J. S. Dennis, M. Heeney, and J. Nelson, Hybridization of Local Exciton and Charge-Transfer States Reduces Nonradiative Voltage Losses in Organic Solar Cells, J. Am. Chem. Soc. 141, 6362 (2019).
  * Burke  _et al._ [2015] T. M. Burke, S. Sweetnam, K. Vandewal, and M. D. McGehee, Beyond Langevin Recombination: How Equilibrium Between Free Carriers and Charge Transfer States Determines the Open-Circuit Voltage of Organic Solar Cells, Adv. Energy Mater. 5, 1500123 (2015).
  * Menke  _et al._ [2018a] S. M. Menke, N. A. Ran, G. C. Bazan, and R. H. Friend, Understanding Energy Loss in Organic Solar Cells: Toward a New Efficiency Regime, Joule 2, 25 (2018a).
  * Liu  _et al._ [2019] X. Liu, Y. Li, K. Ding, and S. Forrest, Energy Loss in Organic Photovoltaics: Nonfullerene Versus Fullerene Acceptors, Phys. Rev. Appl. 11, 024060 (2019).
  * Azzouzi  _et al._ [2019] M. Azzouzi, T. Kirchartz, and J. Nelson, Factors Controlling Open-Circuit Voltage Losses in Organic Solar Cells, Trends Chem. 1, 49 (2019).
  * Troisi [2011] A. Troisi, Charge Transport in High Mobility Molecular Semiconductors: Classical Models and New Theories, Chem. Soc. Rev. 40, 2347 (2011).
  * 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 (2006).
  * Veldman  _et al._ [2009] D. Veldman, S. C. J. Meskers, and R. A. Janssen, The Energy of Charge‐Transfer States in Electron Donor–Acceptor Blends: Insight into the Energy Losses in Organic Solar Cells, Adv. Funct. Mater. 19, 1939 (2009).
  * Hallermann  _et al._ [2009] M. Hallermann, I. Kriegel, E. Da Como, J. M. Berger, E. von Hauff, and J. Feldmann, Charge Transfer Excitons in Polymer/Fullerene Blends: The Role of Morphology and Polymer Chain Conformation, Adv. Funct. Mater. 19, 3662 (2009).
  * Deibel and Dyakonov [2010] C. Deibel and V. Dyakonov, Polymer–fullerene Bulk Heterojunction Solar Cells, Rep. Prog. Phys. 73, 096401 (2010).
  * Coffey  _et al._ [2012] D. C. Coffey, B. W. Larson, A. W. Hains, J. B. Whitaker, N. Kopidakis, O. V. Boltalina, S. H. Strauss, and G. Rumbles, An Optimal Driving Force for Converting Excitons into Free Carriers in Excitonic Solar Cells, J. Phys. Chem. C 116, 8916 (2012).
  * Wang  _et al._ [2018] Y. Wang, D. Qian, Y. Cui, H. Zhang, J. Hou, K. Vandewal, T. Kirchartz, and F. Gao, Optical Gaps of Organic Solar Cells as a Reference for Comparing Voltage Losses, Adv. Energy Mater. 8, 1801352 (2018).
  * Qian  _et al._ [2018] D. Qian, Z. Zheng, H. Yao, W. Tress, T. R. Hopper, S. Chen, S. Li, J. Liu, S. Chen, J. Zhang, X.-K. Liu, B. Gao, L. Ouyang, Y. Jin, G. Pozina, I. A. Buyanova, W. M. Chen, O. Inganäs, V. Coropceanu, J.-L. Brédas, H. Yan, J. Hou, F. Zhang, A. A. Bakulin, and F. Gao, Design Rules for Minimizing Voltage Losses in High-Efficiency Organic Solar Cells, Nat. Mater. 17, 703 (2018).
  * Yang  _et al._ [2018] W. Yang, Y. Yao, P. Guo, H. Sun, and Y. Luo, Optimum Driving Energy for Achieving Balanced Open-Circuit Voltage and Short-Circuit Current Density in Organic Bulk Heterojunction Solar Cells, Phys. Chem. Chem. Phys. 20, 29866 (2018).
  * Zhang  _et al._ [2015] Q. Zhang, B. Kan, F. Liu, G. Long, X. Wan, X. Chen, Y. Zuo, W. Ni, H. Zhang, M. Li, Z. Hu, F. Huang, Y. Cao, Z. Liang, M. Zhang, T. P. Russell, and Y. Chen, Small-Molecule Solar Cells with Efficiency over 9%, Nat. Photonics 9, 35 (2015).
  * Jing  _et al._ [2016] L. Jing, S. Chen, D. Qian, B. Gautam, Y. Guofang, J. Zhao, J. Bergqvist, F. Zhang, W. Ma, H. Ade, O. Inganas, K. Gundogdu, F. Gao, and H. Yan, Fast Charge Separation in a Non-Fullerene Organic Solar Cell with a Small Driving Force, Nat. Energy 1, 16089 (2016).
  * Meng  _et al._ [2018] L. Meng, Y. Zhang, X. Wan, C. Li, X. Zhang, Y. Wang, 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 (2018).
  * Cui  _et al._ [2019] Y. Cui, H. Yao, J. Zhang, T. Zhang, Y. Wang, L. Hong, K. Xian, B. Xu, S. Zhang, J. Peng, Z. Wei, F. Gao, and J. Hou, Over 16% Efficiency Organic Photovoltaic Cells Enabled by a Chlorinated Acceptor with Increased Open-Circuit Voltages, Nat. Commun. 10, 2515 (2019).
  * Liu  _et al._ [2020] S. Liu, J. Yuan, W. Deng, M. Luo, Y. Xie, Q. Liang, Y. Zou, Z. He, H. Wu, and Y. Cao, High-Efficiency Organic Solar Cells with Low Non-Radiative Recombination Loss and Low Energetic Disorder, Nat. Photonics 14, 300 (2020).
  * Zhang  _et al._ [2020] G. Zhang, X.-K. Chen, J. Xiao, P. C. Y. Chow, M. Ren, G. Kupgan, X. Jiao, C. C. S. Chan, X. Du, R. Xia, Z. Chen, J. Yuan, Y. Zhang, S. Zhang, Y. Liu, Y. Zou, H. Yan, K. S. Wong, V. Coropceanu, N. Li, C. J. Brabec, J.-L. Brédas, H.-L. Yip, and Y. Cao, Delocalization of Exciton and Electron Wavefunction in Non-Fullerene Acceptor Molecules Enables Efficient Organic Solar Cells, Nat. Commun. 11, 3943 (2020).
  * Zhong  _et al._ [2020] Y. Zhong, M. Causa’, G. J. Moore, P. Krauspe, B. Xiao, F. Günther, J. Kublitski, R. Shivhare, J. Benduhn, E. BarOr, S. Mukherjee, K. M. Yallum, J. Réhault, S. C. B. Mannsfeld, D. Neher, L. J. Richter, D. M. DeLongchamp, F. Ortmann, K. Vandewal, E. Zhou, and N. Banerji, Sub-Picosecond Charge-Transfer at Near-Zero Driving Force in Polymer:Non-Fullerene Acceptor Blends and Bilayers, Nat. Commun. 11, 833 (2020).
  * Chen  _et al._ [2018a] 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 (2018a).
  * Firdaus  _et al._ [2019] Y. Firdaus, V. M. Le Corre, J. I. Khan, Z. Kan, F. Laquai, P. M. Beaujuge, and T. D. Anthopoulos, Key Parameters Requirements for Non‐Fullerene‐Based Organic Solar Cells with Power Conversion Efficiency >20%, Adv. Sci. 6, 1802028 (2019).
  * Karki  _et al._ [2020] A. Karki, J. Vollbrecht, A. J. Gillett, P. Selter, J. Lee, Z. Peng, N. Schopp, A. L. Dixon, M. Schrock, V. Nádaždy, F. Schauer, H. Ade, B. F. Chmelka, G. C. Bazan, R. H. Friend, and T. Nguyen, Unifying Charge Generation, Recombination, and Extraction in Low‐Offset Non‐Fullerene Acceptor Organic Solar Cells, Adv. Energy Mater. 10, 2001203 (2020).
  * Alvertis  _et al._ [2020] A. M. Alvertis, R. Pandya, L. A. Muscarella, N. Sawhney, M. Nguyen, B. Ehrler, A. Rao, R. H. Friend, A. W. Chin, and B. Monserrat, Impact of Exciton Delocalization on Exciton-Vibration Interactions in Organic Semiconductors, Phys. Rev. B 102, 081122 (2020).
  * Pensack and Asbury [2009] R. D. Pensack and J. B. Asbury, Barrierless Free Carrier Formation in an Organic Photovoltaic Material Measured with Ultrafast Vibrational Spectroscopy, J. Am. Chem. Soc. 131, 15986 (2009).
  * Vandewal  _et al._ [2014] K. Vandewal, S. Albrecht, E. T. Hoke, K. R. Graham, J. Widmer, J. D. Douglas, M. Schubert, W. R. Mateker, J. T. Bloking, G. F. Burkhard, A. Sellinger, J. M. J. Fréchet, A. Amassian, M. K. Riede, M. D. McGehee, D. Neher, and A. Salleo, Efficient Charge Generation by Relaxed Charge-Transfer States at Organic Interfaces, Nat. Mater. 13, 63 (2014).
  * Bässler and Köhler [2015] H. Bässler and A. Köhler, “Hot or cold”: How do Charge Transfer States at the Donor–Acceptor Interface of an Organic Solar Cell Dissociate?, Phys. Chem. Chem. Phys. 17, 28451 (2015).
  * Falke  _et al._ [2014] S. M. Falke, C. A. Rozzi, D. Brida, M. Maiuri, M. Amato, E. Sommer, A. D. Sio, A. Rubio, G. Cerullo, E. Molinari, and C. Lienau, Coherent ultrafast charge transfer in an organic photovoltaic blend, Science 344, 1001 (2014).
  * Song  _et al._ [2014] Y. Song, S. N. Clafton, R. D. Pensack, T. W. Kee, and G. D. Scholes, Vibrational coherence probes the mechanism of ultrafast electron transfer in polymer–fullerene blends, Nat. Commun. 5, 4933 (2014).
  * Rafiq  _et al._ [2015] S. Rafiq, J. C. Dean, and G. D. Scholes, Observing Vibrational Wavepackets during an Ultrafast Electron Transfer Reaction, J. Phys. Chem. A 119, 11837 (2015).
  * Sio  _et al._ [2016] A. D. Sio, F. Troiani, M. Maiuri, J. Réhault, E. Sommer, J. Lim, S. F. Huelga, M. B. Plenio, C. A. Rozzi, G. Cerullo, E. Molinari, and C. Lienau, Tracking the Coherent Generation of Polaron Pairs in Conjugated Polymers, Nat. Commun. 7, 13742 (2016).
  * Brédas  _et al._ [2017] J.-L. Brédas, E. H. Sargent, and G. D. Scholes, Photovoltaic Concepts Inspired by Coherence Effects in Photosynthetic Systems, Nat. Mater. 16, 35 (2017).
  * Sio  _et al._ [2018] A. D. Sio, F. V. d. A. Camargo, K. Winte, E. Sommer, F. Branchi, G. Cerullo, and C. Lienau, Ultrafast Relaxation Dynamics in a Polymer: Fullerene Blend for Organic Photovoltaics Probed by Two-Dimensional Electronic Spectroscopy, Eur. Phys. J. B 91, 236 (2018).
  * Sio  _et al._ [2021] A. D. Sio, E. Sommer, X. T. Nguyen, L. Groß, D. Popović, B. T. Nebgen, S. Fernandez-Alberti, S. Pittalis, C. A. Rozzi, E. Molinari, E. Mena-Osteritz, P. Bäuerle, T. Frauenheim, S. Tretiak, and C. Lienau, Intermolecular Conical Intersections in Molecular Aggregates, Nat. Nanotechnol. 16, 63 (2021).
  * Smith and Chin [2015] S. L. Smith and A. W. Chin, Phonon-Assisted Ultrafast Charge Separation in the PCBM Band Structure, Phys. Rev. B 91, 201302 (2015).
  * Lee  _et al._ [2015] C. K. Lee, J. Moix, and J. Cao, Coherent Quantum Transport in Disordered Systems: A Unified Polaron Treatment of Hopping and Band-like Transport, J. Chem. Phys. 142, 164103 (2015).
  * Kato and Ishizaki [2018] A. Kato and A. Ishizaki, Non-Markovian Quantum-Classical Ratchet for Ultrafast Long-Range Electron-Hole Separation in Condensed Phases, Phys. Rev. Lett. 121, 026001 (2018).
  * Tamura and Burghardt [2013a] H. Tamura and I. Burghardt, Ultrafast Charge Separation in Organic Photovoltaics Enhanced by Charge Delocalization and Vibronically Hot Exciton Dissociation, J. Am. Chem. Soc. 135, 16364 (2013a).
  * Tamura and Burghardt [2013b] H. Tamura and I. Burghardt, Potential Barrier and Excess Energy for Electron–Hole Separation from the Charge-Transfer Exciton at Donor–Acceptor Heterojunctions of Organic Solar Cells, J. Phys. Chem. C 117, 15020 (2013b).
  * Hughes  _et al._ [2014] K. H. Hughes, B. Cahier, R. Martinazzo, H. Tamura, and I. Burghardt, Non-Markovian Reduced Dynamics of Ultrafast Charge Transfer at an Oligothiophene-Fullerene Heterojunction, Chem. Phys. 442, 111 (2014).
  * Chenel  _et al._ [2014] A. Chenel, E. Mangaud, I. Burghardt, C. Meier, and M. Desouter-Lecomte, Exciton Dissociation at Donor-Acceptor Heterojunctions: Dynamics Using the Collective Effective Mode Representation of the Spin-Boson Model, J. Chem. Phys. 140, 044104 (2014).
  * Huix-Rotllant  _et al._ [2015] M. Huix-Rotllant, H. Tamura, and I. Burghardt, Concurrent Effects of Delocalization and Internal Conversion Tune Charge Separation at Regioregular Polythiophene–Fullerene Heterojunctions, J. Phys. Chem. Lett. 6, 1702 (2015).
  * Chaudhuri  _et al._ [2017] S. Chaudhuri, S. Hedström, D. D. Méndez-Hernández, H. P. Hendrickson, K. A. Jung, J. Ho, and V. S. Batista, Electron Transfer Assisted by Vibronic Coupling from Multiple Modes, J. Chem. Theory Comput. 13, 6000 (2017).
  * Polkehn  _et al._ [2018] M. Polkehn, H. Tamura, and I. Burghardt, Impact of Charge-Transfer Excitons in Regioregular Polythiophene on the Charge Separation at Polythiophene-Fullerene Heterojunctions, J. Phys.B 51, 014003 (2018).
  * Bian  _et al._ [2020] Q. Bian, F. Ma, S. Chen, Q. Wei, X. Su, I. A. Buyanova, W. M. Chen, C. S. Ponseca, M. Linares, K. J. Karki, A. Yartsev, and O. Inganäs, Vibronic Coherence Contributes to Photocurrent Generation in Organic Semiconductor Heterojunction Diodes, Nat. Commun. 11, 617 (2020).
  * Yao  _et al._ [2016] Y. Yao, X. Xie, and H. Ma, Ultrafast Long-Range Charge Separation in Organic Photovoltaics: Promotion by Off-Diagonal Vibronic Couplings and Entropy Increase, J. Phys. Chem. Lett. 7, 4830 (2016).
  * Duan  _et al._ [2020] H. G. Duan, A. Jha, V. Tiwari, R. J. Miller, and M. Thorwart, Dissociation and Localization Dynamics of Charge Transfer Excitons at a Donor-Acceptor Interface, Chem. Phys. 528, 110525 (2020).
  * Tamascelli  _et al._ [2018] D. Tamascelli, A. Smirne, S. F. Huelga, and M. B. Plenio, Nonperturbative Treatment of Non-Markovian Dynamics of Open Quantum Systems, Phys. Rev. Lett. 120, 30402 (2018).
  * Somoza  _et al._ [2019] A. D. Somoza, O. Marty, J. Lim, S. F. Huelga, and M. B. Plenio, Dissipation-Assisted Matrix Product Factorization, Phys. Rev. Lett. 123, 100502 (2019).
  * Tamascelli  _et al._ [2019] D. Tamascelli, A. Smirne, J. Lim, S. F. Huelga, and M. B. Plenio, Efficient Simulation of Finite-Temperature Open Quantum Systems, Phys. Rev. Lett. 123, 090402 (2019).
  * Mascherpa  _et al._ [2020] F. Mascherpa, A. Smirne, A. D. Somoza, P. Fernández-Acebal, S. Donadi, D. Tamascelli, S. F. Huelga, and M. B. Plenio, Optimized Auxiliary Oscillators for the Simulation of General Open Quantum Systems, Phys. Rev. A 101, 052108 (2020).
  * Gelinas  _et al._ [2014] S. Gelinas, 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 (2014).
  * Nan  _et al._ [2015] G. Nan, X. Zhang, and G. Lu, Do “Hot” Charge-Transfer Excitons Promote Free Carrier Generation in Organic Photovoltaics?, J. Phys. Chem. C 119, 15028 (2015).
  * Dennler  _et al._ [2009] G. Dennler, M. C. Scharber, and C. J. Brabec, Polymer-Fullerene Bulk-Heterojunction Solar Cells, Adv. Mater. 21, 1323 (2009).
  * Cheung and Troisi [2010] D. L. Cheung and A. Troisi, Theoretical Study of the Organic Photovoltaic Electron Acceptor PCBM: Morphology, Electronic Structure, and Charge Localization, J. Phys. Chem. C 114, 20479 (2010).
  * Idé  _et al._ [2014] J. Idé, D. Fazzi, M. Casalegno, S. V. Meille, and G. Raos, Electron transport in crystalline PCBM-like fullerene derivatives: A comparative computational study, J. Mater. Chem. C 2, 7313 (2014).
  * Chen  _et al._ [2018b] 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 (2018b).
  * Chen and Brédas [2018] X.-K. Chen and J. Brédas, Voltage Losses in Organic Solar Cells: Understanding the Contributions of Intramolecular Vibrations to Nonradiative Recombinations, Adv. Energy Mater. 8, 1702227 (2018).
  * Menke  _et al._ [2018b] S. M. Menke, A. Cheminal, P. Conaghan, N. A. Ran, N. C. Greehnam, G. C. Bazan, T.-Q. Nguyen, A. Rao, and R. H. Friend, Order Enables Efficient Electron-Hole Separation at an Organic Heterojunction with a Small Energy Loss, Nat. Commun. 9, 277 (2018b).
  * Hinrichsen  _et al._ [2020] T. F. Hinrichsen, C. C. S. Chan, C. Ma, D. Paleček, A. Gillett, S. Chen, X. Zou, G. Zhang, H.-L. Yip, K. S. Wong, R. H. Friend, H. Yan, A. Rao, and P. C. Y. Chow, Long-lived and Disorder-free Charge Transfer States Enable Endothermic Charge Separation in Efficient Non-Fullerene Organic Solar Cells, Nat. Commun. 11, 5617 (2020).
  * Gregg [2011] B. A. Gregg, Entropy of Charge Separation in Organic Photovoltaic Cells: The Benefit of Higher Dimensionality, J. Phys. Chem. Lett. 2, 3013 (2011).
  * Ono and Ohno [2016] S. Ono and K. Ohno, Combined Impact of Entropy and Carrier Delocalization on Charge Transfer Exciton Dissociation at the Donor-Acceptor Iolution in multidimensional ultrnterface, Phys. Rev. B 94, 075305 (2016).
  * 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 (1989).
  * Tanimura [2006] Y. Tanimura, Stochastic Liouville, Langevin, Fokker–Planck, and Master Equation Approaches to Quantum Dissipative Systems, J. Phys. Soc. Jpn. 75, 082001 (2006).
  * Rivas and Huelga [2012] Á. Rivas and S. F. Huelga,  _Open Quantum Systems_, SpringerBriefs in Physics (Springer Berlin Heidelberg, 2012).