Exciton Dissociation in a Model Organic Interface: Excitonic State-Based Surface Hopping versus Multiconfigurational Time-Dependent Hartree.

Quantum dynamical simulations are essential for a molecular-level understanding of light-induced processes in optoelectronic materials, but they tend to be computationally demanding. We introduce an efficient mixed quantum-classical nonadiabatic molecular dynamics method termed eXcitonic state-based Surface Hopping (X-SH), which propagates the electronic Schrödinger equation in the space of local excitonic and charge-transfer electronic states, coupled to the thermal motion of the nuclear degrees of freedom. The method is applied to exciton decay in a 1D model of a fullerene-oligothiophene junction, and the results are compared to the ones from a fully quantum dynamical treatment at the level of the Multilayer Multiconfigurational Time-Dependent Hartree (ML-MCTDH) approach. Both methods predict that charge-separated states are formed on the 10-100 fs time scale via multiple "hot-exciton dissociation" pathways. The results demonstrate that X-SH is a promising tool advancing the simulation of photoexcited processes from the molecular to the true nanomaterials scale.

Organic solar cells (OSCs) have attracted much research interest over the last few decades because of a number of beneficial features, such as low manufacturing cost, mechanical flexibility, light weight, and environmentally friendly materials. One of the major parameters determining the performance of a solar cell is the power conversion efficiency (PCE), and the best OSC designs to date have reached over 19%.[1−3] To further optimize OSC materials and devices, an improved understanding of the fundamental optoelectronic processes is essential.

Various spectroscopy approaches have been applied to elucidate the charge separation mechanisms in OSCs.[4−11] However, because of the complexity of such systems, contradictory results have been reported in the literature. One prominent example is the role of the lowest-lying charge transfer (CT1) state, namely, the bound electron–hole state at the donor–acceptor interface. Though many studies suggest it is the gateway state for charge separation independently of the time scale of exciton dissociation,[4−8] other works claim that the charge separation is mainly through ultrafast hot carrier dissociation instead of the CT1 dissociation.[9−11] It has been suggested that these discrepancies are primarily because of differences in the interfacial energetics, which are in turn because of the different molecular constituents, compositions, and diverse fabrication processes adopted that cause distinct morphology for donor and acceptor domains,[12,13] which subsequently alter the parameters that govern the light-induced processes in the OPVs.

Theoretical modeling, however, can be carried out for structurally well-defined and increasingly realistic model systems. This allows us, in principle, to investigate how the mechanism of charge generation and separation depends on structural and energetic parameters.[14−23] Common theoretical models often describe the excited-state processes in OSC materials in terms of rate theories that have been developed for molecular donor–acceptor systems (e.g., Marcus–Levich–Jortner). When applied to materials, they typically come with a number of restrictive assumptions. For instance, they neglect spatial delocalization of excitons and charge carriers, which is expected to have a profound effect on exciton dissociation and charge transport efficiencies.

Direct quantum dynamical simulation of the coupled electron–nuclear motions that account for such phenomena would be highly desirable, but they tend to be computationally expensive when applied to materials.[11,24−29] In this regard, important progress has been made by carrying out Multiconfiguration Time-Dependent Hartree (MCTDH) and related Multilayer MCTDH (ML-MCTDH) calculations,[30,31] for the quantum dynamics of nonadiabatic nuclear motion on vibronic coupling models parametrized to excited-state electronic structure calculations.[32] A variety of complex processes in model systems could be treated in this way, e.g., exciton diffusion in an oligothiophene chain,[33] electron–hole separation in donor–acceptor conjugated co-oligomers,[34] and charge separation in the donor–acceptor junctions.[35,36]

