Revealing generation, migration, and dissociation of electron-hole pairs and current emergence in an organic photovoltaic cell.

Using an innovative quantum mechanical method for an open quantum system, we observe in real time and space the generation, migration, and dissociation of electron-hole pairs, transport of electrons and holes, and current emergence in an organic photovoltaic cell. Ehrenfest dynamics is used to study photoexcitation of thiophene:fullerene stacks coupled with a time-dependent density functional tight-binding method. Our results display the generation of an electron-hole pair in the donor and its subsequent migration to the donor-acceptor interface. At the interface, electrons transfer from the lowest unoccupied molecular orbitals (LUMOs) of thiophenes to the second LUMOs of fullerene. Further migration of electrons and holes leads to the emergence of current. These findings support previous experimental evidence of coherent couplings between electronic and vibrational degrees of freedom and are expected to stimulate further work toward exploring the interplay between electron-hole pair (exciton) binding and vibronic coupling for charge separation and transport.

INTRODUCTION

Organic solar cells (OSCs) have attracted considerable research interest over the past decades owing to their low cost, light weight, and potential use in flexible devices. After more than two decades of development, the power conversion efficiency (PCE) of OSCs has reached more than 16% in single-junction devices (, ). This value is, however, lower than that achieved in silicon-based single-junction solar cells (). One crucial factor that makes the difference in PCEs between the OSCs and the conventional inorganic photovoltaics is the exciton binding energy (). In inorganic materials, free charge carriers are generated directly after light absorption since the binding energy [Δ; e.g., Δ ≈ 25 meV ≈ kBT in silicon ()] can be easily overcome by thermal energy at room temperature. In organic materials, however, spatially localized and strongly Coulomb-bound electron-hole pairs (Frenkel excitons) are generated immediately upon light absorption. In high-performance OSCs, the donor-acceptor (D-A) heterojunction provides an internal electrochemical force that drives exciton dissociation. Various additional mechanisms, such as delocalization () and vibronic couplings (, ), have been proposed, and in particular, the important role of vibronic couplings for the dissociation of excitons and separation of charges has found much experimental and theoretical evidence (–). A detailed analysis of the effect of the vibronic couplings on charge separation across the D-A interface is given in (). These mechanisms roughly refer to two categories, beneficial energetics and advantageous morphology (). It is generally believed that the following processes occur to transform the energy of absorbed sunlight into electrical power: Photoinduced Frenkel excitons diffuse to the D-A interfaces and dissociate into holes on the donor and electrons on the acceptor (, ), and the resulting free carriers migrate through the conducting channels and thus give rise to photocurrent. Efficiencies of all these steps, photoexcitation of Frenkel exciton, exciton diffusion to the D-A interface, exciton dissociation into free electrons and holes at the interface, and free charge carrier migration to the cathode and anode, determine the overall cell performance (). However, the direct observation of the individual steps remains challenging. In this work, we use an innovative computational method, real-time time-dependent density functional tight binding for open systems (TDDFTB-OS), combined with Ehrenfest molecular dynamics (MD) to simulate the electron-hole pair photogeneration in the donor, electron-hole pair migration to and dissociation into free electrons and holes at the D-A interface, and the subsequent free charge transport and occurrence of photocurrent in electron- and hole-conducting materials.

