Quantum simulations of charge-separation at a model donor-acceptor interface: role of delocalization and local packing

We investigate the electronic dynamics of a model organic photovoltaic (OPV) system consisting of polyphenylene vinylene (PPV) oligomers and a [6,6]-phenyl C61-butyric acid methylester (PCBM) blend using a mixed molecular mechanics/quantum mechanics (MM/QM) approach. Using a heuristic model that connects energy gap fluctuations to the average electronic couplings and decoherence times, we provide and estimate of the state-to-state internal conversion rates within the manifold of the lowest few electronic excitations. We show that the electronic dynamics of the OPV are dramatically altered by varying the positions of the molecules simulated at the interface. The lowest few excited states of the model interface rapidly mix allowing low frequency C-C out of plain torsions to modulate the potential energy surface such that the system can sample both intermolecular charge-transfer and charge-separated electronic configurations on sub 100 fs time scales. Our simulations support an emerging picture of carrier generation in OPV systems in which interfacial electronic states can rapidly decay into charge-separated and current producing states via coupling to vibronic degrees of freedom.


I. INTRODUCTION
Advances in both materials and device fabrication have lead to the development of highly efficient organic polymer-based photovoltaic cell (OPV) in which the power conversion efficiency is in excess of 10-11% under standard solar illumination 1 and efficiencies as high as 12% in multi-junction OPVs. This increase in power conversion efficiency indicates that mobile charge carriers can be efficiently generated and collected in well-optimized devices; however, the underlying photophysical mechanism for converting highly-bound molecular (Frenkel) excitons into mobile and asymptotically free photocarriers remains elusive in spite of vigorous, multidisciplinary research activity. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] The generic photophysical pathway that underlay the generation of mobile charge carriers in the OPV begins with an exciton being formed inside the system. The exciton dissociates at the interface into a charge transfer state with a small electron/hole separation. The charge transfer state is bound by Coulombic attraction to the interface, impeding the ability of the particles to migrate away from the interface to form free charge carriers. Ultrafast spectroscopic measurements on OPV systems have reported that charged photoexcitations are generated on ≤ 100-fs 2,12,15,[19][20][21][22] time-scales, despite the strong Coulombic attraction due to the low dielectric constant prevalent in OPV's. Experiments by Gelinas et al., in which Stark-effect signatures in transient absorption spectra were analysed to probe the local electric field as charge separation proceeds, indicate that electrons and holes separate by as much as 40Å over the first 100 fs and evolve further on pico-second time scales to produce unbound and hence freely mobile charge pairs. 4 a) bittner@uh.edu; k2.chem.uh.edu Bittner and Silva recently presented a fully quantum dynamical model of photo-induced charge-fission at a polymeric type-II heterojunction interface. 23 This model supposes that the energy level fluctuations due to bulk atomic motions create the resonant conditions for coherent separation of electrons and holes via long-range tunnelling. Simulations based upon lattice models reveal that such fluctuations lead to strong quantum mechanical coupling between excitonic states produced near the interface and unbound electron/hole scattering states resulting in a strong enhancement of the decay rate of photoexcitations into unbound polaronic states.
A microscopic model of the interface is required to understand the mechanisms that promote ultra-fast charge separation. This requires knowing the morphology and the finest details of the interface and its electronic structure. Currently ab inito methods provide a good tool for exploring the formation of charge transfer states in donor/acceptor pairs, however, it is ill suited for the task of simulating the formation of charge separated states, because it is to computationally cumbersome to expand the interface from a single donor/acceptor pair into the size needed to allow charges separation to occur.
In this paper, we take a molecular mechanics/quantum mechanics (MM/QM) approach to study Poly(pphenylene vinylene) (PPV)/ Phenyl-C61-butyric acid methyl ester (PCBM) heterojunctions. Polymer microstructural probes have revealed general relationships between disorder, aggregation and electronic properties in polymeric semiconductors. 24 The distribution of torsion angles for the PPV molecules at the interface are larger than in the bulk, adding to the structural disorder of the PPV molecules closest to the interface. Moreover, aggregation (ordering) can be perturbed by the blend-ratio and composition of the donor/acceptor polymers. 24 We explore the effect that the positioning of the molecules at the interface has on the electronic properties and estimate the state-state transition rate con-arXiv:1612.06875v1 [cond-mat.mtrl-sci] 20 Dec 2016 stants of the systems.