Widely considered the gold standard for quantum dynamics, MCTDH is computationally expensive. Notwithstanding the impressive progress that has been recently made to increase the number of vibrational modes that can be treated in ML-MCTDH (several hundreds),[33−36] applications to truly nanoscale systems (>10 nm) that require several hundreds if not thousands of molecules to be simulated at quantum level is likely to be out of reach for some time. Computationally more efficient approximate methods, e.g., mixed quantum-classical nonadiabatic molecular dynamics (MQC-NAMD) such as decoherence-corrected trajectory surface hopping,[37] are a viable alternative for large systems.[38−42] Some of us have recently developed surface hopping methodologies (termed Fragment Orbital-Based Surface Hopping (FOB-SH[43−45]) and Frenkel Exciton Surface Hopping (FE-SH[46])) for the quantum dynamical propagation of charge carriers or Frenkel excitons in nanoscale materials, e.g., organic crystals[46−48] and (disordered) thin films.[49] Experimental charge mobilities and exciton diffusion constants could be well reproduced for a range of materials, and the simulations revealed a clear correlation between the spatial delocalization of the charge carrier or exciton and the respective diffusion constants.

In this work, we have implemented an extension of our previously developed FOB-SH and FE-SH methods by combining and coupling the electronic-state spaces for charge transfer and localized (Frenkel) exciton states. This enables us to simulate exciton transport coupled to exciton dissociation, charge generation and charge separation at the level of trajectory surface hopping, which was not possible with either FOB-SH or FE-SH. We dub our extension eXcitonic state-based Surface Hopping (X-SH) to emphasize the applicability of the method to light-induced excited-state processes in molecular materials. For simplicity, we consider a 1D chain composed of N electron donating chromophore molecules and one effective electron-accepting molecule, as shown in Figure a, and only nearest-neighbor couplings. Generalization to larger systems in 2D or 3D and to couplings beyond nearest neighbors is straightforward. In X-SH, the total Hamiltonian is constructed aswhere Tcl is the classical kinetic energy of the nuclei, is the electronic Hamiltonian, and is the electron–phonon coupling term. can be divided into XT and CS blocks and a XT–CS off-diagonal block,whereand

Figure 1

(a) Model system for the fullerene–oligothiophene single junction. (b) Electronic Hamiltonian matrix in the basis of diabatic Frenkel exciton (XT) and charge-separated states (CS). See the main text for definition of the matrix elements. (c) Energies and diabatic-state populations of the adiabatic electronic states at the electronic ground-state minimum geometry (absolute values squared of the eigenvectors, in orange for XT and in blue for CS diabatic states). The sidebar to the right depicts the projection of the diabatic XT1 state, used as the initial state for X-SH and MCTDH simulations, onto the adiabatic eigenstates.

In eq , εXT is the electronic energy offset between diabatic XT and CS1 states, where εXT = EXT – ECS in the minimum energy configuration of the ground state (RGS). εCS is the electronic energy of diabatic states CS, where donor molecule n carries a positive charge and the acceptor molecule a negative charge at configuration RGS. The local electron–phonon part is written as

The term uXT/CS(R) is the electronic energy of state n at a nuclear configuration R, relative to the energy of that state at RGS. Thus, EXT(R) = εXT + uXT(R) is the electronic energy of state XT at nuclear configuration R, i.e., the total potential energy of the system when site n is in the excited state and all other sites are in the (neutral) ground electronic state. Similarly, ECS(R) = εCS + uCS(R) is the energy of state CS at nuclear configuration R, i.e., the total potential energy of the system when the acceptor molecule is negatively charged and the nth donor molecule is positively charged while all other molecules are in the neutral ground electronic state. Here, all electronic and excitonic couplings are fixed; i.e., nonlocal electron–phonon coupling is not included. Moreover, energetically accessible excited states other than the ones considered, e.g., charge transfer excitonic states within the donor manifold, triplet states, etc., are not considered but could, in principle, be included by extension of the electronic-state space.