The accurate modeling of light-to-current conversion processes in OSCs requires a theoretical framework that can adequately describe the evolution of charge carriers with coupled nuclear dynamics. Nonadiabatic MD formalisms are crucial for describing electron-nuclear dynamics in OSCs. Surface hopping () and the Ehrenfest method () are the two mostly used and documented formalisms (). Ehrenfest MD is used in this work. Considering the large system size of OSCs and the time scale of photoinduced current processes, the electronic dynamics is simulated with the real-time TDDFTB method (). Similar to the semiempirical quantum mechanical methods (, ), TDDFTB contains effective electron-electron Coulomb interaction with a 1/r asymptote as the intercharge distance goes to infinity and thus provides correct descriptions of the long-range behavior of the charge carriers (). The DFTB parameter used in this work is obtained by fitting the density functional theory (DFT) results with the Perdew–Burke–Ernzerhof (PBE) functional that belongs to the class of generalized gradient approximation (GGA) type functionals. It provides estimates of the Coulomb correlation between electrons and holes, respectively, but lacks quantitatively correct description of the Coulomb interaction within optically excited electron-hole pairs. As such, the electron-hole pairs modeled in this work cannot be quantitatively considered as Frenkel excitons. Accurate treatment of excitons can be achieved by using methods with extra correlations [e.g., DFTB+U () or LC-DFTB ()] and is currently ongoing. It is well known that the PBE functional does not describe charge transfer excitations well as its exchange interaction does not approach the 1/r limit asymptotically. In the present work, we use the TDDFTB-OS method for the separation of electrons and holes across the D-A interface. As an electron transfers across the D-A interface, another electron would exit the simulation system via the boundary at the acceptor side to keep the charge balance of the acceptor region; in the meantime, an electron enters the acceptor side of the simulation system to keep the charge balance of the donor region. Thus, the charge separation across the D-A interface in the open quantum system is different from the charge transfer state of the isolated system. The correct alignment of the single-electron energy levels across the interface is needed to account for the charge separation process. The effective electron-electron Coulomb interaction with a 1/r asymptote describes well the charge separation across the D-A interface in the TDDFTB-OS calculation. Note that the validity of DFTB parameters has been confirmed by comparing the calculated absorption spectra to the experimental data, as shown in fig. S2 in our previous work ().

The realistic OSCs consist of different functional components. The core processes that we are interested in occur in the active layer composed of donor and acceptor materials. Because of the limitation of computational resources, we can only model a small representative volume of the active layer, which consists of a small section of donor (or hole conducting) materials, a small section of acceptor (or electron conducting) materials, and the D-A interface. To model the electrodes and electron-/hole- conducting channels correctly, the open boundary is introduced, and the representative active layer is modeled as an open quantum system. For this, the time-dependent density functional theory for open system (TDDFT-OS) is used (–). In the calculation, DFT is replaced by DFTB, and thus, our method is termed Ehrenfest TDDFTB-OS.

Previous studies have been reported to analyze the dynamics of electrons and holes, such as the simple yet powerful HOMO (highest occupied molecular orbital)–LUMO (lowest unoccupied molecular orbital) analysis. In addition, the transition density matrix (TDM) provides spatial information about electron-hole pairs for specific excited states (–). The natural transition orbitals technique builds a compact orbital representation for the TDM (). The particle-hole map identifies where the electron is excited and transferred (). Energy and charge transfer in organic functional materials has also been widely investigated theoretically using the multiconfiguration time-dependent Hartree method and its variant (). These approaches are more suitable for isolated systems. Their analysis focuses more on the energy domain and provides relatively coarse-grained spatial information. Thus, extra effort is required to further interpret the dynamics. To obtain a more straightforward description of electron and hole dynamics in OSCs, we use the real-space induced Mulliken charge () distribution method and the fragment molecular orbitals (FMOs) () occupation method to analyze the electronic dynamics obtained with the Ehrenfest TDDFTB-OS formalism. We use Mulliken charge to reveal the charge fluctuations in real space in accordance to that used in DFTB. The different FMOs are used to investigate the charge transfer between fragments through their occupations. These two analysis methods are complementary and focus on spatial and energetic information of excited-state dynamics, respectively. As such, they allow us to study the light-induced energy conversion mechanisms from multiple perspectives. We consider a six-ring oligothiopene (T6) and a C60 molecule (see Fig. 1A) as a minimal representative fragment of the widely investigated blend of the donor poly-3-heylthiophene (P3HT) and acceptor [6,6]-phenyl-C61 butyric acid methyl ester, respectively (, –). The simulation region consists of three T6 and three C60 units and is sandwiched between semi-infinite repeating donor and acceptor molecular units. Here, we provide a first principles method to study the generation and migration of electron-hole pairs and of the separated electrons/holes from both spatial and energetic perspectives in the open system. These calculation and analysis methods of open systems provide an advanced research foundation and tools for future study on photovoltaic materials.

Fig. 1

Atomistic model of the system and analysis approaches.