II. METHODS
Our simulations employ a modified version of the TINKER molecular dynamics (MD) package 25 in which the MM3 26 intra-molecular bonding parameters are allowed to vary with the local π-electronic density as described by a Parisier-Parr-Pople (PPP) semi-empirical Hamiltonian 27,28 . Similar approaches have been described by Rossky 29 and Warshal 30 to include electronic dynamics into an otherwise classical force field description.
At each time-step of the simulation, we compute the Hartree-Fock (HF) ground state for the π system and use configuration interaction (singles) (CI-S) to describe the lowest few π → π * excitations. Intermolecular interactions within the active space are introduced via nonbonding Coulombic coupling terms and static dispersion interactions contained within the MM3 forcefield. All electronic excitations are confined to the π-active orbitals. We used a total of 10 occupied and 10 unoccupied Hartree Fock molecular orbitals to construct the electron/hole configurations for the CI calculations. The excited state bond charge density matrix is constructed by assuming that electron densities are added to the virtual orbitals and hole densities are subtracted from the Hartree Fock ground state. The equilibrium bond orders are modified by where I BO and T BO are the initial equilibrium bond order and the rate of bond length increase with bond order decrease found in the MM3 parameter set, P is the excited state bond order found in the excited state bond charge density matrix. Modifying the equilibrium bond order allows the bond lengths in the MM calculation to change in response to local changes in bond order due to the migration of charge.
During the equilibration steps we assume the system to be in its electronic ground state, after which we excite the system to the first CIS excited state and allow the system to respond to the change in the electronic density within the adiabatic Born-Oppenheimer approximation. It is important to note that the excited state we prepare is not the state which carries the most oscillator strength to the ground nor do we account for non-adiabatic surface hopping-type transitions in our approach. [31][32][33] The dynamics simulations shown reflect the longer-time fate of the lowest-lying excited state populations and sample possible configurations that can be accessed by the system. The combination of a classical MD forcefield with a semi-empirical description of a selected few molecules within the system seems to be a suitable compromise between a fully ab initio approach which would be limited to only a few molecules and short simulation times and a fully classical MD description which would neglect any transient changes in the local electronic density. 22 We specifically chose three separate clusters of molecules to represent model bulk-hetrojunctions in order to study how varying the blend and positioning of the molecules affect the penetration of extended intra molecular electronic states into the bulk region. In each case study we selected a subgroup of PCBM molecules and PPV chains, treating the π-electrons in these units explicitly while the remaining molecules in the simulations are treated using the purely classical MM3 force-field. The number of π active PCBM molecules vary between each simulation, allowing each system to have a different blend ratio inside the π active system. The placement of PPV molecules vary in two of the simulations, changing the number of PPV molecules in direct contact with PCBM molecules, fundamentally changing the hetrojunction. Figure 1 shows representative configurations of the three cases studied. In each, the red and blue coloured spheres represent atoms included in the quantumchemical description. Each corresponds to a periodic simulation cell in the xy plane following equilibration at 100K and 1 atm pressure. In Case A, we selected 2 interfacial PCBMs and 3 nearby π-active PPV oligomers that penetrate into the bulk polymer region, including a total of 230 carbon 2p z orbitals. In Case B, we selected 3 PCBM and 3 nearby PPV oligomers expanding the π active molecules that form the inter-facial hetrojunction, including a total of 288 2p z orbitals. The system has the same boundary conditions as case A and is set up such that the main difference is in the placement of the PPV oligomers. In simulation C we selected 1 PCBM and 3 PPV oligomers consisting of 172 carbon 2p Z that penetrate into the bulk using and employ periodic boundary conditions in xyz. This simulation was set up to be very similar to case A, only adding a single PCBM molecule to the inter-facial region. The case studies are representative of typical interfacial configurations and in no way are a comprehensive sampling.