The time-dependent wave function is expressed as a linear combination of the XT and CS stateswhere and , n = 1, N, and u is the corresponding expansion coefficient. The wave function is propagated according to the time-dependent Schrödinger equation,where are the elements in the electronic Hamiltonian described above. The nonadiabatic coupling elements (NACE) between the quasi-diabatic states, , are typically very small and are neglected. We have shown in our previous study on charge transport that neglecting the NACE gives essentially the same dynamics but accelerates the calculations considerably.[48] The nuclei are propaged on a single adiabatic potential energy surface (the active state) at every time step and are allowed to hop stochastically between different surfaces according to Tully’s fewest switching hopping probablility.[37] The Heisenberg decoherence time for the exponential damping of the coefficients of all except the active electronic states is applied to correct the overcoherence problem in the surface hopping method.[50,51] A state tracking algorithm is used to detect trivial crossings between the adiabatic states.[44]

Here we apply X-SH to simulate the exciton decay dynamics in an one-dimensional (1D) model of a fullerene–oligothiophene interface, which was previously investigated at the level of MCTDH.[35,36] We also report new ML-MCTDH computations significantly extending the simulation time from previously 1 to 10 ps. Both approaches employ the same electronic Hamiltonian (see eqs –5 and equation (S1)), which enables us to compare X-SH with MCTDH dynamics on an equal footing for a complex, application-relevant system. This allows us to assess the trade-offs that the two approaches make: classical treatment but inclusion of all nuclear degrees of freedom (DoF) in X-SH versus full quantum mechanical treatment of a limited number of DoF in MCTDH. In particular, it affords an assessment of the effect of classical vs quantum treatment of nuclear motion on the exciton relaxation mechanism and dynamics. A certain caveat is that the MCTDH calculations, which include zero-point energy, were restricted to T = 0 K whereas X-SH calculations were carried out at a classical temperature of T = 300 K. In this respect, it is worth noting that the classical nuclei at 300 K may actually lead to a narrower distribution of initial conditions in the phase space than that obtained from sampling a Wigner distribution for the ground vibrational state.[52]

The atomistic model of the oligothiophene (OT9)–fullerene (C60) donor–acceptor interface is shown in Figure a. In X-SH, the system is composed of 3 C60 molecules and 26 OT9 molecules with periodic boundary conditions along the direction of the stack. However, as in MCTDH, only one C60 unit and 13 OT9 molecules closest to that fullerene are treated as electronically active. The molecules in the inactive region maintain the structural integrity of the electronically active region during X-SH molecular dynamics, thus effectively replacing the bulk environment. Consequently, 13 (quasi-)diabatic Frenkel excitonic states in the donor phase are considered (eq ). The state where the Frenkel exciton (XT) is located on the nth OT9 molecule is labeled XT, n = 1, ..., 13. Similarly, the charge-separated (CS) state (eq ) with the hole on the nth OT9 molecule and the electron on the C60 molecule is labeled CS. The corresponding electronic Hamiltonian matrix for the 26 electronic states is shown in Figure b. The site energies (diagonal elements) for the XT and CS states were calculated with a fully atomistic force field in X-SH (eq and Supporting Information equations (S6)–(S11)), and with a linear vibronic coupling potential along eight effective modes for each molecule in MCTDH[27,37] (Supporting Information equations (S3)–(S5)). The atomistic force field for XT and CS states was parametrized to reproduce exactly the total reorganization energy for exciton transfer and for charge transfer used in MCTDH; see the Supporting Information for details. Furthermore, in both approaches, XT and CS1 states are offset by 0.1 eV at the Franck–Condon point (ΔEoffset = εXT – ε1CS = 0.1 eV), and the Coulomb energy of the CS states was modeled by an effective Coulomb barrier (i.e., “Coulomb barrier 1” of ref (36) representing an OT13-(C60)7 aggregate; shown in Figure S1 and labeled as εCS in Figure b for site i). The coupling (or off-diagonal) elements in the electronic Hamiltonian between nearest neighbors were fixed to the values used in ref (35); all other coupling matrix elements were set to zero. For XT–XT couplings, j = 0.1 eV, for CS–CS couplings, t = −0.12 eV, and for the XT–CS coupling at the interface, λ = −0.2 eV (note the sign for the latter was incorrectly reported in ref (35)). ΔEoffset, the Coulomb barrier, and coupling values are based on previously reported time-dependent density functional theory (TDDFT) and DFT calculations and a diabatization scheme using adiabatic states from TDDFT as basis set.[53] They are reproduced in the SI for convenience.