(A) Illustration of the system setup, where donor T6 and acceptor C60 molecular units in the simulation region (D) are named T6-1 to C60-3. During the simulation, the incident femtosecond laser pulse is applied only on T6-2 along the y direction. (B) Induced charge distribution projected on real-space grids (vertices of the gray wireframe); the red cloud represents the photoinduced electron distribution and the blue cloud represents the photoinduced hole distribution. (C) σFMO obtained by projecting the reduced single-electron density matrix σD onto the FMO basis. A specific block of σFMO corresponding to the T6-2 unit is enlarged for illustration. is further partitioned according to the FMO: occupied (blue) and virtual (red) FMOs. The diagonal elements correspond to the occupation number of the FMOs while off-diagonal elements correspond to the electronic coherence between them.

RESULTS

The D-A structure depicted in Fig. 1A was simulated with the Ehrenfest TDDFTB-OS formalism for 100 fs. We carry out real-space analysis (Fig. 1B) and energy-resolved analysis (Fig. 1C) based on the induced Mulliken charge and FMO occupations to obtain details of the charge transfer dynamics.

Real-space analysis

We first visualize charge transfer processes in real space by analyzing the real-time electronic dynamics with the induced Mulliken charge distribution method (see Eqs. 9 to 11). The applied laser pulse, depicted in Fig. 2C, peaks around 10 fs with a full width at half maximum (FWHM) of the Gaussian envelope of 4 fs. In Fig. 2A, we plot six snapshots of Δρ(, t) (see movie S1 for the video of the light-induced dynamics). For better illustration, Δρ(, t) is averaged over t ± 0.2 fs to filter out high-frequency fluctuations. Increase (red) and decrease (blue) of Δρ reflect the distribution of photoexcited electrons and holes in real space. The time-dependent electric current passing through the left (SL) and right (SR) surfaces is shown in Fig. 2B, in which the snapshots chosen in Fig. 2A are marked with vertical dashed lines. At 7.6 fs, i.e., still during the laser pulse, we observe the generation of an intramolecular electron-hole pair localized at T6-2. Subsequently, the electron-hole pair hops to the adjacent donor units and delocalizes over the T6-1 and T6-3 units at ~10 fs (see the second dashed line in Fig. 2B). When the electron-hole pair reaches the left surface SL (adjacent to T6-1), the left current is detected at SL. Note that the left current at SL oscillates rapidly with a frequency ~4 eV that corresponds to twice the excitation energy (see Fig. 2B). The oscillation of the left current is caused by the electron and hole spatial oscillation in the pair. Positive current is defined as incoming electrons from the environment, i.e., holes leaving the device; negative current corresponds to the electron leaving the device. The net current at SL is positive, meaning that overall holes leave the device through SL. Note that after ~60 fs, the current through SL oscillates around zero. The right current through SR appears at a much later time as compared to the left current through SL and reaches and remains around its saturation value when the left current oscillates around the zero value. In addition, the right current is always negative, as it is due to electron migration through SR. Therefore, the left current is caused by the electron-hole pair transport, while the right current is purely due to electron transport. More specifically, at ~11 fs, electrons start to appear at C60-1. Meanwhile, the oscillation of the left current displays larger positive values, representing that more holes are flowing out from the left surface SL. After that, electrons gradually build up on C60-2 and C60-3. When electrons reach the right surface at 31.6 fs (fourth dashed line), the negative valued right current (red line) through SR starts to emerge, indicating that electrons are leaving the simulation region through SR. As more electrons transfer to the acceptor units and reach SR (see t = 50.8 fs of Fig. 2A), the right current rapidly increases. At t = 60 fs, the right current reaches its saturation value (~50 nA), which persists beyond t = 100 fs (see Fig. 2B). Figure 2A shows the distribution of induced electrons and holes at t = 96 fs, and the electron-hole pairs reside clearly across the D-A interface.

Fig. 2

Real-space induced charge distribution of light-induced charge carrier dynamics.

(A) Induced charge distribution Δρ() at selected time shots. The amplitude of Δρ() of acceptors in the bottom panel is multiplied by 4 for better illustration. (B) Left (cyan) and right (red) time-dependent electric current. Positive current refers to electrons passing through the surface from the environment into the simulation region and vice versa. (C) Impulsive external electric field used to excite the T6-2 oligomer along the y direction.