III. RESULTS
Over the course of eight 10 ps simulations, the lowest lying CIS excitation samples a variety of adiabatic states ranging from localized excitons to chargeseparated, charge-transfer and de-localized configurations. Figure 2 shows the configurations of four distinct adiabatic states, where the probability of the system to be found in a specific state can be seen in Figure 3. We categorize the states into four types, (EX) represents the exciton states, characterized as having ¿50% of the electron/hole density on a single molecule. The exciton state populates ≈ 58% of the states for simulation A, ≈ 32% for simulation B and ≈ 52% for simulation C. (CT) represents the charge-transfer states, characterized as having ¿50% of the electron/hole density occupying adja- cent molecules. The charge-transfer state populates ≈ 14% of the states for simulation A, ≈ 19% for simulation B and ≈ 12% for simulation C. (CS) represents the charge separated states, characterized as having ¿50% of the electron/hole density occupying a PCBM and a PPV separated by a single molecule. The charge separated state populates ≈ 26% for simulation A, ≈ 49% for simulation B and ≈ 26% for simulation C. (DLOC) represents the de-localized states, characterized as having ¡50% of the electron/hole density on a single PPV or PCBM molecule.
The energies of the lowest CIS state following excitation at t = 0 fs for simulation A and B are shown in Figure 4 (Top). After excitation at t = 0 there is very little energetic relaxation in all of the systems simulated. The simulations appear to cycle through many adiabatic states in a short period of time leaving the impression of a weak electron-phonon coupling. This can be rationalized as the electron/hole density often de-localize over multiple molecules and many conjugated C-C bonds. Another striking affect of the systems is the large number of avoided crossings that occur between the lowest lying states. There is also a 20 fs oscillation in the CI energies, driving the systems excited states into many regions of strong coupling. The 20 fs oscillation also appears in the autocorrelation and and the Fourier transform of the gap energies, contributing the C=C bond stretching modes around ≈ 1600 cm −1 . The oscillation is contributed to small thermally activated fluctuations within the simulation, showing that even at 100 K the thermal fluctuations possess sufficient energy to bring these states into regions of strong electronic coupling.
In Figure 4   of excited states that can easily be brought into strong electronic coupling by small fluctuations in the CI energy of the system. The small average band gap and rapid oscillatory nature of the CI energies facilitate the systems ability to rapidly sample a great many different electronic configurations over the course of the simulation.
We next consider the origins of the energy fluctuations evidenced in Figure 4 (Top). While we only show two 200 fs segments of eight 10 ps simulations over this period, one can see that the CI energies are modulated and cover a small range. The autocorrelation plots of the band gap energies, shown in Figure 5, show that the cor-relation times for the three simulations are very short ≈ 8 fs meaning that the system changes rapidly enough that the oscillations observed are independent of one another. By taking the Fourier transform of the gap the IR active modes that contribute to the modulation of the CI energies inside of the systems are found as shown in Figure 5. The modulation of the CI energy appears to be heavily dependent upon the torsion, C=C, and C-H stretching modes. In each of the plots three distinct regions can be seen, the low frequency torsional modes occur between 200 and 500 cm −1 , the C=C stretching modes occur between 1300 and 1800 cm −1 and the C-H stretching modes occur between 2800 and 3300 cm −1 . We conclude that small-scale vibronic fluctuations in the molecular structures and orientations produce significant energetic overlap between different adiabatic states to drive the system from purely excitonic to purely charge-transfer on a rapid time-scale. This is evidenced in the progression of the CI energies, as small fluctuations in these modes can easily bring the excited states into strong coupling regimes.