Prior to reporting the results of the nonadiabatic dynamics simulations, we analyze the electronic eigenstates (or adiabatic electronic states) of our model interface as obtained by diagonalization of the electronic Hamiltonian at the electronic ground-state minimum geometry. In Figure c, the eigenstates are projected on the diabatic (site) basis of XT (orange) and CS states (blue). The lowest-lying eigenstate in the system is the interfacial charge transfer state, which we label CT1 hereafter. It is mostly composed of the CS1 (54%) and XT1 (31%) diabatic states. The energy of this state, −0.20 eV, is lower than those for CS1 (0 eV) and XT1 (0.1 eV) because of the electronic coupling between these two states, λ = −0.2 eV and, to a lesser extent, between CS1 and CS2 states. (A simple two-state model including CS1 and XT1 states only would give an energy E(CT1) = 1/2(ΔEoffset – (ΔEoffset2 + λ2)0.5) = −0.16 eV.) The XT manifold of eigenstates is centered at ΔEoffset = 0.1 eV, and the bandwidth is ∼4j = 0.4 eV, as expected from a tight binding model. Thus, the lower edge of the XT band is at −0.09 eV, 0.11 eV above CT1. It is mainly composed of XT5–XT9 states and denoted as the dark XT state because an H-aggregate of OT9 molecules (1D packing with positive XT–XT couplings) is considered, where the lowest state is optically inaccessible.[35] The upper edge at 0.29 eV is formed of the bright XT state, which is mainly composed of XT6–XT9. Finally, the CS manifold of eigenstates is centered at 0.24 eV, which is slightly below the maximum of the Coulomb barrier, and the bandwidth is ∼4t = 0.48 eV. The lower CS band edge is formed by a charge-separated state (CSS) composed of diabatic CS states that are furthest away from the interface, CS8–CS13. The energy of this state, 0.00 eV, defines the binding energy of the interfacial charge transfer state CT1, 0.20 eV.

We have carried out nonadiabatic dynamics simulations choosing the fully localized XT1 diabatic state as initial state and report the results in Figure . X-SH and MCTDH predict a similar decay dynamics of the XT1 state at short times (<200 fs, Figure a); within 100 fs the populations of the dark XT state and the CSS accumulate markedly, to 0.13 (0.17) and 0.11 (0.11) at 100 fs, respectively, for X-SH (MCTDH). Beyond 100 fs, some differences are clearly discernible. In MCTDH, the population of the dark XT state keeps increasing, up to 0.25 at 1 ps, whereas in X-SH, it decays after 100 fs by converting to the CT1 state (as indicated by the diabatic CS1 population in Figure ). The population difference between the dark XT state from the two approaches is very close to the difference in the CT1-state population (as indicated by two arrows of similar length in Figure b). At longer times, >1 ps, the populations for the dark XT and CSS states decay in MCTDH, but significant populations remain at 10 ps (Figure c and Figure S1). At this point, the quantum equilibrium population of electronic states (from imaginary time propagation; see Figure S4) has not been reached. By contrast, in X-SH, both dark XT and CSS have fully decayed after a few picoseconds. This establishes the classical equilibrium population for all electronic states in the system to a good approximation (dash-dotted lines in Figure c; see the Supporting Information for their calculation), showing that the newly developed X-SH method satisfies the detailed balance in the long-time limit, similarly to FOB-SH[50] and FE-SH.[46]

Figure 2

Population dynamics of the interfacial state (CS1), the charge-separated state (CSS), and the dark XT state (dark XT) from X-SH (solid lines) and MCTDH (dashed lines) simulations, initialized from the XT1 diabatic state. (a)–(c) show the dynamics up to 200 fs, 1 ps, and 10 ps, respectively. The horizontal lines (dash-dotted) in (b) and (c) indicate the thermal equilibrium (Boltzmann) populations of the electronic states at 300 K. (c) is reproduced in linear scale in Figure S2. The nuclear and electronic time steps are 0.1 and 0.02 fs, respectively. The convergence of the populations with respect to time step can be found in Figure S3.