FMO analysis

In Fig. 3A, we plot the density of states (DOS) of semi-infinite stacks of T6 (blue line) and C60 (orange line) that serve as the respective hole- and electron-conducting layers. The local DOS (LDOS) of the simulation region at steady state is shown in the middle panel. The Fermi level is shifted to zero. Here, we denote the highest occupied and lowest unoccupied states of donor and acceptor as HOS and LUS, respectively. The next higher unoccupied states are denoted as LUS + 1, accordingly. It is shown that the energy gap between the donor’s LUS and the acceptor’s LUS + 1 is only ~0.2 eV. Thus, when an electron-hole pair migrates to the interface, the electron can predominantly move from the donor’s LUS to the acceptor’s LUS + 1. In contrast, the transfer of holes is suppressed owing to a large energy gap (~1.4 eV) between two HOSs of the donor and the acceptor. To analyze the energetics of excited charge carriers, we project the reduced single-electron density matrix (RSDM) onto localized FMOs of T6 or C60 to obtain their occupations using Eqs. 12 to 14. The evolutions of populations of the states of interest, i.e., HOS, LUS, and LUS + 1, are plotted in Fig. 3B. We first analyze the evolutions of charge carriers with a simple partition of the system into the donor and the acceptor. On the donor side (left panel), it is observed that electrons start to be excited from HOS to LUS at ~7.5 fs, which is consistent with the time that an electron-hole pair is formed, as seen in Fig. 2. The population of LUS (HOS) reaches a maximum (minimum) and ΔnLUS ≈ 0.9 (ΔnHOS ≈ −0.9) at 15 fs; the applied laser pulse ends around the same time. On the acceptor side (right panel), a more substantial increase of nLUS + 1 is observed as compared to nLUS, indicating that the dissociation of the electron-hole pair at the D-A interface mainly results in electron transfer from the donor’s LUS to the acceptor’s LUS + 1. Note that the acceptor’s ΔnLUS+1 ≈ 0.1, much less than the donor’s ΔnLUS. A small electron-hole pair is observed between the HOS and LUS of the acceptor at 20 fs and remains beyond 100 fs, which can also be seen in the real-space distributions in Fig. 2A. The spatial resolution can be increased, and more detailed information can be obtained by constructing a localized basis set with the FMOs of each individual molecule. As shown in Fig. 3C, electron-hole pairs transfer from T6-2 to other donor molecules, and electrons transfer subsequently from C60-1 to C60-2, and lastly to C60-3. The time for an electron-hole pair to hop onto the adjacent T6 unit is ~15 fs. The oscillations of population were found to be related to the electronic coupling and C═C stretching mode in our previous work (, ). In addition, the trend of total population and (solid black lines) agrees well with and (dashed lines), showing the consistency of two different FMO partitioning schemes.

Fig. 3

FMO populations of donor and acceptor.

(A) DOS of donor T6 (blue line, left panel) and acceptor C60 (orange line, right panel) coupled to the environment. LDOS (middle panel) of the simulation region at steady state. The Fermi level is shifted to 0 eV (dashed line). (B) The evolution of the induced population Δn for selected MOs in the donor (left) and acceptor region (right). (C) The evolution of the induced population Δn for selected molecular orbitals of three donor units (left) and three acceptor units (right). Each curve corresponds to a single donor/acceptor unit. The color from dark (brown line) to light (light yellow line) represents the arrangement along the transport direction. Black lines, denoted as “T6s sum” (“C60s sum”), represent the change of total population of three T6 (C60) units. The red (brick red) dashed lines are identical to nLUS of the donor (nLUS+1 of the acceptor) shown in (B).

Consistency between spatial and energetic analysis