A. Estimating state to state rates
We can estimate the state → state rates using the model by Bittner 23,34 . Consider a two state system with coupling λ in which the energy gap ∆(t) fluctuates in time around its average∆. In a two state basis the Hamiltonian can be written as whereσ k are Pauli matrices. Note that Eq. 1 transformed such that fluctuations are in the off-diagonal coupling, becomes where ∆ 0 =∆ + λ and δV (t) = 0. The fluctuations in the electronic energy levels are attributed to thermal and bond-vibrational motions of polymer chains which can be related to the spectral density, S(ω) viā Averaging over the environmental noise, we can write the average energy gap ashΩ= ∆ 2 0 +V 2 with eigenstates |Ψ+ = cos θ|0 + sin θ|1 (5) where tan 2θ = |V |/∆ 0 defines the mixing angle between original kets. Consequently, by analyzing energy gap fluctuations, we can obtain an estimate of both the coupling between states as well as transition rates. To estimate the average transition rate between states the equations of motion for the reduced density matrix for τ 1 and τ 2 have been introduced as the lifetimes of each state and T d is the decoherence time for the quantum superposition. The decoherence time can be related to the spectral density via T −1 d =V /h. Taking T d to be short compared to the lifetimes of each state, we can write the population of the initial states as where k is the average state to state transition rate. If we integrate over all time we obtain an equation of the form suggesting a form for the exact solution of Eqs.6. Taking the Laplace transform of the equations of motion (Eqs.6) and assuming that our initial population is in state 1 (ρ 11 (0) = 1) the equations of motion become a series of algebraic equations which after a bit of algebra gives a rate constant of the form.
According to the model outlined above the mean (∆ 0 ) and variance (V) of the energy gap distributions shown in Fig. 4 can be used as input to estimate the state to state transition rate for a two level system. We take T −1 d ≈V /h as an estimate of the decoherence time and we introduce τ as the natural lifetimes of each state. The results are shown in Table 1. The estimated transition rates are consistent with the observations that the systems rapidly sample a wide number of possible configurations over the course of the molecular dynamics simulation. On average, the state to state couplings of 56 meV for simulation A and 48 meV for simulation B are comparable to the average energy gaps between the lowest excited states. The strong electronic coupling allows for rapid transitions; however, larger couplings also imply shorter electronic decoherence times, effectively quenching the ability of charges to separate by tunnelling.

IV. CONCLUSIONS
We present here the results of hybrid QM/MM simulations of the excited states of model PPV/PCBM heterojunction interfaces. Our results indicate that varying the blend ratio and placement of the molecules comprising the heterojunction greatly affect the distribution of states yet have little affect upon the rate constants of the system. We also propose that thermal noise can rapidly change the character of the lowest lying excited state from purely excitonic to charge separated on a time scale of sub 100 fs.
Simulations A and C have a very similar placement of molecules, only differing in that simulation A adds a PCBM molecule to the heterojunction. The addition of the PCBM only slightly changes the distribution of states as seen in Figure 4. The exciton states continue to be the most favored state inside the system, even slightly increasing in probability, while the delocalized states slightly decrease in probability. Simulation B completely changes the heterojunction, placing three PPV and two PCBM molecules at the interface as seen in Figure 1. The states generated by simulation B are radically different from those seen in simulation A and C as seen in Figure 4. The probability of finding the system in the exciton state is dramatically reduced while the charge separated state becomes predominant. This result is quite interesting as it highlights that the complexity of simulating heterojunctions resides not only in the size of the system but on how the donor/accepter interface is chosen.
All three systems start with an exciton localized on the PCBM and dissociating into a charge transfer state with the hole (or electron) delocalized over multiple polymer units before localizing to form charge separated states. There are a wide range of electronic states tightly clustered within a small energy band, allowing small changes in local bond lengths to have a dramatic role in modulating the electronic couplings between excited states. We speculate that the dramatic shift in population seen in simulation B can be caused by disorder in the PPV molecules reducing the band gap by 20 meV. The PPV molecules comprising the interface region undergo large distortions in the C-C torsion angles allowing the molecules to cycle through a larger range of configurations inside of a short time interval. The presence of more π active PPV molecules at the interface also appear to lead to more avoided crossing regions and the ability of the system to more efficiently dissociate excitons into charge transfer and charge separated states to a distance to where their Coulombic attraction is comparable to the thermal energy. While the finite size of our system prevents further dissociation of the charges, the results are suggestive that such interstate crossing events driven by bond-fluctuations can efficiently separate the charges. The results presented here corroborate recent ultra-fast experimental evidence suggesting that free polarons can form on ultra fast time scales (sub 100 fs) and that thermally activated low frequency torsional modes are key in effective electron hole separation in PPV/PCBM heterojunctions.