The fast decay dynamics in X-SH after 100 fs could be an artifact of the classical treatment of nuclear motion. Similar observations have been reported by Freixas et al., who found that (decoherence-corrected) surface hopping exhibited a faster deexcitation dynamics than their reference method, multiconfigurational Ehrenfest with ab initio multiple cloning (whereas Ehrenfest, unsurprisingly, gave a too slow decay).[54] However, the decay dynamics in MCTDH could be somewhat underestimated because of the zero-temperature conditions and the limited number of modes included in the model. During radiationless decay of XT1 to the CT1 state, up to 0.5 eV of electronic energy needs to dissipate via the nuclear degrees of freedom. Because only eight intramolecular modes per molecule (112 in total) are included, in the absence of intermolecular modes, the excess electronic energy cannot efficiently dissipate. This is likely to result in a slight overestimation of excited electronic states (dark XT, CSS) and prevents full relaxation to the equilibrium populations. In X-SH, all nuclear DoF, i.e., all intra- and intermolecular modes, are included (5604 in total), in particular the low-frequency modes that, as we show below, act as a heat bath absorbing the excess electronic energy and establishing thermal equilibrium.

In the following, we trace the energy flow to the nuclear vibrations during the electronic relaxation process in X-SH employing a normal-mode analysis (NMA), similarly as in ref (55). To this end, the total kinetic energy is projected on the 189 intramolecular normal modes of each OT9 molecule in the electronically active region, further details are given in the Supporting Information. The kinetic energy of each mode averaged over all X-SH trajectories is shown in Figure as a function of time. Most of the modes fluctuate stably around the mean thermal energy, kBT/2 = 13 meV at 300 K, at any time (Figure b). These modes do not couple with the electronic transitions and can be interpreted as spectator or bath modes during the de-excitation process. Only six intramolecular high-frequency modes (depicted in Figure S5) exhibit significant excess kinetic energy during the first ∼500 fs of electronic relaxation; see Figure a. This is followed by dissipation of the excess kinetic energy to the spectator modes. After a few picoseconds, the vibrational energy redistribution is essentially complete. Thus, our analysis of the X-SH dynamics agrees with the assumption made in MCTDH simulations that only a few high-frequency modes are important in the relaxation process.[35,36] In fact, the six intramolecular modes identified in Figure a have frequencies (734–3228 cm–1) that are similar to those of the six highest frequency modes selected in MCTDH (902–3238 cm–1). However, our results also highlight the importance of including a sufficient number of bath modes to reach the equilibrium distribution.

Figure 3

Dissipation of electronic excitation energy by vibrational modes in X-SH simulations. Normal modes with kinetic energy exceeding 14 meV in at least one instance within the first 2 ps of relaxation dynamics are shown in (a) and denoted “active” modes. The kinetic energy fluctuations of all other normal modes, denoted “spectator” modes, are shown in (b). The six “active” modes are indicated in Figure S5. The kinetic energy of each mode was averaged over 1000 X-SH trajectories; the thermal average, kBT/2, is 13 meV at 300 K.

Analyzing 1000 X-SH trajectories, we identify seven distinct exciton decay pathways to the energetically lowest CT1 state, which we briefly summarize in the following (see also Table with percentage of occurrence given and Figure for representative X-SH trajectories). (1) As the initial XT1 diabatic state projects to 34% on the CT1 adiabat (see the sidebar in Figure c), it can convert directly to CT1 through decoherence. XT1 is also observed to relax to CT1 (2) via the XT manifold, Figure a, (3) by the CS manifold without formation of the CSS state, Figure b, or (4) by first passing through the XT and then through the CT manifold or (5) vice versa or (6) via XT–CT hybrid states. Finally, (7) XT1 can form the charge-separated state (CSS) by passing through the CS manifold, Figure c.

Table 1

Seven Pathways for Exciton Relaxation from the XT1 State, as Obtained from X-SH Simulationa

	pathways	percentage