In Figs. 2 and 3, we study the photoinduced dynamics by analyzing the real-space induced charge distribution and the population of FMOs via σFMO. Next, we demonstrate in Fig. 4 that these analyses are consistent and give complementary information about the photoinduced charge dynamics in these organic materials. We define the induced charge by separately adding the positive and the negative induced Mulliken charge within a molecular fragment , and the photoexcited electrons/holes, , by summing the change of occupations of FMOs with orbital energies above or below the Fermi level. The oscillations due to rapid electronic motion, ΔqA(t) of atom A, can be averaged within a certain time interval T = Nδtwhere N is the number of time steps used for average and δt is the simulation time step used in the time propagation. A suitable choice of T (N) flattens the induced charge fluctuation and helps to observe the induced electrons and holes in real time and space. Figure 4 (A to D) presents the comparison of Δqpos/neg(t, N) in T6-3 with different averaging times. When N = 1, is plotted at every time step (see Fig. 4A), the overall profile of agrees well with except that there is strong oscillation in Δqpos/neg. This oscillation can be attributed to rapid electronic motion upon photoexcitation. As observed from the light-induced dynamics (see Fig. 5), electron-hole pairs move back and forth within a molecular unit as a coupled harmonic oscillator. When excited electrons and holes coincidently have similar spatial distributions, reaches minimum. On the other hand, at the instant when electrons and holes are localized at different parts of T6-3, reaches a peak (). Such oscillation prevents us from making distinct comparison between spatial and energetic information. To show a clearer picture of the real-time and real-space distributions of the induced electrons and holes, the oscillation of Δqpos/neg is smoothed out by time averaging. Figure 4 (B to D) shows the charge dynamics with different averaging times. For the result obtained with larger time intervals of 0.80 or 1.60 fs (N = 40, 80), the averaged curve can show merely the trend and fail to provide more subtle information. Therefore, 0.40 fs is chosen as the time interval for averaging, which provides a clear real-space dynamics picture of photoexcited electron-hole pairs, electrons, and holes, as shown in Fig. 2A.

Fig. 4

Population analysis of T6-3.

(A to D) Positive/negative induced Mulliken charge Δqpos/neg of T6-3 averaged over different time intervals: 0.02 fs (A), 0.40 fs (B), 0.80 fs (C), and 1.60 fs (D). Sum of occupations of FMOs on T6-3 with the energy below (Δnocc) and above (Δnvir) the Fermi level. (E to G) Comparison among Δqpos/neg (time interval = 0.4 fs), Δnocc/vir, and atomic charge contributed from occupied (Δqocc) and virtual (Δqvir) FMO occupation for the donor region (E), a donor unit T6-3 (F), and a thiophene ring in T6-3 (G).

Fig. 5

Induced Mulliken charge distribution of light-induced charge carrier dynamics without averaging.

(A) Positive and negative induced Mulliken charge Δqpos/neg of T6-3. (B) Induced Mulliken charge distribution in real space at selected time shots.

Figure 4E presents the time evolutions of and for the donor region. The result clearly shows that these two quantities agree well with each other, demonstrating the accuracy of the FMO approach. This allows us to further increase the resolution of our analysis by calculating induced charges on individual molecules and functional groups of the molecules. Figure 4F shows the induced charges () and electron/hole () of a single donor molecule T6-3 at the interface, which again shows similar agreement.

We may analyze the electronic dynamics of functional groups in a molecule, for instance, a thiophene ring of T6-3. We calculate the Mulliken charge on atom A based on FMOs

By adding up the contribution with respect to the FMO index p and then with respect to atoms, we can get that serves as the Δnocc/vir localized within a set of atoms D. () recovers () when contributions of all atoms in the donor region or a donor molecule are summed up, as shown in Fig. 4 (E and F). We illustrate the comparison for a thiophene ring of T6-3 in Fig. 4G. The Δqvir (yellow line) and Δqpos (red line) match each other quite well and so does the Δqocc (magenta line) and Δqneg (blue line), confirming the consistency of Δqpos/neg and Δqvir/occ (Δnvir/occ) for fragments of molecules (e.g., a functional group) and thus the representativeness of Δρ for electrons and holes in real space.

Induced Mulliken charge distribution in real space without coarse graining

We compute the induced Mulliken charge distribution Δρ(, t) in real space without averaging (N = 1, T = 0.02 fs) to illustrate the oscillation of induced Mulliken charge due to the varying distribution of electrons and holes. Selected time shots from oscillation on the T6-3 unit and a thiophene ring of T6-3 between 16.46 and 18.74 fs are shown in Fig. 5, along with the Δqpos/neg(t). At 16.64 fs (first time shot in Fig. 5A), Δqpos (Δqneg) is at its peak (valley), indicating the divergent spatial distribution of electrons and holes. Accordingly, a significant amount of induced charge distribution is observed in T6-3 (Fig. 5B). Then, electrons gradually move down in the direction of the T6 chain, and holes move up simultaneously. When the electron and hole’s distributions become closer at 16.98 fs, Δqpos/neg approaches minimum, yielding less induced charge on atoms and in real space. This is the moment that the electron and hole coincidently have similar spatial distribution. After that, electrons and holes continue to move and separate, and Δqpos/neg proceeds toward its maximum as well. At 17.58 fs, electrons localize at the lower region of the T6, which is different from that of the beginning time spot. This period is ~1.1 fs. Then, the electrons and holes continue to move as coupled harmonic oscillators until they return to the initial distribution. This period is about ~2.2 fs, corresponding to the excitation energy. The Fourier transform shows that the oscillation frequency of Δqpos/neg is around 4 eV (see fig. S1). Considering that the oscillation of excited electrons and holes is localized as described above, a factor of 0.5 should be applied to reveal the true frequency of 2 eV that corresponds to the excitation energy. The slower oscillation of the profile of with a period of ~15 fs is due to the exciton migration to the adjacent T6 unit, as evidenced by the electron-hole pair dynamics illustrated in Fig. 3C.

DISCUSSION

Using the Ehrenfest TDDFTB-OS method, we convincingly simulate the detailed dynamics of the photoexcitation of an electron-hole pair in the donor material, pair dissociation into electron and hole at the D-A interface, and electron and hole transport in respective conducting layers leading to emerging photocurrent. Compared to previous experimental and theoretical results, our findings provide much more details about the underlying mechanism of photocurrent generation, particularly the detailed spatial and temporal information of electron-hole pairs, electrons, holes, and photocurrent. Once it is excited instantaneously upon the absorption of a photon, an electron-hole pair migrates in the donor layer and may reach the D-A interface. The electron within the electron-hole pair that reaches the interface hops onto the selective unoccupied molecular orbitals (UMOs) of the acceptor materials, either because the orbital energies of the selective UMOs match that of the electron or, if they do not exactly match, because the electron may transfer some energy onto phonons and hops onto the energetically matched UMOs. Once this occurs, the transferred electron diffuses further along the electron-conducting layers, purely due to kinetics and entropy. As the transferred electron leaves the D-A interface, the corresponding hole leaves the interface as well, owing to the reducing Coulomb attraction of the leaving electron in addition to kinetics and entropy. These processes occur within a few femtoseconds, once the electron-hole pair arrives at the D-A interface. The leaving electrons and holes result in photocurrent. Detection of the photocurrent depends on how long it takes for the charge carriers to reach the detector. In our system, the photoinduced holes transport faster than the induced electrons, and the hole transport is modulated by the electron-hole pair migration and the electron-hole oscillation within the electron-hole pair, as shown in Fig. 2B. When an electron-hole pair dissociates at the D-A interface, the electron does not necessarily hop onto the acceptor’s LUMO, but rather onto the MO whose orbital energy matches or is closest to that of the electron. In our case, the electrons predominantly move from donor’s LUS to acceptor’s LUS + 1. Efficient electron-hole pair dissociation into free electrons and holes at the D-A interface can be explained by this matching in energy levels. This provides insight into designing donor and acceptor materials to enhance efficiencies of dissociation in organic or hybrid nanomaterials based on the electronic structures of these materials. Furthermore, besides the dissociation of electron-hole pairs at a D-A interface, excitation energy transfer from donor to acceptor is also observed on the same time scale in both Figs. 2A and 3B. Last, carriers leave the simulation region into the environment and are collected as current. This energy transfer at the D-A interface has received comparatively less attention from theoretical modeling. However, experimental work (, ) also shows evidence for ultrafast excitation energy transfer in polymer:fullerene blends. Since this also affects PCEs, it needs to be considered when developing solar cells with higher efficiency. The observed photophysical process is consistent with previous experimental and theoretical results (–, –), which reported that this light-induced process takes place on the scale of tens of femtoseconds.