(1)	XT1 → CT1	33.7
(2)	XT1 → {XTn} → CT1	24.7
(3)	XT1 → {CSn} → CT1	6.3
(4)	XT1 → {XTn} → {CSn} → CT1	0.5
(5)	XT1 → {CSn} → {XTn} → CT1	3.8
(6)	XT1 → {XTn, CSn} → CT1	17.0
(7)	XT1 → {CSn} → CSS → CT1	14.0

The percentage gives the number of trajectories within the ensemble of 1000 X-SH trajectories following the pathway indicated. A trajectory is classified as reaching the {XT}, {CS}, or CSS manifold when the sum of XT1 to XT13, CS1 to CS13, or CS8 to CS13 populations is >95%, >95%, or >80%, respectively, in at least one instance in time.

Figure 4

Three pathways for exciton dissociation as obtained from X-SH: (a) XT1 → {XT} → CT1, (b) XT1 → {CS} → CT1, and (c) XT1 → {CS} → CSS → CT1. Initially, all population is on the XT1 state (black bar). The XT manifold {XT} is color-coded alternately in red and gray, and the CS manifold {CS} in blue and gray. XT1 to XT13 are arranged from center to right, and CS1 to CS13 populations are arranged from center to left. The pathways in (a), (b), and (c) are denoted as pathways (2), (3), and (7) in Table .

The formation of the desired CSS state via pathway (7) may be best described as “hot carrier dissociation”. However, because of the finite system size of our model, virtually all of the generated CSS population converts to CT1 within a few picoseconds. Because the charge carriers in the CSS state are beyond the maximum of the Coulomb barrier, in a larger model, they would become free carriers, swept toward the electrodes by the electric field or recombining bimolecularly at longer times (geminate recombination). Once the system has relaxed to CT1 via pathways (1)–(7), the escape from that state is negligible (<1%) on the time scale of present simulations (10 ps), unsurprisingly because CT1 is >0.2 eV below the next highest states, CSS and the dark XT state in our model. Thus, for the given model parameters employed, our simulations predict that charge carrier generation occurs exclusively via the “hot carrier dissociation” pathway on the 10 ps time scale of present simulations. Though this does not exclude the possibility for dissociation of CT1 to the CSS state (“cold carrier dissociation” pathway), which may well occur on longer, nanosecond time scales and other material parameters.

In summary, we have extended our surface hopping methodology to enable the simulation of exciton dissociation to charge carriers in molecular materials. We found that the developed X-SH approach agrees very well with ML-MCTDH in predicting the decay dynamics and state populations at ultrafast time scales (<100–200 fs) as well as the intramolecular modes to which the excess energy is transferred. Thereafter, X-SH predicts a fairly rapid decay to the equilibrium populations of electronic states in a few picoseconds, not seen in MCTDH, possibly because of the lack of low-frequency bath modes or missing temperature effects in the latter calculations. MCTDH calculations incorporating those effects might clarify this issue in the future. On the whole, nuclear quantum effects such as nuclear tunneling and zero point motion that are missing in X-SH, do not seem to be overly important for the prediction of the population dynamics, at least for the current model at room temperature. This is likely because of the strong electronic coupling between the states in our model facilitating relaxation to the ground state without having to cross high barriers where tunneling and zero point motion are likely to become important. Thus, we conclude that X-SH appears to be a promising method for the simulation of exciton relaxation processes in this regime that is characteristic for organic optoelectronic materials.

A major advantage of X-SH compared to fully quantum dynamical methods is its computational efficiency permitting the simulation of realistic 2D and possibly 3D nanoscale systems (>10 nm) comprising hundreds or more molecular sites. Moreover, it is straightforward to include in a realistic manner many of the effects that have not been included in the present study: charge delocalization in acceptor and donor, thermal fluctuations of excitonic and electronic couplings (i.e., off-diagonal electron–phonon couplings) in donor and acceptor phases and at the interface, recombination to the electronic ground state, influence of interface geometry and static disorder. Work in our laboratory is currently ongoing to include these effects in X-SH.