Here, we report the results of our simulation at room temperature. To study the temperature effects, we used a simplified model as the full TDDFTB-OS calculation is quite time-consuming and found that the temperature does not have a strong effect on the photocurrent mechanism. This is because the energy differences between MOs across the D-A interface are on the order of 0.1 to 1 eV, which are much larger than the thermal energy around room temperature. Another issue is the variation in conformation. Simulations on several structural conformations with varying interfacial structures have been carried out. Although quantitatively different results such as charge separation rate across the D-A interface are obtained, the overall characteristics of electron-hole pair’s generation, migration, and separation across the D-A interface and subsequent photocurrent emergence remain qualitatively unchanged.

The two approaches based on σFMO and induced Mulliken charge provide consistent information of charge carrier dynamics. The coincidence among Δn computed from different partition schemes (Fig. 3C) confirms the FMO occupation analysis. The good agreement between Δn and Δq (Fig. 4, E and F) further confirms our findings on the spatial and temporal properties of photoinduced electronics dynamics. The comparison among Δq confined to various sizes of the spatial area (Fig. 4, E to G) further confirms the representativeness of the electron and hole distribution using induced Mulliken charge, allowing us to directly visualize the temporal evolution of electrons and holes in real space.

Ehrenfest dynamics is known for its deficiencies for isolated systems. However, its applicability has not yet been well studied for open quantum systems. Our results here show qualitatively consistent findings with experimental results (, ). We have investigated further the validity of Ehrenfest dynamics for open quantum systems via a simplified model and found qualitatively consistent results between Ehrenfest dynamics and full quantum mechanical calculations. The Ehrenfest TDDFT-OS or TDDFTB-OS method is a powerful tool to simulate the photophysical processes in photovoltaic devices. It provides detailed spatial and temporal information of electron-hole pairs, free electrons and free holes, and photoinduced current. It is applied to organic photovoltaic devices here and can be readily extended to inorganic optoelectronics as well.

MATERIALS AND METHODS

Ehrenfest TDDFTB-OS formalism

In the second-order self-consistent-charge DFTB (SCC-DFTB), the total energy reads ()where Erep is the pairwise repulsive energy between atoms, is the Fock operator of the referenced neural system, and Ψ is the single-particle wave functions expanded into a slater-type atomic orbitals (AOs). The last term (denoted as E2nd) contains second-order contribution of Coulomb and exchange-correlation energy due to the density fluctuation where γ is the two-electron integral related to Coulomb and exchange-correlation interaction. The corresponding Kohn-Sham Fock operation is given aswhere the off-diagonal matrix elements within an atom (A = B) are neglected (, ). The matrix elements of h0 and γ are parameterized through interpolating the DFT results. The parameter set used in this work is generated according to the GGA/PBE functional. It is worth mentioning that although the SCC-DFTB outperforms its non-SCC version in scenarios including charge transfer, further corrections [e.g., LC-DFTB ()] are needed to obtain more accurate results.

In the TDDFTB-OS method, the Kohn-Sham Fock operator h and RSDM σ are partitioned according to the left environment–simulation region–right environment alignment in real spacewhere D, L, and R represent the central simulation region and the left and right environments, respectively. Typically because of the large spatial distance, the left and right environments are assumed to have no interaction and the corresponding LR and RL blocks are set to zero. The equation of motion (EOM) of electrons in the simulation region is given aswhere SD is the overlap matrix of the simulation region; Qα is the dissipation term induced by the environment, which is evaluated using the nonequilibrium Green’s function (, ). To simplify the EOM and allow the simulation for large systems, the wide band limit approximation is applied to the electrodes’ self-energies (). In the TDDFT-OS formalism, exchange of energy between neighboring sites is included and electrons are allowed to cross the surfaces (SL and SR) between the simulation region and the environment. The time-dependent current at two interfaces SL and SR can be computed from the trace of QαThe (TD)DFTB method for open systems allows efficient simulation for realistic systems at the atomic level and has been reported to be successfully applied in several fields (, –).

The nuclear motions in the simulation region are dealt with classically using Ehrenfest dynamics. The corresponding EOM for the nuclei reads ()where A, B, and C are the nuclear indices; μ and ν are indices of AOs. A is the nuclear displacement of atom A; is the force exerted on nuclear A; is the Fock matrix of the neural system used in DFTB as reference; is the partial derivative of the overlap matrix SD with respect to the nuclear displacement A and ϕμ is AO; ΔqA is the induced Mulliken charge on atom A; γAB is the interaction between two point charges located on nuclear A and B; Erep is the pairwise repulsive energy between two atoms. Both γ and Erep are predetermined in DFTB parameters.

The time-dependent RSDM σD(t) contains overall electronic information about the simulation region. From the RSDM, we can analyze the real-time spatial distribution and the energetics of photoexcited electrons and holes, which allows us to visualize the charge dynamics in these organic materials.

Real-space induced Mulliken charge distribution analysis

The Mulliken charge () qA(t) of atom A at time t is computed from the RSDM σD(t) asand the induced Mulliken charge is defined aswhere is the ground-state Mulliken charge with the system geometry at time t computed by treating the simulation region as an isolated system. This equation characterizes electronic dynamics due to electron excitation and transfer through deducting the change of the ground-state charge (the latter square bracket) that reflects the charge redistribution on atoms owing to the structural variation. The induced Mulliken charge ΔqA(t) can then be projected onto real space by assuming a spherical symmetry distribution for each atom () (red curve in Fig. 1B)where UA is the Hubbard parameter of atom A and is the coordinate of the real space grid.

FMO occupation analysis

The RSDM σD(t) is projected onto the ground-state FMOs () usingHere, the block diagonal matrix CFMO(t) contains a series of FMOs coefficients CDi(t)where CDi(t) corresponds to the ground-state eigenstates of isolated fragment D at time t and is localized on D. S(t) is the overlap matrix for the nonorthogonal AO basis used in DFTB. σFMO(t) is the density matrix in the nonorthogonal FMO representation. We can choose various partition arrangements and apply them on the same system, so as to construct the corresponding CFMO(t) for analysis at a certain spatial resolution.

The projected σFMO(t) yields the local charge distribution contributed from local MOs. As shown in Fig. 1C, σFMO(t) can be divided into blocks. The off-diagonal blocks correspond to the electronic coherence between two fragments. Given a fragment D, the diagonal element represents the occupation number, n, of the pth molecular orbital of this fragment. We define the population of a set of FMOs by summing their occupationswhere Ξ is a preset orbital energy range to filter out the charge contributions from the FMOs of interests. ΔnΞ, the increase (decrease) of population nΞ, can further represent the change of electrons (holes) in the corresponding energy states.

Computation details

The structure within the shaded region (see Fig. 1A) is optimized using DFTB. The convergence criterion for the self-consistent field calculation is 10−10, and the maximum force is reduced to less than 10−4 Hartree·bohr−1. The TDDFTB simulation and all related properties are carried out and computed with the quantum chemistry program LODESTAR (, , , ). The fourth-order Runge-Kutta algorithm is used for the time propagation with a time step of 0.02 fs. In the simulation, the optical excitation is introduced by applying an external electric field (Fig. 1A, red curve) in the y direction that interacts with the local electric dipole moment of T6-2. The external field is modeled as a short Gaussian-enveloped cosine optical pulse with E0 = 0.06 eV·bohr−1; ℏω = 2.1 eV, which matches the strongest absorption peak of donor thiophene molecules; tc = 10 fs; and FWHM = 4 fs. The nuclear motions of the two T6 and two C60 molecular units within the shaded region are treated with the Ehrenfest TDDFTB-OS method, while the other ones are kept fixed. The dynamics in the shaded region are carried out with initial velocities that are set randomly following the Maxwell distribution at room temperature. We have carried out several simulations starting from different initial nuclear velocities, and similar charge transfer processes were observed. A typical trajectory is chosen to illustrate the real-space (spatial) and FMO (energetic) analysis. All simulations were carried out on a single node with a single Intel Core i7-8700K 3.70-GHz processor. Note that we can expand the simulation region to include more donor/acceptor units or use larger structures given modern high-performance computing clusters.